CC BY
26
Вычислительные технологии
Том 15, № 4, 2010
О численном моделировании динамики локального возмущения поля плотности в устойчиво стратифицированной среде*
А. Н. Зудин, Г. Г. Черных Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: [email protected]
Построены линейные и нелинейные численные модели динамики локального возмущения поля плотности в устойчиво стратифицированной среде. Оценено влияние вязкости на процесс генерации и распространения внутренних волн, генерируемых локальным возмущением поля плотности в пикноклине. Рассмотрена задача о динамике локального возмущения поля плотности при наличии волнового фона.
Ключевые слова: устойчиво стратифицированная жидкость, локальное возмущение поля плотности, внутренние волны, линейные и нелинейные модели.
Введение
Рассматривается плоское течение, генерируемое локальным возмущением поля плотности в стратифицированной среде. Задача представляет интерес в связи с изучением ряда геофизических явлений [1-5]. Подробная библиография работ данного направления содержится в [3, 6-19]. Анализ известных исследований показывает, что результаты численного исследования динамики локального возмущения поля плотности в устойчиво стратифицированной жидкости [10-14] недостаточно полны. Так, отсутствуют данные подробного численного анализа применимости линейных и нелинейных математических моделей к расчету характеристик течения. В настоящей работе представлены результаты расчетов, иллюстрирующие эволюцию области частично перемешанной жидкости в устойчиво стратифицированной среде, выполнены численные эксперименты по оценке применимости предложенных математических моделей, рассмотрена задача о динамике локального возмущения поля плотности в жидкости с нелинейной стратификацией при наличии волнового фона. Исследования являются развитием работ [10, 12-14, 19].
1. Постановка задачи с применением уравнений Эйлера в приближении Обербека—Буссинеска
Одной из основных математических моделей динамики нестационарных стратифицированных течений являются уравнения Эйлера в приближении Обербека—Буссинеска. С использованием переменных функция тока—завихренность уравнения Эйлера могут быть записаны в виде
* Работа выполнена при финансовой поддержке РФФИ (грант № 07-01-00363) и гранта Президента РФ № НШ-931.2008.9.
дш дш дш g dpi ,
— +u— + v— =--—, (1)
dt dx dy p0 dx
dx2 dy2 dpi dpi dpi dps
— +u— + v— + v— = 0. (3)
dt dx dy dy
x y y
силы тяжести; pi = p — ps, ps = ps(y) — плотность невозмущенной среды вне зоны смешения; p0 = ps(0); u,v — горизонтальная и вертикальная составляющие скорости жидкости: u = дф/дy, v = — d^/dx; g — ускорение силы тяжести. Стратификация среды предполагается устойчивой, т.е. (l/po)dps/dy < 0,
Граничные и начальные условия в рассматриваемой задаче следующие:
ф = ш = pi = 0, x2 + y2 ^ то, (4)
p = p0(x,y) = po s 1 —
„Л _ „ J 1 Po — Ps(y)
1 _ ее-(1-25г)8
Р0
—то < х < то, —то < у < то, £ = 0, (5)
ф = ш = 0, —то < х < то, —то < у < то, £ = 0. (6)
Здесь е = eonst = 0.75. Условие (6) соответствует и(0,х, у) = 0 ^(0,х,у) = 0.
Приведенное выше начальное распределение плотности соответствует неполному перемешиванию жидкости (известно [6], что при эволюции турбулентных пятен в стратифицированной жидкости в зоне турбулентного смешения происходит неполное перемешивание).
При численном решении задачи система уравнений (1)-(3) предварительно приводится к безразмерному виду путем использования характерных масштабов длины Я (радиуса области смешения) и масштаба времени Т* = 1 /у/ад (а = —(l/po)dps/dy при у = 0); кроме того, используется представление р1 = р0аЯрь Для анализа результатов расчетов удобно ввести масштаб времени Т = 2пТ — период Вяйсала—Брента. Ниже, по возможности, знак обезразмеривания (~) опускается.
Численное интегрирование системы (1)-(3) осуществляется с помощью подхода, основанного на эйлерово-лагранжевых переменных [9, 10, 12, 14]:
£ = х' = £ = х, у' = п : = 0, (7)
п(0,х,у) = п°(х,у), —то < х < то, —то < у < то,
п(£,х,у) ^ (у), х2 + у2 У то, £> 0. (8)
При этом уравнения (1)-(3) переписываются в виде
dy дф dt' =
1 J д дф д
L (!kY +Л ?±
J\\dU / дг,
д f ду дф\ д f ду дф\
дш 1 дф дш dp 1 др ду dt' J дп д£ д£ J дп 9(t,x,y) = ду_ дп
= ш, (11)
(12)
J
др,£,п)
В граничных и начальных условиях (4), (5) вместо функций р3(у), р0(х,у) появятся функции р3(п)-, Р°(С, п)> получаемые из ря(у), Р0(х,у) путем подстановки в них значений У = Я_2(С, п)-, являющихся решением уравнения п(0,х,у) = п0(х,у) (п°(х,у) предполагается монотонной функцией у при фиксированном х), К условиям (4)-(6) добавятся начальное условие
и граничное условие
У(0,С,П) = Q2(i,n)
У = Яз(п), С2 + п2 ^го, t> 0.
(13)
(14)
п=
нии р = const. Тогда уравнение (9) будет выполняться тождественно, так как р = р(п)-В результате перехода к эйлерово-лагранжевой системе координат в уравнении (12) отсутствует конвективный перенос в вертикальном направлении. Это обстоятельство позволяет избавиться от соответствующих погрешностей при конечно-разностной аппроксимации.
Конечно-разностный алгоритм, основанный на идее метода предиктор-корректор, и его тестирование подробно описаны в работах [9, 12, 14].
2. Внутренние волны, генерируемые локальным возмущением поля плотности в жидкости с нелинейным распределением плотности по глубине
С целью выяснения влияния нелинейности стратификации невозмущенной среды на течение, генерируемое при коллапсе области перемешанной жидкости, выполнена серия расчетов для ряда распределений р3(у)- Применялись уравнения Эйлера в приближении Обербека—Буссинеска (1)-(3). Начальное распределение поля плотности задавалось в виде (5).
На рис. 1 приведены результаты расчетов, иллюстрирующие картину течения в случае нулевого начального распределения вихря и функции тока и задания р3 (у) вида
р8(у) = ро - у, (15)
Рз(у) = Ро - в • Ь&пЦу/р) в = в/К =1.0, 0.5, 0.1, (16)
р(У) = { ро - у, -2Ъ - у - 11 (17)
РЛу) \ Ро - (1&пЦ(у - 2)/0.1) +1&пЦ9) + 1.1), 1.1 - у - 3. К >
При распределении р3(у), определяемом (15), процесс возникновения и распространения внутренних волн характеризуется формированием системы вихрей в каждом
Рис. 1. Линии тока, полученные с применением приближения Обербека—Буссинеска для различных профилей плотности ps (y), t/T = 2
квадранте плоскости xy [6, 7]. Число вихрей с ростом времени увеличивается. В частности, для варианта (15) при t/T = 2 таких вихрей четыре в каждом квадранте, для вариантов (16) при /3 = 1, 0.5 и 0.1 — три, два и один соответственно. Это иллюстрирует рис. 1, где на момент времени t/T = 2 показаны линии тока ф = -0T*/R2 = const. Здесь
и ниже линии тока представлены уровнями: 1 --0.0540, 2 --0.0420, 3 --0.0300,
4 --0.0180, 5 --0.0060, 6 --0.0000, 7 --0.0060, 8 --0.0180, 9 --0.0300, 10 -
—0.0420, 11 --0.0540. Рисунки 1, б—г отвечают распределению (16), ß = 1, 0.5 и 0.1.
В вариантах с в = 1 и 0.5 координатные линии n = const совпадают с линиями р = const с переменным интервалом, задаваемым следующим образом. Вначале строится последовательность значений y- = yj-1 + h0qj (yo = 0 Vj > 0 h0/R = 0.064, q = 1.035, j = 1, ..., 25); затем находятся линии равной плотности р = ps(y), которые и принимаются в качестве координатных в начальный момент времени. Для отрицательных значений y линии n = const строятся аналогичным образом (здесь уместно отметить, что в новой эйлерово-лагранжевой системе координат сетка в обоих направлениях выбиралась равномерной). В случае ß = 0.1 лини и n = const в начальный момент вре-
мени строятся следующим образом. Сначала находится линия п = п* такая, что (при y > 0) значение плотности р(п) при п > П* отличается от постоянной не более чем на 0.1 х 10-9. В области с переменной плотностью 0 < п < П* в качестве линий п = const для t/T = 0 выбираются линии p/aRpo = const с интервалом öpj = Pj — Pj_1 (pj = ps(yj), yj = yj-1 + hj, yo = 0, hj = 0.074, j = 1, ..., 10). В области п > п* ПРИ каждом xi = i х h1, i = 1, ..., 71 лини и при t/T = 0 строятся в соответствии с формулами
п(°,х^У^) = п, j = 11,17,
yi,j = yj (Xi) = yi,j_i + h*(Xi)q(j 10), h* (Xi) = (2.5 — yi,io )(1 — q)/(1 — q7), yi,io = У (0,х^п*), q = 1.1.
Для распределения плотности (16), ß = 0.1, рассчитывается также вариант на сетке с числом узлов по вертикали, в два раза большим рассмотренного. Получено, что максимальные значения функции тока при этом различаются не более чем на 7%. Если ввести величину H, определяемую из соотношения
ps(H) — po = °.99(B — po), B = lim ps(y),
то для /3 = 1, 0.5 и 0.1 величина H = H/R равна соответственно 2.64, 1.32, и 0.26. Видно, что переход от распределения (15) к распределению (16) приводит к существенному
I
II
= а
Ё / А
1 / \
1 / \
=
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 111111 м 1 111111111 111111111
10 12 14 16 18 xR
Рис. 2. Сохранение формы вихря (I) и волны (II) в идеальной жидкости для моментов времени t/T = 15 (а), 108 (б)
изменению волновой картины, генерируемой при коллапсе локального возмущения поля плотности в пикноклине: при малых значениях H течение характеризуется лишь одним вихрем в каждом квадранте плоскости xy (этот факт был хорошо известен авторам к моменту написания [10], однако основанные на эйлерово-лагранжевом подходе детальные численные эксперименты приводятся в настоящей работе). Форма вихря и его интенсивность (следовательно, форма и интенсивность волнового пакета) остаются неизменными вплоть до весьма больших значений времени, С целью детального анализа волновой картины расчеты в последнем случае дополнительно проводились в области 200R х 30R на сетке с числом ячеек 2001 х 301, На рис, 2, I представлены линии тока ф = const для ß = 0.1 е = 0.75 t/T = 15 (а), t/T = 108 (б). Соответствующее распределение р = const = р.Д0.0716), характеризующее внутреннюю волну, представлено на рис, 2, II, Видно, что интенсивность и форма основного возмущения практически не меняются, что свидетельствует о достаточно высокой надежности и эффективности предложенной численной модели,
3. Влияние вязкости на динамику локального возмущения поля плотности в пикноклине
Представляет интерес исследование влияния вязкости на характеристики внутренних волн, генерируемых локальным возмущением поля плотности в пикноклине. Расчеты проводились на основе уравнений (9)—(12), при этом (12) заменялось на следующее:
дш 1 дф дш др 1 др ду
dt' J дп J дп
J_}_\d_ дш д_1_ ~ReJ ] д£ +
+1
дш дп
д iду du\ д iду дш\
(18)
Структура численного алгоритма оставалась неизменной (как и в [9]), На этапе предиктора вязкие члены учитывались в соответствии со схемой стабилизирующей поправки [20],
Некоторые результаты расчетов, соответствующие начальному распределению плотности (5), распределению рs(у) вида (16), ß = 0.1, представлены на рис, 3-6, Расчетная сетка бралась такой же, как и для рис, 2, На рис, 3 изображены линии тока ф = фТь/R2 = const, соответствующие момент у времени t/T = t/2nT* = 8, для идеальной и вязкой жидкости со значениями числа Рейнольдса Re = R2/T*u = 104 и 103, здесь v — кинематический коэффициент вязкости. Линии тока представлены уровнями: 6 —
х
помещенные в начальный момент времени в область смешения. Видно, что в случае идеальной жидкости горизонтальный размер зоны смешения и интенсивность вихря являются максимальными; картины течения для Re-1 = 0 и Re = 104 при t/T = 8 различаются незначительно. Существенное проявление вязкости демонстрируется на рис, 3, в.
В дополнение к рис, 3 на рисунке 4 приведены линии р = const = ps (0.0717), р = p/aRpo. Сплошная линия соответетвует Re-1 = 0, пунктирная — Re = 104, штрих-пунктирная — Re = 103, Результаты, представленные на указанных рисунках, согласуются между собой.
Рис. 3. Сопоставление линий тока в случае идеальной (а) и вязкой жидкости: Re = 104 (б), 103 (е), t/T = 8; маркерами помечены частицы, находившиеся в начальный момент времени в области перемешанной жидкости
Рис. 4. Линии равной плотности р = const = ps (0.0717) в случае идеальной (сплошная линия) и вязкой жидкости: Re = 104 (пунктир), 103 (штрихпунктир); t/T = 8
= в
= h 4
/\ 6 Л 8 10
= ' \ J4
мшим шиш! мшим Iiiiiiiii iiiiiiiii мшим мшим мшим Iiiiiiiii мшим Iiiiiiiii мшим Iiiiiiiii Iiiiiiiii мшим
О 1 2 3 4 5 6 7 8 9 10 11 12 13 14 xR
Рис. 5. Сопоставление линий равной плотности р = const = ps(0.0717) в случае идеальной (о) и вязкой жидкости: Re = 104 (б), 103 (о); t/T = 2, 4, ..., 14
Более подробно волновая картина представлена на рис. 5 в виде изолиний p = ps(yo), yo = yo/R = 0.0717 для различных моментов времени. Кривые 2, 4, 6, 8, 10, 12, 14 соответствуют значениям времени t/T = 2, 4, 6, 8, 10, 12, 14. Изолинии на рис. 5, а-в получены при Re-1 = 0, Re = 104 и 103.
Изменение максимальной амплитуды внутренних волн в зависимости от времени и числа Рейнольдса представлено также в таблице. Здесь приведены обезразмеренные максимальные значения амплитуд внутренних волн для Re-1 = 0, Re = 104 и 103. Видно, что в идеальной жидкости для t/T > 4 максимальная амплитуда внутренних волн практически постоянна. Как и на рис. 4, 5, данные таблицы демонстрируют влияние вязкости на генерируемые внутренние волны.
На рис. 6 приведено изменение во времени абсциссы максимума модуля функции тока и горизонтального размера зоны смешения в зависимости от числа Рейнольдса. При небольших значениях времени графики функций на рис. 6 близки; при t/T > 6 заметно отставание горизонтального размера зоны смешения (рис. 6,6), связанное с неполным перемешиванием жидкости в начальный момент времени. При больших значениях
Изменение максимальной амплитуды внутренних волн в зависимости от времени и числа Рейнольдса
t/T Re"1
0 10"4 10"3
2.0 0.3025 0.3017 0.2987
4.0 0.3239 0.3096 0.2358
6.0 0.3277 0.2977 0.1840
8.0 0.3278 0.2869 0.1439
10.0 0.3279 0.2755 0.1124
20.0 0.3286 0.2210 0.0491
30.0 0.3291 0.1738 0.0478
40.0 0.3298 0.1323 0.0458
50.0 0.3306 0.0989 0.0439
60.0 0.3310 0.0720 0.0420
70.0 0.3287 0.0514 0.0407
Рис. 6. Изменение во времени абсциссы максимума функции тока (а) и горизонтального размера области смешения (б) в случае идеальной (сплошная линия) и вязкой жидкости: Re = 104
времени рост горизонтального размера прекращается; в идеальной жидкости соответствующий основной уединенной волне конвективный вихрь продолжает перемещаться с постоянной скоростью.
Остановимся на контроле точности проведения расчетов при Re = 104. Расчеты проводились в области 0 < x/R < 50, 0 < y/R < 30. Получено, в частности, что для моментов времени t/T = 4, 6 и 8 максимальное значение функции тока фт на сетке с числом ячеек 501 х 301 составило 0.458 х 10-1, 0.434 х 10-1 и 0.409 х 10-1, на сетке с числом ячеек 1001 х 601 — 0.423 х 10-1, 0.414 х 10-1 и 0.391 х 10-1 соответственно. Отклонения массивов сеточных значений функций у, ф, и в равномерной норме не превышают отклонений величин фт.
В случае идеальной жидкости скорость перемещения первого (основного) грсбпя внутренних волн в направлении оси x для t/T > 4 составляла ~ 0.150R/T*. Эта скорость вычислялась как скорость изменения абсциссы максимума модуля функции тока (рис.6). Если воспользоваться асимптотическим соотношением (5.36), (5.37) работы [211
\ а а
\ /// ' 1 у S/ \\ \ V . /У у
\
111111111 111111111 111111111 1 1 1 1 1 1 1 1 1 111111111 111111111 111111111 111111111 111111111 111111111
О 2 4 6 8 10 12 14 16 18 xR
Рис. 7. "Широкий" пикноклин — распределение плотности (16), ß = 0.5; изолинии р = ps(0.477) приведены для момента времени t/T = 8; сплошная линия — Re-1 = 0, пунктир — Re = 104, штрихиунктир — Re = 103
то для максимальной амплитуды воли A = A/R £ [0.32; 0.33] скорость перемещения волны с(1) = c(1)T*/R £ [0.144; 0.145]. Величина с(1) вычисляется для первого (основного) гребня волны, аналогичной изображенной сплошной линией на рис. 6. Сравнение расчетных данных с результатами [21] указывает на их удовлетворительное согласие. В случае вязкой жидкости по аналогии с работой [22] (где рассматривались уединенные внутренние волны на границе раздела двух жидкостей разной плотности) скорость с(1) распространения внутренних волн при найденной в расчетах максимальной амплитуде A = A(t) определялась как непосредственно путем обработки расчетных данных, так и из выражения (19). Получено, в частности, что при Re = 104, 103 и t/T = 10 обработка расчетных данных дает величины с(1), равные 0.138 и 0.0953, а соотношение (19) —
t/T
Наряду с "узким" пикноклином (16), ß = 0.1 рассматривался "широкий" пикноклин (16), ß = 0.5. Расчетные данные в виде изолиний р = ps (0.477) приведены для t/T = 8 на рис. 7. Основное воздействие вязкости проявляется в "поглощении" коротковолновой части линии р = const. Расчеты выполнялись в первом квадранте плоскости xy в прямоугольнике размером 30R х 5R с числом ячеек 301 х 51. На рис. 7 приведена лишь часть расчетной области.
4. Линейные численные модели внутренних волн
Поскольку при изучении волновых движений жидкости широкое распространение получила линейная модель внутренних волн [2, 5, 11, 23], представляет интерес численное моделирование течения, генерируемого локальным возмущением поля плотности в стратифицированной среде для распределений р8(у) вида (15), (16), основанное на линейных уравнениях Эйлера [5]. С использованием независимых переменных ж, у, Ь
и зависимых ф, ш, р1 = р — р5(у) линейная модель волновых движений в стратифицированной среде может быть записана стандартным образом (в обезразмеренном виде):
дрг = &ф(1р1
дЬ дх ¿у ' д2ф д2ф
= о;.
(22)
дх2 ду2
Система уравнений (20)-(22) дополняется начальными и краевыми условиями
ш(0,х,у) = 0, ф(0,х,у) = 0, р1(0,х,у) = р0(х,у) — рДу), (23)
ш(£,х,у) = ф(£,х,у) = р1(£,х,у) = 0, г2 = х2 + у2 — то, £ > 0. (24)
На осях симметрии ставились соответствующие условия симметрии.
Задача (20)-(24) решалась с помощью простейшего конечно-разностного алгоритма типа предиктор-корректор:
, n+1/2 о n pn _ n
Ui,j ~ Uí,j = Pli+ij Pli-ij At/2 2hi
pj - PnUj_ C+u - ( dps
At/2 2hi V dy
(25)
(26)
(A»n+1/2=4n+i/2, (27)
. ,n+1 , ,n pn+1/2 pn+1/2
Ai 2hi
^n+1 -n i n+1/2 i n+1/2 x
At 2h1 V dy
>
(28)
(29)
(AV)J1 = oj (30)
В уравнениях (25)-(30) — конечно-разностный аналог оператора Лапласа, примененный к функции индексы (n +1/2) и (n +1) означают, что соответствующие величины отнесены к моментам времени (n + 1/2)At и (n + 1)At, Применялись конечно-разностные сетки, неравномерные по вертикальной координате. Размеры сеточных
областей при использовании линейной модели аналогичны представленным на рис, 1,
y
чений функции y при x/R = 7 в расчетах вариантов для распределения плотности (15), (16) по модели Обербека—Буссинеска. Предварительно алгоритм тестировался на задаче, имеющей точное аналитическое решение:
-0(í, x, y) = c(t) ■ sin ax ■ sin by,
o(t, x, y) = —c(t) (a2 + b2) ■ sin ax ■ sin by, p1 (t, x, y) = d(t) ■ cos ax ■ sin by,
/2 = а2/(а2 + Ь2), с(*) = (1/1) ёш ф) = (а//2) сое
0 < * < 1, а = Ь =1, 0 < х < п, 0 < у < п.
(31)
Расчеты на последовательности сгущающихся равномерных сеток показали сходимость сеточных решений к точному с порядком 0(Д* + Л2) в норме, являющейся сеточным аналогом нормы пространства непрерывных функций (Д* — величина шага по временной переменной, Н = сог^ — шаг тетки вдоль осей х и у).
Результаты расчетов вариантов линейного распределения плотности невозмущенной жидкости вида (15)
на основе модели (20)-(22) представлены на рис, 8, а-г (а — линейное распределение плотности невозмущенной жидкости (15), б-г — распределения (16), в = 1, 0,5, 0,1),
Сравнение рисунков 1 ,о и 8, о указывает на близость волновых картин в случае линейной стратификации при использовании нелинейной и линейной моделей (это обстоятельство в случае областей полностью перемешанной жидкости отмечалось, в частности, в работе [24]), Отличие наблюдается лишь в окрестности локального возмущения поля плотности. Значительно большее различие в распределении (16) прослеживается с уменьшением /5. При /3 = 0.1 (см, рис, 8, г) линейная модель дает физически неправо-подобное решение: конвективные вихри не перемещаются с течением времени.
Более детально результаты расчетов для У = 0.1 с применением линейной модели представлены на рис, 9, а, б. Вычисления проводились в области 0 < х/Я < 28, 0 < у/Я < 10 Нх = 0.1 Ну = 0.05, Д* = 0.01Т (па рисунке приведена только часть расчетной области). Видно, что с ростом времени центр конвективного вихря остается неподвижным (в противовес расчетам в нелинейной постановке, см, рис, 2,1 и 6), Аналогичные данные получены и на более подробной сетке (рис, 9, в, Нх = 0.05, Ну = 0.025, Д* = 0.0025Т), Поскольку этот результат является нетривиальным, то были проведены расчеты на основе алгоритма типа (25)-(27) (без использования этапа корректора) и получены данные (рис, 9, г), близкие к представленным на рис, 9, б, в. Прокомментируем результаты расчетов на основе анализа закона сохранения энергии.
Из линейных уравнений Эйлера в предположении достаточно быстрого убывания подынтегральных функций при г — то может быть получено соотношение (знак обез-размеривания опущен, П : —то < х < то, —то < у < то)
В случае линейной (или, вообще говоря, невырожденной) стратификации, т, е, при а(у) = — (1/р0)(^р8/^у) > 0, для линейной модели справедлив закон сохранения энергии [25]
рз(у) = ро — У
и нелинейного распределения плотности вида (16)
РДу) = Ро — У • 1апЬ(у/в), /5=1, 0.5, 0.1
(32)
п
Рис. 8. Рассчитанные с применением линейных моделей линии тока для стратификации вида (15), (16); t/T = 2
м м | м м | м м | м м | м м | м м | м м
0 1 2 3 4 5 6 х R
Г Т=5
у,Ж-
I I I I I I I I I I I I I I I I I I I I
2 3 4 5 6 х R
—\ г
/ \
Шу
I I I I I I I I I | I I м | I I I I | I I I
0 1 2 3 4 5 6 х R
Рис. 9. Расчеты по линейной модели распределения плотности невозмущенной жидкости (16), ß = 0.1 hx = 0.10, hy = 0.05, At = 0.01T; a — t/T = 3,5- t/T = 5, в — hx = 0.05, hy = 0.025, At = 0.0025T, t/T = 5, г - схема (25)-(27), t/T = 5
В нелинейной модели, записанной в приближении Обербека—Буееипеека, закон сохранения энергии выглядит следующим образом (а(у) > 0):
д_
m
u2 + v2
+ yp1 1 dxdy = 0.
(34)
Правая часть (32) обусловлена линеаризацией уравнений. Она может приводить к "порождению" энергии в пикпоклипе (в отличие от случая линейной стратификации, где правая часть (32) мала в силу специфики генерируемого локальным возмущением течения). Для анализа факта "порождения" в линейной модели соотношение (32) переписывается в виде
E (t) = E (0) + n(t), n(t) = - у I у у pi vdxdy I dt, E (t) =
0 \ Q ) Q
u2 + v2
+ yp1 I dxdy.
В случае нелинейной стратификации при /5 = 0.1 для моментов времени t/T = 0, 0.5, 2.0, 4.0 величина E(t) принимала соответственно значения 0.202 • 10-1, 0.232 • 10-1, 0.589 • 10-1, 0.132 • 100. Эти значения согласуются со сделанным выше предположением о "порождении" энергии в пикпоклипе.
Наряду с математической моделью (20)-(22) и ее конечно-разностным аналогом
(25)-(30) в качестве линейной модели использовалось также уравнение, являющееся следствием (20)-(22):
= <»>
которое в случае линейной стратификации переходит в известное уравнение Соболева
[26], В качестве начальных и краевого условий ставились следующие, согласованные с (20)—(24):
= ^ = & ' = "• <»>
ф = 0, г2 = х2 + у2 ^ то, г > 0. (37)
Из (36) находим ф, дф/дг для г = 0, Использовалась простейшая конечно-разностная аппроксимация уравнения (35)
(д» £1 -2 (д» I + (Д"ф)з 1 ( Лрз ^ ^ - 2ф» + ф»!
дг2 V ^^ з
Величина (Д^ф)1 определялась из второго соотношения (36):
(38)
Значения (ф)»+1 находились с помощью итерационной схемы стабилизирующей поправки, Результаты численных экспериментов, основанные на сравнении с точным аналитическим решением (31), показали сходимость (на последовательности сеток) сеточных решений к точному с порядком 0(Дг2 + к2), соответствующим порядку аппроксимации (38), Результаты расчетов вариантов для распределений плотности (15), (16) по уравнению (35) представлены на рис, 8, д-з. Они практически идентичны приведенным на рис, 8, а-г для линейной модели (20)-(22).
При использовании линейной модели внутренних волн всегда оговаривается область ее применимости. Приведенные здесь данные расчетов являются лишь иллюстрацией применения линейной модели. Вряд ли следовало ожидать положительного результата ее применения в задаче о динамике локального возмущения поля плотности в пикно-клине.
з
5. Полные уравнения Эйлера
ф
С (п = дф/ду, V = -дф/дх, С = д/ду (рп) — д/дх (рь)) полные уравнения Эйлера принимают вид [27, 28]:
дш дш дш др др д {и2 + др д {и2 + у2\
дЬ 4 дх У ду ^ дх дх ду \ 2 / ду дх \ 2 /
д дф д дф л , дхР дх дуР ду
Начальные и граничные условия для задачи с применением полных уравнений Эйлера аналогичны соответствующим условиям для приближения Обербека—Буссинеска
(4)-(6):
ф = ш = 0, р = р3(у), х2 + у2 ^ то, (43)
р = р0{x,y), р°{х,у) = рs(y), х,у / A, г = 0, (44)
ф = ш = 0, —то < х < то, —то < у < то, Ь = 0. (45)
Обезразмеривание, как и выше, выполняется с использованием масштабов длины Я и времени Т* = 1 /у/ад (а = — (1/ро)с1р3/с1у, у = 0); кроме того, используется представление р = р0аЯр] величина у в уравнении (41) заменяется на д = 1/аЯ, Отметим, что рв(0) = 1/аЯ.
Численное решение указанной системы уравнений осуществлялось на основе эйлеро-лагранжева подхода, аналогичного используемому в разделе 1, Тогда (10) (12) записываются следующим образом:
1+ШИ (46>
dy Зф
dt'
dt
(47)
du 1 дф du g f dy dp dy dp\ 1 /dy dp dy dp\ д iu2 + v
dt' ^ J dr] dt J \drj dt dt dr]
J2 drj
1 J
d_
dt
dy d dn dt
dy/ 2 dn
u + v
J2 \ dn dt 24 dy d dt dn
u
dt dn dn ! + v2
+
2
(48)
дф Ж
d_
dn
d f p dy dy дф\
dt \J dq d^ dr] J
= ш.
(49)
д f p dy dy дф\ dr] \J dr] dt dt J
Конечно-разностный аналог системы уравнений (46)-(49) аналогичен алгоритму, приведенному в [12].
При численном интегрировании уравнения (49) с применением итерационной схемы стабилизирующей поправки экспериментально выбранный итерационный параметр значительно отличается от варианта приближения Обербека—Буссинеска. Последнее обусловлено тем, что в уравнении появляется величина р, которая может сильно изменяться.
Первоначально динамика локального возмущения поля плотности с применением полных уравнений Эйлера рассчитывалась для случая линейной стратификации. Расчеты проводились в области 0 < x < 40Я, — 20Я < y < 20Д. Анализируя результаты вычислений, соответствующие 1/аЯ = 50, можно сделать вывод о значительном различии дальних волновых полей. Уравнения Эйлера в приближении Обербека—Буссинеска приводят к симметричной волновой картине. В случае использования полных уравнений Эйлера при 1/аД = 50 симметрия существенно нарушается. Однако детальный анализ ближнего волнового поля (г < 8Д) показал, что полные уравнения Эйлера, уравнения Эйлера в приближении Обербека—Буссинеска и линейная модель дают близкие результаты. Основное отличие расчетов по линейной модели проявляется лишь в
г/Д < 2
2
2
математических моделей весьма близки. При меньшем значении отношения 1 /aR (например, 1/aR = 5) изменение плотности по вертикали в пределах области интегрирования становится существенным. Данное обстоятельство приводит к большим различиям в расчетах по полной модели и по модели Обербека Буссинеска. Это иллюстрируется рис. 10, где представлены линии тока, полученные по полным уравнениям Эйлера и уравнениям Эйлера в приближении Обербека Буссинеска. Видно, что даже вблизи начата координат имеются значительные различия в картинах линий тока. Величина 1/aR = 5 достаточно мала; в реальном океане 1/aR ~ 103^104. Однако в связи с наличием в океане так называемой тонкой микроструктуры гидрофизических полей [1, 4] могут наблюдаться достаточно произвольные значения 1/aR. При больших вели чинах 1/aR (например 1/aR = 103) различия в линиях тока, как показывают расчеты, практически отсутствуют. Детально было проанализировано распределение (16) при ß = 0.1 (узкий гиперболический тангенс). Расчеты проводились для 1/aR = 103, 102, 101, и 10°. Существенное расхождение с расчетами по модели в приближении Обербека Буссинеска (см. рис. 3) наблюдается при 1/aR = 1. Этот вариант для t/T = 8 представлен на рис. 11. Появляющаяся несимметрия объясняется следующим образом. Приближение Обербека—Буссинеска справедливо в условиях, когда величина плотности р мало отли-
Рис. 10. Линии тока при t/T = 2, полученные по полным уравнениям Эйлера (1/aR = 5)(a) и уравнениям Эйлера в приближении Обербека—Буссинеска (б): линейная стратификация
Рис. 11. Линии тока для уравнений Эйлера, 1/aR = 1, t/T = 8
Рис. 12. Линии тока в момент времени t/T = 2 для уравнений Эйлера (1/aR = 2, ß = 1.0) (а) и для приближения Обербека—Буссинеска {ß = 1.0) {б)
чается от своего характерного значения (например ро). Однако с уменьшением параметра 1/аД изменение плотности уже нельзя считать незначительным. Полные уравнения
Эйлера не допускают симметрии (антисимметрии) решения. В случае сильного измене-р
распределение (16) для ß = 1, 2, 5 и различных значений 1/аД. В качестве примера на рис. 12 сопоставляются линии тока, рассчитанные по полным уравнениям Эйлера н модели Обербека—Буссинеска для ß = 1.0 и 1/аД = 2. Видно существенное отличие в распределениях линий тока, обусловленное значительным изменением плотности по глубине.
6. Наличие волнового фона
С применением уравнений Эйлера в приближении Обербека-Буссинеска и в эйлеро-лаг-ранжевых переменных построена численная модель динамики локального возмущения поля плотности при наличии волнового фона. Фоновое волновое возмущение генерировалось следующими начальными данными для поля плотности [10]:
р(Х,у)=' р0 — у, —2-5 - у - 11
р0 — ^апЬ ((Н (X, у) — 2.0) /0.1) + 1апЬ(9)) — 1.1, 1.1 - у - 3.0,
|у — 3А • сов(рх)
\---2 + А - сов (рх) <у< 3.0,
1 — А • сов(рх)
0.9у+1.1А • сов(рх) „ . ,
——-1.1 <у <2 + А-(Шр.т .
0.9 + А • сов(рх)
Здесь А = 0.15, р = 3п/7. Расчетная область принималась размером 0 — х — 63, —20 — у — 5 Число узлов тетки составляло 631 х 251. Область интегрирования для уменьшения влияния граничных условий выбиралась достаточно большой, поэтому на рис. 13 приведена часть расчетной области (0 — х — 21 —10 — у — 5). Рисунок демонстрирует близкое к линейному взаимодействие волнового фона и внутренних волн. Аналогичные результаты получены и для случая А = 0.15, р = 6п/7.
Рис. 13. Линии тока для момента времени t/T = 2; а — взаимодействие волнового фона и внутренних воли; б— внешний волновой фон; е — внутренние волны, генерируемые локальным возмущением поля плотности в жидкости с нелинейной стратификацией; г — разность полей а и б
Заключение
В работе представлены линейные и нелинейные с применением уравнений Навье—Сток-са и Эйлера в приближении Обербека—Буссинеска численные модели динамики локального возмущения поля плотности в устойчиво стратифицированной среде. Построение нелинейных численных моделей основано на переходе к эйлеро-лагранжевым переменным. Результаты расчетов показали существенное влияние распределения плотности жидкости на генерируемую при коллапсе локального возмущения поля плотности картину внутренних волн. В случае динамики локального возмущения поля плотности в узком пикноклине (идеальная жидкость) продемонстрирован известный факт генерации стационарных внутренних волн солитонного типа. Скорость распространения основной волны и ее амплитуда удовлетворяют асимптотическому соотношению Бенджамина.
Выполнены численные эксперименты по оценке влияния вязкости на генерацию и распространение внутренних волн в пикноклине. Предложена простая алгебраическая зависимость, связывающая максимальную амплитуду волны и скорость ее распространения в вязкой жидкости в "узком пикноклине": при известной из расчетов одной из величин вторая находится из соотношения Бенджамина.
Построена линейная численная модель динамики локального возмущения поля плотности в устойчиво стратифицированной среде. Использовались как линейные уравне-
ния Эйлера, так и эквивалентное им уравнение четвертого порядка для функции тока, переходящее в уравнение Соболева для случая линейной стратификации. Для ряда распределений плотности невозмущенной жидкости даны оценки применимости линейной модели.
Представлены результаты расчетов с применением полных уравнений Эйлера (без использования приближения Буссинеска),
На примере волнового возмущения пикноклина рассмотрена задача о взаимодействии внутренних волн, генерируемых локальным возмущением поля плотности, с фоновыми возмущениями. Для рассмотренного класса возмущений пикноклина получено близкое к линейному взаимодействие.
Авторы благодарят О.Ф, Воропаеву за критические замечания и полезное обсуждение работы.
Список литературы
[1] Федоров К.Н. Тонкая термохалинная структура вод океана. Л.: Гидрометеоиздат, 1976.
[2] лайтхилл Дж. Волны в жидкостях: Пер. с англ. М.: Мир, 1981.
[3] Мадерич B.C., Никишов В.И., Стеценко А.Г. Динамика внутреннего перемешивания в стратифицированной среде. Киев: Наук, думка, 1988.
[4] Монин A.C., Озмидов Р.В. Океанская турбулентность. Л.: Гидрометеоиздат, 1981.
[5] Тернер Дж. Эффекты плавучести в жидкости: Пер. с англ. М.: Мир, 1977.
[6] Васильев О.Ф., Кузнецов Б.Г., Лыткин Ю.М., Черных Г.Г. Развитие области тур-булизованной жидкости в стратифицированной среде // Изв. АН СССР. МЖГ. 1974. № 3. С. 45-52.
[7] Лыткин Ю.М., Черных Г.Г. О внутренних волнах, индуцируемых коллапсом зоны смешения в стратифицированной жидкости // Математические вопросы механики (Динамика сплошной среды): Сб. науч. тр. / АН СССР. Сиб. отд-ние. Институт гидродинамики. Новосибирск, 1975. Вып. 22. С. 116-132.
[8] Гущин В.А. Метод расщепления для задач динамики неоднородной жидкости // Журн. вычисл. математики и мат. физики. 1981. Т. 21, № 4. С. 1003-1017.
[9] Зудин А.Н., Черных Г.Г. Об одном алгоритме расчета нестационарных стратифицированных течений // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Вычисл. центр, Институт теор. и прикл. механики. Новосибирск, 1982. Т. 13, № 5. С. 32-47.
[10] Зудин А.Н., Черных Г.Г. Примеры расчета нестационарных стратифицированных течений с применением эйлерово-лагранжевой системы координат. Новосибирск, 1985 ( Препр. АН СССР. Сиб. отд-ние. Ин-т теор. и прикл. механики. № 9-85).
[11] Степанянц Ю.А., Стурова И.В., Теодорович Э.В. Линейная теория генерации поверхностных и внутренних волн // Итоги науки и техники. МЖГ: Сб. науч. тр. / ВИНИТИ. 1987. Т. 21. С. 93-179.
[12] Зудин А.Н., Черных Г.Г. Внутренние волны, генерируемые локальным возмущением поля плотности в жидкости с нелинейной стратификацией // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Вычисл. центр, Институт теор. и прикл. механики. Новосибирск, 1988. Т. 2, № 4. С. 49-74.
[13] Воропаева О.Ф.,Черных Г.Г. Эволюция зоны турбулентного смешения в жидкости с нелинейной стратификацией // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Вычисл. центр, Институт теор. и прикл. механики. Новосибирск, 1989. Т. 3(20), № 5. С. 3-29.
[14] Chernykh G.G., Voropayeva O.F., Zudin A.N. Numerical simulation of internal waves generated by local density perturbation in stratified medium // Mathematical and Numerical Aspects of Wave Propagation. Proc. of the third Intern. Conf. / Ed. Gary Cohen. Philadelphia, PA: Society for Indust. and Appl. Math., 1995. P. 119-128.
[15] Смирнов Г.В., Веденьков B.E., Галковский A.H. Об условиях генерации нелинейных внутренних волн при коллапсе перемешанного пятна в стратифицированной жидкости // Океанология. 1997. Т. 37, № 3. С. 338-344.
[16] 'Гене/ D.E., Кшо О.М. Nunerical simulations of large-amplitude internal solitary waves // J. Fluid Mech. 1998. Vol. 362. P. 53-82.
[17] Maderich V.S., van Heijst G.J.F., Brandt A. Laboratory experiments on intrusive flows and internal waves in a pvcnocline // Ibid. 2001. Vol. 432. P. 285-311.
[18] Зудин A.H., Черных г.г. Численное моделирование динамики локального возмущения поля плотности в сдвиговом потоке линейно стратифицированной среды // Вычисл. технологии. 2002. Т. 7, № 3. С. 18-28.
[19] Chernykh G.G., Zudin A.N. Linear and nonlinear numerical models of local density perturbation dynamics in a stable stratified medium // Russian J. Numerical Analysis and Math. Modelling. 2005. Vol. 20, No. 6. P. 513-534.
[20] Яненко H.H. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.
[21] Benjamin Т.В. Internal waves of permanent form in fluids of great depth //J. Fluid Mech. 1967. Vol. 29, pt. 3. P. 559-592.
[22] Гаврилов H.B. Уединенные внутренние волны в двухслойной жидкости // Журн. прикл. механики и техн. физики. 1986. № 5. С. 49-54.
[23] городцов В.А., Теодорович Э.В. Линейное описание эволюции (коллапса) симметричных распределений возмущений однородно-стратифицированной жидкости // Методы гидрофизических исследований: Материалы I Всесоюз. школы. Горький, Институт прикл. физики АН СССР, 1984. С. 134-147.
[24] Chernykh G.G., Lytkin Yu.M., Sturova I.V. Numerical simulation of internal waves induced by the collapse of turbulent mixed region in stratified medium // Proc. Intern. Svmp. Refined Modelling Flows. Paris, 1982. P. 671-679.
[25] Гавов С.А., Свешников А.Г. Задачи динамики стратифицированных жидкостей. М: Наука, 1986. 288 с.
[26] Соболев С.Л. Об одной новой задаче математической физики // Изв. АН СССР. Серия математическая. 1954. Т. 18, № 1. С. 3-50.
[27] Гранверг И.Г., Дикий Л.А. Стационирование в нелинейной задаче обтекания гор воздушным потоком // Изв. АН СССР. Серия ФАО. 1972. Т. 8, № 3. С. 264-269.
[28] Мл. [лпк Д.Д. Несколько замечаний об уравнениях гидромеханики и их точных решениях // Вестник ЛГУ. Математика, механика, астрономия. Л., ЛГУ. 1967. Вып. 1. С. 98-105.
Поступила в редакцию 6 марта 2008 г., с доработки — 11 февраля 2010 г.
