Численное моделирование динамики локального возмущения поля плотности в сдвиговом потоке линейно стратифицированной среды Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 7, № 3, 2002

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ЛОКАЛЬНОГО ВОЗМУЩЕНИЯ ПОЛЯ ПЛОТНОСТИ В СДВИГОВОМ ПОТОКЕ ЛИНЕЙНО СТРАТИФИЦИРОВАННОЙ СРЕДЫ *

А. Н. Зудин, Г. Г. Черных Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: [email protected], [email protected]

Calculation algorithm for unsteady stratified flows using Euler — Langrangian coordinate system is presented. Numerical modelling of flow generated by local density perturbation in linear shear flow of linearly stratified fluid is done. Calculation results show the significant distortion of wave pattern by shear flow.

Рассматривается течение, генерируемое локальным возмущением поля плотности в стратифицированной среде. Задача представляет интерес в связи с изучением ряда геофизических явлений. Подробная библиография и анализ работ, связанных с этой задачей, содержатся в [1-6]. Анализ известных работ по численному моделированию течения показывает, что результаты численного моделирования динамики локального возмущения поля плотности в сдвиговом потоке стратифицированной жидкости [7] недостаточно полны. В настоящей работе приведены результаты расчетов, демонстрирующие эволюцию области частично перемешанной жидкости в сдвиговом линейно стратифицированном потоке.

1. Постановка задачи

Будем рассматривать сдвиговой поток вида и = (и(у), 0). Для описания течения привлекаются две математические модели. В качестве Модели 1 используются уравнения Эйлера в приближении Обербека — Буссинеска. С применением переменных функция тока — завихренность уравнения Эйлера могут быть записаны в виде

* Работа выполнена при частичной финансовой поддержке СО РАН, интеграционный проект №2000-1, и Российского фонда фундаментальных исследований, грант №01-01-00783. © А. Н. Зудин, Г. Г. Черных, 2002.

Введение

ди ди ди

ди g др\ р0 дх}

(1)

д 2Ф+=и (2)

дх2 ду2 '

дрл др1 дР1 ¿р3

— + и— + V— + V— = 0. (3)

дЬ дх ду ау

Здесь х, у — декартова система координат; ось у направлена вертикально вверх, против силы тяжести; р1 = р — р3, р3 = р3(у) — плотность невозмущенной среды вне зоны смешения, р0 = р«(0); и^ — горизонтальная и вертикальная составляющие скорости жидкости: и = дф/ду, V = —дф/дх; д — ускорение силы тяжести. Стратификация среды предполагается устойчивой, т.е. (1/р0)а,р3/а,у < 0.

В качестве Модели 2 используются линейные уравнения Эйлера, которые для рассматриваемого сдвигового потока имеют вид

ди1 ттди1 а2и дф1 д др1

—- + и—- =--— + — , (4)

дЬ дх ау2 дх р0 дх

д 2ф1 + д^ф! = ^ (5) дх2 ду2 '

др1 + идр±_ 0Е± дф± = 0 (6)

дЬ дх йу дх '

Величины и1, ф1 — отклонения функции тока и завихренности от соответствующих величин набегающего сдвигового потока.

Граничные и начальные условия в рассматриваемой задаче таковы:

ф = и = ф1 = и1 = р1 = 0, х2 + у2 ^ то, (7)

р = р0(х,у) = рм 1 —

.л „ Ь Ро — Ps(v)

1 - £e-(l-25r)s

Ро

— то < x < то, —то < v < то, t = 0, (8)

ф = Со = ф1 = о1 = 0, —то < x < то, —то < v < то, t = 0. (9)

Здесь е = const = 0.75; ф = ф — фs, Со = о — us; ф3, us — распределения функций тока и завихренности, отвечающие невозмущенному сдвиговому потоку U = U(у). Условие (9) соответствует u(0,x,y) = U(у), v(0,x,y) = 0. Стратификация предполагается линейной: Ps(V) = Ро(1 — ау).

Рассмотрение двух математических моделей — уравнений Эйлера в приближении Обер-бека — Буссинеска (1) - (3) и линейных уравнений Эйлера (4) - (6) является традиционным при изучении динамики локального возмущения поля плотности в устойчиво стратифицированной среде (см. [1, 2, 5, 8] и приведенную там библиографию). В указанных работах анализируется применимость линейной модели внутренних волн к расчету характеристик течения в покоящейся среде. Вопрос об использовании модели (4) - (6) в случае сдвигового течения, насколько это известно авторам, оставался открытым. В связи с этим в настоящей работе предпринята попытка восполнить имеющийся пробел. Достоинством линейной математической модели является возможность получения асимптотических решений в практически важных ситуациях [1-3]. Сопоставление конечно-разностных решений с асимптотическими и анализ применимости последних — весьма непростые задачи дальнейших исследований.

При численном решении задачи системы уравнений (1) - (3), (4) - (6) предварительно приводились к безразмерному виду использованием характерных масштабов длины Я (радиуса области смешения) и времени Т* = 1/^ад (а = — (l/рo)dрs/dy); кроме того, использовалось также представление р\ = р0аЯр\. При анализе результатов расчетов удобно использовать масштаб времени (период Вяйсяля — Брента) Т = 2пТ*.

2. Алгоритмы расчета

Численное интегрирование системы уравнений (1) - (3) осуществляется с помощью подхода, основанного на эйлерово-лагранжевых переменных [4, 8]:

t' = t, X = £ = х, y' = 'q : d'q/dt = 0,

rq(0,x,y) = rq°(x,y), — ж < x < ж, — ж <y < то, rq(t,x,y) ^ q1(y), х2 + y2 ^ ж, t> 0. При этом уравнения (1)-(3) переписываются в виде

др 1 дф др dt' + 1дПд£ = ,

дy дф

(10) (11)

at'

1 I д дф д 1}д£ Ж + dq

+ 1

д

д /дy дф\ д /дy дф

д£\д£дп) дП\д£д£

ди 1 дф ди g (др 1 дрду\

at + = Р0 V— Jana=)

(12) (13) и, (14) (15)

J

д(t,x,y) Оу

В граничных и начальных условиях (7), (8) вместо функций р.(у), р0(х,у) появятся функции Ри(п), р°(£,п), получаемые из р.(у), р0(х,у) подстановкой в них значений У = ,п), являющихся решением уравнения п(0,х,у) = п0(х,у) (п0(х,у) предполагается монотонной функцией у при фиксированном х). К условиям (7) - (9) добавятся начальное

а также граничное

У(0,£,П) = Я2(£ /п)

У = Яз(v), £2 + V2 ^ ж, t> 0 .

(16)

(17)

условия. В качестве координатных линий ц = const в рассматриваемом случае выберем линии р = const. При таком выборе координатных линий уравнение (12) будет выполняться тождественно, так как р = р(ц). В результате перехода к эйлерово-лагранжевой системе координат в уравнении (15) отсутствует конвективный перенос в вертикальном направлении. Это обстоятельство позволяет избавиться от соответствующих ошибок при конечно-разностной аппроксимации.

2

Алгоритм на каждом временном слое сводился к последовательному решению системы

разностных уравнений:

v n+l/2,l+l _ v n

At/2

дф

dU ■■

i,j

n+l/2,l

(18)

n+l/2,l+l n /-, n i

At/2 + V J дп

n+l/2,1 f n+l/2,l+l _ c n+l/2,l+l С - С i- l,j

°i,j''' — ci-l,j'' = _/idpdv\

i,j

hl

n+l/2,l

Ро V Jdnd^J

i,j

1 дф J дп) ij

n+l/2,l

> 0,

n+l/2,l+l n / о i \ n+l/2,l n+l/2,l+l n+l/2,l+l

At/2 +V J дп

(I д1ду\

i,j

1 дф

J дп i,j

hl

n+l/2,l

0,

n+l/2,l

Ро \ Jдпд£^

i,j

(19)

^д n+l/2,l фn+l/2,l+l^

С

n+l/2,l

Vij Vi,j At

i,j

дф

ди ■■

i,j

i,j

n+l/2

(20) (21)

(Jc)n+l — (Ju)ntj At

+

д fдф \n+l/2

V дп

С

i др^У)

\дп ди ij

n+l/2

(22)

(An+^n+\3 = uj.

(23)

В разностных уравнениях (18) - (23) (ду/д£, (ду/дп)г^ — центральноразностные аппроксимации производных; А — конечно-разностный аналог оператора Лапласа в переменных £, п. Индексы (п + 1/2), (п + 1) означают, что соответствующие величины (в операторе А — коэффициенты) отнесены к моментам времени (п + 1/2)АЬ, (п + 1)АЬ; I — номер итерации по нелинейности. Итерации проводились вплоть до выполнения условия сходимости

max

i,j

zn+l/2,l+l zn+l/2,l zi,j zi,j

< ео max

i,j

n+l/2,l

i,j

где ео = const > 0; z — любая из величин у, о, ф.

Решения уравнений (20), (23) находятся с помощью итерационной схемы стабилизирующей поправки.

Алгоритм (18) - (23) на классе гладких функций обладает порядком аппроксимации O(At + hl + h2). При использовании линейной Модели 2 применялись стационарные равномерные и неравномерные сетки и метод предиктор-корректор. Подробности тестирования численных алгоритмов изложены в [4, 5, 8].

3. Результаты расчетов

В связи с тем, что симметрия (антисимметрия) течения при рассматриваемых ниже вариантах параметров сдвигового течения нарушается, решение отыскивалось в прямоугольнике —X < x < X, —Y < у < Y. Величины X, Y варьировались в ходе численных экспериментов. При проведении вычислений по Модели 1 обезразмеренный шаг сетки в горизонтальном направлении полагался равным 0.2; X = 14. Координатные линии П = const по вертикали строились следующим образом [4]. При x = X выбирались значения: y0,... ,yj = yj-1 + hj, j = 1,..., 35; h0 = 0.117R, hj = h0qj, j = 1,..., 35; q = 1.025. Линии n = nj при t = 0 отождествлялись с линиями р = Ps(Уj). В нижней полуплоскости сетка строилась аналогичным образом. При расчетах по Модели 2 значения X, Y существенно варьировались, сетка выбиралась равномерной в обоих направлениях.

Течение, генерируемое локальным возмущением поля плотности в линейно стратифицированной жидкости при наличии невозмущенного горизонтально однородного сдвигового потока вида V = 0, U(y) = ay, а = const, иллюстрируется рис. 1-6.

На рис. 1, а представлены рассчитанные с применением Модели 1 линии уровня дефекта функции тока ф1 = ф — ay2/2 = фТ^/R2 — a(y/R)2/2, полученные для t/2nT* = 3 и а = 0.14 (Ri = 1/T'2(dU/dy)2 = 49). На рис. 1, б приведены линии ф = const, соответствующие тому же значению времени t при а = 0. С применением Модели 2 получен рис. 1, в (Ri = 49). Можно видеть, что рис. 1, а и 1, в достаточно близки и существенно отличаются от рис. 1, б количеством конвективных вихрей, характеризующих генерируемые локальным возмущением внутренние волны.

Эволюция зоны смешения иллюстрируется рис. 2. Расчеты выполнялись по Модели 1. Рисункам справа соответствует коллапс в покоящейся жидкости, слева — сдвиговое течение (а = 0.14, Ri = 49). Динамика зоны смешения прослеживалась путем введения в нее в начальный момент времени маркеров. Следует отметить, что при использовании эйлерово-лагранжевых переменных метод маркеров существенно упрощается, так как маркеры могут перемещаться лишь вдоль координатных линий п = const (р = const). Анализируя рис. 2, можно заметить наличие "шапочек". Их появление обусловлено, по-видимому, особенностями начального распределения плотности в области перемешанной жидкости. Если в качестве области смешения считать круг радиуса г = 0.8 (а не r = 1, как на рис. 2), то эти "шапочки" исчезают. Граница области смешения становится несимметричной в связи с воздействием сдвигового потока; частицы жидкости, соответствующие частицам "шапочек", сносятся потоком со скоростью, близкой к U(у).

Внутренние волны, генерируемые локальным возмущением поля плотности в сдвиговом потоке линейно стратифицированной среды, иллюстрируются также рис. 3, на котором представлена линия равной плотности р = ps(4). На рис. 3, в-3, д демонстрируется зависимость от числа Ричардсона, рис. 3, б соответствует расчетам на вдвое большей в направлении x сеточной области. Можно видеть, что рис. 3, б и 3, в практически совпадают в области x/R£ [—14, 14], y/R£ [3.95, 4.05]. Аналогичная ситуация (близость изолиний) наблюдается и при дальнейшем увеличении области интегрирования как в горизонтальном, так и в вертикальном направлении. Факт существенного взаимодействия сдвигового течения с внутренними волнами, генерируемыми локальным возмущением поля плотности, не зависит от размеров сеточной области. При расчетах по Моделям 1, 2 проводилось также сгущение сетки в два раза как в горизонтальном, так и в вертикальном направлении. Получились результаты, близкие к представленным. В частности, на рис. 3, в штрихами изображена линия р = ps(4), полученная на измельченной в два раза сетке.

Рис. 1. Линии тока: а = 0.14 (а), (в); бессдвиговый случай (б); (а), (б) — Модель 1, (в) — Модель 2.

Рис. 2. Эволюция зоны смешения в сдвиговом потоке, Ш = 49, (слева) и в бессдвиговом случае (справа).

Фазовая картина внутренних волн представлена на рис. 4. В закрашенных областях др/дх < 0, в незакрашенных — др/дх > 0. Границы между закрашенными и незакрашенными областями соответствуют линиям постоянной фазы. С применением Модели 2 получены рис. 3, 4. Модель 1 приводит к близким результатам.

Искажение волновой картины сдвиговым течением допускает достаточно простую физическую интерпретацию. С целью ее иллюстрации исследовались траектории частиц, характеризующие волновые движения. В начальный момент времени положение частицы характеризуется значениями абсциссы х0/Я и ординаты у0/Я. Расчеты проводились с применением Модели 1 в области с размерами X = 28, У = 21 в течение периода времени Ь/2пТ* £ [0, 5]. В процессе расчета брались различные значения х0, у0. Наиболее типичные примеры траекторий представлены на рис. 5, 6: рис. 5 относится к динамике локального возмущения в бессдвиговом случае, рис. 6 — к случаю сдвигового потока (Ш = 49). Можно видеть, что сдвиговой поток существенно меняет траектории. В бессдвиговом случае траектории частиц близки к эллиптическим (с осями переменного размера). В левой полуплоскости их большая ось наклонена под тупым углом к горизонтальной оси координат.

Рис. 3. Линия плотности р = ps(4), t/T = 3 при трех значениях а (а - г), бессдвиговое течение

(д).

Рис. 4. Фазовая картина внутренних волн, др/дх = 0, t/T = 3: бессдвиговый случай (а), Ri 196 (б), Ri = 49 (в).

4.03

у/К

У У

/ ¿У

4.03

у/К

4.02

х/К

V V

л

\

N

х/К

Рис. 5. Траектории частиц в бессдвиговом случае (уо/Л = 4) при Хо/Я = 4 (а) и Хо/Я = —4 (б).

4.01

у/К

4.00

3.99

-4

I I I I I I I I

/ \

I I I I I I I I I

I I I I I I I I I

10 12 14

х/К

Рис. 6. Траектории частиц при наличии сдвига (Ш = 49, уо/Я = 4) при Хо/Я = 4 (а) и Хо/Я = —4 (б).

В правой полуплоскости траектории расположены симметричным образом. Воздействие сдвигового течения иллюстрируется рис. 6. Траектории частиц превращаются в кривые, имеющие ярко выраженный волновой характер, который обусловлен коллапсом области перемешанной жидкости.

Основные результаты работы сводятся к следующему. Выполнено численное модели-

рование динамики локального возмущения поля плотности в линейном сдвиговом потоке линейно стратифицированной среды. Для описания течения привлекаются две математические модели: уравнения Эйлера в приближении Обербека — Буссинеска и линейные уравнения Эйлера. Конечно-разностные алгоритмы основаны на применении переменных ф, и и метода предиктор-корректор. Продемонстрировано существенное искажение картины течения сдвиговым потоком. Анализ результатов численных экспериментов позволяет сделать вывод о возможности применения линейной модели к расчету параметров внутренних волн, генерируемых локальным возмущением поля плотности в линейном сдвиговом потоке.

Список литературы

[1] СТЕПАНЯНЦ Ю.А., Стуровл И. В., ТЕОДОРОВИЧ Э. В. Линейная теория генерации поверхностных и внутренних волн // Итоги науки и техники. ВИНИТИ. Сер. Механика жидкости и газа. 1987. Т. 21. C. 93-179.

[2] МАДЕРИЧ В. С., Никишов В. И., СТЕЦЕНКО А. Г. Динамика внутреннего перемешивания в стратифицированной среде. Киев: Наукова Думка, 1988. 239 с.

[3] HARTMAN R. J. The development of a partially mixed region in a stratified shear flow // J. Fluid Mech. 1975. Vol. 71, Pt 2. P. 407-415.

[4] ЗУДИН А. Н., ЧЕРНЫХ Г. Г. Об одном алгоритме расчета нестационарных стратифицированных течений // Числ. методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1982. Т. 13. №5. C. 32-47.

[5] ЗУДИН А. Н., ЧЕРНЫХ Г. Г. Внутренние волны, генерируемые локальным возмущением поля плотности в жидкости с нелинейной стратификацией // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1988. Т. 2. №4. C. 49-74.

[6] Смирнов Г. В., Веденьков В. Е., ГАлковский А. Н. Об условиях генерации нелинейных внутренних волн при коллапсе перемешанного пятна в стратифицированной жидкости // Океанология. 1997. Т. 37, №3. C. 338-344.

[7] Зудин А. Н., Черных Г. Г. Примеры расчета нестационарных стратифицированных течений с применением эйлерово-лагранжевой системы координат. Новосибирск, 1985 (Препр. / АН СССР. Сиб. отд-ние. ИТПМ; №9-85).

[8] CHERNYKH G.G., VüROPAYEVA O.F., ZüDIN A.N. Numerical simulation of internal waves generated by local density perturbation in stratified medium // Cohen, Gary (ed.). Mathematical and Numerical Aspects of Wave Propagation: Proc. Third Intern. Conf. Mandelieu-La Napoule, France, April 24-28, 1995. Philadelphia, PA: Society for Industrial and Appl. Math. P. 119-128.

Поступила в редакцию 21 января 2002 г., в переработанном виде — 14 февраля 2002 г.