Вычислительные технологии
Том 1, № 3, 1996
О ЧИСЛЕННЫХ АЛГОРИТМАХ ДЛЯ НЕЛИНЕЙНО-ДИСПЕРСИОННЫХ МОДЕЛЕЙ МЕЛКОЙ ВОДЫ В ДВУМЕРНОМ
СЛУЧАЕ*
Л. А. Компаниец Вычислительный центр СО РАН, Красноярск, Россия
Рассматриваются разностные схемы для двумерных вариантов нелинейно-дисперсионных моделей мелкой воды. Анализируются диссипативные и дисперсионные свойства разностных схем, приводятся результаты численных расчетов.
1. Описание нелинейно-дисперсионных моделей в двумерном случае
Следуя [1], проанализируем дисперсионные соотношения нелинейно-дисперсионных моделей (н.-д.м.) в двумерном случае. Этот анализ позволяет выделить по крайней мере три класса.
Один класс составляют н.-д.м. Грина — Нагди [2], Пелиновского — Железняка [4], Баз-денкова — Морозова — Погуцце [5], первая модель Перегрина [6], третья модель Дорфма-на — Яговдика [7], имеющие дисперсионное соотношение
2 = Нд(К2 + К2)
" 1 + ¥■ К + к2)'
в котором частота есть вещественная функция волновых чисел К1? К2 и возможно построение устойчивых разностных алгоритмов.
Уравнения модели Грина — Нагди [2, 3], имеют вид
к + У(кУ) = 0,
БУ + дУп = -1/6(-Б2НУ(2п - Н) + Б2пУ(4п + Н) + кУ(2Б2п - Б2Н)),
где х, у — пространственные переменные, У = (и,п) — вектор скорости, п — возвышение свободной поверхности, к — полная глубина, к = п + Н(х,у,£), Н — глубина бассейна, Б = е/т + УУ, У = (д/дх, д/ду).
*© Л. А. Компаниец, 1996.
При построении устойчивых численных алгоритмов необходимо встречающиеся в уравнении движения производные п*, П*х, П*У и т. д. заменить на производные от и, V по пространственным переменным [3] (Н = H(ж, у))
П + V(frV) = 0,
V + (^)У + ^п = 1^[2п - Н][(^)Н]4+ +1^[2п - - 1^[4п + Н][^(йУ) + ^)п]*-
-1^[4п + ^^Н^^) + ^)п] - + ^)п]*-
+ ^)п))] + 1/6^^)^)^ + 1/6^^)^*. (1)
В численных расчетах применяется следующая форма записи:
П + V(hV) = 0,
и - 1/6д/дж[2п - ^^^ + 1/6д/дж(4п + H) + ^)п] +
+1/3йд/дж[^(^) + ^п)] - ^йд/д^^^) = В1, V - 1/6д/ду[2п - ^^^ + 1/6д/ду(4п + ^[^(^ + ^)п]+ +1/3йд/ду[^(^) + (VVn)] - 1/6йд/ду((УV)H) = В2,
В* = Ф(n,u,v,H), (2)
где В = (В 1, В2),
Ф(п,и, V, H) = 1^[2п - ^С^)^^ - 1^[4п + ^^Н^йУН
+^)п] - )) + (V+
+2/3(^(^) + ^)пМ^(^)) + +^)п)(^(^))) - 1/6(У^- ^п -
Несложные выкладки показывают, что уравнения (1) совпадают с уравнениями модели Пелиновского — Железняка
П* + V(hV) = 0, V* + (У V)V + ^п = Е,
где
1 й3 й2 й е = I К+т - vн (2 д+д),
К = + (V^ - (^)2,
д = vívн + (V V) (V vн) (3)
и модели Базденкова — Морозова — Погуцце
П* + V(hV) = 0,
1 9 У + (УУ) V + дУп = т У А - ТУН,
кк
где
к1 А = к2Б[3 УУ +2(У УН)],
к
9 = УУ + (У УН)],
если в последней сделать такую же замену производных п, , , как при преобразовании модели Грина — Нагди.
Для построения численных алгоритмов используется как форма (2), так и форма (3). Модель Пелиновского — Железняка может быть записана в обезразмеренном виде с явным выделением малых параметров нелинейности а и дисперсии в [8]. Так же, как в одномерном случае [1], полагая а = О (в) и отбрасывая в модели Пелиновского — Железняка члены порядка О(ав), получим н.-д.м.
п + У(кУ) = о,
У + (УУ)у + дУп = ЙНу(УНУ) - 1 Н2У(УУ)] (4)
2 6 4
(первую модель Перегрина), оставляя члены О (а), получим модель нелинейной мелкой воды
п + У(кУ) = 0, У + (УУ)У + дУп = 0, (5)
и оставляя только члены порядка О(1), получим модель линейной мелкой воды
п + У(НУ) = о, V + дУп = 0. (6)
Третья модель Дорфмана — Яговдика
д
У + У (У У) + дУп + 1/2У(Н крй) = У(1/3Н2—УУ + 1/2Н УН/ У), кг + У[(Н + п - кр)У] + 1/2У(НУкр) = У[3Я (Укр)2У + Я2УНУУ]
позволяет учесть зависимость изменения положения дна от времени. Здесь функция, задающая дно Н(х,у, ¿), представлена как разность функции Н"(х,у), не зависящей от времени, и функции кр(х, у, ¿), учитывающей зависимость этой функции от времени Н(х, у, ¿) = Я(х,у) - кр(х,у,¿).
Второй класс моделей составляют вторая модель Перегрина
У + У (У У) + дУп = 0,
кг + У(кУ) + 1/2У(-1/3Н3 У(УУ) + Н2У(УНУ)) = 0,
Н = Н (х, у)
и ее обобщение на случай, когда глубина бассейна зависит от времени — вторая модель Дорфмана — Яговдика.
Они имеют дисперсионное соотношение с частотой, которая является мнимой функцией волновых чисел, и построение устойчивых разностных схем невозможно. В третий класс входят первая модель Дорфмана — Яговдика
V* + V + ^п + V(H йр*) = 1^2(Я2^),
Ь + V[V(н + п - йр)] + 1^2(Я2йр*) = 1^3(#У)
с теми же обозначениями, что и в третьей модели Дорфмана — Яговдика H(ж,у, ¿) = //"(ж, у) - йр(ж,у, ¿) и Алешкова [9]
Ь + V(hV) = V[hVH ^н) + (VH (VV) + V(VVH))h2/2 + V(VV)h3/6],
Ví+V(g(h-H )+1 /2| V |2) = V[1/2(VVH )2+^н+VV(VVH )]h+[VVí+VVV-(VV )2]й2/2],
H = H (ж, у)
с дисперсионным соотношением
1 + (К2 + к2)
= ^(К2 + к2 )-
1 + ^ (К2 + к2)'
в котором частота является вещественной функцией волновых чисел.
При численных расчетах используется запись модели Алешкова в виде
ь + v(hv) = д = (д1,д2),
д1 = д/дж[Л^н (VVH) + ^н (VV) + V(VVH ))й2/2 + V(VV )й3/6],
д2 = д/ду[Л^н (VVH) + (VH (VV) + V(VVH ))й2/2 + V(VV )й3/6],
V - V[(hVVH) + 1/2h2VV] = VC, С* = Ф(п,и^,н),
Ф(п, и, V, H) = -й^н - йй^ - - H) - 1 /2| V|2+
+1/2(V VH )2 + V(VVH) + h2(VVVV + ^ )2). (7)
Отметим, что модели Грина — Нагди и Базденкова — Морозова — Погуцце выведены без предположения о потенциальности течения, остальные из перечисленных выше — при этом предположении.
2. Описание численных алгоритмов в двумерном случае
Выпишем формулы двупараметрического семейства схем для нелинейно-дисперсионной модели Грина — Нагди в форме (2), аналогичного рассмотренному в [1] для одномерного случая:
+ Дх(Ь? ■ <,■) + Ду (й? V™.) = 0,
п
Вп+1 - В'
^ Д ^ = Ф(Пс«,ип,.,<, ,нп),
ипп+1 - 1/6Дх[2п„- - Д™ ](иП,+^хЯ™ + ^Ду н?) + 1/6Дх(4п„- + н? )х
хьдлк;1) - Ду (к;1)+«^д*^,+^дА, } + 1/3^, [-для (К+1)--джАу (к;1)+дл«™;^™,+«щ^Ду п™,)] - 1Мг,, Ах«;^:+«^Ду я™) = в V™;1 - 1/6Ду[2%, _ я™.](иГ++1ДхЯ™. + ^™+1Дуя™.) + 1/6Ду(4пг>, + ЯГ,)х х^д^к;1) - Ду (Ц+1) + <+1 ДхП, + ^™+1Ду } + 1/зНг,, [-ДхДу (К+1)--Дуу (к;1)+Ду (и, 1Дхп!:,+г^Ду п™,)] - 1/б^г,, Ду К+^Я:.+^Ду я™,.) = В2™;1,
ПП+1 - п»™.
" д ^ + А*^,<,) + Дуг™,) = 0, (8)
Пег, = ^Пад + (1 - , Нсг, = Пег, + Я*,, Н», = Пг, + Я*,,
«С»,, = к;1 + (1 - ¿к,, V«, = ¿г:;1 + (1 - ¿к™,,
/га _ /га /га _ /га
д /™ = ^г;1', ^»—1', д /™ — ^г',;1 1
г,, = ол^ ' ду
где
2Ах ^ 2Ау
/ га _ 2 + / ™ _ 2 / П + / ™
д <-га = ^ 2/г,' + ^ д = 2fг'j ' + ^г,— 1
Дхх /г,, = дх2 ' Дуу /г,, = Ду2 '
Отбрасывая здесь "лишние"члены, получим разностные схемы для н.-д. модели Перегрина, нелинейной и линейной мелкой воды.
Двупараметрическое семейство для н.-д. модели Грина — Нагди при предположении потенциальности течения отличается от (8) тем, что в них условие иу = гх позволяет при вычислении и™;1^™;1 заменить в дисперсионных членах операторы
Д ?/™+1 + Д г™;1 Д Д ?/™+1 + д г™;1
ДХХ + джу гг, , ду дхиг, + дуу гг,,
на
Д ?/™+1 + Д ?/™+1 Д г™;1 + Д г™;1 (9)
+ дуу , джжгг,, + дуу гг,, (9)
соответственно.
Опуская "лишние" члены, получаем разностный алгоритм типа (9) для модели Перегрина, но для него замена иу = гх в дисперсионных членах справедлива только для Я (х, у) = сопэ^
Двупараметрическое семейство разностных схем для модели Алешкова (7) имеет вид ^ Д^' + А*^- и?,-) + Ду(Н™,г™,) = Дх^(Н™,' и?,-, г™,, Яг,,) + Ду^(Н™,, и™,, г™,, Яг ,,),
с ™+1_с™
^ Д ^ = Ф(Псг,, , гг™,, Яг,),
<+1 - ДхНг,, «+ 1ДхЯг,, + г, Ду Яг,,) - 1/2(Д*<+1 + Ду гт;1 )Дх (Нг,, )2-
_1/2(Нг,, НАх*«™;1+д«у<+1) = А*сг;1'
гг™;1 - ДуНг,,«+1АхЯг,, + г-^ДуЯг,,) - ^(Дхи™;1 + Дугг:;1)Ау(Нг,,)2--1/2(Нг,, )2(Дххгг:;1 + Дуу гг:;1) = Ду ,
nnJ+1 - пГ
д t + Ax(hcij uci,j ) + Ay (hci,j vci,j ) AicQ1(hci,j ) uci,j ) vci,j ) Hi,j ) +
+AyQ2(hci,j, Mcij, vci;j, Hij). (10)
При использовании условия (9) в схеме (8) и всегда в модели Алешкова при нахождении скоростей на следующем шаге по времени решается эллиптическая система уравнений и можно использовать известные методы их решения.
Четкое выделение в моделях Грина — Нагди и Алешкова гиперболической и эллиптической частей позволяет применять для гиперболической части различные известные аппроксимации. Так, разностные алгоритмы типа leap-frog для моделей Грина — Нагди и Алешкова получаются заменой в (8) и (10) разностных аппроксимаций для производных по времени на центральные разности, а алгоритм типа предиктора-корректора для модели Пелиновского — Железняка [10] легко обобщается на двумерный случай через схему вращения, предложенную в [11] для двумерного уравнения переноса.
Определяя значения пГ;, Н— в целых точках разностной сетки, а u™+1/2j, v™j+1/2 в полуцелых, получим разностную схему с разнесенными значениями возвышения свободной поверхности и скорости.
3. Анализ диссипативных и дисперсионных свойств численных алгоритмов в двумерном случае
В таблице выписаны собственные значения Pi(£, ф) матриц перехода линейных аналогов разностных схем для н.-д. моделей первого и третьего классов [12]. Здесь £ = K1 Ax, ф = K2Ax, к1 = At/Ax, к2 = At/Ay. Анализ диссипативных свойств показывает, что для разностных схем leap-frog условия устойчивости схемы становятся более ограничительными, чем в (8) при условии (9) для н.-д.м. первого и в (10) для н.-д.м. третьего класса. Если 8 + и = 1, то |pi| = 1 для любого i и всех двупараметрических семейств схем, что согласуется со свойствами гармонических решений для самих моделей.
На рис. 1 изображены фазовые скорости гармонических решений, определяемые из дисперсионного соотношения для моделей первого класса и дисперсионного соотношения для линейных аналогов разностных схемы (8) и схемы (8) при условии (9). На рис. 2 изображены фазовые скорости гармонических решений, определяемые из дисперсионного соотношения для моделей третьего класса и дисперсионного соотношения для линейных аналогов разностной схемы (10) и разностной схемы с разнесенными разностями для моделей третьего класса. Дисперсионные соотношения для разностных схем вычислялись при 8 + и = 1, g = Но = 1, Ax = Ay = 0.065, At = 0.045.
Анализ рис. 1, 2 показывает, что фазовые скорости для гармонических решений (8) и (8) при условии (9) отличаются только для больших значений волновых чисел K1, K2. Так же, как и в одномерном случае, для модели Алешкова фазовые скорости гармонических решений модели и разностной схемы (10) для больших значений модуля волнового числа K отличаются сильнее, чем чем для модели Грина — Нагди и аппроксимирующих ее уравнений (8) и (8) при условии (9). Применение разностной схемы с разнесенными значениями позволяет улучшить дисперсионные свойства разностной схемы для моделей третьего класса.
Таблица
Схема Р Условие устойчивости
(8) Р1,2 = 1 - E1 + 3) ± ¿^(1 - ^ + 3)2 El), Р3 = 1, _ дН0(к2 sin2 £Су + sin2 фСх) — 2H0gx1x2 sin £ sin |Р1,2| < 1 для К1 = К2 2 < 2 К < + 3)2gH0
1 С С С 2 ' CxCy Cxy Cx = 1+4/3Я2 sin2(£/2)/Дж2, Су = 1 + 4/3Я2 sin2^/2)^y2, Cxy = H д£ д Ф
(8) при условии (10) Р1,2 = 1 — B21 (" + 3) ± ^Bi(1 — Bi), Р3 = 1, B1 = дН0(к2 sin2 £ + sin2 ф)£1; E1 = 1 1 1 + 4/3H02 (sin2 (£/2)/Дж2 + sin2 (ф/2) / Ду2) |Р1,2| < 1 К2 + < (* + 3)2gH0
(11) Р1,2 = 1 — B2 (" + 3) ± i^B (1 — B2), Р3 = 1, B2 = дН0(к2 sin2 £ + sin2 ф)Е2, E 1 + 2/3^2 (sin2 (£/2)/Дж2 + sin2(^/2)/Ду2) 2 1 + 2H02 (sin2 (£/2) / Дж2 + sin2 (ф/2) / Ду2) |Р1,2| < 1 К2 + К22 < („ + 3)2gH0
Leap-frog р2,2 = 1 — 2B ± ¿2^B(1 — B), р3,4 = 1, B = B1 для модели Грина — Нагди, B = B2 для модели Алешкова |Р1,2| = 1 к2 + к22 <
Схемы с разнесенными разностями Р1,2 = 1 — + 3) ± ^F(1 — +4 3)2 F), Р3 = 1, F = 4дЯ0(к2 sin2(£/2) + К sin2(ф/2))Е1 для модели Грина —Нагди, F = 4gH0(K sin2(£/2) + К sin2(ф/2))Е2 для модели Алешкова |рз| = 1 к2 + < (W + 3)2gH0
Схема вращения для модели Грина — Нагди E3 / E4 Р1'2 = 1 — AT±¿VA!' Р3 = 1' E3 = дН/2(к2(1 — cos £)(1 + cos ф)) + xf(1+cos £)(1 — cos ф) E4 = 1/4(к2 sin2 £(1 + cos ф)2 + sin2 ф(1 + cos £)2) |Р1,2| < 1 тах(к1, к2) < 1
Диссипативные свойства разностных алгоритмов иллюстрирует решение задачи о распаде начального возвышения при нулевой начальной скорости в квадратном бассейне, на сторонах которого ставятся условия отражения от твердой стенки. Рис. 3 дает начальное возвышение на сетке т х п:
п(ж,у) = А^г2 - (16.8 - ж)2 - (16.8 - у)2
п = 85, т = 85, г =10, Дж = 0.4, Ду = 0.4, А = 0.05
и решение этой задачи на момент времени 80Д£, 160Д£, Д£ = 0.1 по разностной схеме для модели Перегрина (8) с и + 8 = 2.
Рис. 4, 5 дают волновые картины на момент времени 200At, когда волна отразилась от стенок, для моделей Перегрина, нелинейной и линейной мелкой воды. Влияние сглаживающего параметра ш + 8 для моделей нелинейной и линейной мелкой воды более существенно, чем для нелинейно-дисперсионных моделей, так как в последних |р^| при малых Ax, Ay близки к 1 для всех ш + 8.
Список литературы
[1] KOMPANIETS L. A. Analysis of difference algoritms for nonlinear dispersive shallow water models. Russ. J. Num. Anal, and Math. Model., 11, №3, 1996, 205-222.
[2] GREEN A. E., NAGHDIP.M. A derivation of propagation in water of variable depth. J. Fluid Mech., 71, 1976, 237-246.
[3] ErtekinR. C., WEBSTER W. C., WEHAUSENJ.V. Waves caused by a moving disturbance in a shallow channel of finite width. J. Fluid Mech., 169, 1986, 275-292.
[4] ВольцинГЕР Н. Е., Клеванный К. А., ПЕЛИНОВСКИй Е. Н. Длинноволновая динамика прибрежной зоны. Гидрометеоиздат, Л., 1989.
[5] Базденков С. В., МОРОЗОВ Н. И., ПогуццеО.Р. Дисперсионные эффекты в двумерной гидродинамике. Докл. АН СССР, 293, 1987, 819-822.
[6] PEREGRINE D. H. Long waves on a beach. J. Fluid Mech., 27, №4, 1967, 815-827.
[7] ДОРФМАн А. А., ЯГОВДИК Г. И. Уравнения приближенной нелинейно-дисперсионной теории длинных гравитационных волн, возбуждаемых перемещениями дна и распространяющихся в бассейне переменной глубины. Числ. методы мех. спл. среды, 8, №1, 1977, 36-48.
[8] Железняк М. И., Пелиновский Е. Н. Физико-математические модели наката цунами на берег. В "Накат цунами на берег", (Под ред. Е. Н. Пелиновского), ИПФ АН СССР, Горький, 1985, 8-33.
[9] АлЕШКОВ Ю. З. Теория взаимодействия волн с преградами. Изд-во Ленингр. ун-та, Л., 1990.
[10] Железняк М. И. Воздействие длинных волн на сплошные вертикальные преграды. В "Накат цунами на берег" (Под ред. Е. Н. Пелиновского), ИПФ АН СССР, Горький, 1985, 122-139.
[11] GOURLAYA. R., MORRIS J. L. Finite-difference methods for nonlinear hyperbolic systems. Math. Comput., 22, №104, 1968, 28-39.
[12] РихтМАйЕР Р., МОРТОнК. Разностные методы решения краевых задач. Мир, М., 1972.
Поступила в редакцию 15 сентября 1995 г.
Рис. 1. Фазовые скорости, соответствующие: а — дисперсионному соотношению для моделей первого класса, б — линейному аналогу разностной схемы (8), в — линейному аналогу разностной схемы (8) при условии (10).
Рис. 2. Фазовые скорости, соответствующие дисперсионному соотношению для моделей третьего класса (а), линейному аналогу разностной схемы (11) (б), разностной схеме с разнесенными разностями (в).
Рис. 3. Волновые картины задачи о распаде начального возвышения, модель Перегрина, ш + 6 = 2 при г = 0 (а), 80Дг (б), 160Дг (в).
Рис. 4. Волновые картины на момент времени 200Д£: а модель Перегрина, и + 5 = 2 (а), и + 5 = 1 (б), модель нелинейной мелкой воды, и + 5 = 2 (в).
Рис. 5. Волновые картины на момент времени 200Д£: модель нелинейной мелкой воды, и + 5 = 1 (а), модель линейной мелкой воды, и + 5 = 2 (б), и + 5 = 1 (в).