История развития и анализ численных методов решения нелинейно-дисперсионных уравнений гидродинамики. I. одномерные модели Текст научной статьи по специальности «Математика»

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

Том 20, № 5, 2015

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

З.И. Федотова*, Г. С. ХАкимзянов, О. И. Гусев

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

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

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

Введение

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

К этому времени, т. е. к 70-м годам прошлого столетия, уже были разработаны и активно применялись на практике конечно-разностные (КР), конечно-элементные, псевдоспектральные и другие методы численного решения как эволюционных (гиперболических и параболических), так и эллиптических систем дифференциальных уравнений (ДУ). Были не только развиты практические подходы к решению этих уравнений, но и созданы соответствующие теории (обоснование существования, единственности, сходимости, точности), связанные с именами таких математиков, как П. Лакс, С.К. Годунов, Г.И. Марчук, А.А. Самарский, Н.Н. Яненко и др.

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

© ИВТ СО РАН, 2015

Благоприятным фактором, способствующим разработке численных алгоритмов для НЛД-моделей, было то, что путем введения новых переменных НЛД-уравнения удалось "расщепить" таким образом, что одна составляющая представляла собой систему обыкновенных дифференциальных уравнений (ОДУ) с правой частью, зависящей от пространственных производных гидродинамических функций, или неоднородную систему гиперболических уравнений, а другая составляющая не содержала производных по времени и в простейшем случае выглядела как эллиптическое уравнение. Это открыло широкие возможности для применения как методов решения ОДУ и конечно-разностных/конечно-объемных схем, так и метода конечных элементов, а также разнообразных комбинаций этих и других подходящих методов. Имеющееся разнообразие НЛД-моделей, а также многовариантные способы "расщепления" привели к большому количеству работ, предлагающих те или иные численные методы, обозреть которые в одной работе невозможно. Здесь мы ограничились КР-методом. Именно этот подход, вместе с примыкающим к нему методом конечных объемов, получил наиболее широкое распространение при решении практически значимых задач о поведении поверхностных волн в реальных акваториях.

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

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

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

Так как исходные уравнения и их конечномерные аппроксимации в полной постановке не поддаются аналитическим исследованиям, в дополнение к описанию полного метода привлекаются упрощенные постановки для изучения свойств в частных случаях. Как правило, анализ КР-схем начинается с одномерного случая с последующей цепочкой упрощений: это случаи ровного дна, линейные аналоги, нелинейное скалярное уравнение, описывающее однонаправленное распространение волны. Такие исследования проясняют суть численных методов для НЛД-моделей, их отличие от методов решения бездисперсионных уравнений мелкой воды, выстраивают иерархии по свойствам, характеризующим работоспособность алгоритмов и, наконец, приводят к набору

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

1. Нелинейно-дисперсионные модели

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

иа (х,у,г) = и(х,у,ха (х,у,г),г), (1)

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

V

и = -1 J и(х,у,г) <1г. (2)

-к

1.1. Базовая модель

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

И(ж, у, г, £) = ист(х, у, ¿) + ^2и1(х, у, г, £). (3)

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

Нг + У-

н и + 3)

0, (4)

К)4 + К■ VK + "ГТ = 77^к- 1\(Н3)г + К■ V)(ЯJ) + Н (3 ■ V) ист + (У-ист)Нз], (5)

л л л I

где V = (д/дх, д/ду); р — проинтегрированное по толщине слоя давление НЛД-модели, р0 — давление на дне,

дН2 Н3 Н2 Н2

р =---з"- ^ р° = 9Н - - HQ2, (6)

дг = В (V- ист) - (V- ист)2 ^2 = П2к, (7)

I [ д

J = — и^г, И = — + иа • V.

-н

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

Нг + V • (ниа) = 0, (8)

dHUa дНиа Ua dHva UCT dHJuua dHJv uc

+--^--1--^--+ vp +--^--+

dt

dx

dy

где

U = ua + J, u0

ua Va

dx

J

dy

p0Vh

(9)

Ju Jv

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

Рх Ро

1

Ht + [Н(и + J)]x = 0, ut + иих + Щ: = ^0hx - — (HJ)t + u(HJ)x + 2HJux

H H H

:iq)

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

н2 н3 н2 н2

Р = 3-2---^ (D(ux) - ul) - — D2h, ро = gH - — (D(ux) - и2х) - HD2h,

где D = d/dt + ud/dx. Если исходное 3В-течение является потенциальным, то [5]

j = -

\Н 1

— - (za + h) (hxt + 2ихhx + uhxx) - r L 2 J 6

H2 (za + h)

2 -,

Ua

:i2)

2

1.2. Нелинейно-дисперсионные модели при J = 0

Известные НЛД-модели получаются из базовой модели (8), (9) при соответствующем выборе вектор-функции J [5, 6]. Рассмотрим класс НЛД-моделей, для которых скорость определена по формуле (2) либо изначально предполагается независимость горизонтальной составляющей U вектора скорости трехмерного течения от вертикальной координаты z. При этих предположениях выполняется равенство J = 0 и уравнения выглядят проще, чем в базовой модели. К наиболее известным моделям этого класса относятся модели Перегрина [2], Грина — Нагди [7], Железняка — Пелиновского [8], Федотовой — Хакимзянова [9, 10] и многие другие, вытекающие из указанных моделей при преобразованиях формы дисперсионных членов уравнения движения с использованием уравнения неразрывности. Возможны также и другие модификации, проводимые с целью упрощения дисперсионных членов (например, с использованием ограничений на форму дна или характер течения).

В случае одной пространственной переменной рассматриваемые уравнения имеют следующий вид:

Ht + (Ни)х = 0, (Hu)t + (Ни2 + р)х = Pohx, (13)

где величины с давлением вычисляются по тем же формулам (11), что и в базовой модели. При анализе разностных схем нам зачастую будет достаточно рассматривать НЛД-уравнения, описывающие течения над горизонтальным дном к(х, Ь) = к0. В этом случае уравнения (13) упрощаются до следующих [11]:

Щ + (Ни)х = 0, (Ни)г + (Ни2 + р)х = 0, (14)

где Н = к0 + и дисперсионная составляющая давления имеет кубическую зависимость от полной глубины:

Н 2 н 3

Р = 9---зт{ихл + иихх -иХ). (15)

Уравнения (14), (15) имеют решение в виде уединенной волны:

~н (х, г)

т](х, Ь) = аоБвоЪ2 [к(х — х0 — и Щ, и(х, Ь) = и ^ , ) , (16)

где и — скорость перемещения волны,

и = инлд = к = кнлд = ^ (1 + (17)

ао — ее амплитуда, х = Хо — положение вершины волны при Ь = 0, Со = Vдк0.

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

Н2 к2

Р = д-^т -Нк-(ихЛ + иихх). (18)

Уравнения (14), (18) также имеют решение в виде уединенной волны, при этом формула для возвышения г/(х, Ь) отличается от (16), а скорость и(х, Ь) опять определяется через г/ по второй из формул (16) си = иНЛд.

Уравнения (14), (18) получены в [10] из уравнений полной НЛД-модели (14), (15), записанных в безразмерных переменных, при отбрасывании в осредненном давлении р/Н членов, имеющих порядок 0(ац2), где а и ^ — параметры нелинейности и дисперсии соответственно, а = ао/ко. Слабо нелинейную модель Перегрина [2] (модель типа модели Буссинеска) также можно получить из полной отбрасыванием членов такого же порядка малости, но уже из выражения рх/Н:

к2

Щ + (Ни)х = 0, щ + иих + д Г]х = -Зиххг. (19)

В отличие от НЛД-моделей (14), (15) и (14), (18), модель Перегрина (19) не имеет в качестве точного решения уединенную волну вида (16), не обладает законом сохранения полной энергии [10], и уравнение для полного импульса Ни не удается записать в консервативной форме вида (14). Различаясь в нелинейном случае, после линеаризации все

модели [2, 7, 8, 9] переходят в одну и ту же, уравнения которой в случае одномерных течений с нулевой фоновой скоростью имеют следующий вид:

h2

rit + h0ux = 0, ut + giix = у uxxt. (20)

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

1.3. Модельное BBM-уравнение

Изучение ряда существенных свойств численных методов целесообразно начинать с рассмотрения модельного аналога системы НЛД-уравнений — скалярного уравнения, описывающего волны, движущиеся в одном направлении. Впервые такое уравнение было выведено в работе [12]. Это так называемое RLW-уравнение (Regularised Long Wave equation), или BBM-уравнение (Benjamin — Bona — Mahony equation [13]). Его выводили на основе соотношений для функций, аппроксимирующих главную часть инвариантов Римана бездисперсионных уравнений мелкой воды [14], а также другим способом, опирающимся на длинноволновое разложение [15] дисперсионного соотношения линеаризованной модели потенциальных течений. Здесь мы приведем более прозрачный вывод, основанный на использовании безразмерных переменных и выборе специальной связи между скоростью и возвышением свободной границы.

Для перехода к безразмерным переменным в настоящей работе использованы такие же соотношения, как в [10]:

_ х - tc0 — h _ rq _ u _ p x = —, t =——, h = —, ^ =—, u =—, p = —2, L L ho ao Co c0

при этом далее черта над безразмерными переменными опущена. Перепишем систему уравнений (14), (18) в безразмерном виде

2 2 ß aß rqx

Ht + (Hu)x = 0, ut + uux + arqx —— (uxt + uuxx)x---—£ (uxt + uuxx ) = 0, (21)

3 3 л

где H = 1 + arq, и для получения скалярного уравнения вместо скорости используем некоторое ее приближенное представление через г/. Воспользуемся тем, что уединенная волна является примером возмущения свободной границы, распространяющегося в одном направлении. Поэтому при выводе BBM-уравнения будем считать, что скорость определяется как в уединенной волне, т.е. по второй из формул (16). В безразмерных переменных имеем

u = aUrj/H, U = у/1 + а. (22)

Используя в (22) разложения по степеням малого параметра а, приходим к выражениям

u = атц (1 + а/2 - arj) + 0(а3), Hu = ocq (1 + а/2) + 0(а3).

Отсюда вытекает, что система НЛД-уравнений (21) может быть переписана в виде двух уравнений относительно г/:

2

Vt + + T^j Vx = 0(a2), rit + (l - 2) Vx + 3ащх = уVxxt + О [a(a + ß2)] .

Полусумма этих уравнений

3 а л2

Ш + Лх + — Щх = — г]ххл + О [а(а + л 2 6

после отбрасывания остаточного члена и перехода к размерным переменным дает скалярное нелинейно-дисперсионное ББМ-уравнение

3 с о к 2о

Щ + со 'Г]х + -гт-'П 'Пх =г)ххл, (23)

2к0 6

описывающее распространение волны в одном направлении.

Приведем некоторые соображения в пользу того, что скалярное ББМ-уравнение лучше подходит в качестве модельного для систем НЛД-уравнений, чем классическое уравнение Кортевега — де Вриза [14]

3со сок2о

'Щ + со 'Пх + 7ТГ- Щх =--'Пххх. (24)

2к0 6

Уравнение (24) отличается от (23) формой дисперсионного члена и имеет решение г/(х, Ь) в виде уединенной волны (16) с отличающимися от (17) значениями и и к:

(1 +

V 2ко)

и = иМу = со ( 1 + ^ I > инлд, к = ккау = ^^ > кнлд. (25)

Задача Коши для ББМ-уравнения (23) также имеет решение г/(х, Ь) в виде уединенной волны (16). Впервые оно было получено в работе [16] для уравнения в безразмерных переменных. Формулу (16) получим здесь на основе 1апЪ-метода [17], согласно которому решение ищется в виде г](х, Ь) = г](Т(х)), где Т(х) = 1апЬ(х), X = кх - шЬ + 5. Все производные от функции Т по переменной х являются полиномами от Т:

Т' = 1 -Т2, Т" = -2ТТ' = —2Т + 2Т3, Т'" = -2Т'(1 - 3Т2) = -2 + 8Т2 - 6Т4,

а частные производные от вычисляются по формулам

щ = -ш'П'Т', пх = кп'Т', т]ххь = к2шТ' [ - г]'''(Т')2 + 6г]''ТТ' + 2г/(1 - 3Т2)],

где г/' = А'ц/АТ. Подставляя эти частные производные в уравнение (23) и сокращая на Т , получаем обыкновенное дифференциальное уравнение третьего порядка

(-ш + сок + ^к^ Г]' = к2к2ш[2(1 - 3Т 2Уп' + 6Т (1 - Т 2)Г]'' - (1 - Т 2)2ц'''], (26) 2к0 6

за решение которого можно взять квадратичную функцию

г](Т) = А0 + АгТ + А2Т2. (27)

Подставляя ее в (26), приходим к переопределенной системе из четырех нелинейных уравнений относительно трех искомых коэффициентов А0, Аг, А2:

3 2 к2

АЛ + °°к -ш - к2ш-3) = 0,

3

3сок ,2 . (Зсок . ,2 4к2\

А1 + 2А2 А + с0к -ш - к ш —° = 0, 2ко \ 2ко 3 )

Аг( + к2шк2^ = 0, А^ + 2к2шк2^ = 0.

Система (28) имеет единственное нетривиальное решение

2ho 3

й -1

Ai = 0,

4кшк1

А2 =--0,

3со

подставляя которое в (27) и принимая величину -А2 за амплитуду волны ао, получаем решение вида (16) с параметрами

-1/2

^Л.,, Л / 31 1.Г\ I I / I 1.П \

и = иввм = Со ( 1 + — ) = иК<ж, к = кввш

V 2Н0)

а0

У3ао9 Л + 2к0с0 I 2к0

(29)

при этом

кнлд < кввш < k~Kdv. (30)

Итак, уединенные волны — решения KdV- и BBM-уравнений — распространяются с одинаковой скоростью, большей, чем скорость аналогичной волны полной НЛД-модели, но волна (16), (29) BBM-уравнения несколько длиннее, чем волна (16), (25) KdV-уравнения и ее профиль ближе к форме уединенной волны (16), (17) (см. рис. 1, а, на котором изображены графики функции z = ri(x, 0) для четырех моделей при ao/h0 = 0.4, x^/h0 = 20). Таким образом, можно сказать, что уединенную волну полных НЛД-уравнений модельное BBM-уравнение аппроксимирует точнее, чем KdV-уравнение.

Кроме того, BBM-уравнение имеет преимущество перед KdV-уравнением также и по дисперсионным свойствам. В самом деле, дисперсионные соотношения ш = ш(к) линейного аналога полной модели (20), линеаризованные KdV-уравнения (24)

и BBM-уравнения (23)

Vt + с0'Пх = -

ОоЬ0

~'Пх

h0

Vt + с0 Vx = — Vxxt 6

(31)

(32)

8 1 fU 10

K/h0

Рис. 1. Форма уединенной волны (а) и зависимость фазовой скорости от длины волны (б) в НЛД-модели и для КБЖ-уравнения (кривая 1), для КёУ- и ББМ-уравнений — соответственно кривые 2 и 3

6

записываются в виде

c0k ( k2h0\ с0к

/ к2к2\

ШНЛД = ±0 + ^/з, = Сок{1 , шввм = 1+кк2/б, (33)

где к = 2ж/А — волновое число, А — длина периодической волны

ф, г) = Е0е(34)

Е0 — амплитуда волны, ш — частота. Функцию в = шЬ — кх называют фазой. Для фазовых скоростей ур = ш(к)/к волн, распространяющихся вправо, имеют место выражения

со ( 2тт2к0\ со . ,

УрНЛД + 4т°21г2/3А2, ^ = Ч1 — ~ЗАА2~) , ^'ввм = Т+ЪЩ/3\2 ■ (35)

Графики зависимостей фазовых скоростей (35) от А приведены на рис. 1, б. Видно, что во всем диапазоне длин волн модельное ББМ-уравнение лучше аппроксимирует фазовые скорости полной модели, чем К^У-уравнение. Поэтому мы отдаем предпочтение скалярному ББМ-уравнению как модельному для исследования ряда свойств разностных схем, аппроксимирующих НЛД-уравнения. На этом же рисунке изображен график (пунктирная линия) для фазовой скорости

4-

tanh(2nho/X) CoW 2,ho/X (36)

линейной модели потенциальных течений [14], в которой, в отличие от уравнений мелкой воды, вертикальные перемещения жидкости учитываются. Видно, что скорости движения коротких волн, описываемых в рамках НЛД- и ББМ-моделей, отличаются от эталонных значений (36), однако эти отличия все же меньше, чем в KdV-модели.

1.4. Модельное KRW-уравнение

В работах [18, 19] из модели Перегрина [2] выведено уравнение второго порядка по времени, которое описывает распространение с учетом силы Кориолиса квазиодномерных слабо нелинейных длинных волн по слабо изменяющемуся дну. Эту модель можно рассматривать как промежуточную между скалярной ББМ-моделью и моделями типа Буссинеска с осредненной по глубине жидкости скоростью. В размерных переменных это уравнение имеет вид

1,

Vtt -gV^ (hV г] ) = - f 2г] V- hV (hr]tt) +- gV^V( v2) . (37

3

2

В работе [20] уравнение такого вида выведено с учетом трения. Оказалось, что уравнения вида (37) можно использовать не только как модельные для предварительной апробации численных алгоритмов, предназначенных для решения систем НЛД-уравнений, но и для решения практических задач волновой гидродинамики. Так, например, в [21] успешно решены сложные задачи взаимодействия достаточно длинных волн с различного рода препятствиями.

Видимо, впервые уравнение вида (37) было приведено (без вывода) в [22], поэтому далее будем называть его модельным KRW-уравнением (Kim — Reid — Whitaker equation)

vp

Здесь будет рассматриваться одномерный вариант этого уравнения (при / = 0) в случае переменной глубины к(х)

1 Г 1 3 Ли - д(^Пх)х = з (к'п)х\ха + 2з^2)хх (38)

и постоянной к,:

К2 3

Ли - 4лхх = у 'Пххы + 29(Л2) хх . (39)

Уравнение (39) часто называют модифицированным уравнением Буссинеска [23, 24]. Некоторые свойства схем исследуются на линейном аналоге уравнения (39):

2 к2

Щ — Со1хх = у 'пххЦ. (40)

Поскольку система линейных уравнений (20) и линейное скалярное уравнение (40) эквивалентны, дисперсионные соотношения этих моделей совпадают. Но более поразительным является то, что скалярное нелинейное ККЖ-уравнение (39) имеет в качестве своего решения г/(х,1) в точности ту же уединенную волну (16), (17), как и система НЛД-уравнений (14), (15).

2. Конечно-разностные схемы для модели Перегрина

Учет дисперсии в НЛД-моделях обеспечен членами со смешанными по Ь и х производными третьего порядка. К тому времени, когда эти модели стали применяться для расчетов, уже была создана теория КР-схем для эволюционных систем дифференциальных уравнений первого порядка [3, 25], поэтому построение формальных аппроксимаций НЛД-уравнений не вызывало трудностей. Однако для достижения эффективности вычислений потребовалась разработка новых КР-схем.

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

2.1. Конечно-разностные алгоритмы для BBM-уравнения

Впервые разностные схемы для регуляризованного уравнения длинных волн были исследованы в работе [26], в которой приведены пять неявных разностных схем и дан анализ их точности и устойчивости. Позднее три из них были обобщены на случай систем НЛД-уравнений. Эти три схемы будут рассмотрены здесь применительно к нелинейно-дисперсионному ББМ-уравнению (23), записанному в виде

'Ш + с0 'Пх + ащх = (41)

где введены сокращенные обозначения для коэффициентов: а = 3с0/2к0, V = /6.

Пусть к иг — шаги равномерной сетки в плоскости переменных х и (хз, £п) — узлы сетки, г]п — значения сеточной функции г/ в этих узлах. При написании КР-схем будем использовать стандартные обозначения производных, взятые из [27]:

пп _ у?+1- у! п _ у! - 'пи п _ у?+1- 'пи '1х'3 к ' '1х'3 к ' ''х,з 2к '

пп+1 'П 'П гп-1 пп+1 гп-1

Пп _ 'Ъ - Уп „ _ У3 - Уз „ _ Уз - Уз

_ т ' Щ'3 _ т ' % _ 2т '

П _ У'п+1 - 2УЧ + У]-1 п _ Т,'п+1 - 2^ + ^ 1 п _ Уп+1 - Ухх,з

'1хх,3 к2 ' ' 1Ы>3 ' ' 1ххг,з '

Первая из разностных схем, исследованных в [26], предложена ранее в [12]. Используя вышеприведенные обозначения и опуская нижний индекс ], запишем ее в следующем виде:

< + (Со + аг1п) (2Уп+1 + 2Уп) _ ^хг- (42)

Эта неявная схема имеет второй порядок аппроксимации по пространству, первый — по времени и реализуется скалярными прогонками.

В другой схеме [26] конвективные члены аппроксимируются центральными разностями с п-го слоя:

< + (Со + аг]п) 'П1 _игЦл, (43)

а функции с (п + 1)-го слоя по времени используются только в слагаемых, содержащих производные по времени. Аналог этой схемы лежит в основе одного из алгоритмов метода расщепления для систем НЛД-уравнений [28 - 30].

Еще одна рассмотренная в [26] двухслойная схема может быть интерпретирована как аналог известной схемы Кранка — Никольсон [27], широко используемой для решения эволюционных уравнений:

У'п + 2(со + аг]п+1)^ + 2(Со + аяп)^ _ и^- (44)

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

Замечание 1. Стремление авторов работы [26] использовать во всех своих схемах центральные разности объясняется, видимо, желанием повысить порядок точности. Однако надо иметь в виду, что явные схемы с центральными разностями могут быть абсолютно неустойчивыми [31], а преимущество неявных схем по устойчивости может быть потеряно из-за ограничения на шаг по времени, связанного с обеспечением корректности метода прогонки [3, 31]. ■

Замечание 2. С позиций современной вычислительной гидродинамики существенным недостатком предложенных в [26] схем является их неконсервативность. Известны примеры [27, 32] неконсервативных схем, которые могут расходиться, поэтому консервативность является в настоящее время одним из ключевых требований, предъявляемых к разностным схемам для НЛД-уравнений [33]. С этой точки зрения более предпочтительной, чем, например, неконсервативная схема (44), является консервативная схема

У'п + 1 +1 + ¡хп) _щпяЛ' соу + а^' (45)

аппроксимирующая со вторым порядком дивергентную форму уравнения (41):

^ + ¡х = ЩXXI. (46)

Консервативная схема (45), записанная в дивергентной форме, может быть переписана и в эквивалентной недивергентной форме

'пП + + + + аГУпХх = щХХ*, (47)

которая отличается от (44) нелинейными слагаемыми за счет использования в (47) осредненных значений 'г/п = (гцп+1 + г/?^) /2. ■

Перейдем к анализу разностных схем из [26]. Исследование нацелим на получение следующей информации о схеме: точность, устойчивость, поведение диссипации и дисперсии. Эти свойства изучим на примере схемы с весами

п

$ + * + (1 - 7= (48)

аппроксимирующей линейное BBM-уравнение (32). Применительно к этому уравнению схемы (42) и (44) совпадают и записываются в виде (48) при 7 = 1/2. При 7 = 0 имеем схему (43) в линейном приближении.

Все указанные свойства позволяет исследовать анализ элементарных решений вида (34), для которых легко выписать оператор перехода [4] со слоя времени t на слой (t + т) как для разностной схемы, так и для аппроксимируемого ею дифференциального уравнения. Подставляя решение вида (34) в уравнение (32), получаем для него оператор перехода

П = е-гшт = \П\егф, (49)

где \ П \ = 1, а дисперсионное соотношение ш = wBBM определяется по третьей из формул (33) и введено обозначение Ф = arg П = —шт.

В случае скалярных разностных схем оператор перехода обычно называют множителем перехода и обозначают через р: rqn+1 = рг/п, при этом

Р = е-шт = \р\егф, (50)

где ш = ш(к) — дисперсионное соотношение разностной схемы, ф = arg р.

Множитель перехода р будем рассматривать как функцию от переменной £ = kh (£ Е [0, ). Устойчивость и аппроксимация разностной схемы определяются через свойства функции р($,). Информацию о диссипативных и дисперсионных свойствах КР-схемы дает поведение диссипации % = \ \ — \^^\ КР-схемы и ее фазовой ошибки Аф = Ф — ф [4].

Подставляя гармонику (34) в схему (48), получаем следующее выражение для множителя перехода:

= а — i(1 — 7 )Ь

а + i'jb

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

2 1 ■ t л 2 2 2 -С т

—s , b = со ж sin t, с = 4спж s , s = sin -, ж = —. h2 vi 2 h

а

Легко видеть, что при 7 > 1/2 имеет место неравенство |р| < 1, поэтому при таких 7 схема (48) абсолютно устойчива. В частном случае 7 = 1/2 выполняется равенство |р| = 1, следовательно, схемы (42) и (44) абсолютно устойчивы в линейном приближении и диссипация для них отсутствует.

При 7 = 0 получаем, что |р|2 = 1 + Ь2/а2 > 1 и

г2т2 с2

max (И2) = 1 + -< 1 + г J h2 + 4 v 4v

т.е. схема (43) слабо устойчива в линейном приближении [25].

Замечание 3. Если в BBM-уравнении (32) пренебречь дисперсионным слагаемым и для аппроксимации полученного уравнения переноса использовать линейный вариант схемы (43) с v = 0, то такая схема с центральной разностью будет при ж = const абсолютно неустойчивой [31]. Учет дисперсионного слагаемого улучшает устойчивость схемы, переводя ее в класс слабо устойчивых схем. Этот факт, установленный нами для схемы (43), имеет место, как увидим позднее, и в общем случае, а именно, схемы для бездисперсионной модели мелкой воды устойчивы при более жестких условиях, чем аналогичные схемы для нелинейно-дисперсионных уравнений. ■

Исследуем дисперсию КР-схемы (48) при 7 = 1/2. Преобразуем выражение (51) для к виду

а -ib/2 (а -ib/2)2 1 , 2 2 ,

Р =-т~т = —tWt~ = —т^гг (а2 - Ъ2/4 - гаЬ)

' а + i Ь/2 а2 + Ь2/4 а2 + Ъ2/4К ' '

и рассмотрим малые значения £, соответствующие длинным волнам. Тогда

а2 - b2/4 ^ Со ж о с°ж° о и .о 2£4\

ф=- *«** ^/ц = - **+—е+е++о 2 у.

(53)

С другой стороны, для ББМ-уравнения (32) с дисперсионным соотношением ш = шввм имеем следующую формулу:

ф = -т ^ввм = -1 T7e/h2 = - ^ +СоЖ + о(и 2 hd ■

54)

1 + v i2/h2

При получении этих соотношений предполагалось, что

ж = const, (55)

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

Заметим, что несмотря на то, что в формулах (53) и (54) в знаменателях остаточных членов содержится малая величина h, эти остаточные члены имеют порядок малости О (¡л4), где л — параметр дисперсии. В самом деле, поскольку ^ = kh = 2T\h\i/h0} то

Р4 4-к4

= —л4, (56)

h4 9 и ' v J

т. е. остаточные члены в (53) и (54) имеют тот же порядок малости, что и члены уравнений, отброшенные при получении НЛД-моделей [10].

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

3 3

Д0 = - ^ е3 - ^ + (57)

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

33

^+^ ? « «« -¡¡¡е.

Следовательно, для числа Куранта с0ж должно выполняться неравенство

< Л/22(¥—Т), (58)

где

* = ¡0. (59)

Из ограничения (58) следует, в частности, что шаг к по пространственной переменной должен быть меньше глубины к0 и на достаточно мелких сетках (8 ^ 1) схемная дисперсия не будет подавлять дисперсию модели.

2.2. Разностные схемы для KRW-уравнения

В работе [19] для аппроксимации уравнения (38) предложена следующая разностная схема второго порядка аппроксимации:

Чп - 9 (-Лх)П = 1 [(ь (-л)х )Л + 29 (V2)Пх . (60

Здесь к(х) — переменная глубина и использовано стандартное обозначение [27] для оператора второй разностной производной от произвольной сеточной функции и в случае произвольного переменного коэффициента V:

( _ ) = У3+1/2их,]+1 - У]-1/2их,з

(Ьих>х,] = - ,

где к — шаг равномерной сетки, ь^+1/2 = (ь^ + )/2. Если в модели из работы [20] не учитывать поверхностное натяжение и вязкость, то приведенная в [21] КР-схема также записывается в виде, близком к (60).

Опишем свойства схемы (60), перейдя к ее линейному аналогу

'Пп - с2лПх = V('Пхх)п, V = ¡2/3, (61)

аппроксимирующему уравнение (40). Множитель перехода р удовлетворяет квадратному уравнению

,2 -Л с

- Ч1 - Та)? +Т=°'

в записи которого использованы обозначения (52). При с < 4а для его корней выполнены равенства 1р\^1 = 1. Поэтому разностная схема (61) является условно устойчивой при следующем ограничении на число Куранта [19, 21]:

1

Со Ж <\11 + 4 &2, (62)

где параметр 8, характеризующий относительный размер пространственного шага к, определен соотношением (59). Интересно отметить, что при измельчении сетки условие устойчивости (62) становится все менее и менее ограничительным и при 8 ^ 1 превращается в ограничение на шаг по времени:

т-7То

1.1Ьт0,

(63)

где т0 = к0/с0 — характерное время, необходимое для того, чтобы волна, распространяющаяся со скоростью с0, прошла расстояние, равное характерной глубине к0. При условии (62) изменение фазы множителя перехода выражается формулой

ф = ± arccos (1--)

V 2aJ

±

«а — е+^ (с>2 -е3+°(е4)

24

поэтому для фазовой ошибки Лф = Ф — ф, где

Со

Ф = ±

имеет место выражение

у/1 + vi2jh2

Л ф

С0Ж

~2л

±( соа* — соя + 0(£4))

(с0я2 — 1) е3 + о(с4).

(64)

(65)

(66)

Из формулы (66) следует, что если с0ж = 1, то главная часть фазовой ошибки не содержит "схемную дисперсию". Поскольку при таком выборе с0ж схема (61) устойчива (см. (62)), можно сказать, что она имеет преимущество перед схемой (48) второго порядка аппроксимации (при 7 = 1/2), так как влияние "схемной дисперсии" в рассматриваемой здесь схеме будет минимальным.

Если с0ж = 1, то дисперсия модели превалирует над "схемной дисперсией" при выполнении условия

С0Ж < у/1 + 482. (67)

Замечание 4. При подходе длинных волн к мелководной части акватории возможно появление ондулярных боров [12, 34], которые имеют дисперсионную природу и, в отличие от НЛД-моделей, не описываются в рамках бездисперсионной модели мелкой воды. Как следует из результатов исследований [34], приемлемое описание дисперсионных эффектов на мелководье получается при 8 = 2... 4. Аналогичные исследования [35] показывают, что при решении практических задач о длительном распространении волн в океане с последующим их выходом на шельф достаточно выбрать пространственный шаг сетки так, чтобы он соответствовал значению = 1 на мелководье и = 4 в глубоководной части океана. ■

2.3. Алгоритмы для одномерной модели Перегрина на ровном дне

В работе [12], наряду с КР-схемой (42) для скалярного уравнения (41), предложена КР-схема для одномерной системы уравнений Буссинеска (19) в случае ровного дна. Указывается, что "прямой перенос КР-аппроксимации (42) на случай (19) (h(x) = h0 = const) приводит к неустойчивой схеме", и автор [12] предлагает разностную схему с пересчетом, в которой вначале по явной формуле

+ (ho + ) „? + гГ,{.; = 0

т V 2 ) х х

вычисляется вспомогательное значение г] , затем методом скалярной прогонки находится решение ип+1 разностного уравнения в неявной форме:

<+ип{ 2ип+1+й)+*(? +2 4')

11

-0

п\ — иип и = —

затем по явной формуле

/1 1 \ ип+1 +ип < + (-о + Уп){2ипх+1 + -ип^ + и +и г£ = 0

2 х 2

выполняется пересчет для нахождения свободной границы г]п+1 на новом временном

слое.

Аналог КР-схемы Перегрина для системы линейных уравнений (20) принимает вид

Г]* — Г]п

+ коип = 0, ип + д[\г! * + 1уп) = иип^, чп + ко( 1ип+1 + -ип) = 0.

1 п

Г,-'»

2 х

х 2 х 2 х

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

2

Ч,п + -2о{ип+1 +ип) =0, ип + д г% = 1^ипх + ии:

п , ххЛ,,

(68)

где

п

п

и О и о

п = х,3 + 1 х,3-1

иоо . --.

хх,з 2-

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

Подстановка в разностную схему (68) решения в виде гармоник (34) и и(х, Ь) = и0е-г-кхх) приводит к системе двух уравнений относительно амплитуд гармоник. Условие существования нетривиального решения Е0 и и0 равносильно обращению в нуль определителя этой системы:

Р- 1

• , 1\ • С

ь-^КР+ЧЫЪЯ,

2

%д&вт^ (р — 1)а+ —

0,

(69)

где использованы обозначения (50) и (52). Таким образом, для вычисления р имеем

квадратное уравнение

—К1—I >+1=0-

корни которого вычисляются по формуле

Р1'2 = 1 — 2Ь ■

(70)

(71)

Для устойчивости необходимо, чтобы корни (71) были комплексно-сопряженными, т. е. необходимо выполнение неравенства Ь2 < 4а или

2 2 2/1 2\ ^ 1 , 2 с0а; в (1 — 5 ) < 1 + — в .

к2

Легко проверить, что необходимым и достаточным условием для справедливости неравенства (72) при всех в Е [0, 1] является ограничение

со а - 1 + \/1 + р2, (73)

где параметр 8, введенный по формуле (59), характеризует степень разрешимости сетки по отношению к характерной глубине к0. Для достаточно мелких сеток (8 ^ 1) условие устойчивости (73) можно заменить на (63), т.е. получаем условие устойчивости в виде независимого от к ограничения на шаг по времени .

Замечание 5. Проведенное исследование показывает, что не всегда разностные схемы для систем уравнений обладают теми же свойствами, что и лежащие в их основе схемы для скалярных модельных уравнений. Так, схема (42) в линейном приближении абсолютно устойчива, а построенная на ее основе схема (68) для системы уравнений Перегрина является условно устойчивой. ■

Интересно отметить, что необходимое условие устойчивости (73) оказалось менее жестким, чем условие (62) для модельного ККЖ-уравнения. Оно является и менее жестким по сравнению с условием устойчивости с0ж < 2 разностной схемы

$ + к0 + ч1) =0, и? + дЦ = , (74)

совпадающей с (68) при = 0 и аппроксимирующей бездисперсионные линейные уравнения мелкой воды (уравнения (20) при V = 0):

щ + к0их = 0, щ + д т]х = 0. (75)

Если в разностной схеме (68) не использовать пересчет г/, что эквивалентно применению схемы

к0 2

то вместо (71) имеем выражения

< + Т (UT + Ut) =0, UX +9 Vnx = vul

- 0, ul + qrfl = vuXxt, (76)

у a \16a J

^ = l — £ Ч a — 1 I' (77)

В случае вещественных корней (Ь2 > 16а) для меньшего из них справедлива оценка

b2 b2 ( b2 \ Ь2 16а

Pi = 1 — — \ —--1 - 1 — - 1 — = —3,

4а \/ а \ 16а 4а 4а

поэтому устойчивость возможна только в случае комплексно-сопряженных корней, что реализуется при условии Ь2 < 16а, или

Со а - 2(1 + ^1 + 4$2 ).

При выполнении последнего условия для корней (77) справедлива оценка 1 < max (|р|2) = 1 + 2(-°QT V < 1 + -33, т2,

«е[о,*] J \h + 7h2T4VJ 2т2 '

что означает слабую устойчивость схемы (76). Отметим, что аналогичная схема без пересчета для бездисперсионных уравнений (75), т.е. схема

ко

*п + -Цип+1 +ип) =0, ип + д ^ = 0,

0 . _о . ^о . _г , ,,о (79)

2 х х х

будет абсолютно неустойчивой при законе предельного перехода (55). Таким образом, схема без пересчета по отношению к свойству устойчивости аналогична схеме (43) для ББМ-уравнения.

При выполнении условия (73) корни (71) являются комплексно-сопряженными числами и | р1>21 = 1, при этом изменение фазы множителей перехода выражается формулой

(I)2 \ "

1 — 2^) = ± [со^ — ^ + (с0ж2 — ^ ^ + 0(^4)

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

Д0 = — 4) е3 + о(е4). (81)

Видно, что если с0ж = 2, то дополнительная дисперсия, вносимая разностной схемой (68), будет минимальной, при этом схема будет устойчивой.

Подводя итог, можно сказать, что первые написанные для дисперсионных уравнений разностные схемы характеризуется тем, что они были получены путем прямой аппроксимации всех членов уравнений, включая смешанные производные третьего порядка, при этом разностные уравнения конструировались так, чтобы избежать применения векторных прогонок. В [2] выполнено обобщение схемы из работы [12] для случая плоского откоса, в [36] — для неровного дна, в [37-40] выполнено обобщение КР-схемы Перегрина на разнесенные сетки, на которых, в частности, упрощается постановка разностных граничных условий.

3. Специальные приемы построения численных методов решения НЛД-уравнений

Продолжим рассмотрение НЛД-моделей, для которых скорость определена по формуле (2). Мы начали с анализа КР-схем для ББМ-уравнения и модели Перегрина потому, что они были исторически первыми [12, 2] численными методами, предложенными для решения задач длинноволновой гидродинамики. При этом надо отметить, что слабонелинейная модель Перегрина, полученная в работе [2], может быть выведена из более точных НЛД-моделей, созданных позже другими авторами. Особого внимания заслуживают первые полностью нелинейные НЛД-модели Грина — Нагди [7] и Железняка — Пелиновского [8]. Выведенные независимым путем, эти модели оказались эквивалентными, что показано, например, в [9] путем преобразования выражения для дисперсионных членов уравнения движения с использованием уравнения неразрывности. Однако надо отметить, что "однотипные" формальные КР-аппроксимации этих моделей могут обладать разными (непохожими) свойствами.

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

Мы не будем останавливаться на таких свойствах, как консервативность, сохранение энергии и других инвариантов, а приведем простые примеры.

Рассмотрим линейные аналоги моделей Грина — Нагди и Железняка — Пелиновского для ровного дна к = ко в одномерном случае, которые обычно рассматриваются для того, чтобы выполнить анализ Фурье (гармонический анализ) модели и аппроксимирующей ее разностной схемы. Первая из моделей записывается в виде

к

щ + коих = 0, щ + д 'Пх = —30 ^хи. (82)

Исключив из правой части уравнения движения путем использования уравнения неразрывности, получим линейный аналог модели Железняка — Пелиновского, совпадающий с линейным аналогом модели Перегрина (20). Дисперсионные соотношения этих эквивалентных моделей совпадают и задаются формулой ш = шнлд из (33). Нетрудно видеть, что фазовая и групповая скорости при к ^ то стремятся к нулю, что означает несущественность высокочастотных гармоник (свойство, полезное при конструировании численных алгоритмов).

Приведем "однотипные" разностные схемы

к

< + ко ипх =0, ип + д гр-1 = — -0 ^ (83)

^ + -оипх = 0, ^ + д т,п+1 = ^ипх, (84)

х х 3

для моделей (82) и (20) соответственно, где

V

I /Х1-, 1 — I !П

п =_

itх,j 2к

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

Методом гармонического анализа можно показать, что разностная схема (83) абсолютно неустойчива. В то же время схема (84) является условно устойчивой, поскольку уравнение для р совпадает с (70) и поэтому необходимым условием устойчивости будет неравенство (73). Этот результат оказался неожиданным ввиду эквивалентности моделей (82) и (20) и общего для них дисперсионного соотношения. Таким образом, "однотипные" (похожие) разностные схемы при применении к разным НЛД-моделям могут существенно различаться по своим свойствам.

Нетрудно видеть, что для получения устойчивой разностной схемы (84) надо вместо (83) использовать неявную схему

к

цп + коип+1 = 0, ип + д т,п+1 = — к0 ^ (85)

х 3

Тогда действительно приходим к эквивалентным разностным схемам (84) и (85) с одинаковыми условиями устойчивости. Однако первая КР-схема предпочтительнее второй, так как ее реализация сводится к скалярной трехточечной прогонке, в отличие от схемы (85), для которой требуется векторная прогонка. Из примера вытекает, что, вероятно, модель Железняка — Пелиновского более предпочтительна для проведения численных расчетов, чем эквивалентная ей модель Грина — Нагди.

Надо отметить, что в статье [7], где приведена наиболее известная форма уравнений Грина — Нагди, при построении численного алгоритма для случая ровного дна производные ^ были исключены с использованием уравнения неразрывности. Описанный в этой работе алгоритм является одним из первых для модели Грина — Нагди. Он основан на модифицированном методе Эйлера, который применяется для аппроксимации по времени. Аппроксимация по пространству основана на центральных разностях. Алгоритм расчета является слабо устойчивым, однако авторы отмечают, что неустойчивость в расчетах не была замечена, так как изучались слабо меняющиеся течения. При необходимости авторы предлагают использовать пятиточечный фильтр, который в одномерном случае для произвольной функции принимает следующий вид:

Л- = 1 (-Л--2 + 4Л- +10Л- + 4— ъ+2).

В дальнейшем при построении КР-схем использовалась только модифицированная форма модели Грина — Нагди, по сути совпадающая с моделью Железняка — Пелинов-ского, но которую в англоязычных публикациях по-прежнему называют моделью Грина— Нагди, иногда — моделью Серре — Грина — Нагди [41]. При построении численных алгоритмов используются различные формы записи этих моделей, при этом члены в уравнениях группируются в зависимости от способа аппроксимации в конечномерном пространстве.

Вопросам численной реализации моделей с осредненной по толщине слоя жидкости скоростью посвящено большое количество работ. Особенность численного решения систем уравнений, определяющих эти модели, состоит в том, что группа членов, описывающих дисперсию, зависит от производной скорости и по времени. Наиболее простой способ решить эту проблему состоит в том, чтобы в уравнении движения все члены, содержащие щ, представить в виде выражения Qt, что позволяет свести реализацию соответствующего численного алгоритма к поэтапному решению системы ОДУ и эллиптического уравнения. Этот подход впервые описан в работе [26], где скалярное уравнение щ + их + иих = иххЛ было записано в виде (и — ихх)1 + их + иих = 0, удобном для расщепления на два уравнения

Qt + их + иих = 0, и — ихх = Q, (86)

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

Распространим этот прием на базовую модель (10), предварительно записав ее в дивергентном виде

Н + (Ни)х = 0, (Ни\ + (Нии)х + (НЗи)х + Рх = Р0кх, (87)

где и = и + 3. Далее преобразуем часть уравнения, связанную с "давлением". Используя выражения (11), выделим из рх — р0кх следующие члены:

. Н 3 \ н 2 /Н з \ /Н 2 \ /Н 2 \

— (— ихг) + —ихЛкх = —у — их) г + (Н 2Ни)х + —ихкх)^ — кх)Пх = Qt + Q2,

обозначив оставшиеся через Q1. При этом

, Н 3 \ Н 2 / Н 2 \

Q = — — их) + —ихкх, Q2 = (Н 2Них)х — — кх)Пх.

Теперь уравнение движения из (87) можно переписать в виде

У + (Нии )х + (НЗи)х = —Я1 — Я2, (88)

где

ни + Я = У. (89)

Такое представление позволяет строить КР-методы, основанные на расщеплении уравнения движения на систему ОДУ и уравнение эллиптического типа. Один из наиболее распространенных способов построения численного алгоритма состоит в следующем. Сначала из уравнения неразрывности по двухслойной явной схеме (или неявной относительно Н) вычисляется Нп+1. Затем явным способом с использованием вычисленных Нп+1 находится Уп+1 на основе разностного уравнения, аппроксимирующего (88). После этого решается эллиптическое (в одномерном случае — обыкновенное дифференциальное) уравнение Ни + Q = Уп+1 относительно ип+1 (в одномерном случае — методом скалярной прогонки).

В работах [28, 42] аналогичный подход выделения членов с производными по Ь использован для построения конечно-разностных аналогов для модифицированной модели Грина — Нагди и модели Алешкова, написанных в свободной неконсервативной форме, не придерживающейся никакой структуризации. Как следствие, построенные разностные схемы выглядят очень громоздкими. Все пространственные производные аппроксимированы в "лоб" центральными разностями, для аппроксимации по времени применяется описанный выше метод Перегрина с пересчетом переменной . Эти схемы являются условно устойчивыми и имеют порядок аппроксимации 0(т + к2). В работе [42] предложены также схемы, в которых скорость и определена в полуцелых узлах. Показано, что это приводит к лучшему воспроизведению КР-схемой фазовых скоростей аппроксимируемых моделей.

При 3 = 0 варианты КР-схемы с выделением Qt рассматривались также в работах [29, 39, 40], в которых предложен ряд модификаций, например, использовались разнесенные сетки и КР-схемы с весами. Для случая 3 = 0 в статье [43] предложен метод предиктор-корректор расчета на основе КР-схемы Кранка — Никольсон. В работе [44] для аппроксимации модели Нвогу была разработана КР-схема предиктор-корректор, основанная на известном методе Адамса — Мултона — Бэшфорта четвертого порядка аппроксимации.

Два других приема, которые также можно проиллюстрировать на базовой модели, состоят в том, что с помощью дополнительных переменных р = щ и ф = щ удается разделить производные по времени и пространству, получив для нахождения переменных р и ф эллиптические уравнения, а для и и г/ — систему ОДУ. В случае 3 = 0 для такого расщепления требуется лишь одна переменная р = щ. Для численной реализации полностью нелинейной НЛД-модели, допускающей преобразование Галилея, эффективным приемом оказывается введение переменной d = щ+иих [45]. Все эти варианты рассмотрены в работах [29, 30], где на их основе построены КР-методы для модели Железняка — Пелиновского и НЛД-модели, в которой скорость и(х, Ь) определяется как скорость течения полной модели на расстоянии (2/3)1/2к(х) от донной поверхности. Конечно-разностные схемы [29] имеют порядок 0(т + к2) и условно устойчивы.

Разделение производных по времени и пространству путем использования переменной р = щ осуществлено также в работе [46] применительно к обобщенному ББМ-уравнению. По пространственным переменным использованы аппроксимации высокого

порядка точности, однако для решения уравнения иг = р применена схема Эйлера первого порядка точности. Такая же переменная р = щ использована в работе [47], а полученное при этом эллиптическое уравнение решалось методом конечных элементов.

Недостатком рассмотренных подходов является использование уравнений в недивергентной форме. Также следует отметить, что в двумерном случае получается система двух дополнительных уравнений эллиптического типа относительно компонент вектора скорости, что усложняет вычислительный алгоритм. Кроме того, в вычислительных алгоритмах не используются методы решения гиперболических уравнений, развитые для модели мелкой воды первого гидродинамического приближения. В последнее время появились работы (см., например, [48, 49]), в которых вместо системы ОДУ (88) рассматриваются системы уравнений гиперболического типа, что дает возможность построения КР-схем с ТУВ-ограничителями.

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

Первое применение модели Перегрина для решения инженерных задач описано в работе [37]. Для построения КР-схемы НЛД-уравнения были записаны в квазиконсервативной форме. Для аппроксимации ДУ использована неявная разностная схема типа "ящик" на разнесенной сетке. Проведенный авторами анализ аппроксимационных ошибок показал, что "схемная" дисперсия имеет тот же вид, что и дисперсия НЛД-модели, что при конечных шагах разностной сетки может привести к искажению дисперсионной картины течения, описываемого НЛД-моделью. Для избежания этого эффекта в [37, 38] использован известный прием повышения порядка аппроксимации, когда разностную схему модифицируют, аннулируя члены аппроксимационной ошибки "нежелательного вида" путем добавления их в разностную схему уже в явном виде, но с противоположным знаком. Это позволяет занулить главные (по порядку аппроксимации) члены дисперсии схемного происхождения и улучшить воспроизведение дисперсии аппроксимируемой модели. В итоге авторы [37] построили неявную разностную схему повышенной точности, для реализации которой требуются расщепление по пространственным переменным (АВ1-метод) и векторная прогонка.

Подобный прием повышения точности численных расчетов использован в широко известной работе [43], где основное внимание сосредоточено на выводе НЛД-модели с улучшенным воспроизведением дисперсии 3В-модели. Напомним, что полностью нелинейный вариант модели Нвогу, полученный в работе [50], вытекает из базовой модели при выборе скорости в точках поверхности (х) = —0.531к(х). Здесь для описания численного алгоритма рассмотрим семейство одномерных НЛД-моделей на ровном дне

щ + (Ни)х + (( + 3)к0иххх = 0, иг + иих + д г]х + (Ь0иххг = 0, (90)

где Н = к0 + г],

( _ ^д + ^

5 = 2к20 + к0 ,

га = —0.531к0. Отметим, что, выбирая величину ( другими способами, можно получить известные модели. Например, при ( = —1/3 получаем модель Перегрина (19). Для численной апробации этой модели взята КР-схема, опирающаяся на известную схему Кранка — Никольсон.

Чтобы продемонстрировать суть подхода, примененного в работах [37, 43] для повы шения точности разностной схемы, рассмотрим схему Кранка — Никольсон (схему (48 при 7 = 1/2)

^ +1Со (^Т + - ^ = 0, и = (91)

аппроксимирующую линейное ББМ-уравнение (32) со вторым порядком. Г-форма [4] второго дифференциального приближения схемы (91) записывается в следующем виде:

т к2 т2

Уг + °оУх - VУххг + 2 (У* + °оУх - иУххг)1 - 12иУхххх* + - иУхх)ш+

т2 ^2

+-тСоУхы + — Со 1]ххх = 0(т3, к3). (92)

4 6

Если рассмотреть это уравнение на классе функций, удовлетворяющих дифференциальному представлению, то члены порядка 0(т) удается исключить из (92), при этом получаем

т2

Уг + Со'Пх - Щххг - 12Уш + — Со'Пххх - ('Пхх - а2Уы)хх1 = 0(т3, И3). (93)

Так как аппроксимационная ошибка содержит третьи производные — члены г/ш и т/ххх, "схожие" с дисперсионным членом гцххЛ НЛД-модели, в [37, 43] предложено исходную схему модифицировать добавлением разностного аналога этих членов, взятых с противоположным знаком, и тем самым улучшить дисперсионные свойства схемы. Если следовать работе [43], то "улучшенную" разностную схему надо записать в следующем виде:

^ + 1 со (С1 + ^) - ^ + ^ = 0, (94)

где

т2 Г]п+1 - 3Г]п + 3Г]п-1 - Г]п-2 к2

Вя

- к^со (Ф1 + .

12 \ XXX XXX/

'аррг 12 тз 12'

Полученная схема — неявная, многослойная. Ее дисперсионным соотношением является кубическое уравнение с комплексными коэффициентами, два корня которого — "ложные".

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

,пп+1/2 'п ('п+1/2 ,пп)

У - У , п \У - У )ХХ п п п+1/2 п п

-1/2- + СоУ° - V-1/2- = 0, Уь + СоУо ' - "УххЛ = 0,

и и и П+1

каждый из которых реализуется скалярной прогонкой, а полученное значение ^ используется затем в (94) как для вычисления добавочного члена Даррг, так и для задания начальных значений итерационного процесса, реализующего неявную схему Кранка — Никольсон.

Описанный алгоритм нетрудно распространить как на систему уравнений (90), так и на нелинейную НЛД-модель Нвогу [43] с переменным дном. Отметим, что в [43], в отличие от [37], уравнения рассматриваются в неконсервативной форме. Анализ устойчивости разностной схемы для линейной системы уравнений, вытекающей из (90), проводится без учета члена Даррг (он трактуется как правая часть, вычисленная заранее).

Авторы показали, что в случае системы уравнений, в отличие от скалярного уравнения, разностная схема Кранка — Никольсон является условно устойчивой. Приведено огрубленное (по сравнению с (73)) условие устойчивости в виде неравенства с0ж < 2.

Замечание 6. В работах [37, 43] утверждается, что предложенные схемы имеют порядок аппроксимации 0(т3, к3), однако это только для случая и = 0, поскольку при и = 0 некоторые дисперсионные члены аппроксимированы с порядком 0(т2,к2) (см. уравнение (93)). ■

В качестве альтернативного подхода в работе [44] предложена КР-схема четвертого порядка аппроксимации для НЛД-модели Нвогу. Тогда "схемная дисперсия", вводимая производными третьего порядка, исключается автоматически. Опишем эту схему для Ш-случая, переписав уравнения (90) в специальном виде:

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

где

и = и + (к20ихх. (95)

В [44] предложена КР-схема предиктор-корректор, на первом шаге которой по явной схеме Адамса — Бэшфорта третьего порядка аппроксимации вычисляются значения

пп+1 ип+1:

о ' 1

г!]+1 = + 12 [23Щ — 16Еп~1 + 5Еп-2] ,

Щ+1 = Щ + 12 [23Е? — 1Щ-1 + 5Е™-2] . (96)

Затем с применением метода скалярной прогонки численно решается уравнение (95)

и]+1 + (к0 и1+1 = Щ+1 (97)

и вычисленные величины , иг^+1 используются для получения предварительных значений Еп+ и Е™+1.

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

'п]+1 = Г]1] + 2^ [9Е]+1 + 19Щ — 5Щ-1 + Еп~2\ ,

Щ+1 = Щ + 12 [Щ+1 + 1Щ — 5Е™-1 + Е;1-2] , (98)

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

Замечание 7. В схемах (96), (98) все пространственные производные первого порядка аппроксимируются центральными разностями четвертого порядка. В то же время производные по пространственной переменной, входящие в группу дисперсионных членов (95), аппроксимируются со вторым порядком. Поэтому общий порядок аппроксимации схемы по пространству — второй. ■ В работе [50] НЛД-модель Нвогу распространена на полностью нелинейный случай и для полученной модели выписана КР-схема, являющаяся обобщением вышеописанной.

Для полностью нелинейной базовой модели (87), (88) этот алгоритм построения разностных схем, рассмотренный в п. 3 в более общем случае (консервативная форма

уравнений, переменное дно), сводится к решению системы ОДУ для переменных Н, V и эллиптического уравнения (89) для нахождения и.

В работе [51] на основе такой формулировки построен КР-метод четвертого порядка аппроксимации для полностью нелинейной модели Серре, подправленной дополнительными членами для улучшения аппроксимации фазовой скорости. Для аппроксимации пространственных производных использована конечно-объемная компактная схема четвертого порядка аппроксимации, а для интегрирования системы ОДУ по времени применен стандартный метод Рунге — Кутты такого же порядка точности.

Аналогичная структура численного метода для Ш-модели Серре использована в работе [33]. Для аппроксимации пространственных производных здесь разработан конечно-объемный метод второго порядка аппроксимации. Для решения системы ОДУ применен четырехстадийный метод Рунге — Кутты третьего порядка аппроксимации. Для исследования точности построенного метода применен псевдоспектральный метод, который является высокоточным, но имеет ограниченную область применения. Сравнение методов на тестах показало перспективность применения метода конечных объемов.

5. Алгоритмы, основанные на выделении эллиптического уравнения для дисперсионной составляющей давления

Система уравнений НЛД-модели (8), (9) не является системой типа Коши — Ковалевской в силу того, что уравнение движения (9) содержит смешанные производные третьего порядка по времени и пространственным переменным от компонент вектора скорости. В предыдущих разделах на упрощенных модельных уравнениях мы рассмотрели приемы непосредственной аппроксимации этих производных, а также аппроксимации расширенной системы уравнений, которая не содержит таких смешанных производных и получена выделением эллиптической части для скорости.

Здесь мы рассмотрим другое расщепление, результатом которого является расширенная система уравнений, состоящая из эллиптического уравнения для дисперсионной составляющей р, проинтегрированного по глубине давления р из (6), и гиперболической системы уравнений с производными по времени первого порядка, которая отличается от классической модели мелкой воды наличием в правой части уравнений движения дополнительных слагаемых, связанных с дисперсионными добавками. Такое расщепление оказалось плодотворным для конструирования численных алгоритмов решения задач как в плоской геометрии [52-55], так и в сферической [56].

Недостатком первых работ этого направления [45, 57, 58] является то, что полученная в результате расщепления система уравнений относительно полной глубины и скорости отличалась от классических уравнений мелкой воды и имела кратные собственные значения. Отметим также работу [59], в которой в рамках слабо дисперсионной модели численно исследовалось цунами в Индийском океане в декабре 2004 г. с использованием расширенной системы, одно из уравнений которой имело эллиптический тип и служило для коррекции давления [60]. Недостаток предлагаемого в [59] алгоритма заключается в неприемлемом условии устойчивости, предписывающем, что размеры пространственных шагов не должны быть слишком малыми, в частности, в одномерном случае необходимо выполнение неравенства к > 1.5к0.

5.1. Схема с расщеплением для линеаризованного BBM-уравнения

Вначале суть подхода и свойства схемы расщепления продемонстрируем на примере линейного BBM-уравнения (32), для которого аналогом дисперсионной составляющей давления (18) является величина р = vr¡xt, где v = hh0/6. Расширенная система уравнений [52]

Vt + соrjx = Px, Pxx - i = СоT/xx (99)

состоит из неоднородного уравнения гиперболического типа для r¡ и уравнения "эллиптического" типа относительно р. Последнее уравнение не содержит производных по времени от функции r¡.

В методе предиктор-корректор [53] уравнения (99) решаются поочередно на каждом

из двух вычислительных этапов. Вначале решается уравнение для р:

рП

pXx,3 - — = Соv*x,j. (100)

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

1W+1 + У) , _ п п

rfj+1/2 2(1l1j+1 + 1lfj) . п п пп-о

-¡- + СО Vn¿ = Px,j, (101)

где т* = т(1 + 9™+1/2)/2, +i/2 — схемный параметр, надлежащий выбор которого гарантирует выполнение TVD-свойства для численного решения модельных скалярных уравнений [61, 62]. При теоретическом исследовании рассматриваемой схемы мы будем предполагать, что в = const.

На шаге предиктор уравнение для p решается еще один раз, теперь уже с привлечением вычисленных значений г/*:

P*j+3/2 - Щ+1/2 + P*j-1/2 P*j+1/2 _ - 2Vj+1/2 + Vj-1/2

К2 V К2

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

^ + С - '-1/2 = Р+1/2 - Р*-1/2. (103)

При 9 = О (К) рассматриваемая явная схема имеет второй порядок аппроксимации.

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

Р*+3/2 - Щ+1/2 + Р*-1/2 Р**+1/2

h2

п * 2 п i* п fin/l^

C0Vxx,j+1/2 — Т CoVxxx,j + Т COpxxx,j, (104)

rffj + <ЩIJ = г*с1пПхх,3 - г*СоРпx,j + ^*+1/2 h Р-1/2 , (105)

где

п _ Vxx,j + Vxx,j+1 п _ Vxx,j+1 Vxx,j п _ lpXX,j+1 lpXX,j

Vxx, j+1/2 2 , ^xxx,j h , pxxx,j h '

Подставляя в уравнения (100), (104), (105) гармоники

г]1] = Е0 рпег*, р] = Фо рпег*, р*+1/2 = Ф*0 pnei(j+1/2)è,

получаем следующие выражения для амплитуд Ф0, Ф0 и множителя перехода р:

^ 4с°и 2 F 2со" ( ■ с ■ *4со 2\

Фо =—rirs Ео, Ф0 = —nrs sin £ -г т —s )Ео, an2 an2 V an J

с(1 + в) b , Л

где использованы обозначения (52). Следовательно, необходимое условие устойчивости |р| < 1 будет эквивалентно выполнению неравенства

с2(1 + в)2 + 4а2 [Ъ2 - с(1 + в)] < 0 (107)

или при всех z Е [0, 1] неравенства

/ 4 V \ 2 c20œ2(1 + в)2z - (^1 + n2z) (z + d) < 0,

где z = s2. Достаточным условием для этого являются следующие ограничения:

в> 0, coœ < , 1 . (108)

> ' < VTTû v 7

Исследуем дисперсионные свойства схемы предиктор-корректор. Из выражения (106) для множителя перехода получаем следующую формулу, описывающую изменение фазы гармоники:

-ф = - arg р = -с0ж£ + c^œ^ + с20ж2(3в + 1) - ^£3 + 0(р4). (109)

Сравнивая с (54), видим, что главная часть фазовой ошибки

Со ж / 2 2

Аф = ^œ2(36 + 1) - 1)е3 + 0(р4)

_ --------- . Д10)

6

имеет в общем случае тот же порядок по £, что и "физическая" дисперсия. Если для заданного значения > 0 число Куранта выбрать по формуле

1

С0& = . , (111)

то условие устойчивости (108) будет выполняться и фазовая ошибка как минимум на один порядок по будет меньше "физической" дисперсии. Аналогичная "бездисперсионная" схема рассматривалась ранее в [63] применительно к уравнению переноса с постоянным коэффициентом.

График изменения фазы (109) множителя перехода схемы предиктор-корректор (100)-(103) при значениях 9 = 0, с0ж =1, 5 = 2 приведен на рис. 2 (линия 1). Для сравнения приведен также график (линия 2) аналогичной зависимости (53), характеризующей дисперсионный свойства КР-схемы (48) при = 1/2, аппроксимирующей

Рис. 2. Графики изменения фазы для множителей перехода схем предиктор-корректор (1) и (48) (2) и для операторов перехода, соответствующих ББМ-уравнению (32) (3) и системе НЛД-уравнений (20) (4)

ББМ-уравнение (32). Видим, что при одинаковом втором порядке аппроксимации схема предиктор-корректор лучше воспроизводит изменение фазы (54) оператора перехода модели, изображенное на рис. 2 в виде графика зависимости (линия 3)

фввм = -е (112)

На рис. 2 изображен также график изменения фазы в линеаризованной модели (20) (линия 4), который выражает собой зависимость (65), взятую со знаком минус:

фнлд = —, СоЖ е- (113)

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

Расширенная система уравнений [52] получается из системы линейных уравнений (20), если последнюю дополнить уравнением относительно новой зависимой переменной р = П^итл/3, которая является линейным аналогом дисперсионной составляющей давления (18):

фх ф 2

^ + Коих = 0, щ + д'Пх = -¡—, '-Рхх--= с^хх, (114)

По V

где V = П0/3.

Схема предиктор-корректор (100)-(103) в случае расширенной системы (114) выглядит следующим образом:

2

'£хх,з = Со1хх,1, (115)

и у

4*3+1/2 - 2 ^1+1 + "Л?) . , п п ^+1/2 " 1 (и]+1 + и1) п 1 п , . ---+ ко их,3 = 0, ---+ дпх^ = (116)

Т Т п о

^+3/2 2^*+1/2 + ^-1/2 ф)+1/2 _ 2 4*3+3/2 Щ+1/2 + ^¿-1/2

к2 у о к2

П17)

Uj+l/2 Uj-l/2 __ „ , ^j+l/2 Vj-1/2 _ 1 Pj+l/2 Pj-l/2

П I 7 J + 1/2 J-1/2 _ п „.«I „ V + 1/2 'J-1/2 _ ~ ' J + 1/2 ' J-1/2 ОЛ

vt,3 + ho — -o, U«3 +g — _ - — . (118)

Исключая из уравнений (117), (118) предикторные величины г/* и U и предполагая, что в _ const, получим другую форму последних трех уравнений

Р*3+Ъ/2 - Щ+1/2 + Р*з-!/2 Р*з+1/2 _ 2 /п -T*hUn \ (119)

—2 ^ _ с0 ^4xx,j+l/2 ' -0Uxxx,j) , (119)

+ —oul3 _ Т*с0пПх,з - тз^«х,з, <3 + g_ т*clUnxJ + —- Vj+l/2 - . (120)

Подстановка в уравнения (115) и (119) решения в виде гармоник

^1 _ Еорпег*, un _ Uoрпег^, рп _ Форпег^, р*+1/2 _ Ф*0рпеЪ+1/2К

приводит к выражениям

4clv ~ ж, 2clv / . А . А—о 2тт\ Фо _ ~ah2 S о, Ф0 _ -^8{Ео -г тз—s2Uoj ,

поэтому из (120) при учете обозначений (52) получается следующая система двух уравнений относительно амплитуд Е0 и и0:

Е0(р - 1 + 2с°ж2(1 + в) ^ + и^яЪо ъж^ = 0,

Ео[г--- + Uo[p - H----s J _ 0.

Следовательно, множитель перехода удовлетворяет уравнению

(р - 1)2 + (р - »«±<1 + ^ + Ь1

а 4а2 а

корни которого

Р_1 - ^ ± 4 (121)

2 а л/а

отличаются от (106). Условие устойчивости |р|2 < 1 можно записать в виде квадратного неравенства

с2ж2(1 + в)2 г - + + < 0, (122)

которое должно выполняться при всех г = в2 Е [0, 1], а для этого достаточно потребовать, чтобы выполнялись такие же неравенства (108), как для схемы предиктор-корректор (100)-(103), аппроксимирующей ББМ-уравнение.

Что касается дисперсионных свойств схемы (115)—(118), то они также схожи со свойствами схемы (100)-(103) для скалярного уравнения. В самом деле, изменение фазы множителя перехода (121) выражается формулой

ф _ aarg р _ ± arccos

\ с(1 + в) ч 2а

)/1р1

П23)

-1

Рис. 3. Графики изменения фазы для множителей перехода схем предиктор-корректор (1) и Перегрина (2), для операторов перехода НЛД-модели (3) и линеаризованной модели потенциальных течений (4)

поэтому для длинных волн имеет место зависимость

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

На рис. 3 приведен график изменения фазы (123) (взятой со знаком минус) множителя перехода схемы предиктор-корректор (115)—(118) при значениях 9 = 0, с0ж = 1, 8 = 2 (линия 1). Видно, что схема предиктор-корректор лучше воспроизводит изменение фазы (113) оператора перехода НЛД-модели (линия 3), чем схема Перегрина (линия 2). На рис. 3 изображен также график Ф lpf (£) (линия 4) изменения фазы в линеаризованной модели потенциальных течений:

Из рис. 3 видно, что схема предиктор-корректор удовлетворительно описывает дисперсионные свойства решений НЛД-модели для значений £ = кк < ^/4, которые соответствуют длинам волн А > 8к, что при 5 = 2 приводит к неравенству А > 4к0. Таким образом, дисперсионные свойства волн, длины которых по крайней мере в четыре раза больше глубины акватории, будут воспроизводиться предложенной схемой предиктор-корректор с хорошей точностью.

Заключение

Первоначальная цель настоящей работы состояла в составлении литературного обзора численных методов решения НЛД-уравнений гидродинамики. Однако анализ опубликованных результатов для модельных уравнений с дисперсией высветил множество

(124)

(125)

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

В работе даны уточненные условия устойчивости и получены новые сведения о дисперсионных свойствах схем. Условия устойчивости для нелинейно-дисперсионных уравнений оказались совершенно другими по сравнению с условиями для уравнений гиперболического типа. Отличие состоит в том, что для КР-схем, аппроксимирующих НЛД-уравнения, в условие устойчивости входит новый параметр 5, который характеризует размер горизонтального шага сетки по сравнению с характерной глубиной. В пределе при измельчении сетки условия устойчивости формулируются в виде ограничения только на шаг по времени, а не на отношение шагов ж, как это бывает для гиперболических уравнений, при этом в самом ограничении также появляется новый параметр т0 = h0/c0 — характерное время, необходимое для того, чтобы волна, распространяющаяся со скоростью Со, прошла расстояние, равное характерной глубине h0. Оказалось, что в области устойчивости некоторых схем имеются такие значения числа Куранта с0ж, для которых влияние "схемной дисперсии" будет минимальным.

Как уже было сказано во Введении, многого не удалось рассмотреть из-за ограничений на объем статьи. Например, не освещены проблемы неустойчивости, связанные с нелинейностью уравнений или с неровностью дна [64]. Совсем не упомянуты проблемы, связанные с аппроксимацией краевых условий [36]. За рамками работы остался обзор решенных с помощью НЛД-уравнений задач, интересных с практической точки зрения. Не рассмотрены проблемы реализации плановых (2D-) моделей.

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

Благодарности. Исследование выполнено при поддержке Российского научного фонда (проект № 14-17-00219).

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

[1] Mei, C.C., Le Mehaute, B. Note on the equations of long waves over an uneven bottom // Journal of Geophysical Research. 1966. Vol. 71, No. 2. P. 393-400.

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

[3] Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1968. 592 c.

Rozhdestvenskiy, B.L., Yanenko, N.N. Systems of quasilinear equations and their application to gas dynamics. Moscow: Nauka, 1968. 592 p. (in Russ.)

[4] Шокин Ю.И., Яненко Н.Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, Сиб. отд-ние, 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.)

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

[6] Шокин Ю.И., Федотова З.И., Хакимзянов Г.С. Иерархия моделей гидродинамики длинных поверхностных волн // Доклады Академии наук. 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.

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

[8] Железняк М.И., Пелиновский Е.Н. Физико-математические модели наката цунами на берег // Накат цунами на берег: Сб. науч. тр. Горький: ИПФ АН СССР, 1985. С. 8-33. Zheleznyak, M.I., Pelinovsky, E.N. Physico-mathematical models of the tsunami climbing a beach // Tsunami climbing a beach: Collection of scientific papers. Gorky: Institute Applied Physics Press, 1985. P. 8-33. (in Russ.)

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

[10] Fedotova, Z.I., Khakimzyanov, G.S., Dutykh, D. On the energy equation of approximate models in the long-wave hydrodynamics // Russian Journal of Numerical Analysis and Mathematical Modelling. 2014. Vol. 29, No. 3. P. 167-178.

[11] Su, C.H., Gardner, C.S. Korteweg —de Vries equation and generalizations. III. Derivation of the Korteweg —de Vries equation and Burgers equation // Journal of Mathematical Physics. 1969. Vol. 10, No. 3. P. 536-539.

[12] Peregrine, D.H. Calculations of the development of an undular bore // Journal of Fluid Mechanics. 1966. Vol. 25, pt 2. P. 321-331.

[13] Benjamin, T.B., Bona, J.L., Mahony, J.J. Model equations for long waves in nonlinear dispersive systems // Philosophical Transactions of the Royal Society of London. A. 1972. Vol. 272. P. 47-78.

[14] Уизем Дж. Линейные и нелинейные волны. М.: Мир, 1977. 622 с.

Whitham, G.B. Linear and nonlinear waves. N.Y.: John Wiley & Sons Inc., 1974. 636 p.

[15] Пелиновский Е.Н. Гидродинамика волн цунами. Н. Новгород: ИПФ РАН, 1996. 276 c. Pelinovskiy, E.N. Tsunami wave hydrodynamics. Nizhniy Novgorod: Institute Applied Physics Press, 1996. 276 p. (in Russ.)

[16] Olver, P.J. Euler operators and conservation laws of the BBM equation // Mathematical Proceedings of the Cambridge Philosophical Society. 1979. Vol. 85. P. 143-160.

[17] Hereman, W. Shallow water waves and solitary waves // Mathematics of Complexity and Dynamical Systems (Ed. R.A. Meyers). N.Y.: Springer, 2011. P. 1520-1532.

[18] Новиков В.А., Федотова З.И. Численное моделирование распространения длинных волн в бухтах на основе упрощенной модели Буссинеска // Всесоюз. совещание по численным методам в задачах волновой гидродинамики, 23-25 сент., 1990, Ростов-на-Дону: Сб. научн. тр. / Под ред. Ю.И. Шокина. Красноярск: ВЦ СО АН СССР, 1991. С. 9-14. Novikov, V.A., Fedotova, Z.I. Numerical modelling of long wave propagation in bays on the base of simplified Boussinesq model // Proceedings of All-union Conference on Numerical

Methods in Wave Hudrodynamics Problems 23-25 Sept., 1990, Rostov-na-Dony: Collect. Sientific Papers / Ed. Yu.I. Shokin. Krasnoyarsk: Computer Center of the Siberian Branch of Academy of Sciences of the USSR, 1991. P. 9-14. (in Russ.)

[19] Fedotova, Z.I. On properties of two simplified nonlinearly dispersive models of long wave hydrodynamics // International Journal of Computational Fluid Dynamics. 1998. Vol. 10, No. 2. P. 159-171.

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

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

[21] Литвиненко А.А., Хабахпашев Г.А. Численное моделирование нелинейных достаточно длинных двумерных волн на воде в бассейнах с пологим дном // Вычисл. технологии. 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 // Computational Technologies. 1999. Vol. 4, No. 3. P. 95-105. (in Russ.)

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

[23] Боголюбский И.Л. Модифицированное уравнение нелинейной струны и неупругое взаимодействие солитонов // Письма в Журн. эксперим. и теор. физики. 1976. Т. 24, № 3. С.184-186.

Bogolyubskii, I.L. Modified equation of a nonlinear string and inelastic interaction of solitons // Pis'ma v Zhurnal eksperimental'noy i teoreticheskoy fiziki. 1976. Vol. 24, No. 3. P. 184-186. (in Russ.)

[24] Хабахпашев Г.А. Влияние трения жидкости о дно на динамику гравитационных возмущений // Изв. АН СССР. Механика жидкости и газа. 1987. № 3. С. 119-127. Khabakhpashev, G.A. Effect of bottom friction on the dynamics of gravity perturbations // Fluid Dynamics. 1987. Vol. 22, No. 3. P. 430-437.

[25] Рихтмайер Р., Мортон K. Разностные методы решения краевых задач. М.: Мир, 1972. 420 с.

Richtmyer, R.D., Morton, K.W. Difference methods for initial-value problems. N.Y.: Interscience Publishers, 1967. 405 p.

[26] Eilbeck, J.C., McGuire, G.R. Numerical study of the regularized long-wave equations 1: Numerical methods // Journal of Computational Physics. 1975. Vol. 19. P. 43-57.

[27] Самарский А.А. Теория разностных схем. М.: Наука, 1989. 616 с. Samarskii, A.A. The theory of difference schemes. USA: Marcel Dekker, 2001. 788 p.

[28] Kompaniets, L.A. Analysis of difference algoritms for nonlinear dispersive shallow water models // Russian Journal of Numerical Analysis and Mathematical Modelling. 1996. Vol. 11, No. 3. P. 205-222.

[29] Fedotova, Z.I., Pashkova, V.Yu. Methods of construction and the analysis of difference schemes for nonlinear dispersive models of wave hydrodynamics // Russian Journal of Numerical Analysis and Mathematical Modelling. 1997. Vol. 12, No. 2. P. 127-149.

[30] Chubarov, L.B., Fedotova, Z.I., Shokin, Yu.I., Einarsson, B.G. Comparative analysis of nonlinear dispersive models of shallow water // International Journal of Computational Fluid Dynamics. 2000. Vol. 14, No. 1. P. 55-73.

[31] Годунов С.К., Рябенький В.С. Разностные схемы. М.: Наука, 1977. 440 с. Godunov, S.K., Ryabenkii, V.S. Difference schemes / Translation by E.M. Gelbard. Amsterdam: North-Holland, 1987. 486 p.

[32] LeVeque, R.J. Numerical methods for conservation laws. Berlin: Birkhauser Verlag, 1992. 214 p.

[33] Dutykh, D., Clamond, D., Milewski, P., Mitsotakis, D. Finite volume and pseudospectral schemes for the fully nonlinear 1D Serre equatios // European Journal of Applied Mathematics. 2013. Vol. 24, No. 5. P. 761-787.

[34] Grue, J., Pelinovsky, E.N., Fructus, D., Talipova, T., Kharif, C. Formation of undular bores and solitary waves in the Strait of Malacca caused by the 26 December 2004 Indian Ocean tsunami // Journal of Geophysical Research. 2008. Vol. 113. C05008.

[35] Glimsdal, S., Pedersen, G.K., Atakan, K., Harbitz, C.B., Langtangen, H.P., Lovholt, F. Propagation of the Dec. 26, 2004, Indian Ocean Tsunami: Effects of dispersion and source characteristics // International Journal of Fluid Mechanics Research. 2006. Vol. 33, No. 1. P. 15-43.

[36] Shokin, Yu.I., Chubarov, L.B. The numerical modelling of long wave propagation in the framework of non-linear dispersion models // Computers and Fluids. 1987. Vol. 15, No. 3. P. 229-249.

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

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

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

[40] 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 // Russian Journal of Numerical Analysis and Mathematical Modelling. 1998. Vol. 13, No. 4. P. 289-306.

[41] Bonneton, P., Barthelemy, E., Chazel, F., Cienfuegos, R., Lannes, D., Marche, F., Tisser, M. Recent advances in Serre — Green —Naghdi modelling for wave transformation, breaking and runup processes // European Journal of Mechanics-B/Fluids. 2011. Vol. 30. P. 589-597.

[42] Компаниец Л.А. О разностных схемах для нелинейно-дисперсионных моделей Грина — Нагди и Алешкова с улучшенной аппроксимацией дисперсионных соотношений // Вы-числ. технологии. 1995. Т. 4, № 11. С. 144-153.

Kompaniets, L.A. On difference schemes for nonlinear dispersive Green —Naghdi and Aleshkov models with improved approximation of dispersion relations // Computational Technologies. 1995. Vol. 4, No. 11. P. 144-153. (in Russ.)

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

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

[45] Barakhnin, V.B., Khakimzyanov, G.S. On the algorithm for one nonlinear dispersive shallow-water model // Russian Journal of Numerical Analysis and Mathematical Modelling. 1997. Vol. 12, No. 4. P. 293-317.

[46] Ramos, J.I. Explicit finite difference method for the EW and RLW equation // Applied Mathematics and Computation. 2006. Vol. 179. P. 622-638.

[47] Walkley, M., Berzins, M. A finite element method for the one-dimensional extended Boussinesq equations // International Journal for Numerical Methods in Fluids. 1999. Vol. 29. P. 143-157.

[48] Shi, F., Kirby, J.T., Harris, J.C., Geiman, J.D., Grilli, S.T. A high-order adaptive time-stepping TVD solver for Boussinesq modelling of breaking waves and coastal inundation // Ocean Modelling. 2012. Vol. 43-44. P. 36-51.

[49] Kirby, J.T., Shi, F., Tehranirad, B., Harris, J.C., Grilli, S.T. Dispersive tsunami waves in the ocean: Model equations and sensitivity to dispersion and Coriolis effects // Ocean Modelling. 2013. Vol. 62. P. 39-55.

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

[51] Cienfuegos, R., Barthelemy, E., Bonneton, P. A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations // International Journal for Numerical Methods in Fluids. 2006. Vol. 51. P. 1217-1253.

[52] Гусев О.И., Шокина Н.Ю., Кутергин В.А., Хакимзянов Г.С. Моделирование поверхностных волн, генерируемых подводным оползнем в водохранилище // Вычисл. технологии. 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 // Computational Technologies. 2013. Vol. 18, No. 5. P. 74-90. (in Russ.)

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

[54] Шокин Ю.И., Бейзель С.А., Гусев О.И., Хакимзянов Г.С., Чубаров Л.Б., Шо-кина Н.Ю. Численное исследование дисперсионных волн, возникающих при движении подводного оползня // Вест. ЮУрГУ. 2014. Т. 7, № 1. С. 121-133.

Shokin, Yu.I., Beisel, S.A., Gusev, O.I., Khakimzyanov, G.S., Chubarov, L.B., Shokina, N.Yu. Numerical modelling of dispersive waves generated by landslide motion // Bulletin of the South Ural State University. 2014. Vol. 7, No. 1. P. 121-133. (in Russ.)

[55] Khakimzyanov, G.S., Gusev, O.I., Beisel, S.A., Chubarov, L.B., Shokina, N.Yu.

Simulation of tsunami waves generated by submarine landslides in the Black Sea // Russian Journal of Numerical Analysis and Mathematical Modelling. 2015. Vol. 30, No. 4. P. 227-237.

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

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

[57] Барахнин В.Б., Хакимзянов Г.С. Использование расщепления при решении нелинейно-дисперсионных уравнений мелкой воды // Доклады Академии наук. 1999. Т. 364, № 4. С. 444-446.

Barakhnin, V.B., Khakimzyanov, G.S. The spliting technique as applied to the solution of the nonlinear dispersive shallow-water equations // Doklady Mathematics. 1999. Vol. 59, No. 1. P. 70-72.

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

Khakimzyanov, G.S., Shokin, Yu.I., Barakhnin, V.B., Shokina, N.Yu. Numerical simulation of fluid flows with surface waves. Novosibirsk: FUE "Publishing House SB RAS", 2001. 394 p. (in Russ.)

[59] Dao, M.H., Tkalich, P. Tsunami propagation modelling — a sensitivity study // Natural Hazards and Earth System Science. 2007. Vol. 7. P. 741-754.

[60] Horrillo, J., Kowalik, Z., Shigihara, Y. Wave dispersion study in the Indian Ocean-tsunami of December 26, 2004 // Marine Geodesy. 2006. Vol. 29. P. 149-166.

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

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

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

[63] Shokina, N.Yu. To the problem of construction of difference schemes on movable grids // Russian Journal of Numerical Analysis and Mathematical Modelling. 2012. Vol. 27, No. 6. P. 603-626.

[64] Lovholt, F., Pedersen, G. Instabilities of Boussinesq models in non-uniform depth // International Journal for Numerical Methods in Fluids. 2009. Vol. 61. P. 606-637.

Поступила в редакцию 21 июля 2015 г.

History of the development and analysis of numerical methods for solving nonlinear dispersive equations of hydrodynamics. I. One-dimensional models problems

Fedotova, Zinaida I.*, Khakimzyanov, Gayaz S., Gusev, Oleg I.

Institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia * Corresponding author: Fedotova, Zinaida I., e-mail: [email protected]

The article describes the main stages of development of finite-difference methods for numerical solving of nonlinear dispersive (NLD) hydrodynamic equations and presents new results on the theoretical research. The peculiarity of NLD-equations is a presence

© ICT SB RAS, 2015

of mixed derivatives of the third order in time and space, which bring special features into difference schemes used for their approximation. This affects their theoretical properties (stability conditions, balancing in modeling of numerical and physical dispersion, etc.) and well as the ways for implementation of numerical algorithms.

The paper has analyzed four groups of finite difference algorithms. The first one is the schemes, for which the discrete models are designed by direct approximation of all members of equations, including the mixed derivatives of the third order, with scalar sweeps for realization of numerical algorithm.

The second group includes those numerical algorithms that are based on the decomposition of the NLD-equations into a system of ordinary differential equations and into the differential equations that do not contain time derivatives.

To the third group we refer the algorithms which use special schemes of higher order accuracy in the methods of the first and second groups. Finally, the last group consists of difference schemes based on the splitting of the NLD-equations in such a way that there is the scalar equation of an elliptic type for the dispersive component of the pressure and there is a hyperbolic system of equations, which mimics the nondispersive shallow water model.

Since in the full statement the NLD-equations and their finite-difference approximation are not amenable to analytical study, it is necessary to simplify formulations for the study of properties in particular cases. Here are considered the dispersive scalar equations and systems of linear analogues of the NLD-equations in the cases of both a plane and other profiles of the bottom. Such studies help to clarify the essence of numerical methods for the NLD-models and highlight the advantages and disadvantages of methods. They also show their distinction from methods of solving nondispersive shallow water equations.

In the process of this study the corrected conditions of stability and new knowledge on the dispersion properties of finite-difference schemes are obtained. The conditions of stability of schemes for dispersive equations turned out to be different than for hyperbolic ones. The difference is that the specified conditions of stability includes a new parameter characterizing the fineness of the grid compared to the characteristic depth, and in the limit of grid refinement the new conditions are formulated in the form of restrictions on the time step only. It was shown that in the area of stability of some schemes, there are the values of the Courant number, for which the influence of "scheme dispersion" is minimal. Difference schemes for which the approximation error of dispersion is decreased only by grinding of the computational grid are identified.

Keywords: nonlinear dispersive equations, numerical algorithms, finite-difference methods, accuracy, stability, dissipation, dispersion.

Acknowledgements. The study was supported by the Russian Science Foundation (Project No. 14-17-00219).

Received 21 July 2015