CC BY
28
УДК 534.13+532.582+532.581
A.A. Коробкин, A.A. Папин, К.А. Шишмарев
Аналитическое и численное исследование квазиизотермической задачи взаимодействия ледового покрова канала и поверхностных волн*
A.A. Korobkin, A.A. Papin, K.A. Shishmarev
Analytical and Numerical Study on Quasi-Isothermal Problem of Iteraction between Ice Cover of Channel and Surface Waves
Для уравнений совместного движения воды и ледового покрова канала формулируется начально-краевая квазиизотермическая задача в случае, когда ледовый покров закреплен на стенках канала. Проводится аналитическое и численное исследование задачи.
Ключевые слова: ледовый покров канала, поверхностные волны, дисперсионная зависимость.
Initial-boundary quasi-isothermal problem when ice cover is fixed on the channel walls, equations of combined motion of the water and the channel ice sheet is consider. Analytical and numerical study of the problem is conducted.
Key words: ice sheet channel, surface waves,
dispersion relation.
1. Постановка задачи. В рамках линейного приближения рассматривается безвихревое движение жидкости в канале шириной Ь и высотой Н, покрытом льдом. Потенциал <р{х,у, г,Ь) скоростей течения удовлетворяет задаче
дУ д2ср д2ср _
дх2 ду2 дг2 ’
(-Н < г < 0, 0 < у < Ъ, —со < х < оо);
тг- = °, (у = М);
ду
£ = », (; = -Я);
здесь ьз{х, у, £) - прогиб ледовой поверхности.
Интеграл Бернулли берется в виде
-р + <рг+ ди> = О,
Р
где р - давление; р - плотность воды; д - ускорение свободного падения [1].
Функция ьз{х,у^) удовлетворяет уравнению вида [2, 3]
ти>и + -ОУ4«; - А(С}(х, у)ги) = р(х, у, 0, ¿),
(О < у < 6, —оо < х < оо)
*Работа выполнена при финансовой поддержке аналитической ведомственной целевой программы «Развитие научного потенциала высшей школы (2009—2011 гг.)»(проект №2.2.2.4/4278), а также федеральной целевой программы «Научные и научно-педагогические кадры инновационной России» на 2009—2013 гг. (государственные контракты №14.740.11.0355, №14.740.11.0878).
и условиям
w = wy= 0, (у = 0,6),
где то - масса покрова на единицу площади (то = Pihi; hi - толщина покрова); Q(x, у) - температурные напряжения, связанные с образованием ледового покрова; V - оператор градиента.
Изгибная жесткость D вычисляется по формуле D = Ehf/[ 12(1 — г/2)], где Е - модуль Юнга и v - коэффициент Пуассона для льда. Толщина льда hi и его плотность pi считаются постоянными. Краевые условия для w означают, что ледовое покрытие приморожено к стенкам канала.
Решение сформулированной задачи для w ищется в виде w{x,y,t) = AF(у/Ь) sm(kx -\- u>t), где А = const - амплитуда волны; к = const -волновое число; из = const - частота волны; искомая функция F(y), у = удовлетворяет условиям F = F =0 при у = 0, у = 1. Температурные напряжения Q считаются постоянными, следовательно, для каждой заданной температуры получим собственную функцию w.
В безразмерных переменных (знак тильда опускается) исходная задача для w сводится к следующей задаче относительно F и Л
0_ж0 + (*-4 + вт!,) =
= А(а + A(K))F, (0<S<1); (1)
F = Fy = 0, (у = 0,1), где if = кЪ- 2К= (2if2 + ^);a = ^;A=^;
Учитывая условие = го( и вид функции го, потенциал скоростей 95 будем искать в виде (р = Аф(у, г) соз(/гж + а;£), где функция ф(у, г) удовлетворяет задаче
+ ^ = К^ < °’ ° < У < !)’
дф , .
-т^—0, (у — 0,1),
Й = 0’ (г=“',):
^ = т, = о),
где А(К) < ^ >= ф{у1 0) и к = Н/Ъ.
Требуется, при заданных параметрах К, ее, а найти собственные значения Л и соответствующие функции -Р(у), а также исследовать зависимость со {к) при уменьшении толщины льда и сравнить вычисления с дисперсионным соотношением, полученным для задачи со свободной поверхностью [4]. Полученный результат будет использован для изучения возможности разрушения ледового покрова в канале движущимся судном на воздушной подушке. Ранее подобные задачи о ледовом покрове в канале не рассматривались. Задача без учета температурных напряжений исследовалась в [4]. Задача без учета температурных напряжений с одной вертикальной стенкой исследовалась в [5].
2. Алгоритм численного решения задачи (1). Определим ортонормированную систему функций {'фп(у)}^= 1 на интервале (0, 1) как решения следующей спектральной задачи
Фп = РІ’Фп,
(0<у < 1);
Фп = Фп = 0, (у = о, 1);
1 1
І Фііу)Лу = 1, ¡ФпФт<1у = 0 (п^то).
(2)
Функции фп(у) называются балочными функциями. Они дают формы колебаний балки постоянного поперечного сечения, защемленной на обеих концах; соответствующие собственные значения /3„ связаны с собственными частотами колебаний балки.
Решение задачи (1) будем искать в виде
Р{у) = ^^апФп(у)■
(3)
П= 1
Подставляя (3) в (1) и используя (2), получим систему линейных однородных уравнений относительно коэффициентов ап и параметра Л. Система имеет бесконечное число уравнений. При ее приближенном решении ограничимся конечным числом N уравнений и находим собственные числа
полученной матрицы. Проводя расчеты при различных 14, получаем последовательности вида
\-iVo ^N0 + 1 \Лг0 + 2
Л1 I Л1 I Л1 I
лЛГ0 л N0 +1
Л2 , Л2 ,
\ N о \ N0 +1
к 1 к > •">
сходимость которых при фиксированном к и растущем N проверяется численно.
Функции фп(у) вычисляются аналитически, а собственные значения /3„ - численно.
Подставляя (3) в (1) и умножая обе части полученного равенства на фт(у), после интегрирования у с учетом (2) приходим к равенству:
Рт + а,пСПт + (К4" + ае)
П= 1
(4)
где
1
Ьт = j фт(у)Ф(у,0)Лу',
о
1 1 1
Спт = - I ФпФтЛу = J Ф'пФ'т<1у = ~ J фпф'^<1у. 0 0 о
Для вычисления Ът рассмотрим вспомогательную задачу относительно ф(у, г):
У2ф = К2ф, (0 < у < 1, —к < 2 < 0); дф
= у = М;
ду
дФ п ,
— =0, г=-к;
ОХ
(г = 0,0 < у < 1),
где -Р(у) определен в (3).
Функцию ф(у, г) представим в виде
Ф{у,г) = ^2акфк(у,г)
к=1
где новые искомые функции фк(у, %) удовлетворяют следующей задаче:
V фк = К фк, (0 < у < 1, -к < г < 0);
дфк
1р=0, у = 0,1; ду
дфк
дг
0, г = —к:
^- = Фк{у), = 0,0 < у < 1).
Тогда для Ът имеем
1 1
Ът = / 'фтфйу = ТХ=\ ак I ФтФк(у, 0)<1у =
О О
СО 1
= £<**/ ЧгФ^У-
к= 1 О
Повторяя анализ, проведенный в [4], находим: 1
I Г<^Т^фкЛу = I (V фк V фт + К2фкфт)с1усЬ.
О П
Следовательно,
[ Э(^т Л, /7 (\ 8г! - М
/ ~о^Фк У = / Фт-г^-ау = Мтк.
о о
Поэтому приходим к равенству
СО
= ^ ^ М-ткак ? Мтк = Мкт • й=1
Используя введенные элементы С*пт и Мт&, равенство (4) представим в виде
__г\ со
(/^т + -^4 + 3е)ат + 2.К" ^ Стпап
и 1
со
= А(а;ат -|- АГтп(2п),
п= 1
или в матричной форме
(Ъ + 2¥2С7) а = А(а/ + М)а, (6)
где а = (ах, а2, аз, ...)т - вектор; В = (Иад{/3^ + ТГ4 + ее, /?| + + 5!,...} - диагональная матрица;
С = {Спт)п,т=1 и М = {М„т}~т=1 - симметричные матрицы.
Умножим обе части (6) на (а/ + М)-1 и обозначим
А = (а1 + М)-1 ■ (О + 2К2.
Приходим к задаче вида
[А — Х1)а = 0. (7)
Элементы матрицы А вычисляются численно по указанным выше формулам. После этого находятся собственные значения Л„ и вычисляется функция Рп.
3. Численное решение задачи (7). Будем решать численно задачу (7) при следующих значениях параметров: К = К = 100; аё = 144.452; а = 0.0183.
Решение в вычислительной среде МаЛаЬ строится в несколько этапов:
1. Матрица С не зависит от параметров задачи, вычисляется заранее и сохраняется в файле.
Это же справедливо для последовательности /3„ и функций фп(у).
2. Вычисляются матрицы М и Б, зависящие от параметров задачи.
3. Вычисляется матрица А заданной размерности.
4. Вычисляются собственные значения Л„ и собственные вектора ап.
5. По найденным значениям Л„ находятся частоты изп.
В таблице 1 приведены 4 первых значения изп, соответствующие полученным собственным значениям Л„ для разного числа уравнений в (7)
Таблица 1
Значения изп в зависимости от числа уравнений
5 уравнений 20 уравнений 50 уравнений
113,9748 205,7769 328,8829 419,2954 105,6593 193,6093 305,2569 397,1653 105,1798 192,8032 302,798 394,6687
Заметим, что разность чисел, полученных для матриц размерности 50 х 50 и 20 х 20, на порядок меньше от разности чисел, полученных для 50 х 50 и 5 х 5 . Разность первых двух частот для матриц 50 х 50 и 20 х 20 по модулю не превышает 1, когда разности 3 и 4 чисел больше разности первого числа в 6 раз. Для более высокой точности полученных значений из требуется построение матриц с большим числом уравнений в (7).
На рисунках 1-2 представлены графики изменения прогиба ледового канала го, соответствующие частоте волны изп = 105, 6593, при t = const =
1, ж Є [—10,10] и х = const = 0, і Є [0,0.1] соответственно.
Будем увеличивать волновое число к от 1 до 4 с шагом по к равным 0,1. Получим при каждом к свои собственные частоты изп. Начнем уменьшать толщину льда h и сравним полученные последовательности изп с частотами, полученными для задачи со свободной поверхностью [4]
из2 = дл/к2 + /¿2 tanh[л/к2 + /л^Н],
где к - волновое число; из - частота волны; Н -высота канала; ¡лп = тт/Ъ; b - ширина канала.
Результаты вычислений представлены на рисунке 3.
Рис. 3. График функции ш(к) в зависимости от толщины канала /?
Библиографический список
1. Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. — М., 1963. - Ч. I, II.
2. Тимошенко С.П., Гудьер Дж. Теория упругости. — М., 1975.
3. Nowacki W. Thermoelasticity. - 2nd Edition.
— Pergamon Press Ltd., 1986.
4. Коробкин A.A., Папин A.A., Шишма-
рев К.А. Поведение ледового покрова канала под действием поверхностных волн // Известия Алт-ГУ. - 2012. - № 1/1.
5. Brocklehurst P., Korobkin A.A., Parau E.I. Hydroelastic wave difraction by a vertical cylinder. Philos Transact A Math Phys Eng Sci. — 2011. — Vol. 369, №19.
