CC BY
33
УДК 519+531
DOI: 10.14529/mmph160206
КОНЕЧНО-РАЗНОСТНАЯ СХЕМА ПОЛУРАСЩЕПЛЕНИЯ МАТРИЧНЫХ КОЭФФИЦИЕНТОВ И ЕЕ ПРИМЕНЕНИЕ ДЛЯ ОПИСАНИЯ РАСПРОСТРАНЕНИЯ ЦУНАМИ
В.А. Симоненко, Н.А. Скоркин, А.С. Углов
ФГУП РФЯЦ-ВНИИТФ, г. Снежинск, Российская Федерация E-mail: a.s.uglov@vniitf.ru
Представлены конечно-разностная схема полурасщепления матричных коэффициентов одномерных уравнений газовой динамики и результаты расчетов по одномерному вычислительному коду, составленному на основе этой схемы, для решения задачи о распространении мегацунами в океане и выхода волны на береговую зону океана. Показана пригодность предлагаемого в данной работе вычислительного кода для решения задачи о накате морской волны на сушу. Дальнейшее развитие работы видится в разработке двухмерной программы расчета распространения поверхностной волны.
Ключевые слова: конечно-разностная схема полурасщепления; теория мелкой воды; цунами; морское побережье; космическое тело.
Введение
Процесс столкновения астрономического объекта с поверхностью океана и последствий этого явления можно разделить на три этапа: 1 - расчет начальной фазы удара и формирование поверхностной волны в океане, 2 - решение задачи о движении поверхностной волны в океане и выходе ее на шельф и 3 - распространение воды по суше. Различие пространственных и временных интервалов этого сложного явления вызывает значительные трудности при проведении сквозного расчета всех этапов такого процесса с помощью единой физической модели и одного вычислительного кода. Процедура решения 1-го и 2-го этапов задачи о распространении мегацунами приведена в работе [6]. В ней отсутствует описание конечно-разностной схемы, используемой при решении задачи о распространении мегацунами, но отмечено ее преимущество по сравнению с другими численными схемами: обладание минимальным сглаживающим эффектом. Для решения 2-го и 3-го этапов задачи предложенная в данной работе конечно-разностная схема в своей основе содержит метод неполного расщепления матричных коэффициентов - модификация метода расщепления матричных коэффициентов. Авторам не известны работы, в которых бы решалась поставленная выше задача о распространении мегацунами. Так, в работе [7] приводится основательный обзор численных методов решения задачи лишь для непродолжительного по времени процесса распространения волны по мелководью и ее наката на берег. Следует отметить, что в последнее время для решения задачи о накате волны широко используются уравнения На-вье-Стокса и SPH-метод их решения [8], требующий, однако, использования мощных суперкомпьютеров.
Постановка задачи
Рассмотрим две формы представления одномерного движения идеального газа в эйлеровых координатах. Пусть первой формой является представление в виде
Г 1
ut + иих =--px,
р
' Pt + ирх +рих = 0, (1)
Et + иЕх =- PUx, p = (y-1)рE,
I р
здесь и - компонент вектора скорости, р - объемная плотность, p - давление, Е - внутренняя энергия в единице массы газа, показатель адиабаты у= 1,4.
Представим систему (1) в другом виде. Для этого введем соотношения, известные из курса термодинамики, для скорости звука с и внутренней энергии Е, а именно:
ut + uux +--ccx = 0,
Симоненко В.А., Скоркин Н.А., Конечно-разностная схема полурасщепления
Углов А.С. матричных коэффициентов и ее применение ...
с' =YP' Е =-тЦ с2.
Р 7(7-1)
Выражая из данных формул значения давления и энергии через скорость звука и подставляя в систему уравнений (1), получим уравнения в виде [9, с. 325-326, (33.2)-(33.4)].
' 2_ Y-1
Pt + UPx +PUx = 0 (2)
Y- 1
ct + ucx + 2 cux = 0.
Доказательство. Третье уравнение для внутренней энергии системы (1) с помощью второго уравнения системы (1) представимо в виде полных дифференциалов в следующем виде
dE + pdV = 0, (а1)
в котором V = 1/ р - удельный объем. Рассмотрим 1-е уравнение системы (1) одномерного движения из нашей статьи
dp ,
ut + uux =--, (а2)
рдх
или в эквивалентной форме
p д ln p c2 д ln p .
ut + uux =---^— =---^—, (а3)
р дх у дх
так как c2 =ypjр. Обратимся к уравнению (а1), используя в нем соотношение V = 1/р . В результате будем иметь
dE =-p- dр. (а4)
р
1 p
Из уравнения состояния p = (у- 1)рЕ следует Е =--. Отсюда дифференцированием по-
Y-1 р
следнего соотношения получаем
dE = р (а5)
Y-1 р
Подставляем соотношение (а5) в выражение (а4):
1 р - рар=р_ dр (а6)
7-1 р р
После несложных преобразований сокращения на общий множитель 1/ р2 и приведения подобных слагаемых, сокращения на общий множитель —1— получим — = у<~р. Отсюда следует ре-
7-1 Р Р
шение в виде
1п р = у(1п р + 1п С0), (а7)
в котором Со - произвольная постоянная величина. Из соотношения (а7) легко получить выражение для давления
Р = Ср, (а8)
где С1 = С1 - константа.
Из выражения с2 = ур/р = уС1р^-1 определяем плотность
, 1 ^ —
Y-1
cY-1 =р. (а9)
( i >
Логарифмируем выражение (а9), введя для удобства новую постоянную дем иметь
Y-1
In C2 +--In c = In p.
Y-1
Из равенства (а7) с учетом (а 10) следует
In p = Y(ln C2 + --- In c + In Co),
(Y-1)
= ^ бу-
(а10)
(а11)
отсюда немедленно следует соотношение для частных производных (т. к. 1п , 1п С0 - постоянные величины)
дln p = 2 Y<d Inc = 2y cx
дх у-1 дх у-1 с Таким образом, уравнение (а3) с учетом соотношения (а 12) примет надлежащий вид
2ссх
(а12)
ut + uux +
(Y-1)
= 0.
Конец доказательства.
Отметим важное обстоятельство. В системе уравнений (2) уравнение состояния р = (у- 1)рЕ включено в уравнения системы, в отличие от системы (1), в которой уравнение
состояния рассматривается отдельно. Систему (2) можно переписать в матричном виде, введя обозначения
u u 0 2 Y-1
U = P , A = P u 0
c Y-1 c 2 0 u
Тогда систему уравнений (2) можно записать в виде
и, + лих = 0. (3)
Для численного решения уравнений вида (3) существует много различных методов. В работе [1] приведен перечень этих методов, в том числе метод полного расщепления матричных коэффициентов. Схема полного расщепления матричных коэффициентов предполагает представление системы уравнений движения идеального газа в эйлеровых координатах в виде (2). Однако это не всегда возможно, например, в случае табличного уравнения состояния или двучленного уравнения состояния Ми-Грюнайзена. Метод неполного расщепления матричных коэффициентов применим как для уравнений в форме (1), так и к уравнениям в форме (2).
Важное замечание. Дальнейшее изложение не использует систему уравнений в форме (2).
Результаты
Опишем метод неполного расщепления матричных коэффициентов. Прежде всего, отметим, что производная их аппроксимируется следующим образом [1]
и+ ~ и—и~1, если и^ + с^ > 0 и и+ ~ и+1—— в остальных случаях; (4) Х1 - Х1 -1 Х1+1 — Х1
их ~ ——и=1, если и1 - сг- > 0 и их ~ и——— в остальных случаях. (5) Х1 - Х1 -1 Х1+1 - Х1 Далее запишем систему (1) в матричном виде и{ + ЛиХ = В, где
и и 0 0"
и = р , А = 0 и 0
Е 0 0 и
В =
1
--Рх
р
-рих р..
р
Перепишем матрицу А в другом виде. Для этого представим элементы и матрицы А в виде
и=2 (и+с)+2 (и -с)
тогда
А =
К ч 1/ ч -(и + с) + -(и - с)
0
0
о
К ч Ь ч 2 (и+с)+2 (и - с)
0
0 0
1( ч 1. ч - (и+с)+-(и - с)
Конечно-разностные уравнения в этом случае записываются в виде
и, + А+иХ+ + А-и-= В,
(6)
где А = ТЛ+Т-1 + ТЛ-Т-1 = А+ + А- .
Здесь Т-1 - матрица, строки которой левые собственные векторы матрицы А, матрица Т -обратная матрица к матрице Т-1. Л+ и Л- - диагональные матрицы положительных и отрицательных собственных чисел матрицы А соответственно. Производные и+ и и- аппроксимируются конечными разностями аналогично конечным разностям (4) и (5).
Формулы частичного расщепления (6) для расчета сильных разрывов были модифицированы путем замены вектора В на вектор
В =
1 - р
--Рх, -рих ,--их
р р
в котором давление р вычисляется по формуле
Г р - на волне сжатия, р = 1
[ р - на волне разрежения, где р - давление, вычисляемое по форме, предложенной В.Ф. Куропатенко в работе [2].
Для уравнения состояния р = (у- 1)рЕ, рассматриваемого в данной работе, в узле с номером / + 0,5, согласно работе [2], выражение для давления принимает вид
р$5 = рП+0,5 + Ь(Ли)2 ^(д(Ли)2)2 +(аП+0,5 )2,
, 7+1 п / <2 Л 2 р
где Ь =-р +0 5, а = (рс) , Ли = и,+1 - и,, с = у—, п - номер временного слоя.
4 р
Т
Распространение поверхностной волны в мелкой воде
Для описания расходящейся кругами по воде волны и ее наката на побережье с учетом профиля дна и шельфа мирового океана был разработан вычислительный код ТМВ, в котором реализована теория мелкой воды.
Уравнения теории мелкой воды, описывающие распространение одномерной круговой поверхностной волны, имеют вид
du du дп
— + u— = g —,
dt dR SdR (7)
дВ дВ J du u Л n
— + u — + В\ — + — I = 0,
dt dR 4 dR R J
где R - расстояние вдоль невозмущенной поверхности океана от места входа тела в воду, t -время, g - ускорение свободного падения, u = u (R, t) - скорость частиц на поверхности воды,
?7 = n(R,t) - конфигурация волны, h(R) - глубина мирового океана, В = h(R) + n(R,t).
Уравнения (7) суть уравнения гиперболического типа. Особенностью рассматриваемой задачи является то, что распространение поверхностной волны осуществляется на очень большие расстояния 1000 км. Время движения волны до побережья океана - это десятки тысяч секунд [6]. Численное решение уравнений потребует осуществления большого числа циклов по времени. Поэтому конечно-разностные схемы, аппроксимирующие уравнения (7), должны иметь минимальный эффект численной диссипации [6].
Для решения данной задачи были протестированы схемы МакКормака, Лакса-Вендроффа и Неймана. Все они оказались неудовлетворительными по причине наличия сильного сглаживания схемной вязкостью [6]. Наилучшей оказалась Я-схема Моретти [1], она формирует монотонные профили, и в ней практически нет схемной вязкости. Применение Я-схема Моретти к уравнениям (7) давало физически неправдоподобные решения [6]. По этой причине была осуществлена модификация Я-схема сохранившая положительные свойства схемы Моретти. Эта модификация представляет собой описанную выше схему неполного расщепления матричных коэффициентов. Для уравнения (7) была написана одномерная вычислительная программа ТМВ. Для ее верификации было проведено сравнение результатов расчетов одной и той же задачи о распространении волны по программе ТМВ и аналогичной программе SWAN, представленной в работе [3]. В работе [3] представлено численное решение задачи о движении одиночной волны в океане глубиной 4,55 км. На удалении в 460,0 км от берега, то есть левым краевым условием счетной области являлось условие для скорости частиц на поверхности воды в виде
[0,04666 • sin2 (0,004713 • t)м/с, если t < 660 с; u = \ (8)
[0 м/с, если t > 660 с.
Данное условие описывает формирование уединенной волны высотой 1 м и шириной 140 км, которая движется со скоростью 210 м/с. В качестве правого граничного условия счетной области задавалось равенство нулю скорости частиц воды, то есть правый берег - отвесная скала.
Рассмотрен был и случай, когда левое граничное условие имело вид
0,04666 • sin2 (0,004713 • ^м/с, (9)
что означало формирование последовательности волн, распространяющихся со скоростью 210 м/с с амплитудой от -1 м до 1 м и шириной 280 км. В работе [3] расчеты были проведены при помощи двух программ: SWAN и ZUNI (SOLA). Геометрия счетной области имела вид, представленный на рис. 1 (левая часть). Глубина океана - Y = 4,55 км. На расстоянии X = 283,50 км от источника волны начинается скальный уступ с уклоном 1:15. На расстоянии 344,25 км начинается шельф постоянной глубины 500 м [3]. На рис. 1 (правая часть), который взят из статьи [3], приведена рассчитанная с помощью вычислительных кодов SWAN и ZUNI поверхность воды, соответствующая двум моментам времени. При этом в качестве источника волн на левой границе принято условие (9).Та же задача была рассчитана и с помощью кода ТМВ. Результаты расчетов вполне согласуются с данными из работы [6]. В этом можно убедиться, сравнив рис. 1 (правая часть), на котором представлены результаты работы программ SWAN и ZUNI, и 2, на котором представлен расчет по программе ТМВ.
Обратимся к результатам расчетов, приведенным в статье [3]. Согласно им волна (граничное условие (8)), взаимодействуя с подводным уступом с уклоном 1:15, увеличила свою высоту с 0,96 м до 1,50 м. Дальнейшее ее движение над шельфом до момента встречи со скальной стеной отмечено постоянством высоты волны (см. на рис. 1 пунктирную огибающую линию в верхнем правом углу).
Та же задача (граничное условие (8)) была рассчитана с помощью вычислительного кода ТМВ для теории мелкой воды [6]. Так же получено увеличение высоты волны с 1 м до 1,5 м и дальнейшее движение волны без изменения ее высоты [6]. На рис. 3 показано увеличение высоты волны в случае шельфа переменной глубины: глубина шельфа при X = 459,0 км равнялась нулю
[6]. г; кт
0. №35 О.ООЗО 0.0025 0.0020 0,0015 0.0010 0.0045 0.0000 -0.0005 -0.0010 -0.0015 -0,0020
100 200 300 400 ft,,,
Рис. 2. Волны на воде в моменты времени 3 000 и 3 840 с
1 1
г 3000s 3840 S
---time 1
1
Л r\ 1
1 - T" t \i
/ \ I N i ' Г
f / ' r V
/ \ t n / * T (
У -__ 4, / /
Ч У \
1 ton
0.003
time = 1440B
__ ____ _ time time = 18005 = 5200s
- - --- -- time = 2400s f.
time = 3000s 1 \
1
/ \ if 1
/ У I
4 / 1
/ V t 1 1
' \ - /i — f—
■ J
Я, fo«
Рис. 3. Профили волн на различные моменты времени для шельфа с уклоном
Итак, показано, что результаты расчетов аналогичных задач по вычислительной программе ТМВ и уже проверенной программе SWAN хорошо согласуются друг с другом.
Накат волны на сушу. Задача считается по уравнениям (11) до прихода волны к береговой линии Rq . Как в работе [4], для области R > Rq , т. е. на суше, вводится «глубина» h (R)< 0,
представляющая собой возвышение суши над невозмущенной поверхностью океана, взятое с отрицательным знаком. Для области R > Rq высота невозмущенной поверхности «воды» (до прихода волны к береговой линии) равна возвышению суши над невозмущенной поверхностью океана, взятому с положительным знаком n(R ) = -h (R )>0. На линии, разделяющей
воду от суши, принято n(R) = h (R) = 0. Поэтому в промежутке 0<R<Rq функция -з.оз; H (R, t = 0) ф 0, а для R > R0 функция H (R, t = 0) = 0. В результате получается задача о распространении сильного разрыва по воде и
по суше, которая решается по однородной модифицированной конечно-разностной схеме неполного расщепления матричных коэффициентов. (Описание схемы приведено выше).
Метод расчета наката волны на берег протестирован сравнением результатов численных расчетов с аналитическими решениями (см. рисунок 4), которые взяты из работы [4] со ссылкой на работу [5]. На рисунке представлен профиль волны Z (R, t ) на разные моменты времени. На этом
рисунке, по данным работы [4], приведено сравнение численного расчета наката волны (цунами) на берег, который имеет наклон 32°, с аналитическим решением [5]. Пунктирная линия - начальное положение волны, штриховая линия - положение волны на последующий момент времени, сплошная линия - момент максимального наката волны на берег, который обозначен линией с маркерами в виде кружков. Аналитическое решение представлено линией с маркерами в виде треугольников. Из анализа рисунка видно, что сравниваемые величины хорошо согласованы.
Выводы
Были рассмотрены метод и конечно-разностная схема полурасщепления матричных коэффициентов уравнений газовой динамики. Решена прикладная задача о движении поверхностной волны в океане и было осуществлено сравнение полученных результатов с результатами решения аналогичной задачи из работы К. Мейдера [3]. Проведена верификация вычислительной программы, в основе которой лежит рассматриваемая схема, путем сравнения результатов расчетов по предлагаемой схеме полурасщепления матричных коэффициентов с аналитическим решением задачи о накате волны на сушу. Было отмечено, что предложенная процедура вычислений работает правильно. Схема полурасщепления матричных коэффициентов является вполне пригодной для численных расчетов, где имеют место уравнения теории мелкой воды. Указание на применение конечно-разностной схемы неполного расщепления матричных коэффициентов приведено в работе [б]. Показана пригодность предлагаемого в данной работе метода и для решения задачи о накате морской волны на сушу. Дальнейшее развитие данной работы видится в разработке двухмерной программы расчета распространения поверхностной волны.
Литература
1. Флетчер, К. Вычислительные методы в динамике жидкостей. Т. 2 / К. Флетчер. - М.: Мир, 1991. - 552 с.
2. Куропатенко, В.Ф. Методы расчета ударных волн / В.Ф. Куропатенко // Дальневосточный математический журнал. - 2001. - Т. 2, № 2. - С. 45-59.
3. Mader, Charles L. Numerical Simulation of Tsunami / L. Charles Mader // Journal of Physical Oceanography. - 1974. - V. 4. - P. 74-82.
4. Marchuk, A.G.A. Method for Numerical Modeling of Tsunami Run-up on the Coast of Arbitrary Profile / A.G.A. Marchuk, A.A. Anisimov // ITS 2001 Processing, Session 7, Number 7-27, P. 933-940.
5. Мазова, Р.Х. Линейная теория набегания волн цунами на берег / Р.Х. Мазова, Е.Н. Пелиновский // Изв. АН СССР. Физика атмосферы и океана. - 1982. - Т. 18, Ш 2. - С. 1бб-171. (Препринт. ИПФ АН СССР № 25, Горький, 1981, 15 е.)
6. Симоненко, В.А. Математическое моделирование мегацунами / В.А. Симоненко, Н.А. Скоркин, В.П. Елсуков, А.С. Углов // Вестник ЮУрГУ. Серия «Математика, физика, химия». - 2008. - № 22, вып. 11. - С. б2-б9.
7. Бейзель, С.А. О некоторых численных алгоритмах расчета наката волны цунами в рамках модели мелкой воды / С.А. Бейзель, Н. Ю. Шокина, Г.С. Хакимзянов и др. // Вычислительные технологии. - 2014. - Т. 19, № 1. - С. 40-б2.
8. Ketabdari, M.J. Simulation of waves group propagation and breaking in coastal zone using a Navier-Stokes solver with an improved VOF free surface treatment / M.J. Ketabdaria, M.R.H. Nobarib, M. Moradi Larmaei // Appl. Ocean Res. - 2008. - Vol. 30, № 2. - P. 130-143
9. Кочин, Н.Е. Теоретическая гидромеханика. Часть II / Н.Е. Кочин, И.А. Кибель, Н.В. Розе. - М.: Физико-математическая литература, 19бЗ. - 728 с.
Поступила в редакцию 16 декабря 2014 г.
DOI: 10.14529/mmph160206
A SEMI-SPLITTING FINITE DIFFERENCE SCHEME OF MATRIX COEFFICIENTS AND ITS APPLICATION TO DESCRIBE THE TSUNAMI PROPAGATION
V.A. Simonenko, N.A. Skorkin, A.S. Uglov
Russian Federal Nuclear Center - Zababakhin All-Russia Research Institute of Technical Physics, Snezhinsk, Russian Federation E-mail: a.s.uglov@vniitf.ru
The well-known finite-difference scheme of Moretti of splitting the matrix coefficients of the system of gas dynamics equations involves writing equations in a special form - pressure and internal energy are excluded from the equation using the equations of state for an ideal gas. In this paper the author proposes a modification of Moretti scheme as a finite-difference scheme of semi-splitting of matrix coefficients which do not intend to constitute a system of equations in a special form. The semi-splitting scheme allows solving equations of hyperbolic equations of state of any type, for example, even those in tabular form. For one-dimensional equations of propagation of circular wave on the water surface, that are the equations of hyperbolic type, the results of the calculations of the problem of the propagation of a surface wave in the ocean and the output of wave on the shore of the ocean area are given according to one-dimensional computational code of the shallow water theory. Verification of semi-splitting finite-difference scheme is performed by comparing the calculation results for the problem of the propagation of a single surface wave in the ocean and the problem of propagation of a wave train on the ocean surface with the results of calculations of the same problems cited in the work by C. Mader. To calculate the wave setup on the shore the approach is used in the computational code of the shallow water theory, which is described in the work by An. G. Marchuk, A.A. Anisimov. By comparing the calculation results with analytical solutions, the suitability of the computational code proposed in this work for the problem solution of the sea wave setup on land is demonstrated. Further development work is seen in the development of a two-dimensional program for calculating the surface wave propagation.
Keywords: semi-splitting finite-difference scheme; shallow water theory; tsunami; seashore; astronomical body.
References
1. Fletcher K. Vychislitel'nye metody v dinamike zhidkostey. Tom 2 [Computational Techniques for Fluid Dynamics. Vol. 2], Moscow, Mir Publ., 1991, 552 p. [Fletcher C.A.J. Computational Techniques for Fluid Dynamics, Vol. 2, Springer Verlag Ed., 1988.]
2. Kuropatenko V.F. Methods of shock wave calculations. Dal'nevostochnyy matematicheskiy zhur-nal, 2001, Vol. 2, no. 2, pp. 45-59.
3. Mader Charles L. Numerical Simulation of Tsunami. Journal of Physical Oceanography, 1974, Vol. 4, pp. 74-82.
4. Marchuk Andrey G.A., Anisimov A.A. Method for Numerical Modeling of Tsunami Run-up on the Coast of Arbitrary Profile. ITS 2001 Processing, Session 7, no. 7-27, pp. 933-940.
5. Mazova R.Kh., Pelinovskiy E.N. Lineynaya teoriya nabeganiya voln tsunami na bereg [Linear Theory of Wave Run Up on a Beach]. Izv. AN SSSR. Fizika atmosfery i okeana, 1982, Vol. 18, Sh 2, pp. 166-171. (in Russ.).
6. Simonenko V.A., Skorkin N.A., Elsukov V.P., Uglov A.S. Matematicheskoe modelirovanie megatsunami [Mathematical Modeling of Magatsunamis]. Vestnik YuUrGU, ser. "Matematika, fizika, khi-miya", 2008, no. 22, issue 11, pp. 62-69. (in Russ.).
7. Beizel S.A., Shokina N.Y., Khakimzyanov G.S., Chubarov L.B., Kovyrkina O.A., Ostapenko V.V. On some numerical algorithms for computation of tsunami runup in the framework of shallow water model.I. Computational technologies, 2014, Vol. 19, no. 1, pp. 40-62. (in Russ.).
8. Ketabdari M.J., Nobari M.R.H., Larmaei Moradi M. Simulation of waves group propagation and breaking in coastal zone using a Navier-Stokes solver with an improved VOF free surface treatment. Applied Ocean Research, 2008, Vol. 30, Issue 2, pp. 130-143. DOI: 10.1016/j.apor.2008.08.005
9. Kochin N.E., Kibel' I.A., Roze N.V. Teoreticheskaya gidromekhanika Chast' II. [Theoretical Hydromechanics Theoretical Hydromechanics. Part II]. Moscow, Fiz-matlit Publ., 1963, 728 p. (in Russ.).
Received December 16, 2014
