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

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

Том 1, № 3, 1996

ОБ АЛГОРИТМЕ ЧИСЛЕННОГО РЕШЕНИЯ

УРАВНЕНИЙ ОДНОЙ НЕЛИНЕЙНО-ДИСПЕРСИОННОЙ МОДЕЛИ

_ __и __I

МЕЛКОЙ ВОДЫ**

В. Б. Барахнин, Г. С. Хлкимзянов Институт вычислительных технологий СО РАН Новосибирск, Россия

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

1. Введение

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

В численных исследованиях наиболее часто используемой приближенной моделью является модель мелкой воды первого приближения. Недостатком модели является то, что она дает достоверные результаты лишь для волн малой амплитуды. Работы [9, 12, 14, 16, 17, 19, 21, 24] посвящены исследованию нелинейно-дисперсионных моделей и алгоритмов,

*© В. Б. Барахнин, Г. С. Хакимзянов, 1996.

^ Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №97-01-00819.

основанных на этих моделях. Из множества нелинейно-дисперсионных моделей (см., напр., [1, 2, 6, 18, 20, 22]) особый интерес вызывает модель Железняка — Пелиновского, так как при ее получении не используется предположение о малости амплитуды волн. Эта модель представляет собой систему нелинейных уравнений относительно возвышения свободной поверхности и скорости. Получается она путем разложения основных гидродинамических функций в ряд по параметру дисперсии в = (к0/Ь)2 и параметру нелинейности а = а0/Н0, подстановкой этого разложения в трехмерные уравнения, описывающие потенциальные течения жидкости со свободной границей, и последующим отбрасыванием членов порядка О (в2) с сохранением членов порядка 0(атв) (см., напр., [10]). Здесь Ь — характерный размер по горизонтали, к0, а0 — соответственно характерные глубина и амплитуда. Вывод, изучение и применение этой модели даны в [5, 7, 8]. Ее одномерный аналог получен также в работе [23].

Численная реализация многих нелинейно-дисперсионных моделей осложняется наличием в соответствующих уравнениях производных высокого порядка искомых функций, например смешанных производных по времени и пространственным переменным. Поэтому в настоящее время при построении разностных схем для нелинейно-дисперсионных уравнений предпочтение отдается алгоритмам, основанным на решении эллиптических и гиперболических уравнений, получаемых тем или иным способом из исходных дифференциальных уравнений. Эта идея применялась, например, в работах [9, 12, 16, 21]. Численная реализация таких алгоритмов сводится к тому, что на каждом временном шаге сначала решают эллиптические уравнения, а потом — неоднородную гиперболическую систему с правой частью, зависящей только от производных по пространственным переменным. Такой подход позволяет строить эффективные алгоритмы, так как для решения эллиптических и гиперболических уравнений имеются хорошо зарекомендовавшие себя численные методы.

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

2. Математическая постановка

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

Н + а1у(Я и) = 0, и + (и -У)и + Уп = Б, (1.1)

где

1 / Н з Н2 \ (Н

в = Ну ( — ^ + — яЛ -уы у+ п2

Я1 = (ё1у и)4 + (и ■ у)а1у и - (а1у и)2, Я2 = щ -УН + (и ■ у)(и ■ уН),

П = П(г) — меняющаяся, вообще говоря, со временем односвязная ограниченная область в плоскости декартовых координат х1, х2, иа(х1, х2, г) — декартовы компоненты вектора скорости и (а = 1, 2), п(х1, х2, г) — возвышение поверхности над невозмущенным уровнем, Н = п + Н — полная глубина. Обезразмеривание проводилось по формулам

ха = хако, п = пНо, Н = ННо, иа = иал/ф0, г = гл/ко/д,

где Н0 — характерная глубина, д — ускорение свободного падения, символом "~"обозначены размерные величины.

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

что

ди

&у[(и ■ у)и] = (и ■ у)а1у и + (а1у и2 - 2 ае^д"),

где

а / ди\ ди1 ди2 ди1 ди2

-дх) дх1 дх2 дх2дх1 С учетом этого

д и

я1 = а1уа - 2(а1у и)2 + 2 ае^д^), я2 = а -ун + и ■ ((и -у)ун). (1.2)

Здесь а — ускорение частиц жидкости, а = и4 + (и ■ у)и.

Переписав уравнения движения системы (1.1) в виде а = Ю - уп, получим

1 / Н3 Н 2 \ (Н \

а = Ну ( — Я1 + — яП - ум -Я1 + яП - Уп.

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

имеем

где

Тем самым

1 Н2

-уп = -уН + уН = -НУ (— ) + уН

На = Уф -Унф, (1.3)

Н3 Н2 Н2 Н2

Ф = —Я1 + —Я2--, ф = —Я1 + НЯ2 - Н.

3 1 2 2 2 2 1 2

12ф - 6Нф -6ф + 4Нф + Н2 Я1 = --, Я2 =

Н3 Н2

Используя равенства (1.2) и (1.3), мы получим систему уравнений для неизвестных ф и ф:

' уф УНф \ 2 /ди \ 12ф - 6Нф

-2(а-уи)2 + 2ае<1х) * ■

(^ - т \ .уН + и. ((и .у)уН) = -6ф + %Ф + Н 2. (1.4)

Из второго уравнения системы (1.4) следует, что

Н УН -Уф + 6ф + Н 2и ■ ((и ■ у)уН) - Н2

Ф

Нг

где г = (уН)2 + 4. Подставляя выражение для ф в первое уравнение системы (1.4), после очевидных преобразований получим уравнение для ф:

• уф_улул^уЛ - / /у.

Н Н г \ \Н2г

12 г -

+ Н-з — ]Ф

(1.5)

где

^ = а1у

(и ■ ((и ■ у)уН) - 1)уН\ 6 и ■ ((и ■ у)уН) - 1

+2<а1у и)2 - 2 <4 I)-

Н

Предполагая, что Н(х1,х2,Ь) > Н0 > 0, нетрудно показать, что уравнение (1.5) является эллиптическим.

С учетом (1.3) уравнения системы (1.1) запишутся в виде

Н + а1у(Н и) = о,

и + (и ■ у)и = —(Уф - уНф). Н

(1-6)

Итак, решение системы (1.1) свелось к решению эллиптического уравнения (1.5) и системы (1.6), правая часть которой не содержит производных по времени. Уравнения (1.6) можно записать и в дивергентной форме:

д V ЗЕ1 д Е2 дЬ дх1 дх2

О,

х

(х1,х2) е П.

(1.7)

Здесь

V

Н

Ни1 Ни 2

Е1

Ни1 Ни\

Ни1и2

Е2

Ни2 Ни1и2

и22

О

Ни2

0

фх 1 - фНх 1

фх2 - фНх2

Рассматриваемые уравнения дополняются начальными и краевыми условиями. Например, на непроницаемых неподвижных участках Г0 имеем

и • п

0,

(1.8)

Го

где п — единичный вектор внешней нормали к границе.

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

в • п

к (и ■ и),

(1.9)

Го

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

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

Из формул (1.3) и (1.9) имеем следующее краевое условие на функцию ф:

дф дк Н У к •Уф + 6ф дп дп Нг

дк Ни • ((и • у)ук) — Н

Го

дп

+ кН (и • и).

(1.10)

В частном случае, когда функция к, описывающая рельеф дна, удовлетворяет на Го равенству У к = — |ук| п, для уравнения (1.5) возникает третья краевая задача с граничным условием

Ж + ЗНф) = Н (|уК|( — (и • у)(и • уК) + 1) + 4к (и • и)) .

' Го

Если же дно плоское и горизонтальное, к = 1, то краевое условие (1.10) запишется в виде условия Неймана

дф = кН (и • и). (1.11)

Го

дп

Уравнение (1.5) примет в этом случае более простой вид:

д ( 1 дф) д ( 1 дф ) 3

где

дхН Ндх1) + дх2\ Ндх2) Н3 ф '

^ = 2(^ и)2 — 2 Ц £) + А.

(1.12)

3. Постановка задачи в криволинейных координатах

Пусть преобразование координат

г = т' ха = ха(д\д2,т) а = 1, 2 (2.1)

является достаточно гладким, взаимно-однозначным и невырожденным (для определенности, с положительным якобианом), причем область П(г) при этом преобразовании является образом единичного квадрата Q = (0; 1) х (0; 1). В новых координатах система (1.7) имеет вид

дУ дЕ1 д Е2 . 1 2. „

ж + ац + ^ = О д= (^2) е Q'

где

У

Н.1

Ни.

ре

/ Ньа3 Нщиа. \Ни2иа.];

О

0

(2.2) \

дф ,дк\ дх2 ) дф \дд2 дф2) дд1

дх1

( дф ,дк \ дх2

удд1 дд1) дд2

/ дф ^дк \ дх1 + / дф ^дк\ удд1 дд1) дд2 \дд2 дд2) дд1)

.] — якобиан преобразования (2.1), уа = — контравариантные компоненты вектора скорости:

дх1^ дда ( дх2\ дда щ--— —г + щ —

дг дх1

дда ( —1)а+в дх3-в

дг дх2'

дхв

. дд3-

а

Нам потребуется также недивергентная форма записи уравнений (2.2)

д V л1 д\ л2 д V

--+ А1—+ А —

дЬ дд1 дд2

д V дv 42 д V .

+ А1— + А2^ = f, (2.3)

где

V = | «1 | , А

П2

( уа Едда Едда \

дх1 дх2

0 уа 0

V 0 0 уа /

1

ЕЗ

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

д ( др др \ д ( др др \

\к11 + к12 дф) + к12 дд! + к22 Эф) — ^ = 3 (2^

здесь

= 3 /_ав _а7 дв6 дк дк\ д / 3_ав 8к\ 123 г - 3

кав = г дд^ д^у' ко = 6 д^у ЯГ д^у + я3 ' ^ = 1А (£ З_ав А V А + А Аз (V - дда V 2 det( ди \

З дда у г ддв) Ег З2 дд« \ дЬ ) 3 \дд)'

Г = _ дд" ддв + 4 £ = ^ — дЬ ) дд*{дду Эх?) — 1 _ = 32 9з-а3-в'

9а@ — ковариантные компоненты метрического тензора; во всех формулах по повторяющимся индексам ведется суммирование.

В криволинейных координатах условие непротекания (1.8) примет вид

Отсюда сразу следует, что

уа I = уа I =0.

\да=0 1да = 1

™ I = ™ I =о

\да=0 \да = 1

4. Конечно-разностный алгоритм

Для численного решения уравнений (2.2), (2.4) замыкание области О (квадрат О) покрывается прямоугольной сеткой с количеством узлов Ыа по направлению оси Ода и шагами ка = 1/(Д^ — 1).

Алгоритм расчета на произвольном шаге по времени с номером п состоит в следующем. Сначала решается конечно-разностное уравнение для нахождения функции р, аппроксимирующее уравнение эллиптического типа (2.4), причем р вычисляется в целочисленных узлах сетки. Для решения системы (2.2) нами использовалась явная конечно-разностная схема предиктор—корректор с автоматически настраиваемой аппроксимационной вязкостью. На шаге предиктор этой схемы аппроксимируется недивергентная форма системы (2.2), то есть система (2.3), правая часть которой содержит уже найденные значения р и ф. После этого вновь численно решается уравнение (2.4) с использованием величин Е, и1} и2,

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

Схема предиктор—корректор описана в нашей работе [15], поэтому здесь мы опишем лишь схему, использованную для численного решения эллиптического уравнения (2.4). Конечно-разностное уравнение получено на основе интегроинтерполяционного метода.

Рис. 1. Контур интегрирования и шаблон разностной схемы.

Проинтегрируем уравнение (2.4) по прямоугольнику ABCD (рис. 1), стороны которого параллельны координатным осям вычислительной области Q и делят расстояния между соседними узлами пополам. Имеем

Í(kii 1т + к12й) dq2 - í(к11 TTI + к12й) dq2 + í(k2l+ k22^2

J \ dql dq2) J \ dql dq2) J \ dql dq2 /

(BC) (AD) (DC)

-/ (к21 di + к22д^ dq1 - JJко pdqldq2 = jjFJdqldq2. (3.1)

(AB) ABCD ABCD

Коэффициенты кар вычисляются в центрах ячеек сетки — в точках qil+l/2,i2+l/2. Для подсчета интегралов в (3.1) по сторонам прямоугольника ABCD применяется формула трапеций. Тогда из (3.1) получим следующее разностное уравнение во внутреннем узле

qil,i2 • _ _

= (Dqi Fi + Dq2F2 - ко)&1,ъ = (FJ)ii,i2, (3.2)

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

Dq1 f (qii + l/2,i2 + l/2)

fil + l,i2 + l + fil + l,i2 fil,i2 + l fil,i2

2h

l

гл rr \ fil + l,i2 + l + fil,i2 + l fil + l,i2 fil,i2

Dq2 J (qil + l/2,i2 + l/2) = -2h¡-

Flf (qil+l/2,i2+l/2 ) = (кllDql f + h2 Dq2 f) (qil+l/2,i2+l/2),

F2f (qil+l/2,i2+l/2 ) = (кl2Dql f + к22 Dq2 f) (fe+l/2, i2+l/2) ,

y F ( ) Dql Fl (N)+ Dql Fl(S) D D q2 F2(E)+ Dq2 F2(W)

Y Fl(qil,i2) =-2-' Dq2 F2 (qil,i2) =-2-'

Далее будем рассматривать случай плоского горизонтального дна, к = 1.

Рис. 2. Нумерация типов узлов. Рис. 3. Контур интегрирования и шаблон на

нижней границе.

Опишем аппроксимацию граничного условия Неймана (1.11). Заметим, что на нижней

границе вычислительной области п = —■— (х21' — хМ. На этой границе выполняется

л/Зп

равенство

дф дф 1 ( 2 (дф х\г дф ) 1 ( дф дф х^1 -

, дф . , дф _ 1 ( 2 ( дф хЯ2 дфхд1 ) . 1 ( дфхд2 дфхч1 )Л

к12 ОО1 + к22 "до2 = Н\ — ^ 1№. — №. + ^ ^ — +

1 ( 2 дф . 1 У3^ дф .- / л

н I—^ дх1+^ охо = = ^1К1 (и •и)' (3.3)

где, согласно формуле подсчета кривизны [11],

2 1 1 2 ГУ*^ ГУ* _ ГУ* ГУ

К1 = -

3/2

911

На верхней границе области Q также имеет место равенство вида (3.3). Наконец, для левой и правой границы области Q получим следующую формулу:

к11 дф + к12ОФ = ^к2 (и • и)' (3.4)

где

1 2 2 1 гу -1- гр^ _ гр^ гу 1-

¿Ь О О ¿Ь О О О Ж О

Я2 Я2 Я2 Я2 Я2 Я2

К2 = 3/2 .

9 22

Приведем конечно-разностные уравнения в граничных узлах, разбив множество узлов на типы в зависимости от положения узла на сетке (рис. 2). Шаблон и контур интегрирования для узлов, принадлежащих нижней границе (тип 2), изображены на рис. 3. Рассматривая интегралы вида (3.1) и используя равенства (3.3) и (3.4), получим следующие разностные уравнения (черта над оператором Еа означает, что берется среднее арифметическое значений оператора в двух соседних ячейках; к0 берется из центров ближайших ячеек). Тип 1:

2_ 2

Лф1,12 = — Е1ф(Е) + Д*¥2ф(Е) — коф1,г2 = + (—^к2(и • и)) ^. (3.5)

Тип 2:

2_ , 2 )

Лфги1 = Бд1 Лф(#) + —) — кофч,1 = (^.)гь1 + (—^9Ик1 (и • и))г1 1. (3.6)

к2 к2 1'

Тип 3:

2

ЛфМ1,г2 = —Т~ Р1ф(Ш) + Дд2 ) — кофМ1 ,г 2

к1

(Р.)^,г2 (к л/~922К2 (и • и)) ^1,г2 .

(3.7)

Тип 4:

2__( 2 )

Лфг 1,М2 = Д л) — — Р2ф(Б) — кофг= )г1,^2 — (т-К1 (и • и))

к2 к2

(3.8)

Тип 5:

Тип 6:

Тип 7:

Тип 8:

22 Лф1,1 = т Р1ф(°) + Т Р2ф(С) — кофц = к1 к2

2 2 = )1,1 + (—-^922к2(и • и)) м + (h^v/911к1(u • и))IX

22

Лфм1,1 = — т- ^ф(Д) + — ^2ф(Д) — коф^д = к1 к2

2 2 — (и • и))+ (—\/9пК1(и • и))

22

ЛфмъМ2 = — Г" ^\ф(А) — — ^ф(А) — кофм1,М2 = к1 к2

2 2 — (к^\/922к2(и • и)) — (к^\/9^к1(и • и))

22

Л ф 1, N2 = т- ^1ф(£) — — Р2ф(Б) — коф1,м2 = к1 к2

2 2

)1,М2 + (и • и)) ^ — (и • и)) 1,м2.

(3.9)

(3.10)

(3.11)

(3.12)

Таким образом, формулы (3.2), (3.5)—(3.12) полностью определяют оператор Л и систему конечно-разностных уравнений

Лфг 1,г2 рг 1,г2' ^а 1'...' Nа.

(3.13)

Коэффициенты этих уравнений приведены в [3]. Нетрудно убедиться, что для оператора Лапласа, аппроксимируемого на квадратной сетке, описанная схема переходит в схему "косой крест".

Докажем самосопряженность и положительную определенность оператора А = —Л. Оператор А действует из М-пространства сеточных функций, заданных на сетке Qh, в М. Введем в М скалярное произведение:

N1-1 N2-1

1

N2-1

(ф/Ф) = к1к2 ^ ^ (фф) г1, г 2 + 2 к1к2^2 (фф)1,г2 +

г1=2 г2=2

г2=2

1 N1-1 1 N2-1 1 N1 -1

¿1=2 ¿2=2 ¿1 =2 + 1 -1-2(фф)1,1 + 4-1-2(фф)^,1 + 4-1-2(фф)^,^ + 1 -1-2(фф)1^2 •

Пусть f Е М. Определим f Е М так: /г1гг2 = /г1гг2 (ко)п,г2• Имеет место следующее утверждение.

Лемма 3.1. Для любых функций ф, ф из пространства М имеет место равенство

(Аф,ф) = -1-2^ ^ (^фДд ф

+1/2) + (Ф,Ф). (3.14)

¿1 = 1 ¿2 = 1

Доказательство. Преобразуем слагаемые, входящие в выражение (Аф,ф). Преобразование зависит от типа узла, к которому относится соответствующее слагаемое.

Рассматривая выражение (3.2), используемое в узлах типа 0 и взятое со знаком " —", подробно запишем преобразование тех слагаемых скалярного произведения, в которые входит сомножитель Рд1 Р1ф:

N1-1 N2-1 , , N1 -1 N2-1

-1-2

-1-2X X (РЯ1 Р1ф ф)^2 = т^ X X] ( — Р1ф(^1+ 1 / 2 ,¿2 + 1 / 2)Ф ¿1 ,¿2 +

д1

г1=2 г2=2 1 г1 =2 г2=2

+ ^ 1ф(Ц%1-1/2,12 + 1/2)ф ¿1,12 — Р1ф(Я%1+1/2,12-1/2 У^М + ^1ф(Ч%1-1 / 2,%2-1 / 2)ф %1,%2

= — X X р1 ф^+Ф^+Ф^Ч^ + X X РМЧ+Ф^+1/2^+1,2— 1 ¿1=2 ¿2=2 ¿1 = 1 ¿2=2

N1-1 N2-2 N1-2 N2-2

— X X р1ф(^1+1/2,02+1/2^1,02+1 + X X FlФ(q¿,1+lA¿,2+l/2)ф¿,1+M,2+l) =

¿1=2 ¿2 = 1 ¿1 = 1 %'2 = 1

N1 -2 Щ-2 — — N^2-2

= -1-2^ ^ {Р1фОд1ф) + РМ^+Ф^^^^^ + ф1+1,¿2)+

¿1 =2 ¿2=2 1 ¿2=2 N1-2 N2-2

+ X ^^¿1+1/2,1+1/2)^¿1+1,1+1— ^1,1+1) +^2 Р^ш-Фъ+Фк — фNl-l,¿2+l— фт-1,ъ) +

¿1=2 ¿2=2

N^2

+ X FlФ(Q¿l+1/2,N2-1/2) ^¿1 + 1^2-1 — 2-1) + Р1ф(д1+1/2,1+1/2)ф1+1,1+1—

¿1=2

— Р1ф(Я^-1/2,1+1/2)ф^-1,1+1 — Р1ф(Я^-1/2,^-1/2)ф^-1^2-1 + Р1ф(Я1+1/2^2-1 /2)ф 1+1^2-^ • Осуществив аналогичные преобразования для выражения

-1-2 Е X (Р

д2

¿1=2 ¿2=2

полностью определим вклад узлов типа 0 в скалярное произведение (Аф,ф).

Тип 0:

N1-1 N2-1 N1-2 N2-2

к1к2 X ^(АфФ)»1,г2 = к1 к2^ X (Р1фДд1 ф + ^фДд2ф)(дг 1+1/2^+1/2) +

г1=2 г2=2 г1=2 г2=2

и и N2-2 N1-2

+ кф2 ■

( X] Р1ф(д1+1/2М+1/2){ф1+1,г2+1 + ф1+1,г2) + X Р1ф(дг1+1/2,1+1/2) (^1+1,1+1 — фгъ1+1) + г2=2 г1=2

N2-2

+ X ^1ф(д^-1/2,г2+1/2К — ф^-М2+1 — ф^-1,г2) + 2=2 N1 —2

+ X Р1ф(дг1+1/2^2-1/2){фг1+1^2-1 — фг ь^-1) + 1 =2

ф(д1+1/2,1+1/2)ф1+1,1+1 — р1ф(д^-1/2,1+1/2)Ф^-1,1+1 — — Р1 ф(д^-1/2^2-1/2)ф^-1^2-1 + Р1ф(д1+1/2^2-1/2)ф1+1^2-^ +

к1к2 ( N2 2

+^1~2( X! ^2ф(д1+1/2,г2+1/2){ф1+1,г2+1— ф1+1,г2) + X ^ф(дг 1+1/2,1+1/2^фг 1+1 ,1+1 +Ф 1, 1+1 +

2 г2=2 г1=2

N2-2

-1/2, 2 +1/2К фNl-l, г 2+1 — фNl-l, г 2) +

2 =2 N^2

+ У] ^2ф(дг 1+1/2^2-1/2К — фг 1+1^2-1 — фг-1) +

1=2

+ ^2ф(д1+1/2,1+1/2)ф1+1,1+1 + ^2ф(д^-1/2,1+1/2)ф^-1,1+1 — ^2ф(д^-1/2^2-1/2)ф^-1^2-1—

N1 —1 N2-1

— Р2ф(д1+1/2^2-1/2)ф1+1^2-1^ + к1к2^ X (ф ф)г1,г2. (3.15)

1=2 2=2

Для других типов узлов будем иметь следующие равенства. Тип 1:

Ъ Ъ N2 1 Ъ Ъ ■ N2 2 к1к2 ^^ , . , ч к1к2

Х(Аф Ф)1,г2 = "¿"К X ^ф^Аг2+1^ ( — ФМ2 + 1 — ^2) — ^1ф(д1+1/2,1+1/2)ф1,1+1-

__О 1 __о

2

г2=2 1 г2=2

) к1к2 ( N2-2

—Г1ф(д1+1/2,щ-1/2)ф1,щ-1) + "¿""Ч X ^2ф(д1+1/2,г2+1/2Кф1,г2+1 — фМ2) +

\ кгк2 N2 1

+ Р2ф(д1+1/2,1+1/2)ф1,1+1 — ^2ф(д1+1/2^2-1/2)ф1^2-1 ) +--— 2 £ (ф ф)^ . (3.16)

2 2=2

Тип 2:

Ъ Ъ ^ Ъ Ъ ■ N1-2

к1к2 , . к1к2

X (Аф Ф)г 1,1 = -2г1(£ ^ф(д +1/2,1+1/2)(фг 1+1,1 — фг 1,1) + ^1ф(д1+1/2,1+1/2)ф1+1,1-

2 ^ч-г г,^1 2к1

г1=2 1 г1=2

16 В. Б. Барахнин, Г. С. Хакимзянов

-F1<<(qNl-1/2,1+1/2^Nl-1,1) + Ц2 ( X F2<(qil+1/2,1+1/2)( - фil+1,1 - Фчл)-

2<(qil+1/2,1+1/2) \ — Wil+1,1 — Will,

2

i1=2

Nl - 1

— F2<(ql+1/2,1+1/2)ф 1+1,1 - F2<(qNl-1/2,1+1/2)ф^-1,1 (< ф)и,1• (3.17)

hlh2 J

il=2

Тип 3:

h h N2-1

(A< ^Nl,i2 =

2

i2=2

7 7 N2-2 h1h2

( X F1<(qNl-1/2,i2+1/2)(ф^^+1 + + F1 <(qNl-1/2,1+1/2)ф^,1+1 +

2h1 1 i2 =2

+ ^1<^1-1/2^2-1/2)ф^^2-1) + ¡¡¡^ ^ F2<(qNl-1/2,i2+1/2)(фт,Ь + 1 — +

2 i2=2

Л hlh2 N2 1

+ F2<(qNl-1/2l1+1/2УЫ11+1 — F2<(qNl-1/2,N2-1/2Yhl,N2-l) + X (< ф^2 • (3-18)

hlh2 Y

i2=2 Тип 4:

h h N1-1

1 2 X (A< =

2

il=2

7 7 N1-2 hlh2

X F1<(qil+1/2,N2-1/2 ){фп + 1^2 — Фп ,N2) — F1<(qNl-1/2^2-1/2)ф^-1^2 +

1

2h

il=2

7 7 . N1-2 h1h2

+ ^1<^1+1/2^2-1/2)ф1+1,^) + ¡¡¡~( X F2<(qil + 1/2,N2-l/L^il+lN + ^,N2) +

2 il =2

Л hlh2 Nl-1

+ F2<(qNl-1/2,N2-1/2)фNl-1,N2 + ^2<^1+1/2,^-1/2)ф1+1^2 J + ^ (< Ф)il,N2 • (3Л9)

Тип 5:

¡a (A< ф)1,1 = — ¡¡iF1Ф1+1/2,1+1/2)ф1,1 — ¡¡Ц2F2Ф1+1/2,1+1/2)ф1,1 + ¡i¡2 (< ф)1,1•

1 2 (3.20)

Тип 6:

(A< ^N1,1 = F1<(qNl-1/2,1+1/2)ф^,1 — ¡¡L Р"2<<^1-1/2,1+1/2)ф^,1 + ¡А № ф)^,1.

1 2 (3.21)

Тип Т:

(A< ^NlN =

= ¡¡1 ^<<^1-1/2^2-1/2)ф^,Щ + ¡¡¡2 ^<<Ы1-1/2,Щ-1/2)ф^,Щ + ((Р ^NlN . (3.22)

Тип 8: ^^

(аф Ф)1,м2 =

= - Г1'АЯ1+1/2,М2-1/2)Ф1,М2 + Е2ф(0-1+1/2,Щ-1/2 )Ф1,М2 + (ф Ф)1,М2 • (3-23)

Складывая равенства (3.15)-(3.23) и приводя подобные слагаемые в правой части, получим требуемое равенство (3.14).

Лемма 3.2. Для любых функций ф, ф из пространства М имеет место равенство

(Р1фБд1 Ф + Е^фВ д2 ф) (дг1+1/2,г2 + 1/2)= (^ф^! ф + Е2фОд2 ф) ^ + 1/2^+1/2)-

Доказательство. Из определения операторов Еа имеем

^фБд! ф + Е2 фОд2 ф) (^+1/2,г2+1/2) =

= ((к 11Од1 ф + киОд2 ф)Од1 ф + (кмО д1 ф + к22Од2 ф)Од2

ф) + 1/2,12 +1/2) =

= {(киВд1 ф + ^Од2 ф)Од1 ф + (к^Дд ф + к22 Од2 ф) О д2 ф) (^+1/2^2 + 1/2) =

= (Е1фБд1 ф + Е2 гфВд2 ф) (^+1/2,г2+1/2)-

Теорема 3.1. Оператор А является самосопряженным и положительно определенным в пространстве М, причем имеет место оценка

^ (ф,ф) < (Аф,ф) < (С21 + ^)(ф,ф),

где

. 3/ 3/ /44

Й1 = Ш1П ——, ^ = шах ——■, I = шах ,

дед Н3 дед Н3 ^ ^2

кц + к22 + \/(кц - к22)2 + 4к1

2 12

С2 = шах-г

дед 2

Доказательство. Согласно леммам 3.1 и 3.2

(Аф,ф) = ^2 ^ {Е1ф°д1ф + Е2фОд2ф) (Я11 + 1/2^2+1/2) + (ф, ф) =

¿1 = 1 ¿2 = 1

= ^2 X X {Е1фОд1 ф + Е2фОд2ф)(5г1+1/2^2+1/2) + (ф/',ф) = (Aф,ф),

¿1 = 1 ¿2 = 1

что доказывает самосопряженность оператора.

Рассмотрим в области ф дифференциальный оператор

2 5 / 5

ь = X к

у •

а,в=1 4 4 ^ 7 Нетрудно показать, что для произвольных (2 имеет место оценка

2

*(<? + С22) < x кав(а(в < С2(С? + С22),

а,в=1

. kn + k22 - л/D kn + k22 + VD —

где c1 = min---, c2 = max---, D = (kn — k22)2 + 4k22. Если в Q

q&Q 2 qgQ 2

выполняются неравенства kn > 0, k11k22 > k^, то L является равномерно эллиптическим оператором и выполняется неравенство С\ > 0. Используя определения операторов Fa, имеем

Ni — l N2 — 1

(Аф, p) = hl h2^ X {F1^Dql V + F2^Dq2 tp) (qii + 1/2,i2 + 1/2) + (p ф) =

il = 1 i2 = 1

Ni — l N2 — 1

hih2 £ Ё ((k 11Dql P + k12 Dq2 p)Dql p + (ki2D ql p + k22Dq2 p)Dq2 P (Qii + 1/2,i2 +1/2) + (p,P),

il = 1 i2 = 1 откуда следует, что

Nl — 1 N2 — 1

Clhlh2^ Ё ((Dql P)2 + (Dq2 P)2) (q il + 1/2,i2 +i/2) + di(p,p) < (Ap,p) <

il = 1 i2 = 1

Ni — 1 N2 — 1

< C2hi h2^ ((Dql p)2 + (Dq2 p)2) (q il + 1/2,i2 +1/2) + d2(p,p). (3.24)

il = 1 i2 = 1

Рассмотрим в области Q следующую задачу Неймана:

—-p=о, dp

= 0.

г

Аппроксимируем эту задачу на прямоугольной сетке с шагами Н1 и к2 с помощью конечно-разностного оператора В, получаемого из оператора А путем замены коэффициентов на постоянные: к11 = к22 = 1, к12 = 0, к0 = 0. Очевидно,

Nl — 1 N2 — 1

(Bp,p) = hih2 £ £ ((D

qi

p)2 + (Dq2p)2^j (qii+l/2,i2+l/2).

il = 1 i2 = 1

Собственные значения оператора B имеют вид

Л-nl,n2 = -Щ (1 — COS nhiUi)(l + cos nh2U2) + JJ (1 + COS nhiUi)(l — cos nh2U2), где na = 0, ..., Na — 1, при этом

( 4 4

Amin = 0, Amax = l = max ( -jü , -j2

то есть

0 < (Bp, p) < l (p,p).

Таким образом, мы можем записать неравенство (3.24) в виде

Ci(Bp, p) + di(p, p) < (Ap, p) < C2(Bp, p) + d2(p, p)

и, окончательно,

( (ф, ф) < (Аф, ф) < (С21 + (12)(ф, ф),

что и требовалось доказать.

Разработанный алгоритм был применен для численного исследования процесса взаимодействия уединенной волны с плоской вертикальной стенкой, расположенной под некоторым углом к фронту волны, а также для моделирование течения, возникающего при разрушении дамбы, перегораживающей водоем [3]. Полученные результаты сравнивались с результатами расчетов по модели мелкой воды первого приближения и по модели потенциальных течений [4, 13, 15]. На основании проведенного сравнения можно сделать вывод о том, что разработанный алгоритм решения нелинейно-дисперсионных уравнений в достаточной мере продемонстрировал свою работоспособность, надежность и экономичность.

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

[1] АлЕшков Ю. З. Теория взаимодействия волн с преградами. ЛГУ, Л., 1990.

[2] Базденков С. В., Морозов Н. И., Погуцце О. Р. Дисперсионные эффекты в двумерной гидродинамике. Докл. АН СССР, 293, 1987, 819-822.

[3] Барахнин В. Б., ХАкимзянов Г. С. Численная реализация условий непротекания для одной нелинейно-дисперсионной модели мелкой воды. Актуальные проблемы современной математики, НИИ МИОО НГУ, Новосибирск, 3, 1997, 3-13.

[4] Барахнин В. Б., ХАкимзянов Г. С., Чубаров Л. Б., Шкуропацкий Д. А. Некоторые проблемы численного моделирования волновых режимов в огражденных акваториях. В "Вычислительные технологии", ИВТ СО РАН, Новосибирск, 1, №2, 1996, 3-25.

[5] Вольцингер Н.Е., Клеванный К. А., ПЕлиновский Е. Н. Длинноволновая динамика прибрежной зоны. Гидрометеоиздат, Л., 1989.

[6] Дорфман А. А., Яговдик Г. И. Уравнения приближенной нелинейно-дисперсионной теории длинных гравитационных волн, возбуждаемых перемещениями дна и распространяющихся в бассейне переменной глубины. Числен, методы мех. сплошной среды, 8, №1, 1977, 36-48.

[7] Железняк М. И. Воздействие длинных волн на сплошные вертикальные преграды. В "Накат цунами на берег", ИПФ АН СССР, Горький, 1985, 122-139.

[8] Железняк М. И., ПЕлиновский Е. Н. Физико-математические модели наката цунами на берег. Там же, 8-33.

[9] Марчук Ан. Г., Чубаров Л. Б., Шокин Ю. И. Численное моделирование волн цунами. Наука, Новосибирск, 1983.

[10] Овсянников Л. В., Макаренко Н. И., Налимов В. И. и др. Нелинейные проблемы теории поверхностных и внутренних волн. Наука, Новосибирск, 1985.

[11] ПОГОРЕЛОВ А. В. Дифференциальная геометрия. Наука, М., 1974.

[12] ФЕДОТОВА З. И. О свойствах разностных схем для длинноволновых приближений уравнений гидродинамики. В "Вычислительные технологии", ИВТ СО РАН, Новосибирск, 2, №7, 1993, 237-249.

[13] ХАКИМЗЯНОВ Г. С. О численном моделировании на адаптивных сетках трехмерных течений жидкости с поверхностными волнами. В "Тр. Всесоюзн. совещ. по числ. мет. волн, гидродин", Ростов-на-Дону, 1990, ВЦ СО АН СССР, Красноярск, 1991, 103-108.

[14] Шокин Ю. И., Чубаров Л. Б., МАРЧУК Ан. Г., Симонов К. В. Вычислительный эксперимент в проблеме цунами. Наука, Новосибирск, 1989.

[15] Barakhnin V. B., Khakimzyanov G. S. Adaptive-grid numerical solution of one-dimensional and two-dimensional problems for the shallow-water equations. In "Advanced Mathematics: Comput. and Appl. Proc. of AMCA-95", Novosibirsk, 1995, 144-153.

[16] Barakhnin V. B., Khakimzyanov G. S. On the application of adaptive grids to the numerical solution of one-dimensional problems in the shallow-water theory. Russ. J. Numer. Anal. Math. Modelling, 10, №5, 1995, 373-391.

[17] Eilbek J. C., McGuire G. R. Numerical study of the regularized long-wave equations, I. Numerical methods. J. Comput. Phys. 19, №1, 1975, 43-57.

[18] Ertekin R. C., Webster W. C., WEHAUSEN J.V. Waves caused by a moving disturbance in a shallow channal of finite width. J. Fluid Mech., 169, 1986, 275-292.

[19] Fedotova Z.I., Pashkova V. Yu. On the numerical modelling of the dynamics of weakly nonlinear waves with dispersion. Russ. J. Numer. Anal. Math. Modelling, 10, №5, 1995, 407-424.

[20] Green A. E., Naghdi P. M. A derivation of propagation in water of variable depth J. Fluid Mech, 71, 1976, 237-246.

[21] KOMPANIETS L. A. Analysis of difference algorithms for nonlinear dispersive shallow water models. Russ. J. Numer. Anal. Math. Modelling, 11, №3, 1996, 205-221.

[22] Peregrine D. H. Long waves on a beach. J. Fluid Mech., 27, pt. 4, 1967, 815-827.

[23] SEABRA-SANTOS F.T., Renouard D. P., Temperville A.M. Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle. J. Fluid Mech., 176, 1987, 117-134.

[24] Shokin Yu. I., Khakimzyanov G.S., Chubarov L.B. New potentialities of computational experiment in tsunami problem. In "Proc. of the Int. Tsunami Symp., TSUNAMI'93", Wakayama, Japan, August 23-25, 1993, 277-284.

Поступила в редакцию 15 сентября 1996 г.