Аналитическое решение одной модели ветрового движения вязкой жидкости Текст научной статьи по специальности «Физика»

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

Том 14, № 4, 2009

Аналитическое решение одной модели ветрового движения

вязкой жидкости

Л. А. КОМПАНИЕЦ, Т. В. ЯкувАйлик, О.С. ПиТАЛЪСКАЯ Учреждение Российской академии наук Институт вычислительного моделирования СО РАН, Красноярск, Россия e-mail: kla@icm.krasn.ru, tanyakub@torins.ru, olga.pitalskaya@gmail.com

Найдено аналитическое решение стационарной модели ветрового движения вязкой однородной жидкости в замкнутом водоеме (случай движения в вертикальной плоскости). Предполагается, что дно бассейна ровное, коэффициент вертикального турбулентного обмена постоянен или изменяется по глубине либо по линейному, либо по экспоненциальному закону. Полученное решение сравнивается с известным решением для модели экмановского типа (без учета члена горизонтальной вязкости), что позволяет определить область применимости этой более простой модели. Полученное решение может быть использовано для тестирования вычислительных алгоритмов, предназначенных для расчета ветрового движения жидкости.

Ключевые слова: вязкая жидкость, движение в вертикальной плоскости, аналитическое решение.

Введение

Уравнения нестационарного трехмерного ветрового движения вязкой однородной жидкости [1, 2] с соответствующими граничными условиями, учитывающими влияние ветра и твердых стенок, часто используются для определения движения в конкретных водоемах. Однако аналитических решений для таких задач в литературе не встречается ни для трехмерного, ни для двумерного течения. Отметим, что аналитические решения для моделей экмановского типа (без учета горизонтальной вязкости) в случае трехмерного или двумерного стационарного течения известны для постоянного и переменного коэффициентов вертикального турбулентного обмена и широко используются для анализа решения [2-7] и расчета конкретных задач. Например, в работе [7] по аналитическому решению для этой модели оценивается эффект возникновения компенсационного противотечения в водохранилищах равнинного типа.

Сделаем ряд упрощающих предположений, позволяющих найти аналитическое решение хорошо известной системы уравнений, выписанной в [1, 2], для горизонтальной составляющей скорости как при постоянном, так и при переменном коэффициентах вертикального турбулентного обмена.

Предположение 1. Движение является однородным и двумерным в вертикальной плоскости, коэффициент горизонтального турбулентного обмена — величина постоянная.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 07-01-00153).

© ИВТ СО РАН, 2009.

Предположение 2. В уравнении движения и в кинематическом краевом условии на свободной поверхности можно пренебречь нелинейными членами.

Предположение 3. Отклонение свободной поверхности от невозмущенного положения мало, и влияние ветра можно рассматривать на невозмущенной поверхности бассейна [2].

Предположение 4- Дно ровное; на дне задано условие, линейно связывающее касательное напряжение и горизонтальную скорость течения.

Предположение 5. На вертикальных стенках бассейна ставятся условия прилипания.

Предположение 6. Существует стационарное решение уравнений с указанными граничными условиями.

Тогда получаем следующую краевую задачу, которая решается в области 0 < х < Ь, -Н < г < ((х):

д( к д2и д iк ди

g дх x дх2 dz \ z dz

ди д'ш о

дх дz '

с граничными условиями

где в(х) = т/ро,

Kz -

z дz

ди

= в (х),

z=0

Kz дz

kbu\z=-

H

z=-H

U\x=0 = U\x=L = 0 W\z=0 = W\z=-H = 0

(1) (2)

(3)

(4)

(5)

(6)

Здесь H > 0 — глубина бассейна; и = и(х,z),w = w^,z) — соответственно горизонтальная и вертикальная компоненты вектора скорости течения; g — ускорение свободного падения; Kx > 0 — коэффициент горизонтального турбулентного обмена; Kz (z) > 0 — коэффициент вертикального турбулентного обмена; ( = ((х) — отклонение поверхности жидкости от равновесного положения; т — касательное напряжение ветра на водной поверхности; р0 — плотность воды; kb = const — коэффициент придонного касательного напряжения. Отметим, что вариант kb = 0 отвечает условию скольжения без трения, а kb = то — условию прилипания. Ось z направлена вертикально вверх.

1. Точное решение при постоянном коэффициенте вертикального турбулентного обмена (Kz = const)

д(

Найдем решение системы (1)-(2) при условиях (3)-(6). Исключим величину —, для чего проинтегрируем уравнения (1)-(2) по z £ [-H, 0]:

НК

дх

Kx

дх2

ts- ди udz + Kz — дz

-H

(7)

-H

0

0

о

н = —— ивг. |-н дх У -н

С учетом граничных условий (5), (6) при любых х £ (0,Ь) справедливо:

о

У ивг = 0.

-н

Таким образом, соотношение (7) преобразуется к виду

гг дС ту- ди °

дх дг -н

а уравнение (1) —

д2и д2и в — кьи1х=-н

Кх дХ2 + К д^ =-Н-• (8)

Отметим, что для модели экмановского типа (когда членом с горизонтальной вязкостью можно пренебречь) уравнение для горизонтальной составляющей скорости можно записать следующим образом:

д 2и в — кьи1х=-н

К* а* =-Н-• (9)

Решение уравнения (8) ищем в виде

^ те

/ Ч ✓ N • [пих\ Г . пш. * ^ * 1 . /пих\ . .

и(х, г) = ип(г) яп ^—) =1^ [Апв ь * + Бпв- ь * + Бп яп ^—) • (10)

п=1 п=1

С учетом граничных условий (3), (4) Ап, Бп и Оп определяются из решения следующей системы уравнений:

К-.пи

^ ~^(Ап — Бп) = вп,

—z ПП ( . —nn^H f —nn^H nn^H \

[Ane L - Bne L J = kb \^Ane L + Bne L + DnJ

/ППЛ 2 ( —nnpH nn jj.H \

-H—x{—J Dn = pn - kb i^Ane L + Bne L + DnJ ,

где ^ =л —X, а Pn — коэффициенты разложения функции — в ряд по синусам:

V —z Ро

те

Т

Р0 х>ч^). Д.=!/fs^d, (1Ц

Если ветер постоянен вдоль всего бассейна, то

в (—1 + cos (ПП))

Pn = -2.0

ПП

L

и ненулевые вп для больших значений п имеют порядок О — . Докажем в этом случае

\п)

сходимость ряда (10).

г

Коэффициенты пп(г) имеют вид -, где

д

г = -вп Ь (-к, , кь Ь еЬ (^) + К, , кь Ь еЬ (Ь + Н-

—п2 n2Kz ц H Kx ch

ппц (z + H)

L

+ L^2KZ 2nn sh

пnцH

L

— kb LH Kx nn sh пnцH

ппц (z + H) ^ +

q = пnцKz [ Kz ц kb L2 sh

+ Kz ц kb L2 ch пnцH

L

L

— Kz ц kh L

L

— n2n2 Kz ц H Kx sh

+ kb LH Kx nn ch

пnцH

L

пnцH\

L

+

i H

Разделим г и д на е ь п2п2. Так как -Н < г < 0, то все члены, содержащие

п п(г-Н)^ п пН^ п п(г+2Н)^ п п2Н^ п пг^

е ь , е ь , е ь , е ь , стремятся к нулю при п ^ те, а величины е ь ,

п п(г+Н)д

е ь для данного диапазона изменения переменной г остаются ограниченными, и с учетом оценки для вп величины пп убывают с возрастанием п не медленнее, чем 1/п2, что гарантирует абсолютную сходимость ряда (10).

В двумерном случае можно сравнить известное стационарное течение экмановского типа и стационарное течение для модели с учетом горизонтальной вязкости.

Особенность течения с учетом горизонтальной вязкости в том, что решение обращается в нуль на концах отрезка х £ [0, Ь], в то время как решение уравнения (9) удовлетворяет на границах условию обращения в нуль полного потока. На рис. 1 и 2 сплошной линией изображено решение (10), а ромбиками — решение для модели Эк-мана, которое не зависит от точки х. Рассматривался водоем длиной 1000 м, при этом

Рис. 1. Распределение скоростей при Kz = const, Kx = 20 м2/с, L = 1000 м

Рис. 2. Распределение скоростей при Kz = const, Kx = 200 м2/с, L = 1000 м

т/р0 = 10-4 м2/е2, коэффициент придонного трения кь = 2 м/с, коэффициент вертикального турбулентного обмена Кг = 0.02 м2/с, глубина Н = 20 м. Расчеты проводились для различного числа членов ряда (10), и было найдено то значение п0, начиная с которого решения уже практически не отличаются друг от друга (в данном случае по = 40).

Видно, что в области, достаточно удаленной от берегов, решения будут близки. При увеличении коэффициента горизонтального турбулентного обмена Кх различие между решением (10) модели с учетом горизонтальной вязкости и решением для модели экма-новского типа становятся все более существенными. Например, на расстоянии 20 м от берега при Кх = 20 м2/с решение (10) в 5 раз меньше, чем экмановское решение, при Кх = 200 м2/с решения различаются в 10 раз.

2. Точное решение при линейном распределении Кг по глубине

Известно, что распределение скорости по глубине существенно зависит от выбора Кг. Найдем решение для модели с учетом горизонтальной вязкости при Кг = аг + Ь. В этом случае уравнение (8) примет вид

д2и д ( дп\ Р- кьп\г=-и

КхдХ2 + д~г КаХ + Ь) д~г) =-Н-■ (12)

Будем искать решение в виде ряда

то

nn

u(x,z)^Фn(z)sin — х. (13)

Ь

п=1

Подставим это выражение в уравнение (12):

д2фпг ,1Л , дФп (пп)2 рп - кЬфп\,=-и

(аг + Ь) + - МТ) Ф(г) =-Н-■

Сделаем замену переменных:

£ = аг + Ь, Фп(г) = Хп(£).

Получим уравнение

>2 д2Хп , 2 дХп _ (пи\2 вп — кЪХп\^=-аН+Ъ п

*а~д? + — КЛт) Хп =-Н-• (14)

Поделив обе части в (14) на а2 и умножив на *, получим

*2д2 Хп + *дХп *К (пи\2 х = *вп — къХп\^=-ан+Ъ (15)

* Ж+— *КЛТа) Хп = * на2 • (15)

Это уравнение представляет собой неоднородное модифицированное уравнение Бесселя с граничными условиями:

дХ

п д*

= вп, (16)

И=ъ

дХп д*

= кЪ Хп\^=Ъ-аН • (17)

$=Ъ-аН

С учетом этого обстоятельства решение имеет следующий вид [8]:

Фп(*) = Ап1° (^^КЩО) + БпК° + Пп, (18)

где ), К°(£) — обобщенные функции Бесселя. Коэффициенты Ап, Бп, Оп находятся из уравнения (15) и граничных условий (16), (17) и зависят от и и функций Бесселя К1, К°, 11, 1°. Для больших значений аргумента эти функции ведут себя подобно экспоненциальным функциям соответствующего вещественного отрицательного (К1 и Ко) и вещественного положительного (11 и 1°) аргументов [9] и доказательство сходимости ряда (13) вполне аналогично доказательству сходимости ряда (10).

На рис. 3 сплошной линией изображено решение (18) при К* = аг + Ь, ромбиками — решение для модели экмановского типа. Решение найдено в водоеме длиной 1000 м, при этом т/р° = 10-4 м2/с2, коэффициент придонного трения кЪ = 2 м/с, для коэффициента вертикального турбулентного обмена Ь = 0.02, а = 0.0009, глубина Н = 20 м. Рассмотрен случай, когда и° = 40.

Рис. 3. Распределение скоростей при а = 0.0009, Кх = 200 м2/с, Ь = 1000 м

При сравнении ветрового движения жидкости с постоянным коэффициентом турбулентного обмена Kz = const (см. рис. 2) и в случае его линейного изменения по глубине Kz = az + b (см. рис. 3), при прочих одинаковых параметрах, скорости движения жидкости на поверхности и ближе ко дну различаются почти в 2 раза. Для модели с учетом горизонтальной вязкости скорость движения жидкости на поверхности и ближе ко дну увеличилась примерно в 1.5 раза. Также изменилось положение точки, в которой скорость движения становится нулевой — она сместилась вниз как для модели экмановского типа, так и для модели с учетом горизонтальной вязкости. Достичь такого эффекта при Kz = const невозможно. Легко показать, что для модели экмановского типа при Kz = const эта точка постоянна (например, при условии прилипания на дне, Ho = -H/3).

3. Аналитическое решение при экспоненциальном изменении Кг по глубине

Для модели с учетом горизонтальной вязкости при Кг = deXz получим уравнение

К, ^ + dXe^z ^ + = в - къ<=-н. (19)

дх2 дг дг2 Н

т

Ищем решение в виде (13), — представляем как разложение в ряд (11). Подставим

Ро

выражения (13), (11) в (19):

-К, (Т)2 Ф„ + ^ + ^ = вп - .

Сделаем замену переменной £ = е- *2:, после чего приходим к уравнению ^ (ПП )2Ж ЛД2 1 д Фп , ,А2 д2Фп Рп - къфп\ е * н

х\т) фп- ^{ИТ + ^=-Н--(20)

с соответствующим образом преобразованными граничными условиями:

d\ дФ

d\ дФп

2С дС

С=1 2С сдС

п

х = кь Фп1 4н ■ (21)

=4H

d\2

Умножим (20) на £2 и разделим на -4-:

К (™ )2 4 ё2Ф ёд Фп + ё2 д 2Фп 4 . 2 вп - къФп1 = *н

-КЛт) Ж2ё Фп - +ё ^ = Ж2ё-н-• (22)

Это уравнение представляет собой неоднородное модифицированное уравнение Бесселя, решение которого выписывается в виде

Фп (С) = Anthl ^JKЧ + BntKi( JjVKiС I + Dn. (23)

Рис. 4. Распределение скоростей при Кх = 200 м2/ с, d = 0.02, Л = 0.05, Ь = 1000 м

Рис. 5. Распределение скоростей при Кх = 200 м2/ с, d = 0.01, Л = 0.01, Ь = 1000 м

Рис. 6. Сравнение экспоненциального распределения (штриховая линия) и линейного распределения (непрерывная линия) К* по глубине

Подставляя данное решение (23) в уравнение (22) и граничные условия (21), находим коэффициенты Ап, Вп и Оп. Сходимость ряда (13) с Фп(£), найденными по формулам (23), доказывается аналогично рассмотренным ранее случаям.

На рис. 4 и 5 сплошной линией изображено решение (23) при К* = ¿ех*, ромбиками — решение для модели экмановского типа. Рассматривался водоем длиной 1000 м, при этом т/р0 = 10-4 м2/с2, коэффициент придонного трения кь = 2 м/с, глубина Н = 20 м.

При сравнении решения, соответствующего экспоненциальному распределению К* (см. рис. 4), с решением, соответствующим линейному распределению К* (см. рис. 3), видно, что на поверхности их максимальные скорости совпадают, а ближе ко дну — различаются. Это можно объяснить тем, что при данных параметрах значение коэффициента турбулентного обмена К* = ¿ех* на поверхности практически совпадает со значением линейного коэффициента турбулентного обмена К* = ах + Ь, а в нижней части они уже существенно различны (рис. 6).

4. Определение области применения модели Экмана в случае безразмерных переменных (Kz = const)

Введение безразмерных переменных помогает оценить количественные характеристики применимости моделей экмановского типа. Проведем обезразмеривание уравнений (8) и (9) следующим образом: x* = x/L, z* = z/H, w* = w/u0, u* = u/u0, (* = gC/Mo, K* = Kx/u0L, K* = Kz/u0H, (т/ро)* = (т/р0) /u°, k* = kb/uo, где u0 — характерная скорость задачи.

Тогда краевая задача для уравнений с вязкостью запишется следующим образом:

H В2ь* d2u*

H * д u , TS* д u d* 1* *l

K

LKxdxx*o + k dz*o = в - kbu iz*=-i du

*

*

z dz*

= в *(x), в *(x) = (т/ро)*

z* =0

du*

K

*

dz *

_ 7 * * I

= kbu |z* = -1

z* = -1

u* |x*=0 = u*|x* = 1 = 0

w* |z*=0 = w*|z*=-1 = 0

lz*=0

Уравнение экмановского типа запишется в виде

д { 8у*\

а? (К'а?) = в'- к«

Один безразмерный параметр Н/Ь определяет область применимости модели экма-новского типа.

Данная задача возникла при анализе ветрового течения в оз. Шира, потому в этом разделе выбор основных параметров обусловлен геометрией оз. Шира и масштабами происходящих в нем процессов.

Рис. 7. Оценка области применимости модели Рис. 8. Поведение решения вблизи границы экмановского типа

При этом бралось в = 0.0001 м2/с2, что приблизительно соответствует скорости ветра 7 м/с и и0 = 0.02 м/с. Коэффициенты Kz = 0.02 м2/с, къ = 2.0 м/с и Кх = 10.0 м2/с взяты из литературы [10] (отметим, однако, что эти величины в различных источниках различаются на порядок). Это соответствует безразмерным коэффициентам К* = 0.05, к* = 100.0, в* = 0.25. На рис. 7 и 8 представлены результаты расчетов в безразмерных переменных. На рис. 7 приведены зависимости скорости от глубины в середине бассейна единичной длины. Штриховая линия — модель с вязкостью при Н/Т = 0.05 (К, = 1.25); точечная линия (крестики) — при Н/Т = 0.03 (К, = 0.75); сплошная линия соответствует Н/Т = 0.02 (К, = 0.5) и совпадает с решением по модели экмановского типа. Уже при Н/Т < 0.05 различие между решениями модели с вязкостью и модели экмановского типа составит не более 10%. На рис. 8 исследовано поведение решения при приближении к границе области единичной длины. Здесь рассмотрен случай Н/Т = 0.004 (К, = 0.1). Сплошной линией обозначено решение модели с вязкостью в центре области, которое совпадает с расчетами по модели экмановского типа. Штриховой линией обозначено решение модели с вязкостью на расстоянии 0.06 от границы, штрихпунктирной линией — на расстоянии 0.02 от границы. Видно, что на расстояниях, меньших 0.02 от границы, различие решений составит более 15%.

Заключение

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

Изменение коэффициента турбулентного обмена по глубине согласно линейному закону Kz = az + b в сторону уменьшения (с увеличением глубины) приводит к возрастанию скорости движения жидкости на поверхности и к изменению положения линии, на которой скорость движения становится нулевой, по сравнению со случаем Kz = const = b. Такие же выводы можно сделать для случая экспоненциального распределения Kz = deXz по глубине.

Авторы выражают глубокую благодарность проф. В.М. Белолипецкому за содержательное обсуждение полученных результатов.

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

[1] Саркисян А.С. Моделирование динамики океана. СПб.: Гидрометеоиздат, 1991. 296 с.

[2] КочЕргин В.П. Теория и методы расчета океанических течений. М.: Наука, 1978. 128 с.

[3] Квон В.И. Гидротермический расчет водоемов — охладителей // Изв. АН СССР. Энергетика и транспорт. 1979. № 5. С. 129-137.

[4] Волкова Г.Б., Квон В.И., Филатова Т.Н. Численное моделирование ветровых течений в Чудском озере // Водные ресурсы. 1981. № 3. С. 91-99.

[5] Компаниец Л.А., ЯкувАйлик Т.В. Аналитическое решение одной модели ветрового движения жидкости // Вычисл. технологии. 2003. T. 8, № 5. C. 78-83.

[6] Модели экмановского типа в задачах гидродинамики. / Л.А. Компаниец, Т.В. Якубай-лик, Л.В. Гаврилова, К.Ю. Гуревич. Новосибирск: Наука, 2007.

[7] Зырянов В.Н., Фролов А.П. Природные компенсационные противотечения в водохранилищах равнинного типа // Водные ресурсы. 2006. T. 33, № 1. C. 5-13.

[8] Камке Э. Справочник по обыкновенным дифференциальным уравнениям. М.: Наука, 1972.

[9] Коренев Б.Г. Введение в теорию бесселевых функций. М.: Наука, 1971.

[10] Wang Y., Kolumban H. Methods of substructuring in lake circulation dynamics // Advances in Water Resources. 2000. N 23. P. 399-425.

Поступила в редакцию 16 декабря 2008 г., в переработанном виде — 1 апреля 2009 г.