ИСТОРИЯ РАЗВИТИЯ И АНАЛИЗ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ НЕЛИНЕЙНО-ДИСПЕРСИОННЫХ УРАВНЕНИЙ ГИДРОДИНАМИКИ. II. ДВУМЕРНЫЕ МОДЕЛИ Текст научной статьи по специальности «Математика»

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

Том 22, № 5, 2017

История развития и анализ численных методов решения нелинейно-дисперсионных уравнений гидродинамики. II. Двумерные модели

З.И. Федотова*, Г. С. Хлкимзянов, О. И. Гусев, Н. Ю. Шокинл Институт вычислительных технологий СО РАН, Новосибирск, Россия *Контактный e-mail: zf@ict.nsc.ru

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

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

Введение

В работе [1], посвященной истории развития численных методов решения нелинейно-дисперсионных (НЛД) уравнений гидродинамики, дан анализ основных свойств популярных разностных методов для одномерных моделей. Следующий шаг — исследование численных алгоритмов для случая двумерных в горизонтальной плоскости моделей (в русскоязычной литературе такие модели часто называют плановыми). Исторически стимулом для развития плановых НЛД-моделей послужила потребность в решении таких задач инженерной направленности, как описание волновых режимов в бухтах, взаимодействие волн с защитными сооружениями, а также плавучими или полупогруженными телами, и многих других задач, которые представляют интерес именно в двумерной постановке и с учетом "реального" рельефа дна. Такие модели для слабонелинейных длинноволновых течений впервые предложены в [2].

Следует отметить, что разработка оригинальных численных методов, их обоснование, исследование и конкретная реализация на ЭВМ сильно отстали от уровня разработки самих моделей. В упомянутой работе [2] построены НЛД-модели с двумя пространственными переменными, но их приложение ограничено рассмотрением одномерной задачи об отражении уединенной волны от однородного уступа, и конечно-разностная схема выписана только для этого случая. Однако плановые задачи имеют свою специфику, достаточно сравнить картины взаимодействия уединенной волны с прямой стенкой, расположенной перпендикулярно направлению движения волны, и с этой же стенкой в случае нахождения ее под углом к движению [3]. Численные алгоритмы привносят

© ИВТ СО РАН, 2017

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

Первые приложения НЛД-моделей к решению инженерных существенно двумерных задач были разработаны в Делфте (Нидерланды) и Хёрсхольме (Дания). В статье Abbott, Petersen & Skovgaard [4] описаны основные подходы к созданию численного метода для расчетов генерации и распространения "коротких волн" (как периодических, так и произвольной формы) в реальных акваториях с произвольной батиметрией. Однако и в этой работе конечно-разностная схема (вернее, ее наиболее важные фрагменты) приведена лишь в случае одной пространственной переменной, а ее описание, как и в работе [2], вынесено в Приложение.

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

Иной подход предложен и подробно исследован в книге [7], где численный алгоритм решения полных 2В-НЛД-уравнений Железняка — Пелиновского опирается на расщепление исходных уравнений на гиперболическую систему с правой частью и скалярное уравнение эллиптического типа.

Ввиду того, что НЛД-модели не поддаются в полной постановке аналитическим исследованиям, именно численный алгоритм отвечает за то, насколько адекватно будет решена исходная задача. Например, улучшая НЛД-модель повышением порядка аппроксимации дисперсионных характеристик, можно эту модификацию свести "на нет" влиянием "схемной" дисперсии, вносимой численным методом. Следует отметить, что исследование двумерного дискретного метода — трудная задача, даже в упрощенных формулировках. Чаще всего статьи, предлагающие метод для решения задач в двумерной области, дают описание свойств только для одномерного случая.

Опубликовано очень мало работ, которые содержат подробное описание и тестирование численных методов для расчета плановых моделей. Так, в недавно опубликованном обширном обзоре [8], посвященном НЛД-моделям с точки зрения "взаимодействия физики, математики и численных методов", непосредственно обзору численных методов отводится весьма скромное место. Что же касается конечно-разностных методов, упоминается лишь один метод, построенный в [4].

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

чением времени. Ввиду сложности рассматриваемых двумерных моделей даже метод фон Неймана, простейший из способов изучения устойчивости, позволяющий определить условия невозрастания компонент Фурье в разложении решения линеаризованных уравнений, вызывает определенные затруднения. Мы провели эти исследования для нескольких известных разностных схем, аппроксимирующих 2В-уравнения Перегрина. При этом были получены результаты, отличные от тех, что имеют место для аналогичных разностных схем, аппроксимирующих бездисперсионные уравнения мелкой воды.

1. Базовая нелинейно-дисперсионная модель и ее частные случаи

Чтобы охватить разнообразные подходы к построению численных методов, будем отталкиваться от базовой НЛД-модели [9, 10], которая описывает течение идеальной несжимаемой жидкости в слое, ограниченном снизу подвижным дном г = -к(х,у,Ь), а сверху — свободной границей г = Г)(х,у,Ь), где £ — время; х, у, г — координаты точки в декартовой системе координат Охуг, ось Ох которой направлена вертикально вверх, а координатная плоскость Оху совпадает с невозмущенной поверхностью воды. НЛД-уравнения выводятся из трехмерных (3В-) уравнений Эйлера, описывающих это течение, в предположении, что оно имеет длинноволновой характер, т. е. отношение ^ = к0/Ь, где Ь, к0 — характерные размеры по горизонтали и глубине соответственно, мало, так что при выводе величинами порядка 0(^4) пренебрегают. В НЛД-модели искомыми величинами являются полная глубина Н = к + ^ и некоторая приближающая горизонтальную 3В-скорость вектор-функция и(х,у,Ь), которую рассматривают как скорость приближенной модели. Обычно принимается, что и = (и, у) — это либо горизонтальная составляющая и вектора скорости уравнений Эйлера на некоторой поверхности г = (х, у, ¿), т. е.

и(х,у,г) = и(х,у,га (x,y,t),t), (1)

либо средняя по толщине слоя воды скорость

V

и = J И(%,У,г) (2)

-к

При выводе базовой НЛД-модели считается, что и приближает горизонтальную составляющую И скорости 3В-модели с точностью до членов порядка 0(^2):

И(ж, у, г, г) = и(ж, у, г) + ^2Иа(х, у, г, г). (3)

При этом условии и удержании членов до 0(^2) в [9] получена базовая НЛД-модель, уравнения которой имеют следующий вид:

Нь + V• (Ни) = -V- (ни), (4)

и4 + (и • V)u + ^ = ^Vh - 1 \(Ни)4 + (и • V)(HU) + Н (Ы •V) и + (V • и)НЫ Н Н Н -

(5)

где р — проинтегрированное по толщине слоя давление НЛД-модели, р0 — давление на дне,

ц2 Н 3 Н 2 Н 2

Р = У~Н~ -р, Р0 = дН - ф, р = —Пг + —Я2, ф = — Кг + НК2, (6) Кг = В (V • и) - (V • и)2 К = 02к, (7)

V

Л Г Л

и = - / и ж, В = - + и •V.

-н

Отметим, что уравнение движения (5) допускает дивергентную форму записи с ис-точниковым членом p0Vh, обусловленным неровностью поверхности дна:

( Ни)4 + V • ( Ни 0 и) + Vp = poVh - Г(Ни) + V • ( Ни 0 и + НЫ 0 и)

Вид вектор-функции Ы задает конкретную модель [9, 10]. Если скорость определена по формуле (2) либо изначально предполагается независимость горизонтальной составляющей вектора скорости трехмерного течения от вертикальной координаты г, то и = 0. В этом случае квазидивергентная форма базовой модели принимает следующий вид:

Н + V • (Ни) = 0, (9)

( Ни)4 + V- ( Ни 0 и) + Vp = poVh. (10)

Термин "квазидивергентная форма" используется в том смысле, что левая часть уравнений записана в дивергентной форме, а правая обращается в нуль на ровном дне, т. е. на ровном дне уравнение (10) принимает вид закона сохранения полного импульса.

В форме (9), (10) записываются известные НЛД-модели Грина — Нагди [11], Железняка— Пелиновского [12], Федотовой — Хакимзянова [13, 14] и др. Из этой модели можно получить целый ряд более простых моделей. Например, упрощая выражение для "давления" путем отбрасывания в р/Н, р0/Н членов, имеющих порядок 0(ац2), где а = а0/^ (а0 — характерная амплитуда) и ^ — параметры нелинейности и дисперсии соответственно, приходим к уравнениям типа уравнений Буссинеска [14], которые имеют тот же вид, что и уравнения (6), (9), (10) полностью нелинейной дисперсионной модели. Отличаются лишь формулы для вычисления величин р и ф — дисперсионных составляющих давления, которые в случае слабонелинейной модели из [14], являясь нелинейными функциями, зависят от полной глубины Н линейным образом:

р = и) + Ф = н(^ и) + В2^. (11)

Первые конечно-разностные схемы для расчета двумерных длинноволновых течений с дисперсией были основаны на аппроксимации слабонелинейной дисперсионной модели Перегрина, выведенной в работе [2] для случая стационарного дна (^ = 0). Уравнение неразрывности модели Перегрина имеет тот же вид (9), что и у полной НЛД-модели, а уравнение движения получается из недивергентной формы уравнения (10) при отбрасывании в нем всех нелинейных слагаемых, входящих в дисперсионные члены (V р)/Н, ф/Н .С учетом формул (11) получаем

и + (и • V)u + gVГ]

hhv(v• ^и)) - ^V (V- и)

= ^Рег. (12)

В работе [15] из модели Перегрина (9), (12) выведено уравнение второго порядка, описывающее распространение с учетом силы Кориолиса квазиодномерных слабонелинейных длинных волн над слабо изменяющимся дном:

щ - gV • (hVrj) = -f 2V + 1V

l 3

hv (Нщц + ^gV • V (r]2).

:i3)

Впервые это уравнение (без вывода) приведено в работе [16], поэтому его называют уравнением Кима — Райда — Витакера. В работе [17] уравнение такого вида выведено с учетом трения. Оказалось, что уравнение вида (13), несмотря на его простоту, можно применять для решения практических задач волновой гидродинамики. Так, например, в [18] успешно решены сложные задачи взаимодействия достаточно длинных волн с различного рода препятствиями. В статье [19] описано применение модификации модели Кима — Райда — Витакера для расчета трансокеанических цунами.

Для изучения свойств разностных схем для НЛД-уравнений необходимо привлекать их упрощенные формулировки. В двумерном случае самой простой моделью является скалярное линейное уравнение Кима — Райда — Витакера

h2

% - c20A<q = -3△щt,

:i4)

где Со = у/дк0; А = V2 — оператор Лапласа. Следующей по простоте моделью, которая содержит уже три независимых переменных, является система линейных уравнений

Перегрина

rjt + h0V

u

0, ut + gV'q

h0V (V • u)t,

:i5)

к которой в случае ровного дна приводятся при линеаризации в окрестности нулевого решения все упомянутые выше системы НЛД-уравнений со скоростью (2).

В том случае когда скорость определяется по формуле (1), известные выводы НЛД-уравнений существенно опираются на предположение о потенциальности 3В-течения [6, 20, 21]. При таком предположении базовая модель (4), (8) замыкается соотношением

=

у - (za + h) (VDh +(V • u)Vh) -

. 6

(Za + h)

2 -,

2

V(V • u).

:i6)

Применяя к уравнениям (4), (5), (16) условие Буссинеска и отбрасывая из дисперсионных членов все нелинейные (относительно ^ и и) слагаемые, получаем слабонелинейную дисперсионную модель для стационарного дна (к = 0):

Ht + V • (Нu) = -V • (hU),

:i7)

ut + (u • V)u + gV'q = TPer - Uu (i8)

где

U ={2 + V(V • (hu)) + (f - hpj V(V • u).

Для горизонтального дна h = h0, горизонтальной поверхности za = const и течения с нулевой фоновой скоростью линеаризованную базовую модель можно записать в следующем виде:

Vt + hoV • u = -hl($ + u), ut + gV'q = -ffh0V (V • u)t, (i9)

где при 3 = /Ьо + г2/2Ь0 имеем уравнения Нвогу [22], а при 3 = — 1/3 — линеаризованный вариант (15) моделей с усредненной скоростью.

Для исследования свойств выписанных уравнений, а также свойств аппроксимирующих их разностных схем будем анализировать поведение гармоник

ф, у, г) = Е0е—(ш4-кх), и(х, у, г) = ^е-гИ-кх), ь(х, у, г) = Уое-г(шЬ-кх), (20)

где Е0, и0, Уо — амплитуды гармоник; к = (к\, к2) — волновой вектор; х = (х1 ,х2) = (х, у); ш — частота. Подставляя гармоники в систему (19), получаем следующие выражения для частот:

Ш1 = 0, Ш2,з = ±Со| ки 1 —

1

ЬЦ к|2/3

1 — 3Ь0| к|2

(21)

Дисперсионные соотношения модели Перегрина (15) получаются из (21) при 3 = —1/3. Нетрудно показать, что уравнение Кима — Райда — Витакера (14) имеет две моды и они совпадают с ш2,3 из (21) при 3 = —1/3.

Для оценки дисперсионных свойств моделей мелкой воды фазовые скорости ур = ш(к)/к, соответствующие этим моделям, можно сравнить с "эталонными" значениями

АапЬ(2тгЬо/ А) \ 1/2 СЧ 2ттко/ А )

22)

2жко/А

модели потенциальных течений [7], в которой, в отличие от уравнений мелкой воды вертикальные перемещения жидкости учитываются. Здесь А = 2^/| к| — длина периодической волны (20). Далее, не ограничивая общности, будем рассматривать волны соответствующие дисперсионным соотношениям (21) со знаком "+". Тогда для величин ур в НЛД-модели с усредненной скоростью и в модели Нвогу [22] со скоростью заданной на поверхности = — 0.531Ьо, получаем следующие выражения:

б

1.0

0.6

0.4

0.2

0.0-

уу / лг / / ——— 1 ^—1—

■1— ' / / у / / / / ' / » » « 3 -----4

' г ' / |( . . . 1 ■ ■ ■ ■ .... ■ ■ ■ ■

АУйд

Рис. 1. Дисперсионные свойства гидродинамических моделей: а — зависимость фазовой скорости от длины волны для моделей Кима — Райда — Витакера и нелинейно-дисперсионной (1), Нвогу (2), бездисперсионной (3) и потенциальных течений (4); б — относительные отклонения (24) фазовых скоростей в моделях мелкой воды (кривые 1-3) от "эталонных" значений (22)

а

( 1 У/2 _ ( Л)2/3 V'2

VpNbD = СЧ1 + (2nho/\)2/ъ] ' ^Nw = Ч1 - 1 - ß(2nh0/X)2 ) , (23)

где ß = -0.390. Графики зависимостей фазовых скоростей от Л изображены на рис. 1, а. На рис. 1, б показаны графики отклонений

S = 1 Vp - Vp'PFl 100% (24)

Vp, PF

фазовых скоростей в моделях мелкой воды от "эталонных" значений (22). По графикам отклонений видно, что различие не превышает 5 % для волн с длиной Л ^ 4ho в НЛД-модели и Л ^ 2h0 в модели Нвогу. Таким образом, последняя модель воспроизводит дисперсионные свойства 30-модели в более широком диапазоне длин волн, чем модель Перегрина и тем более, чем бездисперсионная модель мелкой воды (кривая 3 на рис. 1). Вместе с тем надо отметить, что в численном решении, полученном на мелких сетках, могут присутствовать гармоники с очень малыми длинами волн, для которых, как это видно из рис. 1, б, предпочтительнее модель Перегрина, поскольку для нее, в отличие от модели Нвогу, lim i»pNLd( Л) = 0.

2. Численные алгоритмы, основанные на аппроксимации НЛД-модели Перегрина

При расчетах двумерных длинноволновых течений с дисперсией, базирующихся на НЛД-модели Перегрина (9), (12), часто используется форма уравнений, записанных для потоков Н, uH, vH, так что левая часть этих уравнений, представляющая собой классическую модель мелкой воды, имеет дивергентный вид. Такая форма выбрана, в частности, для конечно-разностного алгоритма, включенного в компьютерную систему SYSTEM 21, предназначенную для решения широкого круга инженерных задач в прибрежных водах [4]. Эта система получила широкое распространение во многих странах, имеющих освоенный выход к прибрежным морям. В [23] отмечено, что в 1970-х гг. было заключено 44 контракта на использование вариантов системы в 13 странах.

Изначально в этой системе использовались бездисперсионные уравнения мелкой воды, для решения которых применялись известные в гидрологии схемы Абботта — Ионес-ку и Прейсмана [24, 25]. Первая из схем характеризуется применением разнесенной сетки в плоскости пространственных переменных и при подходящей аппроксимации по времени использует только центральные разности, что исключает численную диффузию. Вторая схема, имеющая компактный шаблон, в русскоязычной литературе известна как схема "ящик". Она допускает очень удобную реализацию. В [4] использованы модификации этих схем применительно к уравнениям с дисперсионными членами, что позволило моделировать волновые режимы в расширенном диапазоне длин волн.

Так как численное моделирование в рамках системы SYSTEM 21 предполагало решение практических задач о взаимодействии волн и гидротехнических сооружений, то основные уравнения были расширены учетом пористости волноломов. Выпишем эти уравнения, взяв за основу модель из [23]:

nHt + V • (Нu) = 0, (25)

Н 2

n (Нu)t + V • (Нu 0 u) + n2V^ = n2gHVh + Нu$ + nHTPer . (26)

Здесь п — объем пор (пористость),

3Р0и , п - 1^1

* = (п - 1ГЪг + — Т1|и|;

v — кинематическая вязкость; 0О, — коэффициенты сопротивления ламинарного и турбулентного потоков для пористого волнолома соответственно. При п =1 получаем НЛД-уравнения Перегрина (9), (12), записанные для потоков в случае неподвижного дна произвольной формы [4].

При создании численного алгоритма, привязанного к практическому назначению системы SYSTEM 21, разработчики сформулировали ряд требований, обеспечивающих нужную точность расчетов на достаточно крупных сетках, использование которых продиктовано сравнительно невысокой производительностью компьютеров в начале 1980-х гг. Прежде всего было отмечено, что для создания работоспособной численной модели необходимо знать свойства аппроксимируемых НЛД-моделей и условия их вывода, с тем чтобы осуществлять по мере необходимости грамотную модификацию дисперсионных членов, позволяющую улучшить математические свойства модели без изменения порядка длинноволновой аппроксимации (впоследствии аналогичные модификации применялись по отношению к другим НЛД-моделям).

Далее, отталкиваясь от опыта пробных расчетов в рамках НЛД-моделей, авторы отметили сильную зависимость численного решения от вида аппроксимационной ("схемной") ошибки главного порядка. В частности, отмечено, что для достижения требуемой в расчетах точности необходимо использовать схемы третьего порядка аппроксимации, так как для описания реального волнового режима в гаванях на компьютерах 1970-х гг. приходилось на длину волны брать всего несколько (в среднем 6) узлов расчетной сетки.

Важное место в разработке компьютерной модели отведено постановке фильтрующих и неотражающих условий на открытых границах, а также описанию системы тестовых и демонстрационных 2В-задач, на которых были показаны возможности вычислительной модели.

Что касается непосредственного описания конечно-разностного метода, то авторы [4] ограничились описанием основных подходов применительно к (25), (26) в случае ровного дна, одной пространственной переменной и п = 1, выбрав в качестве отправной точки схему Прейсмана. В построенной разностной схеме используется трехслойный по времени и двухточечный по пространству шаблон, а для аппроксимации по времени и пространству применяются центральные разности относительно центра шаблонного "ящика". Такой выбор обусловлен тем, что использование центральных разностей устраняет главные члены аппроксимационной вязкости, а третий порядок аппроксимации достигается методом исключения "дисперсионных слагаемых" схемного происхождения путем добавления в разностную схему тех же слагаемых, но с противоположным знаком (метод исключения). Аналогичный подход, как указано в [4], применен к построению разностной схемы повышенного порядка для системы НЛД-уравнений на основе метода Аббота — Ионеску, использующего разнесенную сетку по пространственным переменным.

Однако следует отметить, что метод исключения просто применять для случая одномерных линейных систем уравнений гиперболического типа. В случае двумерных нелинейных НЛД-уравнений, как отмечено в [23], приходится исключать свыше 100 членов в главных членах аппроксимационной ошибки.

Впоследствии было опубликовано несколько работ, посвященных модификации конечно-разностных методов из [4]. В работе [26] рассмотрена разностная схема из [4], использующая разнесенную сетку. Манипулируя производными третьего порядка, имеющими "схемное" происхождение, а также введя дополнительные дробные шаги, авторам [26] удалось повысить фактическую точность разностной схемы. В частности, устранен эффект численной гравитации.

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

3. Численные алгоритмы для моделей Буссинеска

с улучшенной аппроксимацией фазовых характеристик

Развитие подходов к построению разностных схем для НЛД-уравнений с уточненным описанием фазовых характеристик было продолжено при создании компьютерной системы MIKE 21BW (Danish Hydraulic Institute, Danish Technical University), где за основу взята НЛД-модель с улучшенной аппроксимацией дисперсионного соотношения. Улучшение достигнуто путем модификации слабонелинейной модели с усредненной скоростью (2). Новшество заключалось во введении дополнительного параметра В и добавлении новых дисперсионных членов в правую часть уравнения движения. Модифицированные в [27] НЛД-уравнения, выведенные для случая ровного дна, можно записать в следующем виде:

щ + V ■ ( Нu) = 0, (27)

.Н2

~2

где

(Нu)t + V- ( Нu 0 u)+ gV — = TMad, (28)

TMad = [В + 3) h2v(v ■ (Нu)t) + Bgh3V [V ■ Vr,]. (29)

Дисперсионное соотношение линеаризованного варианта этой модели имеет вид (21) при [ = —В — 1/3. Таким образом, при некоторых значениях В получаются известные модели. При В = 1/15 имеем модель, у которой аппроксимация дисперсионного уравнения достигает четвертого порядка (в одномерном случае). При В = 0 система уравнений (27)-(29) совпадает с рассмотренной в работе [26], имеющей в линейном случае то же самое дисперсионное соотношение, что и модель Перегрина. Аналогичная модель, но уже в сферических координатах использована в работе [28].

Приведем разностную схему [27], допускающую реализацию (27)-(29) одномерными прогонками (конвективные члены опущены). На первом шаге

,пп+1/2 _ 'п 1 1

1 т/2 1 +1 {РП+1 + РП) +1 {ОП+1/2 + ЯП-1/2) = 0, (30)

~ + з) ~ [(Рхх — Рхх) + / — ^п- / )] +

+дН*г]-П+1/2 — Вдк3 {г--- + Ч*хуу) = 0 (31)

неизвестными являются г]п+1/2 и Рп+1, где Р = иН, Q = ьН. На втором шаге определяются г]п+1 и Qn+3/2:

г]п+1 - т]п+1/2 1 1

+ Ö (PT1 + PX) + Ö {QT3/2 + Qny+1/2) = 0, (32)

r/2 2 v X XJ 2

Qn+3/2 — Qn+1/2 fn i 1 ^ h2

3

— {в +1) h22 [(P-+1 — PX,) + (Q-3/2 — Q-1/2)] +

+gH**— Bgh3 (^ + CJ = 0. (33)

Величины, помеченные одной или двумя звездочками, являются оценочными (estimated) значениями, вытекающими из вычислений по явной схеме, аппроксимирующей уравнение неразрывности на временных слоях (п + 1/2) и (п + 1) соответственно. Этот прием впервые использован в [26] для того, чтобы нелинейный гравитационный член gH'W'q вычислять без применения итераций. Производные по пространству аппроксимированы центральными разностями на разнесенной сетке. Благодаря указанной структуре шаблона по времени (разнесенные узлы) расчет удается выполнять поочередно в направлениях х и у, используя прогонки. В работе [29] метод распространен на случай переменного дна.

Указанный подход к построению разностной схемы был использован для расширения рамок компьютерной системы MIKE 21, основанной на бездисперсионной модели, включением в нее НЛД-моделей. Полученная система MIKE 21BW широко применяется для исследования распространения волн в прибрежных морях. Наиболее эффективно система работает в совмещенном режиме с программой SWASH [30], обеспечивающей расчет взаимодействия волны, проникающей в гавань, с берегом по нелинейной бездисперсионной модели мелкой воды [31, 32].

В работе [22] сделан новый шаг в построении моделей Буссинеска с повышенной аппроксимацией дисперсионного соотношения не только для очень длинных, но и умеренно длинных волн. Этого удалось достичь путем выбора скорости модели по формуле (1), в которой вектор u определяется как горизонтальная составляющая скорости 30-течения жидкости на поверхности z = za(х, у) = —0.531h(x, у). Эта модель имеет преимущество перед другими моделями, полученными с той же целью формального расширения области применимости в сторону умеренно длинных (коротких) волн. Оно заключается в том, что модель выведена для двумерного случая и произвольной донной поверхности h = h( х, ) путем асимптотических разложений основных функций 30-модели без использования дополнительных искусственных приемов. Впоследствии модель Нвогу претерпела несколько обобщений и для них были разработаны конечно-разностные методы высокого порядка аппроксимации. Что касается исходной статьи [22], то в ней описан лишь одномерный вариант разностной схемы типа Кран-ка — Николсон с исключением членов "схемной" дисперсии по тому же алгоритму, как это предложено в [4]. Разностная схема, аппроксимирующая систему уравнений Нвогу с четвертым порядком аппроксимации [33], будет выписана в другом параграфе.

4. Численные алгоритмы для моделей Буссинеска,

использующие выделение эллиптической части уравнений

Работа [34] является, по-видимому, первой публикацией, где подробно описан конечно-разностный метод решения двумерных НЛД-уравнений и приведено условие устойчи-

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

^ х у(у ■ (м) - ПVП XV (V ■ и)

ш

П.

(34)

справедливого с точностью 0(ц2 + а) для течений, потенциальных в 3В-постановке. Здесь ш = ух — иу и для двумерных векторов использована операция х, обозначенная так же, как операция векторного произведения для трехмерных векторов. Результатом действия этой операции на векторы а, Ь € И2 (а = (а 1, а2), Ь = (Ь1, Ь2)) будет не вектор, как в И3, а скаляр, определяемый по формуле а х Ь = а-\Ъ2 — а2Ь1.

Ввиду соотношения (34) и равенства ш = 0(ц2) вместо исходного уравнения движения (12) можно с сохранением порядка длинноволновой аппроксимации использовать [34] следующее:

и +

(П и )хх + (П иу )у + (кхУ)

П2 ,гг

--А и,

6 '

(35)

V +

йг+ дг])

(П Ух)х + (П V )уу + (Пуи )а

П2

— 6 А V, 6

где

и = и

(36)

(37)

( и,у).

Таким образом, вместо системы уравнений (9), (12) в работе [34] рассмотрена система (9), (35)-(37), где (35) и (36) являются уравнениями второго порядка, причем старшие производные от и содержатся только в уравнении (35) и соответственно в (36) входят старшие производные только от V. Формулы (37) можно трактовать как замену переменных, позволяющую записать уравнения движения в виде эллиптических уравнений относительно новых зависимых переменных и и V.

Конечно-разностная схема выписана на разнесенной сетке Аракавы-С, которая предполагает не только разнесение переменных по пространственным сеточным узлам, но и чередование целых и полуцелых узлов сетки по времени (схема "чехарда"). На каждом шаге по времени сначала с помощью метода предиктор-корректор решается уравнение неразрывности (9). Затем методом последовательной верхней релаксации находятся решения и и V эллиптических уравнений (35), (36). Наконец, путем решения обыкновенных дифференциальных уравнений (37) вычисляются компоненты скорости и и V.

4.1. Линеаризованная разностная схема Ригга

Для горизонтального дна П = П0 и течения с нулевой фоновой скоростью модель (9) (35)-(37) можно записать в следующем виде:

Т]1 + ПoV

и

0, И + gVr]

ПО АИ,

3 '

и = И.

(38)

Выпишем разностную схему Ригга для этой системы уравнений.

В соответствии с выбранной разнесенной сеткой функция определена в узлах (х^, Уj2) прямоугольной равномерной сетки с шагами Ах и Ау, а компоненты скорости — в центрах сторон ячеек: и — в центрах (х^1+1/2, ) горизонтальных сторон ячеек, V — в центрах (х^1, У]2+1/2) вертикальных сторон. Считается, что компоненты вектора И

у

X

У

определены в тех же точках, что и компоненты вектора скорости. Сеточные функции разнесены также по времени: значения г/ определены на временных слоях Ь = пт, тогда как и, V, и, V — на слоях Ь = (п + 1/2)т, п = 0, 1, ... С учетом областей определения будем использовать сокращения для значений сеточных функций: = ^,

п+1/2 _ п+1/2 п+1/2 _ п+1/2

и 31+1/2 = иЦ+1/2,32 , +1/2 = У31,32 + 1/2 и т.д.

Аппроксимируем систему уравнений (38) следующей схемой второго порядка:

+ Ь>( иП+1/2 + уП+1/2) = 0,

^ V х,з У,з )

ип+1/2 + апп+1 = (тт + ит )п+1/2 (40)

Т 31 + 1/2 + УГ1°х^1 + 1/2 = 3 (и-х + иУУ)31 + 1/2 , (40)

(39)

(41)

h2

3

„п+1/2 = ттП+1/2 п+1/2 = Vn+1/2 (42)

Ut,h+1/2 = U 31 + 1/2, %j2 + 1/2 = Vh + 1/2. (42)

Vn+1/2 + а'пп+1 = h0 (V- + V- )п+1/2 Vh+1/2 + У\h+1/2 = 3 (Vxx + VW)h+1/2

При записи разностных уравнений (39)-(42) использованы следующие обозначения:

^п+1 <п ип+3/2 _ ип+1/2 п+3/2 _ П+1/2

п = Щ_..п+1/2 = 31+1/"2 и 31 + 1/2 п+1/2 = % + 1/2 %+1/2

т , \п + 1/2 т , и1,п + 1/2 т ,

п+1/2 п+1/2 п+1/2 п+1/2

.п+1/2 = ип + 1/2 — ип-1/2 уп+1/2 = +1/2 — -1/2

х,з Ах ' у,з Ау '

71п+1 _ П+1 П+1 _ П+1

п+1 = >31 + 1,32 Чз Г1П+1 = 21+ 1 3

-,31+1/2 Ах , У ,32+1/2 Ау ,

ип+1/2 — ип+1/2 ип+1/2 _ 2ип+1/2 + ип+1/2

и п+1/2 = -,31+1,32 х,з и п+1/2 = и 31+1/2,32 + 1 31 + 1/2,32 ^ и 31+1/2,32-1

и-х,31+1/2 = АХ , УУ,31+1/2 = Ау 2 ,

уп+1/2 _2Лгп+1/2 +^+1/2 — Уоп+1/2

Vп+1/2 = у 31 + 1,32+1/2 ^У 31,32+1/2 * 31-1,32+1/2 уп+1/2 — У,31,32+1 У,3

хх,32+1/2 Ах2 ' УУ^2+1/2 Ау

Для исследования устойчивости схемы (39)-(42) в разностные уравнения подставим гармоники

= Ео рпег(-Э1 , и^Ц = То рп+1/2е Ш+1/2)*+М, = у0 рп+1/2е гт+(32+1/2Ю,

ип+1/2 = и* п+1/2 г((.31+1/2)г+Ы) уп+1/2 = у * п+1/2 Л(з1^+(п + 1/2)0

и31 + 1/2 = и0Р е , Ул+1/2 = У0 Р е ,

где ^ = ^Ах, ( = к2Ау. Из уравнений (40), (41) получаем, что величины ио*, У* выражаются через Е0:

jh/f^2xsrn(£/2)^ g yß 2«sin( (/2) Vh + Ь2

где использованы обозначения

4 r2 • 2

-Oi sin -, „. „2 „„ , * , л

3 1 2 2 3 2 2 1 Ах' 2 Ау

u0 = -1 + h + b2 Ах E0 = -

4 2 2 4 2 2 h0 h0 b1 =-ö 2 sin2 -, b2 = - 62 Sin2-, = , ö2 = . (43)

Поэтому из уравнений (39), (42) получается однородная линейная система алгебраических уравнений относительно амплитуд Е0, и0 и V. Для существования нетривиального решения определитель этой системы должен обращаться в нуль:

р — 1 2г ^ph0ai sin(£/2) 2i у/р h0a2 sin(£/2)

2 г^p дац sin(£/2) (p — 1)(1 + bi + b2) 0

2 г^p ga2 sin((/2) 0 (p — 1)(1 + bi + b2)

(44)

где

а1

(45)

——, ж2 = ——. А х А

Равенство (44) дает кубическое уравнение для множителя перехода р. Один из корней этого уравнения равен единице, а два других р1,2 являются корнями следующего квадратного уравнения:

(46)

р2 — 2р (1 — R) + 1 = 0,

где

R

2 с0

2 2 2 2 , , а 1 sin —+ а2 sin -1 + bi + &2V 1 2 2 2

(47)

Необходимое условие устойчивости | р1,2 | < 1 рассматриваемой схемы эквивалентно условию К < 2 неположительности дискриминанта уравнения (46), а это эквивалентно (см. Приложение А) выполнению неравенства

1

4

с0 а

<\ 1 + Ö

где

При этом имеем

г ho

о = "Т-,

Ао'

а =

А0 ■ Ао

1 Pi,2 1 = 1.

Ах • Ау у/Ах2 + Ау2'

(48)

(49)

(50)

Далее предполагаем, что условие устойчивости (48) выполнено. Информацию о дисперсионных свойствах схемы дает поведение ее фазовой ошибки Аф = Ф — ф [35]. Здесь Ф = шт, ш — дисперсионное соотношение математической модели, ф = arg р. Изменение фазы Ф оператора перехода математической модели и изменение фазы ф множителя перехода разностной схемы будем рассматривать как функции от переменных и (£,С Е [0, Согласно (21), за один шаг по времени изменение фазы в решении линеаризованной НЛД-модели (38) определяется по формуле

^ld = ™nld

со | k | т

оа

где

л/1 + h21 k|2/3 у/1 + 62я2/3'

к|Ао.

(51)

(52)

Из формулы (51) следует, что для модели (38) изменение фазы решения можно записать в виде функции от одной переменной € [0, ж]. Кроме того, изменение фазы не зависит от направления волнового вектора к. Отметим также, что изменение фазы решения, выражаемое формулой (51), одинаково для моделей (15) и (38).

0

Взяв множитель перехода схемы Ригга со знаком "+"

p=1 -R + i^R(2 - R) (53)

и используя равенство (50), получаем (см. Приложение А) следующую формулу для определения изменения фазы = arg р в решении разностной схемы:

0Kg = 2 arcsin y/RJ2. (54)

Рассмотрим теперь малые значения (, соответствующие длинным волнам. Тогда (см. Приложение А)

**=О - Ч- - +^+^

Из формул (51) и (55) следует, что фазовая ошибка определяется выражением

Д0 = £ (ж1$4 + ж2(4 - + 0(я4). (56)

24

Таким образом, дополнительная дисперсия, вносимая разностной схемой, имеет тот же порядок <^3, что и второй член разложения по я "физической" дисперсии (дисперсии модели)

Фкьв = Оо^ - ^С3 + 0(^4). (57)

6

При конечных шагах разностной сетки это может привести к искажению дисперсионной картины течения, описываемого НЛД-моделью.

В [1] показано, что в одномерном случае в области устойчивости некоторых схем имеются такие значения числа Куранта, для которых влияние схемной дисперсии будет минимальным. Для двумерных схем это свойство, к сожалению, нарушается. Например, для произвольного волнового вектора к фазовая ошибка (56) имеет порядок 0(я3) независимо от числа Куранта с0ж. Исключение составляют волновые векторы к = (к\, 0) и к = (0, к2), для которых фазовая ошибка имеет порядок 0(<^4), если задать с0Ж1 = 1 и с0ж2 = 1 соответственно. Порядок малости 0(<^4) достигается также для волнового вектора к, компоненты которого соответствуют шагам сетки согласно формуле к1Дх = к2Ду. Для того чтобы влияние схемной дисперсии было минимальным, необходимо задать

с0ж = 1. (58)

В этом случае фазовая ошибка как минимум на один порядок по будет меньше "физической" дисперсии.

4.2. Вариант модели Перегрина с выделением эллиптического уравнения

В работе [36] рассмотрена аналогичная модификация уравнений Перегрина, однако расставлены другие акценты при построении численного алгоритма. Модифицированное уравнение движения в этой работе рассмотрено в виде

Ъ

и + (и ■ У)и + дЩ = Ъ - , (59)

где

Ь ( д д\

Е = А(Ьи) — -Аи, = —, ~ ] , С = дУЬ хУг],

3 д д х

при этом в случае слабо меняющейся поверхности дна величиной С можно пренебречь. Если ввести новую переменную (37), то из (59) следует уравнение

h ( h \ h

- ¡A(hU) - 3AU) - U = (u ■ V)u + gVV + ^ Vя* G, (60)

которое при замороженной правой части является эллиптическим относительно U. Предложена следующая структура вычислительного алгоритма для системы уравнений (9), (37), (60): из уравнения неразрывности находится г]п+1 с использованием вычисленного на предыдущем временном слое вектора un. Затем методом верхней релаксации решаются два независимых эллиптических уравнения (60) для нахождения компонент вектора U. Наконец, решение ОДУ (37) дает вектор скорости на (п + 1)-м слое по времени.

В работе [36] система уравнений (9), (37), (60) аппроксимирована на разнесенной в 3В-пространстве (х,у, t) сетке с использованием центральных разностей. Полученная разностная схема имеет второй порядок аппроксимации. В [36] подробно описан численный алгоритм, отличающийся в нелинейном случае от [34]. Для линеаризованных уравнений и горизонтального дна (h = h0) рассматриваемая разностная схема совпадает с (39)-(42). Следовательно, анализ Неймана приводит к тому же необходимому условию устойчивости (48). Очевидно, что и все другие свойства разностной схемы, полученные на основе анализа линеаризованных уравнений, также совпадают. В силу (50) разностная схема является недиссипативной, и для расчета нелинейных задач рекомендовано вводить искусственную диссипацию.

4.3. О численных алгоритмах для уравнения Кима — Райда — Витакера

В [15] построение разностного метода для решения уравнения Кима — Райда — Витакера (13) основано на его расщеплении на двумерное волновое уравнение с "правой частью" и эллиптическое уравнение. Разностная схема для волнового уравнения является явной трехслойной и использует пятиточечный шаблон для аппроксимации вторых производных по пространственным переменным. Эллиптическое уравнение решается методом верхней релаксации. Условие устойчивости и дисперсионные характеристики в [15] даны для случая одной пространственной переменной (одномерный вариант этого численного алгоритма обсуждался в [1]). Для линеаризованного варианта двумерной модели (14) эта разностная схема имеет следующий вид:

п h2 п

Vtt — о ( Vxx + Vyy) = "3 ( r!xx + rrlyy)tf (61)

Нетрудно показать, что ее множитель перехода удовлетворяет квадратному уравнению (46). Следовательно, все свойства, полученные выше для разностной схемы Ригга, выполняются и для разностной схемы (61).

В работе [18] уравнение Кима — Райда — Витакера, дополненное учетом поверхностного натяжения и придонной диссипации, аппроксимировано разностной схемой, аналогичной [15]. Приведено без вывода условие устойчивости, зависящее от числа Бонда,

характеризующего капиллярные эффекты. В случае ровного дна, линеаризации на нулевом фоне, а также при отсутствии поверхностного натяжения и трения это условие совпадает с (48).

В работе [19] предложен метод расщепления уравнения Кима — Райда — Витакера, аналогичный [15], и построена разностная схема на разнесенной по всем переменным сетке, дополненная специальным выражением, введение которого позволило уменьшить "схемную дисперсию". Эта модификация схемы основана на упомянутом выше методе исключения. Исследование устойчивости и дисперсионных свойств полученного конечно-разностного метода проведено в указанной работе лишь для одномерного случая.

5. Численные алгоритмы повышенного порядка аппроксимации

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

Один из таких методов получил широкое применение для решения НЛД-уравнений, слабонелинейный вариант которых, известный как модель Нвогу, впервые предложен в [22] и может быть записан в виде уравнений (17), (18) при г = (х, у) = -0.531Ъ(х, у). Как уже отмечалось выше, эта модель имеет расширенную область применимости благодаря повышенному порядку аппроксимации дисперсионного соотношения ее линейного аналога. Последовавшие позже обобщения модели Нвогу позволили получить полную нелинейно-дисперсионную модель [6] с учетом нестационарности донной поверхности [21].

В [33] для уравнений (17), (18) конечно-разностный метод использует прямоугольную сетку с целочисленными узлами. Построенная двухшаговая (двухстадийная) многослойная разностная схема имеет четвертый порядок аппроксимации по временному шагу. Такой же порядок обеспечивается подходящей аппроксимацией (в виде центральных разностей) пространственных производных первого порядка. Встречающиеся пространственные производные второго порядка (в дисперсионных членах) аппроксимированы на трехточечном шаблоне со вторым порядком. Такое понижение порядка предпринято для упрощения реализации разностной схемы, о чем будет сказано ниже.

Продемонстрируем конечно-разностный метод, переписав уравнения (17), (18) в следующем виде [33]:

щ = Е(г/, и, V), (62)

иг = ^ (п ,и, +

(63)

^ = С(г,,и, ^) + [С1(и)]4

(64)

где

и (и) = и + Ъ Ь1Ъихх + Ь2(Ъи)а

(65)

V(у) = V + Н ЬхкУуу + Ъ2(Ну)

УУ

Здесь

Е(г/, и, у) = -(Ни)х — (Ну)у — а^^х + ^у) + (Ни)хх + (Ну)Ху^

(66)

ахН (Ьуу + иху) + а2к ( (Ну)уу + (Ни)

)ху

Е(г/, и, у) = —дг]х — иих — уиу, 0(г], и, у) = —дщ — Vуу — иух,

(у) = —Н [Ь хНУху + Ъ2(Ну )ху ] , О (и) = —Н [Ь \НПху + Ъ2(Ни)ху ].

Выбор констант а\, а2,Ъ\ и Ь2 позволяет задавать разные модели. Так, при

сц = р2/2 — 1/6, а,2 = р + 1 /2, Ъх = р2/2, Ь2 = р, р = /Н = —0.531 имеем модель Нвогу, а при

= 0,

2

0, Ьх = 1/6, Ь2 = —1/2

— модель Перегрина. Возвращаясь к уравнениям (62)-(66), заметим, что выбор и и V в виде (65) и (66) обеспечит разностной схеме реализацию с помощью трехточечных прогонок, если применить для аппроксимации слагаемых в квадратных скобках (65) и (66) центральные разности второго порядка аппроксимации (для достижения более высокого порядка надо расширять шаблон).

Для решения уравнений (62)-(66) авторы [33] предлагают метод предиктор-корректор. Приведем здесь основные детали разностной схемы, использующей равномерную прямоугольную сетку с целочисленными узлами ] = (]\, ]2). Для записи сеточных функций , где / = (г/, и, у, и, V), применим мультииндекс, т.е. будем писать вместо

£П

3 31,32 .

В качестве "предиктора" применяется явная разностная схема Адамса — Бэшфорта третьего порядка аппроксимации

'Пп + ^ (23Еп — 16Еп-1 + 5Еп-2)

щ+1

ип + — (23Еп — 16 Еп-1 + 5Еп-2) + 2Е1 — ЗЕ^-1 + Е{

п-2

Уп + ^ (230п — 160п-1 + 50п-2) + 2Оп — ЗО-1 + О

п— 2

(67)

При этом все производные, входящие в Е, Е, О, Е1, 01, аппроксимированы с помощью центральных разностей. Явным счетом по схеме (67) получаем значения и^1 и , которые далее используются как известные левые части дискретных аналогов уравнений (65) и (66) относительно неизвестных ипп+1 и уп+1. Каждое из этих уравнений можно рассматривать как одномерное в соответствующем направлении, поэтому переменные ип+1 и у'п+1 находятся методом скалярных прогонок.

Таким образом, на этапе "предиктор" вычисляются "предварительные" значения всех гидродинамических величин на (п + 1)-м слое. Значения Е'п+1, Е'п+1, , (Е\)п+1, (01)п+1 вычисляются с помощью дискретных аналогов формул для Е, Е, О, Е1, 01 соответственно.

х

у

3

3

3

На втором шаге значения и™+1 и У-+1 пересчитываются заново с помощью

неявной разностной схемы Адамса — Мултона четвертого порядка аппроксимации:

г1 =

r]n + ^ (9Еп+1 + 19En - 5En-1 + En-2)

иn+1

Un + ^ (9Fn+1 + 19Fn - 5Fn-1 + Fn-2) + F^1 - F\

(68)

V^n+1 = j

Un + ^ (9Gn+1 + 19Gn - 5Gn-1 + Gn-2) + G^1 - Gn,

для решения которой относительно величин на (п + 1)-м слое применяется итерационный процесс, при этом начальным приближением для него выступают вычисленные на шаге "предиктор" значения соответствующих величин. Для расчета окончательных

значений и- и V- повторно решаются разностные аналоги уравнений (65), (66) с новыми значениями и""+1 и У-+1, полученными при пересчете (68). Этим завершается

n+1 Mn+1

и V.

n+1

на полном

процесс вычисления искомых гидродинамических величин r/j шаге разностной схемы.

В [6] этот численный метод распространен на случай полных НЛД-уравнений. Соответствующая система НЛД-уравнений в покомпонентном виде, который используют авторы работы, очень громоздкая, и здесь мы ее приводить не будем.

Описанный алгоритм распространен также на случай нестационарного дна [21]. На основе разработанного конечно-разностного метода создана компьютерная программа FUNWAVE (Delaware school), находящаяся в свободном доступе.

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

j

э

j

6. Численные алгоритмы для полных нелинейно-дисперсионных уравнений с усредненной скоростью

Для полных НЛД-моделей с усредненной скоростью очень мало работ, в которых был бы подробно изложен конечно-разностный метод решения 2В-уравнений. Можно упомянуть статью [5], где для модели Грина — Нагди приведена формальная аппроксимация, а также указать на способ, рассмотренный в предыдущем параграфе, позволяющий построить схему второго — четвертого порядка для общей базовой НЛД-модели, в том числе и для уравнений с усредненной скоростью, о чем упоминается в работах [6, 33].

Недостаток первых работ [3, 7], в которых для построения численного алгоритма решения системы полных 2В-НЛД-уравнений использовалось ее расщепление на гиперболическую и эллиптическую части, заключается в том, что полученная в результате расщепления система уравнений относительно полной глубины и усредненной скорости отличается от классических уравнений мелкой воды и имеет кратные собственные значения. Кроме того, расщепление системы уравнений в [7] выполнено только для неподвижного дна и плоского случая. Эти недостатки устранены в работе [37], в которой представлен численный алгоритм для моделирования процесса распространения длинных поверхностных волн в рамках полной нелинейно-дисперсионной модели, учитывающей сферичность Земли и ее вращение. Алгоритм основан на расщеплении исходной системы НЛД-уравнений и решении расширенной системы уравнений, состоящей из скалярного эллиптического уравнения для дисперсионной составляющей р давления р

вида (6) и гиперболической системы уравнений первого порядка, которые отличаются от уравнений классической модели мелкой воды лишь наличием в правой части уравнений движения дополнительных слагаемых, связанных с дисперсионными добавками. Алгоритм реализован в виде явной двухшаговой схемы предиктор-корректор, на каждом шаге которой поочередно решаются задачи, полученные в результате расщепления.

Поскольку структура уравнений полной нелинейно-дисперсионной модели сохраняется для слабодисперсионной модели на сфере, а также модели слабодисперсионных течений над слабодеформируемым дном [14, 38], то алгоритм численного решения полных НЛД-уравнений переносится без изменений и на эти приближенные модели. Такое расщепление оказалось плодотворным для конструирования численных алгоритмов решения задач как в плоской геометрии [39, 40], так и в сферической [41, 42].

Суть предлагаемого подхода и свойства схемы расщепления продемонстрируем на примере линейной системы уравнений (15). Расширенная система уравнений [39] получается из системы линейных уравнений (15), если последнюю дополнить уравнением относительно новой зависимой переменной

V = ^V • и, (69)

которая является линейным аналогом дисперсионной составляющей (6) давления в НЛД-модели. Из второго уравнения системы (15) получаем следствие

V • и + gV • Vr] = 1V • V р.

П0

Используя обозначение (69), приходим к расширенной системе

щ + hoV• и = 0, и + gVr] =-1 V р, А р - Р = с0Аг], (70)

До У

где V = Д0/3.

Предполагается, что на п-м слое по времени вектор скорости ип и свободная граница г]п известны, причем эти величины определены в целочисленных узлах х, где ] — мультииндекс, ] = ]0). Численный алгоритм решения расширенной системы уравнений (70) основывается на схеме предиктор-корректор, предназначенной для решения бездисперсионных уравнений мелкой воды, и состоит в поочередном решении задач для уравнений (70) как на шаге "предиктор", так и на шаге "корректор". Вначале решается уравнение для р:

рпх,з + рпу,3 ~ = Со(^пх,] + Щу,^) . (71)

Здесь г]ХПх а , Г]УПУ а — вторые разностные производные (см. ниже формулу (98)). Получен-

7п тП Хх,]1 Чуу,]

ные значения рп используются для вычисления предикторных величин г/*, и*, определенных в центрах ячеек:

^+1/0 ^+1/0 + г*. ... 1=0, (72)

•,^'+1/0 У,з+1/0,

и*+1/0 - иП+1/0 + = Д-VРп+1/о, (73)

До

*

где для сокращения записи использованы мультииндексы j + 1/2 = (j1 + 1/2, j2 + 1/2) и обозначения

ri'j+l/2 = 1 (v'j1,j2 + 'rfj1+l,j2 + rl1jl,j2 + l + rfjl + l,j2+l) ,

Uj+\/2 — Д U'iJ2 + U'l + lj2 + Ujl,j2 + 1 + U'l + lj2 + l

Un = 1 (nn + Un — Un — Un

UX,j+1/2 — 2Ax\jl + 1,j2 + 1 + ajl + 1,j2 ajl,jl + 1 Ujl,j2 tП = 1 ftП + tП — tП — tП

y, j+1/2 — 2Ay V jl+1,j2+1 + v jl,n+1 Vjl+1^ V

31,32

V$+1/2 = (^З+Ф' Чз+J , ^^+1/2 = (ih+W -

В уравнениях (72), (73) т* = т(1 + в™+1/2)/2, &f+1/2 — схемный параметр, надлежащий выбор которого гарантирует выполнение TVD-свойства для численного решения модельных скалярных уравнений [43, 44]. При теоретическом исследовании рассматриваемой схемы будем предполагать, что в = const > 0.

Кроме уравнений (72), (73) на шаге "предиктор" решается и уравнение относительно сеточной функции p**+i/2, определенной в центрах ячеек, при этом используются вычисленные ранее значения г/*+1/2:

*

* I * _ Vj+1/2 _ 2 ( * I * А

Vxx,j+1/2 + Vyy,j+1/2 у — Vxx,j+1/2 + Vyy,j+1/2j .

Vj+1/2 _ j2 ,

Полученные значения г/*, U*, р* используются на шаге "корректор":

(74)

V-+1- Vj

+ ho(u* +v* ) —0, (75)

V x,j y,jJ

■.■.П+1 n 1

U • — U, о 1 о

j j + gv v! — i- v р*, (76)

т * h

Q

где

1

A 1

* * * * *

UXj — ~2Ax\jl + 1/2,h + 1/2 - Ujl-1/2,h+1/2 + Ujl + 1/2,h-1/2 - Ujl-1/2,n-1/2 *

* —

(* * I * * \ __Vjl+1/2,h+1/2 - Vh + 1/2,j2-1/2 + Vh-1/2,h + 1/2 - Vj!-1/2,h-1/2 ) ,

VV** = (V* ., V* ) , VP* = (<P*° ., <P*) .

J \ x,3 y,]J J \ x,3 y,]J

Необходимое условие устойчивости схемы (71)-(76) при малых значениях параметра в можно записать (см. Приложение Б) в виде ограничения на число Куранта

min (Ах, Ау) /Ао + 25у/в/3

Соя < — , (77)

где использованы обозначения (49). Условие (77) является менее ограничительным по сравнению с условием устойчивости

Со. (78)

схемы предиктор-корректор (схемы (72), (73), (75), (76) с нулевыми правыми частями в разностных уравнениях движения (73), (76)) для бездисперсионных уравнений мелкой воды (уравнений (70) с р = 0).

На рис. 2 изображено распределение значений модуля множителя перехода |р(£,С)| схемы предиктор-корректор (71)-(76) для следующего набора параметров:

В = 0; с0ж 1 = 1/2, с0ж2 = 1, ¿1 = 1, ¿2 = 2. (79)

Видно, что низкочастотные гармоники подавляются слабо, однако и высокочастотные гармоники диссипируют незначительно, особенно в окрестности угловой точки (£, £) = (ж, ж), что отличает рассматриваемую 2В-схему от одномерной схемы предиктор-корректор, при использовании которой высокочастотные гармоники подавляются в гораздо большей степени, чем низкочастотные.

Рис. 2. Модуль множителя перехода схемы предиктор-корректор для параметров (79)

а б

ф 1.2

0.8-

0.4-

0.0

— - --1 1

-■3 -4 \ X 1 1 1 1 |

1 1 1 1

¡V *Ч | 1 1 1 1

1 Ч\ 1 ........... 1 1 1 1 1 1 .'.-"г-!-Г . ,

Рис. 3. Дисперсионные свойства гидродинамических моделей и разностных схем: а — графики изменения фазы для операторов перехода моделей потенциальных течений (1) и Перегрина (2), для множителей перехода схем Ригга (3) и предиктор-корректор (4); б — сечения поверхностей изменения фазы диагональной плоскостью £ = (

График изменения фазы (см. ниже формулу (135)) множителя перехода схемы предиктор-корректор (71)-(76) при значениях параметров (79) приведен на рис. 3, а в виде поверхности 4, а на рис. 3, б с горизонтальной осью 5 = \/2£ — в виде кривой 4. Для сравнения приведен также график 3 аналогичной зависимости (54), характеризующей дисперсионные свойства схемы Ригга. Видно, что при одинаковом втором порядке аппроксимации схема предиктор-корректор грубее, чем схема Ригга, воспроизводит изменение фазы (51) оператора перехода линеаризованной модели Перегрина. Оно показано на рис. 3 в виде графика 2, хотя на длинных волнах получаются сравнимые приближения. На рис. 3 также изображен "эталонный" график 1 изменения фазы за один шаг по времени

ФРР = ^, (80)

в линеаризованной модели потенциальных течений. Из рис. 3, б видно, что схема предиктор-корректор удовлетворительно описывает дисперсионные свойства НЛД-модели для значений з < 0.7. Учитывая (49) и справедливую при £ = £ цепочку равенств

5 = = /2, = /2| к|До = /2^ До = /2^ = /2^

А 0 А 6 А '

получаем, что для указанных в (79) значений параметров неравенство в < 0.7 будет выполняться при А > 5.7к0, т. е. дисперсионные свойства достаточно длинных волн будут воспроизводиться 2В-схемой предиктор-корректор с удовлетворительной точностью.

Заключение

Дан обзор конечно-разностных схем для решения двумерных НЛД-уравнений. Выявлено, что практически отсутствуют работы, посвященные теоретическому исследованию применяемых разностных методов. Мы взяли наиболее популярные схемы и для них выполнили исследование устойчивости и дисперсионных свойств. И хотя для такого анализа использован простой математический аппарат, в статье приведены подробные выкладки для того, чтобы облегчить в будущем исследовательскую работу по изучению свойств аналогичных конечно-разностных методов.

Один из главных выводов работы состоит в том, что условия устойчивости для разностных схем, аппроксимирующих 2В-уравнения с дисперсией, являются более слабыми по сравнению с аналогичными условиями для схем, аппроксимирующих бездисперсионные 2В-уравнения мелкой воды (для наглядности можно сравнить формулу (77) с (78)). Таким образом, это свойство, ранее обнаруженное для одномерных схем в [1], наследуется при переходе к двумерным. Вид условий устойчивости двумерных схем также оказался аналогичным одномерному случаю (см. для сравнения, например, формулы (77) и (116)). Однако такое наследование не выполняется по свойству дисперсии, поскольку в Ш-случае, как правило, существует значение числа Куранта (см., например, формулу (121)), при котором фазовая ошибка схемы становится минимальной, а для рассмотренных двумерных схем подобное свойство не обнаружено (это видно из сравнения формул (120) и (136)).

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

Сегодня не вызывает сомнения, что численное моделирование на основе уравнений мелкой воды является эффективным инструментом решения задач волновой гидродинамики, имеющих практические приложения. При решении крупномасштабных задач (обширные акватории, значительные времена распространения волн) представляется разумным использовать стратегию мультиподобластей [45]. При этом важно при построении методов расчета придерживаться иерархического подхода, позволяющего в рамках заданных масштабов движения отслеживать выполнение основных балансовых соотношений, не нарушая при этом структуры и основных свойств применяемых методов. Для расчетов вплоть до береговой линии перспективной может стать комбинация НЛД-модели и модели классической нелинейной мелкой воды [28].

Приложение А. Об условии устойчивости и дисперсионных свойствах схемы Ригга

Необходимое условие устойчивости схемы (39)-(42) эквивалентно условию К < 2 неположительности дискриминанта уравнения (46), где величина К задана формулой (47). Неравенство К < 2 эквивалентно выполнению условия

^(г, в) < 0 У(г,в) е я = [0,1] х [0,1] (81)

для функции двух переменных

Р(г,.) = с0 (ж?г + - 3- 3- 1 = (с0®2 - 452)(г ^ + . Д 2) - 1,

где

г = яп2 8 = вт2 ^ (82)

2 2 к ;

и использованы обозначения (49). В силу линейности функции Р выполнение условия (81) на всем квадрате Я эквивалентно выполнению этого условия только в его вершинах, а это с учетом того, что Р(0,0) = -1 < 0, эквивалентно одновременному выполнению следующих трех неравенств:

^ ^ 1) = - 3 ¿2) Д02 -1 < o,

¿Ж2- 4Л-Д^- 1<0, (83)

* (1,0)=(с>2 - 3 #

= (С0Ж2 - 3^ Д^ - 1 < 0, (84)

^(1,1)= ^ж2 - 4+ ) - 1 < 0. (85)

3 / \Дж2 Ду2

Из определения (49) среднего шага Д0 пространственной сетки следует, что

/ Д0 \2 / Д0 \2

(Д) + (Д) = 1 (86)

поэтому неравенство (85) эквивалентно условию

1

С0Ж <\\ 1 + 4¿2. (87)

При выполнении неравенства (87) ограничения (83), (84) выполняются автоматически, поэтому неравенство (87) является необходимым условием устойчивости двумерной схемы (39)-(42) на разнесенной сетке.

Из равенства (50) следует, что изменение фазы 0Rg = arg р множителя перехода (53) разностной схемы (39)-(42) определяется из выражения

cos ф^ = 1 — R

или

■ fag sin ■

2

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

i ГФ г

__I I у + § I

ф^ = 2arcsi^v R/2 = 2arcsin I c0W——-| . (88)

\ V 1 + Ol + »2

Для длинных волн (малых £, £) имеют место выражения

ь = |е2 + о(е4), ъ2 = |с2 + о(с4), л * . =1 -^ + о(е4 + с4), (89)

3 3 у 1 + Ьх + &2 6

в записи которых использованы обозначения (49), (52) и равенство

¿2 + ¿2 = ¿2. (90)

Следовательно,

ф^ = 2 arcsin

2—1+o(s5))2+я2( 2—48+o(i 5))2 О—?i2+о(<4))

Здесь при записи остаточных членов 0(£4 + £4) и 0(£5), 0((5) используются степени переменной (52), что возможно, если учесть равенства £ = 0(<г) и £ = 0(я), которые являются следствием естественного предположения о том, что шаги сетки Дж, Ду имеют одинаковый порядок малости, т. е. при измельчении сетки выполняются неравенства

тДу < Дж < МДу (91)

с некоторыми положительными константами т, М (0 < т < М < то). Из (91) следует, что

т М т2 М2

-, Дж < До < -, Дж, -, Ду < До < -, Ду,

М^ 1 + т2 ~ ~ ту/1 + М2 Му/ 1 + т2 ~ ту/1 + М2

поэтому

f2 = А-?Д;2 < | к|2Д;;2 < (M^^)2| к|2Д2 = M2(l + -),2,

<2 = £2ДУ2 < | к|2Ду2 < ()2| к|2Д2 = M2(- + -V2

- 2 - - 2

т.е. f = O(0,C = O(0.

С учетом равенства

ж Я

(92)

запишем

Фщ

2 агсвт

со

У

ж2 я2

48

2

+ 0(я6) (1 — -+ 0(^4))

2 агсвт

ж? Г ж2е + ж2с4,

^^У1 — 12ж2я2 +

0(^4) (1 — у^2 + 0(я4))

2 агсвт = 2 агсвт

= СоЖЯ

- ж (1—Р) (1

ж2е4 + ж2С4 24а? я2

ж12 4 + ж22 4

2

ож ~2

—

ож 2 12

-Я — Со-

ож 2 "6"

ж2 ¿4 + Ж2 (~4 „3 ж3

3 + ж2<, Сож з

— Со~

24ж

48 жя

+

24

+ 0(я4)) + 0(я4) я3 + 0(я4).

Сравнивая полученное выражение с формулой (57) для дисперсии модели (38), видим, что фазовая ошибка схемы (39)-(42) определяется формулой

ж2 ¿4 + ж2 /-4 „3 ж3

Аф = Со ^ + — с<°ж-я3 + 0(я4) = 0(я3),

24ж

24

(93)

из которой следует, что в общем случае "физическая" дисперсия имеет тот же порядок по я, что и фазовая ошибка. Однако для некоторых волновых векторов к фазовая ошибка будет меньше, чем в общем случае. Например, для к = ( к1, 0) в силу равенства (92)

имеем

ДФ ж?£4 Аф = со-

24ж1

33

с0ж1_^3 + 0(^4) _ с0ж1

24

24

(1 — с2ж2)е3 + 0(я4),

поэтому при условии с0ж1 = 1 выполняется равенство Дф = 0(<^4). Это равенство выполняется и для волнового вектора к = (0, к2), если положить с0ж2 = 1.

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

ЬДх = к2Д у. (94)

Тогда = = и в силу равенства

2 2 2

Ж _L СУ"» - СУ"»

1 + ж2 = ж

(95)

выражение (93) принимает вид

Дф--

ож ~24

(1 — с^ж2) + 0(я4),

поэтому при условии

ож = 1 (96)

фазовая ошибка как минимум на один порядок по будет меньше "физической" дисперсии.

4

Приложение Б. Исследование устойчивости и дисперсии схемы предиктор-корректор

Здесь выписано необходимое условие устойчивости и приведено выражение для фазовой ошибки схемы предиктор-корректор (71)-(76). Для упрощения изложения вначале рассмотрен одномерный случай, для которого приведен подробный вывод уточненного (по сравнению с [1]) условия устойчивости.

Схема с расщеплением для линеаризованной системы НЛД-уравнений в одномерном случае

В одномерном случае расширенная система уравнений (70) принимает вид

Px Р 2

Vt + hoUx = 0, ut + g Г)х = —, pxx--= с0 rqxx, (97)

ho v

а схема предиктор-корректор (71)-(76) для системы (97) выглядит следующим образом:

P?+l - Щ + Р]-1 PL = г2 ~ = 22 -2VL + nu (98)

Ах2 v = C°Vxx'i - С° Ах2 , (98)

у*+1/2- от+1 + V?) + hUmzUL = * (99)

^ +ho Ах =0 (99)

u*+1/2 - °.4ul+1 + uL + VL+I - v1? = 1 PL+1 - PL (100)

T* 9 Ах ho Ах , ()

Pj+3/2 Щ+1/2 + Pj-l/2 P*j+l/2 = c2 Vj+3/2 2Vj+l/2 + V*j-l/2 А х2 - o А х2

vT1 - 'П? + h Uj+l/2 - Uj-l/2 =0

--+ h0-7- = 0,

А х

UT1 -U? + V'j+1/2 - Vj-l/2 _ 1 Pj+l/2 - Pj-l/2

+ 9

:ioi)

:io2)

т Ах к0 Ах

Для исследования свойств разностной схемы (98)-(102) будем использовать гармо-

ники

Г \ ( Ео\ / V* \ iE*

Ио I рпег*; I и* I = I и** | рпег(з+1/2)^ 3 \фо] \р*; 3+1/2 Из уравнений (98) и (99) следует, что

Ф=Е 4С /3 (103)

Фо = Е 1 + 45 2т/3, (103)

Е* = Е со^ - Шокоф + в) яп | , (104)

где 8 = к0/Ах, ж = т/Ах, г = 8т2(£/2). С учетом равенства (103) из уравнения (100) получаем

и* = — Е0 9ж(1+в} вт^ + Щ соз (105)

0 01 + 48 т/3 2 2' К 1

а из (101) следует выражение

=Е* 1г4Йз • (106)

Подставляя гармоники в уравнения шага "корректор" (102) и учитывая выражения (104)—(106), получаем следующие уравнения относительно амплитуд Е0 и U0:

е {? -1+If^f )+uo о«*. sin e)=o,

1 + 45 2 г/3

Л gasin^ \ ,гг/ 1 , 2с0ж2(1 + ^)г\ Ч4 + 452г/J + ЧР - 1+ 1 + 452г/3 J = 0

Следовательно, множитель перехода является решением квадратного уравнения

_ 1 + 2с2а2(1 + А)г\2 + с2a2 sin2^ =0 Р + 1 + 452г/3 у + 1 + 452г/3 0' корнями которого являются комплексно-сопряженные числа

2 с2а2(1 + А)г . Coa sinC Пп7ч

Р1^2 = 1 - 1 + 452г/3 ± \;1 . ,.2 , . з • (107)

Условие устойчивости | р12| < 1 можно записать в виде неравенства

/(г) > 0 Уге [0,1], (108)

где введена квадратичная функция переменной г:

/(г) = 452г2 - (с>2(1 + 0)2 - 1 - 452¿A г + В. (109)

Из условия (108) с необходимостью следует выполнение неравенств /(0) > 0 и /(1) > 0 или

В > 0, (110)

<с' = Г-++г (111»

Графиком квадратичной функции (109) для г Е (-то, то) является парабола с ветвями, направленными вверх, и вершиной, имеющей абсциссу

= с2Ж2(1 + А)2 - 1 - 4520/3 Г* = 852/3 , (112)

при этом

4

/(Г*) = 0 - 452г2.

Рассмотрим все возможные случаи расположения вершины параболы. Пусть 0 < г* < 1. Тогда совместное выполнение неравенств 0 < г* < 1 и /(г*) > 0 эквивалентно условию

0 < г* < min

-Г

"T^J

или

С2 < с0ж < Сэ при 0 < 452/3, С2 < с0ж < С4 при 0 > 452/3,

:ii3)

где С2 < Сэ, С2 < С4,

а/1 + 4 5 20/3

С2 =

1 + 0

Сэ =

1 + 2 5 ^0/3

1 +

С4 =

а/1 + 4 5 2(0 + 2)/3

1 + 0

Если г* < 0 (т.е. СоЖ < С2), то (108) эквивалентно условию /(0) > 0 или ограничению (110), которое предполагается выполненным. В этом случае подходят все числа Куранта, удовлетворяющие условию

с0ж < с2.

Г114)

Наконец, в последнем из возможных случаев г* > 1 (т. е. СоЖ > С4) выполнение условия (108) на отрезке [0,1] эквивалентно неравенству /(1) > 0 или ограничению (111), поэтому решением будет множество

С4 < СоЖ < Ci,

Г115)

при этом 0 > 4 52/3.

Объединяя множества решений (113)—(115), получаем следующее, зависящее от параметров 8 и 0 необходимое условие устойчивости схемы предиктор-корректор:

с0ж < С(5, 0),

П16)

где

С( , )

С( , )

1 + 2 5

1 +

1 + 4 5 2/3 1 + 0

4

при 9 < 4 52,

4

при 9 > - 52.

П17)

:и8)

Графики зависимостей (117), (118) изображены на рисунке. По ним можно судить об области устойчивости схемы предиктор-корректор. Поскольку с практической точки зрения интерес представляет только случай малых значений параметра , то можно ограничиться рассмотрением лишь первой зависимости (117). Для нее С(5, 0) > 1 при 0 < 452/3, поэтому условие устойчивости допускает проведение расчетов с числами Куранта, большими единицы. Интересно отметить, что с увеличением значений параметра верхняя граница допустимых чисел Куранта также растет.

На рисунке изображена также верхняя граница области устойчивости схемы предиктор-корректор для бездисперсионной модели мелкой воды (уравнения (97) с р = 0). В этой схеме отсутствуют этапы (98), (101) вычисления дисперсионной добавки к давлению, а в разностных уравнениях движения (100), (102) правые части полагаются равными нулю. Условие устойчивости такой схемы имеет вид

СоЖ <

1

VT + F

:119)

V................................................... \ \ .............Ч.......i................... - 1 -------- 3 ------- 4

▼......................... ............................................. rs ^^ ----_

, , , , <<<<<> -Да-- - —

К условию устойчивости схемы предиктор-корректор: графики зависимостей при 5 =1, 2 (117) — кривые 1, 3 соответственно, (118) — кривые 2, 4 соответственно; кривая 5 — верхняя граница области устойчивости (119) схемы предиктор-корректор для бездисперсионной

модели мелкой воды

откуда следует, что при 9 > 0 допустимые числа Куранта строго меньше единицы. Легко проверить, что при любых в > 0 имеет место неравенство

С(6,9) >

1

vr+e'

означающее, что условие устойчивости схемы предиктор-корректор для бездисперсионных уравнений мелкой воды является более жестким, чем для уравнений, учитывающих дисперсию волн. При в = 0 эти условия устойчивости совпадают:

С (6,в)

=1=-тт=

0=0 Vl+e

0=0

Из формулы (107) следует, что для длинных волн изменение фазы фрс = argp множителя перехода разностной схемы предиктор-корректор (98)—(102) определяется выражением

фрс = ±

Со -

С0 ж<52 3 с0 ж ( 2 2

^ + с2 ж2 (3 0 + 1) -1) е3 + о(е4)

(120)

поэтому при условии

с0 ж =

^Т+30

(121)

фазовая ошибка имеет порядок 0() и "физическая" дисперсия преобладает над численной.

6

1

Об устойчивости и дисперсионных свойствах 2Ю-схемы предиктор-корректор

Для исследования устойчивости разностной схемы (71)-(76) подставим в разностные уравнения гармоники

( пп\

и

\рП /

( Ео\

Uo V*

\ф0 J

( V* \

рПeí(jiC+hO,

и*

Из уравнений (71) и (72) следует, что

Фо =

V Ч? ) +1¡2

Е0£, h + h

( Щ \

Uo* Vo*

\Ф*о J

рп е%((ji+l/2)a+{j2 + 1/2)0 .

Щ

o

'1 + bi + b2

E0 cos - cos- — ih0(1 + вц U0^1 sin - cos - + V0w2 cos - sin -

:122)

:123)

2 2 04 ' \ 0 1 2 2 0 2 2 2/

где использованы обозначения (43). С учетом равенства (122) из уравнений (73) и (74) получаем

Т* _ „• И< И""1^ 1 " / ^ „„„ Ъ | ТТ „„„ „„„ (124)

U* = — %E0---— sin - cos —+ U0 cos - cos-,

0 01 + h + b2 2 2 0 2 2'

T/* 9^2(1 + 0) £ . с T/ e С V* = —* Eo i + bl + b2 cos 2sin2+ Vo cos 2 cos 2,

Ф*

E*

2

1+ 2

:125)

(126)

1 + Ъг + Ъ2~

Подставляя гармоники в уравнения шага "корректор" (75), (76) и учитывая выражения (123)-(126), получаем следующие уравнения относительно амплитуд Е0, и, Ц):

(р — 1 + ai + ü2)E0 + ihodiUo + ihod2Vo = 0, igdi

1+ 1+ 2 ígd2

1 + bi + b2

Щ + (p — 1 + ai)U0 + y/aia2 V0 = 0, Щ + /aa Uq + (p — 1 + a2)V0 = 0,

П27)

где

ai =

2cWi(1+e) . 2 С 2 С 2c20^2(1+e) . 2 С 2 С

-1-;- sin - cos - , a2 = -2-;- sin - cos -,

1 + h + b2 2 2 2 1 + h + b2 2 2'

d1 = ж1 sinocos2 y d2 = ж2 sinocos2

Для существования нетривиального решения системы уравнений (127) ее определитель должен обращаться в нуль:

- 1 + i + a2

igdi 1 + bi + Ь2

igd2 1 + bi + Ь2

ih0di ih0d2 p — 1 + ai /a ia2

/aia2 p — 1 + a2

*

п

3

Отсюда для множителя перехода р получаем кубическое уравнение

(Р - 1)

¿2 + ¿2

(р - 1 + аг + а2)2 + с0 1 2

1 + Ьг + Ь2

0, (128)

один из корней которого равен единице, а два других являются комплексно-сопряженными числами:

1

Д2 + л2

Р1.2 = 1 - (аг + 02) ± гон/1 % + \ . (129)

1 + Ьг + &2

Условие |р1)2| < 1 эквивалентно неравенству

й 2 + й 2

1 - 2( а г + 02) + (аг + 02)2 + с^ Д < 1,

1 + Ьг + 62

или

('1|Л\|С0(1+^) ( 2 • 2 £ 2 С , 2 • 2 С 2 , 2 £ 2 С ^ п

— (1 + 6') +--;-Г" «2 81П' - СОв2--+ «2 В1П - сое - + сое- сое - < 0,

1 ; 1 + Ьг + Ь2\ г 2 2 2 2 2) 2 2

или

^(Г, в) < 0 У(г, в) еЯ = [0,1] х [0,1], (130)

где

^(г, в) = -(1 + 0) (1 + 4^г + 4&) + с2(1 + 0)2 («?г(1 - в) + «2«(1 - г)) +

/ 4 4 \

+(1 - 0(1 - я)(1 + 4+ 3¿2 в )

4 4

3:2г + 4 ¿2 ^), (131)

а величины г, в определены формулами (82).

Далее будем рассматривать только случай малых значений параметра 0, удовлетворяющих условию из (117) по обоим направлениям х и у:

4

0< 4Ш1П(¿2; . (132)

Для таких 0 шаг по времени будем выбирать исходя из аналога условия устойчивости (116) одномерной схемы

/1 + 2^073 д 1 + 2 52 У0/3 . \

с0т < ш1п I --Дх; --Ау I . (133)

Тогда для функции (131) имеем оценки

^ (Г, 8) <-(1 + 3 ^ +4 Ф )(Г + Я-Г 5 + 0)+ (1 + ^^ )г (1 - 0+ (1 + ^^ (1 - г) =

3 г 3 2 Г ' V у/3 ) у ' \ у/3

( 44 \ / 44 \

- (1 + 3 ^2Г + 3 ^) (г +8 -г 8 ) -в (1 + 3 ^2Г + 3 ^) + т + в - 2 г 5+

П"ТТ + -^)г(1 - 5 ) + ("ТТ + Чп5(1 - г) =

— (з ^2г + 4 ^2S) (r +8 —r s) —r s —

4 8iV

4 8 ff? 482V

, _ r(1 — s)--гs + _

v^ ( ) 3 у/3

(V )2 + 2Ve (V r (1 — s ) + V s(1 — o) — 4(6ir (! — s ) + 5*s (! — r)) 2+

s(l — r) —

4 2

+ 4 (Si r(1 — s) + S2S(1 — rj) 2 — (3tfr + 4Öfs) (r + s — rs) — (l + U(S2 + ¿2))rS <

4

< -

3

Ö1r (1 — S )+ Ö2 S (1 — r)}

+

4

4

8

4

+-81 (гв — г — 1)г8 + -8%(гв — в — 1)г8 + -8\82Гз(1 — г)(1 — в ) < ~тгС(г, в), 3 3 3 3

где билинейная функция С определена формулой

г, з) = 5^(гз — Г — 1) + 5^(гз — 8 — 1) +28г82 (1 — г)(1 — 8).

Поскольку на единичном квадрате Q функция С — неположительна, то неравенство (130) выполняется, т.е. при условиях (132), (133) выполняется необходимое условие |Р\,21 < 1 устойчивости рассматриваемой схемы предиктор-корректор.

Отметим, что формулу (133) можно переписать в виде 2В-аналога условия устойчивости (116) одномерной схемы:

с2ж <

min (Ах/Ао, Ау/Ао) + 28у/в/3

1+

:134)

Для множителя перехода (129) схемы предиктор-корректор, взятого со знаком "+", получаем следующую формулу для определения изменения фазы за один шаг по времени:

фрс = arg р = 2 arcsin [ —Д= л/\р[—Г+а1+~

л/2Й

2

:i35)

Для длинных волн с учетом разложений (89) и условия (91) получаются следующие выражения:

1 + Ь-1 + = 1 + + 0(?4),

3

аг = 2с%а2(1 + в)( ^-т — ^: —

" 0ж2(

+ Ü2 = с20(1 + в)

е4 е2с2 s2e

48 16 12

С4 £2С2 82(2

48 16 12

■■ + ж2С ж2ее

24

я2 + 0(я6)

я2 + 0(я6) ж2 ¿2

я4 + 0(я6)

6

d? + d.2

1 + d2 _ „2„22 + Ж2(4 С2 Ж?82 ^ + o(^6)

1 + bi +b2 2 - 3-е

= ж я —

3

И = 1 — f( ^V +

1

23 ж282я4

44

ж4 я4

2/^21 v +*2<4)+V" ж2ее—о —£2(1+20) —+

4

0( я6)).

1

2

8

Поэтому

фрс = 2 arcsin

1 - 6 ж2?2 - --^ + (1 + 4^)—g— +

Ж2М + ж2л4 ¿2^2 г3ж3

^с^" 3 + Ж2Ч S S , /Л , оП\ЬСЖ 3 | гл< 4\

ссжя---— С - сс----Ссж ---+ (1 + Щ^—я3 + 0(я4).

с - „- „„^ + (1 + 3в)

6 6«<^ 4<^ 6

Сравнивая полученное выражение для фрс с формулой (57) для Фкьб, видим, что фазовая ошибка схемы (71)-(76) определяется как

«2М + «2М ¿2^2 г3«3

Аф = со«г4 + «2С + со«- (1 + 30)+ 0(я4) = О(0), (136) 6« 4 6

т. е. имеет тот же порядок малости, что и "физическая" дисперсия модели. Только для волновых векторов к = (кг, 0) и к = (0, к2) она может быть меньше. Например, для к = ( кг, 0) имеем

Аф = ^ [1 - с0«?(1 + 3 0)1 £3 + 0(<т4),

6

поэтому при выполнении условия

С0Ж1 = , , (137)

0 1 л/Т+30 v 7

аналогичного (121), имеем равенство А ф = 0(<^4)

Благодарности. Исследование выполнено при финансовой поддержке РНФ (проект № 14-17-00219) и гранта Президента РФ для государственной поддержки ведущей научной школы РФ (№ НШ-7214.2016.9).

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

[1] Федотова З.И., Хакимзянов Г.С., Гусев О.И. История развития и анализ численных методов решения нелинейно-дисперсионных уравнений гидродинамики. I. Одномерные модели // Вычисл. технологии. 2015. Т. 20, № 5. C. 120-156.

Fedotova, Z.I., Khakimzyanov, G.S., Gusev, O.I. History of the development and analysis of numerical methods for solving nonlinear dispersive equations of hydrodynamics. I. One-dimensional models problems // Comput. Technologies. 2015. Vol. 20, No. 5. P. 120-156. (In Russ.)

[2] Peregrine, D.H. Long waves on a beach // J. of Fluid Mech. 1967. Vol. 27. P. 815-827.

[3] Барахнин В.Б., Хакимзянов Г.С. Численное моделирование косого наката уединенной волны // Прикладная механика и техническая физика. 1999. Т. 40, № 6. С. 17-25. Barakhnin, V.B., Khakimzyanov, G.S. Numerical simulation of an obliquely incident solitary wave // J. of Appl. Mech. and Techn. Physics. 1999. Vol. 40, No. 6. P. 1008-1015.

[4] Abbott, M.B., Petersen, H.M., Skovgaard, O. On the numerical modelling of short waves in shallow water // J. of Hydraulic Res. 1978. Vol. 16, No. 3. P. 173-204.

[5] Компаниец Л.А. О численных алгоритмах для нелинейно-дисперсионных моделей мелкой воды в двумерном случае // Вычисл. технологии. 1996. Т. 1, № 3. C. 44-56. Kompaniets, L.A. On numerical algorithms for nonlinear dispersive shallow water models in two-dimensional case // Comput. Technologies. 1996. Vol. 1, No. 3. P. 44-56. (In Russ.)

[6] Wei, G., Kirby, J.T., Grilli, S.T., Subramanya, R. A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves // J. of Fluid Mech. 1995. Vol. 294. P. 71-92.

[7] Хакимзянов Г.С., Шокин Ю.И., Барахнин В.Б., Шокина Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001. 394 с.

Khakimzyanov, G.S., Shokin, Yu.I., Barakhnin, V.B., Shokina, N.Yu. Numerical simulation of fluid flows with surface waves. Novosibirsk: Izd-vo SO RAN, 2001. 394 p. (In Russ.)

[8] Brocchini, M. A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics // Proc. of Royal Soc. of London. A. 2013. Vol. 469(2160):20130496.

[9] Федотова З.И., Хакимзянов Г.С. Базовая нелинейно-дисперсионная модель гидродинамики длинных поверхностных волн // Вычисл. технологии. 2014. Т. 19, № 6. С. 77-94. Fedotova, Z.I., Khakimzyanov, G.S. The basic nonlinear-dispersive hydrodynamic model of long surface waves // Comput. Technologies. 2014. Vol. 19, No. 6. P. 77-94. (In Russ.)

[10] Шокин Ю.И., Федотова З.И., Хакимзянов Г.С. Иерархия моделей гидродинамики длинных поверхностных волн // ДАН. 2015. Т. 462, № 2. С. 168-172.

Shokin, Yu.I., Fedotova, Z.I., Khakimzyanov, G.S. Hierarchy of nonlinear models of the hydrodynamics of long surface waves // Doklady Physics. 2015. Vol. 60, No. 5. P. 224-228.

[11] Ertekin, R.C., Webster, W.C., Wehausen, J.V. Waves caused by a moving disturbance in a shallow channel of finite width // J. of Fluid Mech. 1986. Vol. 169. P. 275-292.

[12] Железняк М.И., Пелиновский Е.Н. Физико-математические модели наката цунами на берег // Накат цунами на берег: Сб. науч. тр. Горький: ИПФ АН СССР, 1985. С. 8-33. Zheleznyak, M.I., Pelinovskiy, E.N. Physical and mathematical models of the tsunami run upon a beach // Tsunami Climbing a Beach: Sb. nauch. tr. Gor'kiy: IPF AN SSSR, 1985. P. 8-33. (In Russ.)

[13] Федотова З.И., Хакимзянов Г.С. Нелинейно-дисперсионные уравнения мелкой воды на нестационарном дне // Вычисл. технологии. 2008. Т. 13, № 4. C. 114-126. Fedotova, Z.I., Khakimzyanov, G.S. Nonlinear dispersive shallow water equations for a non-stationary bottom // Comput. Technologies. 2008. Vol. 13, No. 4. P. 114-126. (In Russ.)

[14] Fedotova, Z.I., Khakimzyanov, G.S., Dutykh, D. On the energy equation of approximate models in the long-wave hydrodynamics // Russ. J. of Numer. Anal. and Math. Modelling. 2014. Vol. 29, No. 3. P. 167-178.

[15] Fedotova, Z.I. On properties of two simplified nonlinearly dispersive models of long wave hydrodynamics // Intern. J. of Comput. Fluid Dynamics. 1998. Vol. 10, No. 2. P. 159-171.

[16] Kim, K.Y., Reid, R.O., Whitaker, R.E. On an open radiational boundary condition for weakly dispersive tsunami waves // J. of Comput. Physics. 1988. Vol. 76. P. 327-348.

[17] Хабахпашев Г.А. Нелинейное эволюционное уравнение для достаточно длинных двумерных волн на свободной поверхности вязкой жидкости // Вычисл. технологии. 1997. Т. 2, № 2. C. 94-101.

Khabakhpashev, G.A. Nonlinear evolution equation for moderately long two-dimensional waves on a free surface of viscous liquid // Comput. Technologies. 1997. Vol. 2, No. 2. P. 94-101. (In Russ.)

[18] Литвиненко А.А., Хабахпашев Г.А. Численное моделирование нелинейных достаточно длинных двумерных волн на воде в бассейнах с пологим дном // Вычисл. технологии. 1999. Т. 4, № 3. C. 95-105.

Litvinenko, A.A., Khabakhpashev, G.A. Numerical modeling of sufficiently long nonlinear two-dimensional waves on water surface in a basin with a gently sloping bottom // Comput. Technologies. 1999. Vol. 4, No. 3. P. 95-105. (In Russ.)

[19] Yoon, S.B., Lim, C.H., Choi, J. Dispersion-correction finite difference model for simulation of transoceanic tsunamis // Terrestrial, Atmospheric and Oceanic Sci. 2007. Vol. 18, No. 1. P. 31-53.

[20] Liu, P.L.-F. Model equations for wave propagations from deep to shallow water // Advances in Coastal and Ocean Engineering / P.L.-F. Liu (Ed.). Singapore: World Sci. Publ. Co., 1994. Vol. 1. P. 125-157.

[21] Lynett, P.J., Liu, P.L.-F. A numerical study of submarine-landslide-generated waves and run-up // Proc. of Royal Soc. of London. A. 2002. Vol. 458. P. 2885-2910.

[22] Nwogu, O. Alternative form of Boussinesq equations for nearshore wave propagation // J. of Waterway, Port, Coastal, and Ocean Eng. 1993. Vol. 119, No. 6. P. 618-638.

[23] Abbott, M., Petersen, H., Skovgaard, O. Computations of short waves in shallow water // Coastal Eng. Proc. 1978. Vol. 16. P. 414-433.

[24] Abbott, M.B., Ionescu, F. On the numerical computation of nearly-horizontal flows // J. of Hydraulic Res. 1967. Vol. 5. P. 97-117.

[25] Novak, P., Guinot, V., Jeffrey, A., Reeve, D.E. Hydraulic modelling — an introduction. Principles, methods and application. London, New York: Spoon Press, 2010. 602 p.

[26] Abbott, M.B., McCowan, A.D., Warren, I.R. Accuracy of short-wave numerical models // J. of Hydraulic Eng. 1984. Vol. 110, No. 10. P. 1287-1301.

[27] Madsen, P.A., Murray, R., Sorensen, O.R. A new form of the Boussinesq equations with improved linear dispersion characteristics // Coastal Eng. 1991. Vol. 15. P. 371-388.

[28] Lovholt, F., Pedersen, G.K., Glimsdal, S. Coupling of dispersive tsunami propagation and shallow water coastal response // The Open Oceanography J. 2010. Vol. 4. P. 71-82.

[29] Madsen, P.A., Sorensen, O.R. A new form of the Boussinesq equations with improved linear dispersion characteristics. Pt 2. A slowly-varying bathymetry // Coastal Eng. 1992. Vol. 18. P. 183-204.

[30] Zijlema, M., Stelling, G.S., Smit, P. SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters // Coastal Eng. 2011. Vol. 58. P. 992-1012.

[31] Zijlema, M., Stelling, G.S. Efficient computation of surf zone waves using the nonlinear shallow water equations with non-hydrostatic pressure // Coastal Eng. 2008. Vol. 55. P. 780-790.

[32] Yamazaki, Н., Kowalik, Z., Cheung, K.F. Depth-integrated, non-hydrostatic model for wave breaking and run-up // Intern. J. for Numer. Methods in Fluids. 2009. Vol. 61, No. 5. P. 473-497.

[33] Wei, G., Kirby, J.T. Time-dependent numerical code for extended Boussinesq equations // J. of Waterway, Port, Coastal, and Ocean Eng. 1995. Vol. 121, No. 5. P. 251-261.

[34] Rygg, O.B. Nonlinear refraction-diffraction of surface waves in intermediate and shallow water // Coastal Eng. 1988. Vol. 12, No. 3. P. 191-211.

[35] Шокин Ю.И., Яненко Н.Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, Сибирское отделение. 1985. 364 с.

Shokin, Yu.I., Yanenko, N.N. Method of differential approximation. Application to gas dynamics. Novosibirsk: Nauka, Sibirskoe Otdelenie, 1985. 364 p. (In Russ.)

[36] Chubarov, L.B., Fedotova, Z.I., Shkuropatskii, D.A. Investigation of computational models of long surface waves in the problem of interaction of a solitary wave with a conic island // Russ. J. of Numer. Anal. and Math. Modelling. 1998. Vol. 13, No. 4. P. 289-306.

[37] Гусев О.И., Хакимзянов Г.С. Численное моделирование распространения длинных поверхностных волн по вращающейся сфере в рамках полной нелинейно-дисперсионной модели // Вычисл. технологии. 2015. Т. 20, № 3. С. 3-32.

Gusev, O.I., Khakimzyanov, G.S. Numer. simulation of long surface waves on a rotating sphere within the framework of the full nonlinear dispersive model // Comput. Technologies. 2015. Vol. 20, No. 3. P. 3-32. (In Russ.)

[38] Федотова З.И., Хакимзянов Г.С. Уравнения нелинейно-дисперсионной модели мелкой воды на вращающейся сфере и выполнение законов сохранения // Прикладная механика и техническая физика. 2014. Т. 55, № 3. С. 37-50.

Fedotova, Z.I., Khakimzyanov, G.S. Nonlinear-dispersive shallow water equations on a rotating sphere and conservation laws // J. of Appl. Mech. and Techn. Physics. 2014. Vol. 55, No. 3. P. 404-416.

[39] Гусев О.И., Шокина Н.Ю., Кутергин В.А., Хакимзянов Г.С. Моделирование поверхностных волн, генерируемых подводным оползнем в водохранилище // Вычисл. технологии. 2013. Т. 18, № 5. С. 74-90.

Gusev, O.I., Shokina, N.Yu., Kutergin, V.A., Khakimzyanov, G.S. Numerical modelling of surface waves generated by underwater landslide in a reservoir // Comput. Technologies. 2013. Vol. 18, No. 5. P. 74-90. (In Russ.)

[40] Гусев О.И. Алгоритм расчета поверхностных волн над подвижным дном в рамках плановой нелинейно-дисперсионной модели // Вычисл. технологии. 2014. Т. 19, № 6. С. 19-41. Gusev, O.I. Algorithm for surface waves calculation above a movable bottom within the frame of plane nonlinear dispersive model // Comput. Technologies. 2014. Vol. 19, No. 6. P. 19-41. (In Russ.)

[41] Fedotova, Z.I., Gusev, O.I., Khakimzyanov, G.S. New algorithm for numerical simulation of surface waves within the framework of the full nonlinear dispersive model // Proc. of VII Europ. Congress on Comput. Methods in Appl. Sci. and Eng. (ECCOMAS Congress 2016), 5-10 June 2016 / M. Papadrakakis, V. Papadopoulos, G. Stefanou, V. Plevris (Eds). Crete Island, Greece, 2016. Vol. 1. P. 1093-1103.

[42] Gusev, O.I., Beisel, S.A. Tsunami dispersion sensitivity to seismic source parameters // Sci. of Tsunami Hazards. 2016. Vol. 35, No. 2. P. 84-105.

[43] Хакимзянов Г.С., Шокина Н.Ю. Метод адаптивных сеток для одномерных уравнений мелкой воды // Вычисл. технологии. 2013. Т. 18, № 3. С. 54-79.

Khakimzyanov, G.S., Shokina, N.Yu. Adaptive grid method for one-dimensional shallow water equations // Comput. Technologies. 2013. Vol. 18, No. 3. P. 54-79. (In Russ.)

[44] Шокина Н.Ю., Хакимзянов Г.С. Усовершенствованный метод адаптивных сеток для одномерных уравнений мелкой воды // Вычисл. технологии. 2015. Т. 20, № 4. С. 83-106. Shokina, N.Yu., Khakimzyanov, G.S. An improved adaptive grid method for one-dimensional shallow water equations // Comput. Technologies. 2015. Vol. 20, No. 4. P. 83-106. (In Russ.)

[45] Cai, X., Pedersen, G.K., Langtangen, H.P. A parallel multi-subdomain strategy for solving Boussinesq water wave equations // Advances in Water Resources. 2005. Vol. 28. P. 215-233.

Поступила в 'редакцию 17 апреля 2017 г.

History of the development and analysis of numerical methods for solving nonlinear dispersive equations of hydrodynamics. II. Two-dimensional models

FEDOTOVA, ZlNAIDA I.*, KHAKIMZYANOV, GAYAZ S., GUSEV, OLEG I., SHOKINA,

Nina Yu.

Institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia Corresponding author: Fedotova, Zinaida I., e-mail: zf@ict.nsc.ru

The review of the finite-difference methods for solving two-dimensional NLD-equati-ons is presented. It is found that, despite of the successful application of these methods to the problems of wave hydrodynamics, the works on the theoretical investigation of the used difference schemes are nearly absent.

In order to fill the gap in this knowledge, the investigation of stability and dispersion properties is done for the series of difference schemes. Although the simple mathematical tools were used for the analysis, the detailed computations are presented in order to make future investigations of the properties of similar finite-difference schemes easier.

One of the main conclusions of this work is that the stability conditions for the difference schemes of the 2D-equations with dispersion are weaker than the similar conditions for the schemes which approximate the nondispersive shallow water equations. Therefore, this property, which was earlier discovered by the authors for one-dimensional schemes, is inherited by the two-dimensional schemes. The form of the stability conditions for the two-dimensional schemes also proved to be similar to the one-dimensional case. Nevertheless, such inheritance is not true for the dispersion property. In 1D-case, the value of the Courant number typically exists, for which the phase error of the scheme becomes minimal. But for the discussed two-dimensional the schemes similar property was not discovered.

The conducted investigations set a number of new problems. In particular, it becomes clear that it is necessary to develop the schemes, which are invariant with respect to rotation in order to improve the description of the dispersion properties of the model, in particular, to minimize the dependency of dispersion on the direction of the wave vector.

Keywords: nonlinear dispersive equations, two-dimensional models, numerical algorithms, finite-difference methods, stability, dispersion.

Acknowledgements. The study was supported by the RSF (project No. 14-17-00219) and grant of the President of the Russian Federation for State Support Leading Scientific School (grant No. NSh-7214.2016.9).

Received 17 Aprill 2017

© ICT SB RAS, 2017