CC BY
122
Вычислительные технологии
Том 11, часть 2, Специальный выпуск, 2006
О ПРИМЕНЕНИИ РАЗНОСТНОЙ СХЕМЫ
МАК-КОРМАКА ДЛЯ ЗАДАЧ ДЛИННОВОЛНОВОЙ ГИДРОДИНАМИКИ*
3. И. Федотова
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: zf@ict.nsc.ru
The paper deals with an application of MacCormack finite-difference scheme for simulation of the dynamics of long wave in the framework of nonlinear shallow water theory. Conclusions on the effectiveness of this scheme for a number of problems have been made.
Введение
При решении ряда задач длинноволновой гидродинамики в рамках нелинейной теории мелкой воды [1-3], воспроизводящих трансформацию волн цунами при их распространении к побережью и взаимодействии с прибрежными структурами, использован численный алгоритм, базирующийся на конечно-разностной схеме Мак-Кормака. Эта двухшаговая схема, впервые описанная Робертом Мак-Кормаком в работе AIAA paper, 1969, N 69-354 (недавно полностью перепечатана в [4]), часто применяется для численного решения гиперболических уравнений ввиду хороших диссипативных и дисперсионных свойств, а также компактного шаблона и удобной реализации краевых условий. За период, прошедший после первых применений схемы Мак-Кормака, она претерпела много модификаций и обобщений. Так, существует неявный вариант этой схемы. Схема успешно применена для построения аппроксимации в методе конечных объемов. Предложены варианты расщепления как по геометрии, так и по физическим процессам. Для исследования свойств схемы, в том числе обусловленных нелинейностью аппроксимируемых уравнений, применен метод дифференциального приближения [7, 9]. Выполнено исследование точности разностной схемы, включая граничные условия, путем расчетов на последовательностях вложенных сеток, проведен сравнительный анализ с другими конечно-разностными схемами сквозного счета [1, 10]. Задача о моделировании динамики волн цунами в реальных обширных океанических акваториях имеет ряд особенностей (с точки зрения математика-вычислителя):
1. Расчетная область неодносвязная, неоднородная по глубине.
2. Внешняя и внутренние границы области могут быть сильно изрезанными.
*Работа выполнена при поддержке Президента РФ (грант № НШ-9886.2006.9), Российского фонда фундаментальных исследований (гранты № 5-05-64460, № 06-05-64869), Программы междисциплинарных интеграционных исследований СО РАН (грант № 2006-113).
@ Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
3, Расчет требуется выполнять на очень большие промежутки времени, волны проходят расстояния, в десять раз и более превосходящие их длину, в некоторых подобластях за это время происходит многократное взаимодействие волн,
4, В реальных задачах шаг сетки (сеточная батиметрия) является довольно крупным, причем этим параметром в расчетах вычислитель не управляет,
5, Положение (или выбор) границы расчетной области не является строго определенным из-за масштабов акватории,
6, Точная форма начального возмущения неизвестна.
Для преодоления этих проблем построенный на основе схемы Мак-Кормака алгоритм оказался вполне работоспособным, В работах Е, Елецкого, Л, Чубарова и автора данной статьи схема Мак-Кормака использована при создании компьютерной модели цунами [5],
1. Различные формы записи схемы Мак-Кормака
Впервые схема Мак-Кормака предложена для численного моделирования задач газовой динамики. Благодаря ряду свойств схемы (консервативность, умеренное проявление численной дисперсии и достаточность численной диссипации для подавления высокочастотных колебаний) мы стали применять ее для решения уравнений мелкой воды в областях, где профили волн становятся крутыми, В длинноволновой гидродинамике такие задачи возникают при описании трансформации волн на мелководье и их взаимодействии с берегом и прибрежными конструкциями.
Для вывода нужной модификации схемы Мак-Кормака рассмотрим интегральную форму уравнений мелкой воды. Если сделать допущение, составляющее сущность теории мелкой воды — малость вертикальной составляющей ускорения по сравнению с ускорением свободного падения д, то из общих законов сохранения импульса и массы несжимаемой жидкости для произвольно движущегося объема жидкости вытекают законы сохранения для уравнений мелкой воды в интегральной форме:
Wdxdy + Fdtdy + Gdtdx
<г{ш)
"(t)
Qdxdydt.
(1)
Здесь
W
Hu
Hv,
F
( Hu \
Hu2 + ^H2 Huv
G
Hv Huv
Hv2 + 9-H2 2
Ф = [ gHhx — ru(u2 + v2)1/2 KgHhy — rv(u2 + v2)1/2y
где x,y — точки на горизонтальной плоскости Oxy, t — время; u = u(x,y,t),v = v(x,y,t) —
x, y
ственно; H = H(x,y,t) — полная глубина жидкости; h = h(x,y,t) — глубина акватории,
rg
свободного падения. Произвольный объем жидкости u(t) ограничен поверхностью а(и).
Указанный способ вывода уравнений (1) реализован в работе [6]. Эти же уравнения можно получить исходя из дифференциальных уравнений мелкой воды
д д д + <2>
если проинтегрировать их по произвольному объему, а затем применить формулу Остроградского и перейти к элементам площади и объема в декартовой прямоугольной системе координат [1].
Рассмотрим конкретную задачу. Пусть область течения жидкости в плоскости Оху стационарна, Рассмотрим прямой цилиндр и = О х [0, £], где О — произвольная (вообще гово-
х, у
екцып водной акватории на горизонтальную плоскость. Предположение о неодносвязности соответствует тому, что мы учитываем наличие в акватории островов. Пусть в начальный момент времени £ = 0 заданы скорость жидкости и возмущение свободной поверхности, а границы области (или ее части) являются либо вертикальными непроницаемыми стенками, на которых задано условие полного отражения волны, либо "открытыми" границами, пропускающими без отражения подходящие к ним волны.
Для построения численного алгоритма перейдем к сеточному пространству. Для этого
0
дельными осям Ох Оу, и в последнем построим прямоугольную сетку с узлами (хт,уь). Множество узлов, попавших внутрь области О, обозначим Они образуют внутрен-
ность сеточной области. Выпишем аппроксимацию интеграла (1) по поверхности, ограничивающей контрольный объем — параллелепипед П^, = 0тЬ х [пД£, (п + 1)Д£], где Д£ —
шаг по времени, О^ — прямоугольные ячейки размером -(Джто +Джто+1) х (Ду& +Ду^),
Дхт = хт — хт-1, Дуь = уь — уь-1 и с центром в узле (хт, уь), целиком лежащие в области О:
х ,1 у, ,1 л У ,1 л х ,1
к+-2 , .^Л. t+At к+-2 t+At
ш
roДt ~ ± 1
х I у, I * г/, 1 1x1 к 2 П
™ — — ^ к— тт /г—— -¿.
X . 1 т+7
(Ыу + / в*
Уь+Ъ
(Ых = Ф* (Ых(у.
Уъ-1
(и+1)Ы
(х(у + I I Г *
X 1
т—
-5 т-5
(3)
Здесь хт+1 = (хт + хт+\)/2, ук+1 = (у& + ук+\)/2. Будем считать, что все стоящие под знаком интегралов функции постоянны на гранях прямоугольника Пть, и их значения будем связывать с центральным узлом соответствующих граней, используя общепринятые индексные обозначения. Перепишем интегральное соотношение (3) в следующей форме:
х 1 у, 1 ^ У, 1
т— -гг
л X ,1 л X ,1 У, 1
í+Дí í+Дí к+£
+ У У (СшкЦ~Сшк-ф<1Ых= / / / фтк<й<1х<1у-
^ X 1 I X 1 у 1
т— ^ т— ^ к—т^
Выбор величин, помеченных звездочками, определяет конкретную разностную схему.
(4)
В данной работе взяты следующие аппроксимации:
1 т+\к т+1к
С
1
тк-\-т,
— _ (ГШ+1 I Г<п А ~~ о \^гтк+1 Т ^тк))
Щ1 = ^(Шт+1), = о^х1), Фтк1 = ф(Штк1)
(5)
Величины, помеченные тильдой, вычисляются на шаге предиктор также путем аппроксимации интеграла (1) при использовании ячеек Отк размер ом Ахт х Аук:
Ук
Ук
(Шпк1 - wnk )йхйу+
(^тк - Рт-1к )(й(1у+
Хт-1 Ук-1 + / /
* Ук-1
Ук
(Стк 0тгк-1)(^(х
Фтк (Ых(у.
(6)
* хт—1 * хт—1 Ук— 1
На равномерной сетке (Ахт = Ах, Аук = Ау для всех к, т) из интегральных соотношений (4), (6) вытекает известная схема Мак-Кормака:
где
(Шп+1 - Ш^к)АхАу + (^ - )АгАу + (С;^ - С^—1)А*Ах = = Ф^к А1АхАу,
№ - + - Р^_1к)А1Ау + (С*тк+, - с*тк_,)АгАх
т—т^к'
Фт к АгАхАу,
(7)
Фтк
\
9Нт —1/2к
тк—1/2
0
А" - кп
тк т—1к п п
д^ Гт-1/2кит-1/2к
ьп - кп
птк птк-1 _ п п (, п Л2, (п
1 тк-1/2итк-1/2\\итк-1/2) Т ) у
—1/2к) — 1/2к) ) 2 ) 1/2
Ау
(|п+1 = к
\
дн п+1
¿> 1
к
+1к
- к
к
п+1
п+1
т+1/2к Дж
кщк+1 Ь"тк
тк+1/2 Д^
т+1/2к т+1/2к
(
-т+1 т+1 /2к ^
- т /2к ит + ,2к( (ит +..)2 + (С1^
)2)1/2
дн т+1
¿> 1
тк+1 /2 тк+1/2
)2 +(Дп+1 )2)1/2
тк+1/2) + (СУтк+1/2) ^
Здесь индексы т ± 1/2 к ± 1/2 означают, что помеченные ими величины представляют собой полусуммы соответствующих величин из соседних узлов, ктк+1 = ктк+1 (£ + А£/2); величины со звездочками определены по формулам (5),
Ф
решение дифференциальных уравнений, описывающее состояние покоя жидкости, является также точным решением разностной схемы,
В расчетах обычно используют более компактный вид схемы Мак-Кормака:
Шп+1
_ _ р™ ^ - — к Ах к —1к Ау
жпк -1
(отк - отк—1)+а^Фп
к
ТД/га+1 _ 4- ЛМп _ ^ ( _ _ 1 _ \ ,
^ял — тк 2Аж +1?с к 2Ау к+1 к 2 к '
*
х
гп
0
Из способа построения разностной схемы следует, что она консервативна. Особенно наглядна форма записи
^тк1 ~ ^тк = _ ^тк+§ ^тк-\
Д* Дж Ду
Это по существу формулировка метода предиктор-корректор, приведенная в [8], когда корректор обеспечивает выполнение законов сохранения,
В схеме (8) на шаге предиктор используются разности "назад", а на шаге корректор — "вперед". Разности "назад" и "вперед" можно циклически чередовать. Таким образом устраняется рассогласование, обусловленное аппроксимацией односторонними разностями.
Основные свойства разностной схемы Мак-Кормака хорошо известны. Она имеет второй порядок аппроксимации, условием устойчивости в одномерном случае является обычное условие Куранта, Изучению устойчивости схемы в двумерном случае посвящено много работ, что связано главным образом с неперестановочностью матриц А = (Е/(Ш и В = (С/(Ш. Подробный обзор по исследованию устойчивости этой схемы дан в [7, 10], К свойствам разностной схемы (8) мы обратимся позже, а сейчас перейдем к описанию границы и постановке граничных условий.
2. Аппроксимация граничных условий
Разностная схема Мак-Кормака имеет удобный шаблон с точки зрения реализации краевых условий. Выше мы определили множество внутренних узлов, образующих внутренность сеточной области. Определим граничные узлы сеточной области = О°гМ иГш- Так как в нашем случае область решения задачи стационарна, набор граничных точек также не зависит от временного слоя. Будем говорить, что (т, к) € Гп^, если он не является внутренним, но нужен для вычислений во внутренних точках. Последнее обстоятельство определяется шаблоном разностной схемы. Отметим (и это важно), что те узлы, через которые прошла граница дифференциальной задачи, автоматически входят во множество граничных узлов расчетной сетки, что нетрудно видеть из строения шаблона схемы Мак-Кормака, Внутренние узлы, в которых для вычислений по схеме потребуются граничные условия, будем называть приграничными узлами.
Учитывая специфику шаблона схемы (8), опишем для нее постановку граничных условий, применяя локально-одномерный подход, на примере задачи, когда граница представляется отражающей стенкой. Тогда единственным краевым условием, обеспечивающим корректность системы гиперболических уравнений (2), является условие ип = 0 на границе, где п — нормальная к границе компонента скорости. Для разностных схем, не являющихся схемами с направленными разностями, требуется по крайней мере еще одно дополнительное граничное условие, которое, как правило, берут в виде д(Н — Н)/дп = 0, Это равенство легко вытекает из уравнений движения мелкой воды, если предположить,
ип = 0
Пусть необходимо вычислить значение Ш в узле (т, к) и пусть (т — 1, к) € ГёГ;а. Возможны только два случая,
1, Граница области О пересекает отрезок с концами хт, хт-1 (рис, 1,о),
2, Отрезок границы области О частично налагается на [хт,хт-1 ] (рис, 1,6),
В обоих случаях мысленно выстраиваем проходящую через узел (т — 1, к) € Гп^ локальную границу, перпендикулярную рассматриваемому отрезку. Для осуществления рас-
Рис. 1. Построение локально-одномерных границ.
чета на (п + 1)-м временном слое на границе надо знать ит—1к, уп —1к, Ип —1к. Предполагая, что задача локально-одномерная, на границе принимаем
<—1к =0, ип—1к = ипк- нтк + нт—1к. (9)
Так как в вектор Г скорость уп+11к входит только умноженная на ит — 1к = 0, ее значение определять не нужно. Аналогичным образом рассматриваются остальные точки шаблона. Это правило действует для всех внутренних точек области независимо от того, вычисляется Ш или Ш.
На шаге предиктор приходится вычислять по схеме значения на границе (чтобы затем их использовать на шаге корректор). Пусть узел (М, К) — граничный (т. е. требуемая для вычисления точка шаблона ((М -1, К) или (М, К -1)) принадлежит ^м)- Пусть это для определенности узел (М -1, К). Тогда, если узел (М, К -1) — внутренний, то вычисления ведутся непосредственно по разностной схеме, если же он окажется граничным, то так как в этом узле надо знать вектор —1, вычисляются все переменные: ипмк—1) = 0,
И(М,К—1) = И(М,К) - нм,к + к(м,к—1)> У(М,К—1) = У(М——1)- Последняя формула (условие
(М - 1, К - 1)
или граничной.
Было исследовано много способов задания дополнительных граничных условий для схемы Мак-Кормака. В формуле (9) для одномерного направления по оси Ох использована формула (И - к)х = 0, аппроксимированная с первым порядком точности. Можно также использовать формулы второго порядка точности. Однако расчеты на последовательности измельчающихся сеток показали, что независимо от порядка аппроксимации разностная схема сохраняет второй порядок точности (этот факт теоретически известен).
Алгоритм вычислений в приграничных точках сильно упрощается, если на шаге предиктор искомые функции вычислять только во внутренних точках (как на шаге корректор), а недостающие на границе значения Ш вычислять по формулам (9) для величин с тильдой. Тестовые расчеты показали, что получаемая при этом погрешность вполне приемлема.
3. Свойства схемы
В случае одномерной системы линейных уравнений схема Мак-Кормака после исключения вспомогательного слоя сводится к схеме Лакса — Вендроффа, Здесь схему Мак-Кормака можно рассматривать как один из способов представления в виде схемы типа предиктор-корректор трехточечной по пространству двухслойной схемы второго порядка аппроксимации, В двумерном случае и в случае нелинейных уравнений это разные схемы. Выпишем первое дифференциальное приближение схемы Мак-Кормака для нелинейного уравнения переноса (уравнение Хопфа) щ + иих = 0:
ди д [и2\ А 2
+ 7Г I — ) = Аж
dt дх
д3и ди д2и (ди kiu^—r + + к3 [ —
дх3 дх дх2 \ дх
(10)
Здесь
k2 = ^(2сг2 + аа- 1), h = j(2a + a),
а = и At/Ах (число Куранта), а = 1 или а = 0 в зависимости от того, какая разность (левая или правая) применяется на шаге предиктор. Первое слагаемое в правой части отвечает за дисперсию. Для пологих волн этот член в правой части оказывается наиболее влиятельным. Следовательно, на грубой сетке численная дисперсия может конкурировать с физической дисперсией, присутствующей в реальном течении, но не учитываемой уравнениями мелкой воды. Заметим, что численная дисперсия может быть уменьшена за счет близости числа Куранта к единице. При уменьшении шага по времени дисперсия увеличивается, Правильная оценка дисперсионных эффектов разностной схемы очень важна для интерпретации численных решений задач динамики волн цунами.
Второе слагаемое вносит вклад в диссипацию, которая для устойчивости течения должна быть положительной. Однако, каким бы ни оказался знак выражения k2, знак диссипа-
тивного члена будет зависеть от знака Поэтому иногда приходится в решение вводить
дх
"искусственную" диссипацию, мажорирующую "схемную". Отметим, что другие известные схемы второго порядка аппроксимации (схема Лакса — Вендроффа, трехслойная схема leap-frog) также не лишены этого недостатка.
Второе дифференциальное приближение схемы Мак-Кормака для нелинейного уравнения выглядит довольно громоздко. Отметим лишь наличие в этом уравнении диссипа-тивного члена
Ах3 2 д4и
-—шг( (11)
который подавляет нелинейную неустойчивость (отсутствие аналогичного члена в недис-сипативной (см. [7]) схеме leap-frog в нелинейном случае приводит к разболтке решения).
Заметим также, что согласно методу дифференциальных приближений из выражения (10) вытекает условие устойчивости Куранта — Фридрихса — Леви разностной схемы Мак-Кормака [7]
< 1. (12)
3
В случае двух переменных для уравнения щ + иих + ууу = 0 условие устойчивости схемы, полученное методом дифференциальных приближений, имеет вид
+ 02 |< 1, (13)
где 01 = иАЬ/ Ах а2 = уАЬ/Ау. Формально это верно при о1о2 > 0, Однако в случае противоположного неравенства на практике неустойчивость не была обнаружена (подробнее
в И).
Если вместо скалярных уравнений рассматривать векторные, то в условиях устойчивости (12), (13) величины и и V надо заменить на максимумы собственных чисел соответствующих матриц.
4. Численные эксперименты
Как основа вычислительного алгоритма в работах [1-3, 10, 11] схема Мак-Кормака использовалась для моделирования сложного физического процесса трансформации и наката на берег длинных поверхностных волн в прибрежных акваториях, образовавшихся под воздействием крупномасштабного гравитационно-оползневого процесса. На примерах этих задач выполнены численные эксперименты для сравнения алгоритмов двух типов, один из которых реализует схему Мак-Кормака, а второй построен на основе схемы четвертого порядка точности во внутренних точках.
Численный алгоритм, основанный на схеме Мак-Кормака, применительно к решению задач о выходе волн на берег подробно описан в [1], Для сравнения взята неявная схема Адамса — Мултона — Бэшфорта, Рассмотрим одномерные уравнения мелкой воды:
+ Ях = Ф; (14)
Ъ = Я = (я„2 Нн 2/2) • ф= (Д )■ <16>
Схема Адамса — Мултона — Бэшфорта является двухшаговой. Первый шаг (предиктор) — разностная схема Адамса — Бэшфорта:
Ъп+1 - Ъп 1
тМп т + ~ + 5(Эт ) = 0, я = рх-Ф, д: = « (16)
Второй шаг — схема Адамса — Мултона, использующая четыре слоя по времени: Ъп+1 - 1
тМп т + + 19<& " б^Г1 + Япт~2) = 0. (17)
Производная по пространству аппроксимирована с четвертым порядком аппроксимации:
= Рт~2 ~ 8Рт~^А8хРт+1 " ^ + 0(Ах)\
Для решения неявных уравнений применяется метод простой итерации, причем для нулевого приближения Ът+1'0 используется найденное то схеме Адамса — Бэшфорта Далее схему Адамса — Мултона — Бэшфорта будем называть для краткости схемой Адамса, Сравнение численных решений обеих схем с гладким аналитическим решением уравнений мелкой воды, когда дополнительные граничные условия берутся из аналитического
решения, показывает бесспорное преимущество схемы Адамса, требующей всего одну-две итерации. Однако при усложнении задачи картина меняется. Так, выполнена серия расчетов течений в канале с ровным дном с целью определить глобальный порядок точности в расчетной области с учетом граничных условий. Рассмотрены границы двух типов: "отражающая" и "прозрачная". Как и следовало ожидать, в схеме Адамса, требующей из-за расширенного шаблона дополнительных граничных условий, глобальный порядок сходимости может существенно снизиться, если граничные условия невысокого порядка аппроксимации, тогда как схема Мак-Кормака имеет порядок точности, близкий ко второму, при граничных условиях как первого, так и второго порядка. Отметим, что при расчете течений в акваториях сложной формы с учетом островов крайне трудно применять граничные условия высокого порядка аппроксимации, особенно на "открытых" границах.
Для расчета подвижной границы (линии уреза) рассмотрены два алгоритма — "сквозного счета" для схемы Мак-Кормака и "экстраполяции" для схемы Адамса (другого алгоритма пока не придумано). Рассмотрена задача о накате на берег волны, вызванной погружением в воду "оползневой массы", имеющая аналитическое решение [10]. Схема Адамса хорошо приближает это решение на удалении от берега, однако вблизи подвижной границы погрешность, связанная с низкой точностью экстраполяции, хорошо заметна. Отметим, что при увеличении угла наклона берега погрешность уменьшается. Детальный анализ поведения линии уреза показывает, что на грубой сетке метод сквозного счета по схеме Мак-Кормака по точности уступает методу, основанному на схеме четвертого порядка аппроксимации, но с измельчением сетки численные траектории линии уреза, полученные разными методами, сближаются, что говорит о том, что метод экстраполяции понижает глобальный порядок точности схемы Адамса,
Обе схемы являются схемами четного порядка аппроксимации, что в случае значительной нелинейности или негладкой правой части уравнений h = h(x,t) приводит к появлению выраженной немонотонности. Поэтому в случае необходимости использовались численные фильтры, В частности, если при моделировании оползня донная поверхность задается функцией, имеющей разрывы производных (например, оползень — это "половинка" эллипса), то в схеме Адамса развиваются неустойчивые колебания, особенно в окрестности точки уреза, К таким изломам дна оказались чувствительны обе используемые схемы, но схема Адамса — в гораздо большей степени, так что ее решение становилось неустойчивым и разваливалось даже при использовании процедур сглаживания. Анализ результатов проведенных расчетов показывает, что по совокупности свойств метод, основанный на консервативной схеме Мак-Кормака, в котором линия уреза трактуется как внутренняя граница области вода — суша, более предпочтителен, так как он наименее чувствителен к возможному нарушению гладкости течения и неоднородностям дна.
Алгоритм, построенный на основе схемы Мак-Кормака, использован для моделирования ряда известных цунами, В работе [11] описаны результаты моделирования реального цунами, случившегоя 12 июля 1993 года у берегов Японии, Данное цунами выбрано для математического моделирования благодаря наличию подробной батиметрии морского дна и топографии прилегающей суши, а также достаточно полных данных о возникновении и проявлениях этого события. Индонезийское цунами произошло 26 декабря 2004 года в акватории Индийского океана в результате сильнейшего землетрясения. Возникшая волна цунами повлекла катастрофические последствия для стран этого региона, сопровождавшиеся человеческими жертвами. На рис, 2 показаны фрагменты динамики волн, рассчитанные с использованием алгоритма на основе схемы Мак-Кормака на протяжении 180 минут от начала зарождения. Следует отметить, в этих расчетах использована простей-
Рис. 2. Волновые картины при распространении Индонезийского цунами.
шая модель очага, заданного как часть эллипсоида, форма и координаты которого определялись в соответствии с имеющейся в распоряжении информацией. Расчеты выполнены C.B. Елецким в рамках программной системы Nereus, разработанной с использованием средств языка программирования Фортран (стандарт Fortran 95) и библиотек Win API (авторы: C.B. Елецкий, Л.Б. Чубаров, З.И. Федотова [5]).
Автор выражает благодарность д-ру физ.-мат. наук Л.Б. Чубарову, аспиранту НГУ C.B. Елецкому, студентке С.А. Бейзель за плодотворное сотрудничество при постановке и решении задач динамики поверхностных волн.
Список литературы
[1] федотова З.И. Обоснование численного метода для моделирования наката волн на берег // Вычисл. технологии. 2002. Т. 7, № 5. С. 58-76.
[2] Chubarov L.B., Fedotova Z.I. An effective high accuracy method for tsunami runup numerical modeling // Proc. of the NATO Advanced Research Workshop, Istanbul, Turkey, May 23-26, 2001. N.Y.: Kluwer Acad. Publ., 2003. P. 203-216.
[3] Елецкий С.В., Майоров Ю.Б., Максимов В.В. и др. Моделирование генерации поверхностных волн перемещением фрагмента дна по береговому склону // Вычисл. технологии. 2004. Т. 9, спецвыпуск. Ч. II. С. 194 206.
[4] MacCormack R.W. The effect of viscosity in hypervelocity impact cratering //J. Spacecraft and Rockets. 2003. Vol. 40, N 5. P. 757-763.
[5] Елецкий C.B., Федотова З.И., Чубаров Л.Б. Компьютерная модель волн цунами // Тр. X Байкальской Всерос. конф. "Информационные и математические технологии в науке, технике и образовании". Ч. I. Иркутск: ИСЭМ СО РАН, 2005. С. 138-146.
[6] Лятхер И.М., Милитеев А.Н., Школьников С.Я. Расчет наката волн цунами на берега // Изучение цунами в открытом океане: Сб. науч. тр. М.: Наука, 1978. С. 48-55.
[7] Шокин Ю.И., Яненко H.H. Метод дифференциального приближения. Применение к газовой динамике. М.: Наука, 1985. 365 с.
[8] Рождественский Б.Л., Яненко H.H. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1968. 687 с.
[9] Федотова З.И. о конструктивном подходе к исследованию устойчивости разностных схем // Вычисл. технологии. 2003. Т. 8, спецвыпуск: Тр. сов. рос.-казах, раб. группы по вы-числ. и информ. технологиям. С. 123-133.
[10] Федотова З.И., Чубаров Л.В., Бейзель С.А. Моделирование наката длинных волн в условиях динамически изменяющегося дна // Вычисл. технологии. 2004. Т. 9, спецвыпуск. Ч. III. С. 141-149.
[11] Федотова З.И., Чубаров Л.Б. Численное моделирование наката цунами // Вычисл. технологии. 2001. Т. 6, спецвыпуск. Ч. II. Тр. RDAMM-2001. С. 380-396.
Поступила в редакцию 4 апреля 2006 г.
