CC BY
22
УДК 519.6
DOI:10.22449/0233-7584-2019-6-600-620
Новый численный алгоритм для уравнений многослойной мелкой воды на основе гиперболической декомпозиции и схемы КАБАРЕ
В. М. Головизнин*, Павел А. Майоров, Петр А. Майоров,
А. В. Соловьев
Институт проблем безопасного развития атомной энергетики РАН, Москва, Россия
E-mail: gol@ibrae.ac.ru
Поступила в редакцию 02.08.2019 г., после доработки - 08.10.2019 г.
Цель. Описание новой методики численного решения уравнений гидродинамики несжимаемой жидкости со свободной границей и переменной плотностью в гидростатическом приближении - цель настоящей работы.
Методы и результаты. Алгоритм основан на методе гиперболической декомпозиции - представлении многослойной среды в виде отдельных слоев, взаимодействующих через границы раздела. Силы, действующие на верхнюю и нижнюю границы каждого слоя, трактуются как внешние, не нарушающие свойства гиперболичности системы уравнений для каждого слоя. Для решения системы гиперболических уравнений с переменной плотностью в каждом слое используется явная схема КАБАРЕ. Схема имеет второй порядок аппроксимации и обратима по времени. Ее особенностью является повышенное число степеней свободы - наряду с консервативными переменными, определенными в центрах расчетных ячеек, используются потоковые переменные, отнесенные к серединам граней. Система уравнений многослойной мелкой воды не является безусловно гиперболической и при потере гиперболичности становится некорректной. Гиперболическая декомпозиция не устраняет некорректности исходной системы. Для регуляризации численного решения предлагается использовать следующий набор средств: фильтрацию на каждом временном шаге потоковых переменных скорости, плотности и толщины слоя; сверхнеявную аппроксимацию градиента давления; линейную искусственную вязкость; переход к эйлерово-лагранжевым (СЭЛ) переменным, приводящий к обмену между слоями массой и импульсом. Основным средством, стабилизирующим численное решение на больших временах, является переход к СЭЛ-переменным. Остальные приемы вспомогательные и используются для тонкой настройки.
Выводы. Показано, что для обеспечения регуляризации и гарантированной устойчивости задач необходимо не только перестраивать расчетную сетку на каждом временном шаге, но также использовать фильтрацию потоковых переменных и искусственную вязкость, моделирующую турбулентное перемешивание.
Ключевые слова: нестационарная гидродинамика, свободная поверхность, гидростатическое приближение, некорректные задачи, переменная плотность, численный алгоритм, гиперболическая декомпозиция, схема КАБАРЕ.
Благодарности: работа выполнена при финансовой поддержке гранта РНФ № 18-11-00163. Авторы выражают благодарность В. Б. Залесному и Е. В. Семенову за плодотворные обсуждения и конструктивные замечания.
Для цитирования: Новый численный алгоритм для уравнений многослойной мелкой воды на основе гиперболической декомпозиции и схемы КАБАРЕ / В. М. Головизнин [-и др.] // Морской гидрофизический журнал. 2019. Т. 35, № 6. С. 600-620. doi:10.22449/0233-7584-2019-6-600-620
© Головизнин В. М., Майоров Павел А., Майоров Петр А., Соловьев А. В., 2019
New Numerical Algorithm for the Multi-Layer Shallow Water Equations Based on the Hyperbolic Decomposition and the CABARET
Scheme
V. M. Goloviznin*, Pavel A. Maiorov, Petr A. Maiorov, A. V. Solovjov
Nuclear Safety Institute, Russian Academy of Sciences, Moscow, Russia *e-mail: gol@ibrae.ac.ru
Purpose. The present article is devoted to describing a new method of numerical solution for hydrostatic approximation of incompressible hydrodynamic problems with free surfaces and variable density. Methods and Results. The algorithm is based on the hyperbolic decomposition method, i. e. representation of a multilayer model as a sum of the one-layer models interacting by means of the reaction forces through the layers' interfaces. The forces acting on the upper and lower interfaces of each layer are interpreted as the external ones which do not break hyperbolicity of the equations system for each layer. The explicit CABARET scheme is used to solve a system of hyperbolic equations with variable density in each layer. The scheme is of the second approximation order and the time reversibility. Its feature consists in the increased number of freedom degrees: along with the conservative-type variables referred to the centers of the calculated cells, applied are the flux-type variables related to the middle of the vertical edges of these cells. The system of the multilayer shallow water equations is not unconditionally hyperbolic, and in case hyperbolicity is lost, it becomes ill-posed. Hyperbolic decomposition does not remove incorrectness of the original system of the multilayer shallow water equations. To regularize the numerical solution, the following set of tools is propose: filtration of the flow variables at each time step; super-implicit approximation of the pressure gradient; linear artificial viscosity and transition to the Euler-Lagrangian (SEL) variables that leads to the mass and momentum exchange between the layers. Such transition to the SEL variables is the basic tool for stabilizing numerical solution at large times. The rest of the tricks are the auxiliary ones and used for fine tuning. Conclusion. It is shown that regularizing and guaranteeing the problems' stability requires not only reconstruction of the computational grid at each time step, but also application of the flow-type variables' filtering and the artificial viscosity simulating turbulent mixing.
Keywords: non-stationary hydrodynamics, free surface, hydrostatic approximation, ill-posed problems, variable density, numerical algorithm, hyperbolic decomposition, CABARET scheme.
Acknowledgments: The work was supported by the Russian Science Foundation, project No 18-1100163. The authors are grateful to Zalesny V.B. and Semenov E.V. for productive discussions and constructive comments.
For citation: Goloviznin, V.M., Maiorov, Pavel A., Maiorov, Petr A. and Solovjov A.V., 2019. New Numerical Algorithm for the Multi-Layer Shallow Water Equations Based on the Hyperbolic Decomposition and the CABARET Scheme. Physical Oceanography, [e-journal] 26(6), pp. 528-546. doi: 10.22449/1573-160X-2019-6-528-546
Введение
Достижения последних десятилетий в области моделирования и прогноза морских течений связаны в первую очередь с бурным развитием высокопроизводительных вычислительных систем и новыми математическими алгоритмами [1-6]. Современные задачи гидродинамики океана описываются нелинейными уравнениями с частными производными, решаются в сложных областях с высоким пространственным разрешением [3, 4]. Во многих случаях эти задачи не являются классическими и требуют разработки специальных методов их пространственной аппроксимации, а также методов решения по времени, включая явные, неявные и схемы расщепления [1, 2, 5].
Одним из главных преимуществ использования явных схем является их эффективное распараллеливание на многоядерных вычислительных платформах. Для решения прогностических задач выигрыш за счет эффективности распараллеливания часто превышает потери за счет уменьшения временного шага.
Для моделирования гидродинамики морей и океанов широко используется многослойное гидростатическое приближение. Несмотря на богатую историю развития и большое число численных реализаций, задача дальнейшего совершенствования вычислительных алгоритмов для уравнений многослойной гидродинамики со свободной поверхностью остается актуальной [7, 8].
Уравнения многослойной мелкой воды не являются безусловно гиперболическими [9, 10]. Это приводит к тому, что начально-краевая задача в процессе ее решения становится некорректной и алгоритмы решения безусловно гиперболических уравнений теряют устойчивость*. Считается, что потеря гиперболичности происходит при развитии неустойчивости Кельвина -Гельмгольца на границах раздела слоев [11]. В этих случаях должен возникать интенсивный обмен массой и импульсом между слоями, запрещенный в классическом многослойном приближении [8]. Этот недостаток можно устранить включением в многослойные уравнения так называемой турбулентной вязкости. Она часто зависит от параметров рассчитываемого течения и содержит эмпирические параметры, настраиваемые на различные типы течений [12, 13]. Недостатком такого подхода является критическая зависимость результатов расчета от опыта вычислителя.
Другой подход к решению указанной проблемы заключается в модификации классических уравнений многослойной мелкой воды - снятии запрета обмена массой и импульсом между слоями [14, 15]. Фактически это означает переход от лагранжевого описания вертикального движения границ слоев к смешанному эйлеро-лагранжевому описанию [16]. В этом случае между слоями возникают потоки массы и импульса, регуляризирующие задачу.
При реализации такого подхода мы используем метод расщепления по физическим процессам. Вначале вычисляются значения физических переменных на новом временном слое по классическим уравнениям многослойной мелкой воды с переменной в каждом слое плотностью. Затем задаются предписанные вертикальные координаты границ слоев и находятся возникающие между ними обмены потоками массы и импульса. При построении предписанных вертикальных границ слоев предполагается, что верхний слой со свободной границей остается лагранжевым, нижний - неподвижным. Если координаты промежуточных границ для каждой вертикали пропорциональны друг другу, то такие расчетные сетки можно назвать сигма-сетками [1, 17]. Другими, часто используемыми в океанологии, являются зет-сетки, когда границы всех слоев, за исключением верхнего, остаются неподвижными [1, 17].
Для численного решения классических уравнений многослойной мелкой воды с различной плотностью слоев разработано много алгоритмов, имеющих свои достоинства и недостатки [18, 19]. Как правило, вычислительные алгоритмы строятся методами конечных разностей или конечных объемов,
* Кабанихин С. И. Обратные и некорректные задачи. Новосибирск : Сибирское научное издательство, 2009. 457 с.
а для аппроксимации конвективных потоков применяется какой-либо из вариантов сноса значений вниз по потоку [20].
Следует отметить, что классические методы, основанные на решении однослойной задачи о распаде произвольного разрыва [21], в многослойном случае не применимы. Это связано с тем, что задача о распаде произвольного разрыва в многослойном случае не имеет простого решения.
В предлагаемой новой численной методике на первом этапе расщепления методом экстраполяции локальных римановых инвариантов вычисляются горизонтальные потоки, как это делается в схеме КАБАРЕ [6]. Для этого многослойная модель представляется в виде суммы однослойных (гиперболическая декомпозиция), взаимодействующих посредством сил реакции, приложенных к границам раздела. Для каждой из однослойных моделей записывается схема КАБАРЕ с нелинейной коррекцией потоков на основе принципа максимума [22]. При этом выполняются условие несжимаемости и законы сохранения массы и импульса. В результате (на новом временном слое) находятся потоковые и консервативные переменные без учета обмена массой и импульсом между слоями. Поскольку эта задача является некорректной, новые потоковые переменные подвергаются фильтрации, не нарушающей законов сохранения. Результирующий алгоритм имеет второй порядок аппроксимации и обладает свойством сбалансированности - сохранять при постоянной плотности состояние покоя над любым донным рельефом [23].
На втором этапе расщепления по физическим процессам задаются вертикальные границы слоев и вычисляются потоки между слоями.
Для расширения области устойчивости алгоритма применяются две процедуры, понижающие порядок аппроксимации по времени до первого. Это включение линейной вязкости, приводящей к наибыстрейшему затуханию самой высокочастотной гармоники [24], и сверхнеявная аппроксимация градиента давления [25].
Уравнения однослойной мелкой воды с переменной плотностью, учетом внешнего давления и донного рельефа
В случае одного пространственного измерения выпишем интегральные законы сохранения (баланса) для произвольной подобласти х е[х1, х2 ] (рис. 1). Введем обозначения:
V = | hdx, М = | phdx, П = | puhdx, к = Н - В,
где V - площадь, занятая жидкостью; M - масса; П - импульс; H(x, 0 - свободная поверхность; B(x) - рельеф дна; h(x, 0 - толщина слоя; 0 - плотность; u(x, 0 - горизонтальная скорость жидкости.
Уравнение изменения площади, занимаемой жидкостью, имеет вид
—V = — X ксХ = - ик\ + ик =-Х——с1х. дt дt I 2 * I —х
Р и с. 1. Однослойная модель Fig. 1. One-layer model
Изменение массы и импульса описывается соответственно уравнениями
1 —phu
^^ = — j phdx = - puhl + puh| =-f
—t —tJ x2 x
x J —x
xi
-dx,
—П — X2
-= — [puhdx = -(pu2h + pgh2 / 2 + PTh) + (pu2h + pgh2 / 2 + PTh)
—t — tJ v "ъ v "x
x2 — в I
j PB —dx + j P
—H
—x
dx.
Здесь g - ускорение свободного радения; Pт - давление на свободной поверхности; Pв = Pт + pgh - давление на дне.
Из приведенных интегральных уравнений следуют соответствующие законы сохранения в дифференциальной форме:
—h —hu
— +-= 0,
—t —x
—ph —phu
—t
—x
—puh —phu2 g —ph2 —PTh
- +
■ +
+ ■
= F = P
—H
- P
—B
д1 дх 2 дх дх дх дх
После некоторых преобразований эти уравнения можно представить в виде
^ + А х ^ = В, ф = (И,р,и)т, й = (ОДР/рй),
где
A =
Л
u 0
0 h u0 gh —p
p (ph Г gh—p
B yy ' 2p —x
u
F=P —H - P —B - hdp
—x
—x
—x
Все собственные числа матрицы А действительны:
^ = и + с, Х2 = и - с, Х3 = и, с = ^РВ /р,
откуда следует так называемая характеристическая форма исходных уравнений:
' ди с дк gh др
\
- + —
+ -
д к дt 2рс дt ди с дк gh др
+(и+с)
/ди с дк gh др
Л
■ + —
+ -
дt к дt 2рс дt
др др _ — + и — = 0. дt дх
+ (и - с)
дх к дх 2рс дх
/ди с дк gh др дх к дх 2рс дх
= ( Р/ рк ), = ( Р/ рк ),
Уравнения многослойной мелкой воды с переменной плотностью, учетом внешнего давления и донного рельефа
Уравнения многослойной мелкой воды с переменной плотностью при отсутствии обмена массой и импульсом между слоями можно представить как сумму однослойных уравнений, связанных между собою силами давления и реакции, действующими на границы раздела (рис. 2):
дК , ки )
■ +
= 0,
^рА , а( рки )к
+
дt дх дt дх
д(Рки)к , д(Рки2\ , g д(РкI , д(ккРк)
= 0,
дt
+
дх
+
2 дх
+
дх
■ + Р
дг,
дх
— Р
дх
= 0, (1)
к = 1,...,ы, 2Х = н, = Б, кк = гк -2к+!, Рк+1 = Рк + gPkhk.
Здесь к - номер слоя, отсчитываемый от свободной поверхности; Хк - координаты верхней границы слоя; Р1 = Рт - внешнее давление на свободную поверхность 2 . Так называемая простая форма системы уравнений (1) будет иметь вид:
- ^ , . чг
аГ + ~дх = = ,
где С - матрица размерности N х Л^и I) - некоторая правая часть. Хорошо известно, что уже при N = 2 матрица С может иметь комплексные корни и система не будет безусловно гиперболической [9], что порождает известные вычислительные трудности. Для большего числа слоев ситуация усугубляется. Уже при наличии двух слоев вычисление собственных значений матрицы (4 х 4) - достаточно сложная задача. При большом числе слоев она становится практически неразрешимой и непосредственное использование ба-лансно-характеристических методов [23] оказывается невозможным. Для преодоления этого затруднения мы используем прием, который будем называть гиперболической декомпозицией задачи.
Р и с. 2. Многослойная модель Fig. 2. Multilayer model
Если считать силы, действующие на границы раздела, внешними, то (1) можно представить в виде:
где
A =
V
uk 0 hk
0 uk 0
(Ph )_1 ghk 2Pk
F = P
Fk = Pk
8Z}
__p 8Zk±i _ h P
1 k±l - ' lk -8x 8x 8x
j
что приводит к системе независимых характеристических уравнений:
^ ^ \(8К ghk 8Pk^ k ± ck 1 ±'
(8uk+c18h^+.shk 8Pk
v
8t hk 8t 2Pkck 8t j
(8uk ck 8hk ghk 8P ^
8t hk 8t 2ptck 8t
P ± uk 8Ъ = 0.
8t 8x
±(uk ±ck) ±(uk _ck)
■±-
8x hk 8x 2ptck 8x
(8uk ck 8h ghk 8P ^
8x hk 8x 2ptck 8x
= ( Fkl Pkhk), = (FklPkhk ), (2)
Это дает возможность использовать для численного решения полной системы многослойных уравнений алгоритмы, хорошо зарекомендовавшие себя в однослойном случае. Некорректность полной системы при этом не
j
\
j
исчезнет и будет проявляться в высокочастотном искажении границ. С помощью специальных фильтров искажения можно существенно ослабить, сохраняя порядок аппроксимации и не нарушая законы сохранения.
Схема КАБАРЕ для послойного решения уравнений многослойной
мелкой воды
Область решения задачи покроем неравномерной расчетной сеткой с координатами узлов ж,, 1 = 1, №. В узлах в начальный момент времени зададим вертикальные координаты слоев цк;:
X) = %, (X0)= H, ^ )= B¡, I=1,...,N.
Неизменные координаты узлов по оси х и зависящие от времени вертикальные координаты слоев ХС{ определяют расчетную сетку на плоскости (х, ¿).
В схеме КАБАРЕ используются два типа переменных: потоковые и консервативные [22]. Потоковые переменные относятся к серединам вертикальных граней расчетных ячеек, консервативные - к их центрам. Потоковые переменные, относящиеся к слою с номером к, будем обозначать как
РП,I> < I > ¿С, I = ХС,I - ХГ+1,I , консервативные - как РП1+1/2 , К, 1+1/2 > К 1+1/2 . В начальный момент задаются потоковые переменные. Консервативные переменные определяются по потоковым следующим образом:
Р1,г +1/2 = 0,5 (Р01 + р1,г+1 ) , ик,, + 1/2 = 0,5 (и0,/ + и0,/ +1 ) ,
К, +1/2 = 0,5 ( + И1 +! ).
Вычислительный алгоритм состоит из следующих элементов: вычисление консервативных переменных на промежуточном временном слое по консервативным конечно-объемным схемам (фаза 1); вычисление потоковых переменных на новом временном слое экстраполяцией локальных инвариантов Римана (фаза 2); вычисление консервативных переменных на новом слое по конечно-объемным схемам (фаза 3).
Фаза 1. Для вычисления промежуточных консервативных переменных используем:
СГ - ^ ,(¿и )1к -(¿и )1=0
т/2 Лх
(РЪГГ -(РЪУс к , (РЫ )1к -(РЫ )
= 0,
т/2 Лг
(Р^Г+Г-(РНи)"сл (рИЫ %-(Р^и ^ (КРк+1121 -(КРк+112)"ь
>С,к V !к,к V )ь,к
т/2 Лх Лх
РК,к+1 + РЬ,к+1 ХК,к+1 - ХЬ,к+1 РК,к + РЬ,к ХК,к - ХЬ,к _ д
2 Лх 2 Лх '
где Р+1/2 =(Р + . Данные уравнения обладают первым порядком
аппроксимации по времени и вторым - по пространственной переменной.
Фаза 2. Для нахождения новых потоковых переменных используется линеаризованная система характеристических уравнений (2), которую можно записать в следующем виде:
и к ь Ск п+1/2 СКк + ёКк п+1/2
а V _ К _ ,+1/2 ^ _ 2ркСк _ 1+1/2
Фк а
Л
+[(«* +Ск)]
1+1/2
Сик , Ск п+1/2
Сх V _ К _ 1+1/2 СХ
2РЛ
п+1/2 ^
дРк сх
1+1/2
= [( Рк/ ркК )]П2,
Сик Ск п+1/2 СКк ёКк п+1/2
ат V _ К _ .+1/2 _ 2ркСк _ 1+1/2
срк а
\
+[(и - Ск)]
п+1/2 1+1/2
Сик Ск п+1/2 ёК п+1/2 Л Фк
Сх V _ К ] 1+1/2 СХ _ 2ркСк _ сх 1+1/2 у
= [( Р/ ркнк)]
п+1/2 1+1/2 :
+ Ги г+1/2 = о
сТ + К ]'+1/2 сХ = 0'
что можно записать как
4)
а '2 )к
к,1+1/2 /л \ п+1/2 Л )к,¡+1/2 \п+1/2
+(Л1 л+1/2 —— - (и,к) 1+1/2,
■ + (X 2 )
■ + ( " )
п +1/2 1+1/2
дх д(12 )к.
С(Ь )к„+1/2
д? 4 "1+1/2 дх
п +1/2 1+1/2
дх
С( ¡3 )к+1
■ = (в2к )
п +1/2 1+1/2
(3)
= (&,к )
п +1/2 1+1/2
где
/Л \П+1/2 / \П+1/2 /Л \П+1/2 / \П+1/2 /Л \ П+1/ 2 / \.
(" )п+1/2 =К + Ск )п+1/2 , ("2 )г+Ш = (ик - Ск )п+1/2 , (" )г+1/2 =(ик )
\п+1/2 / \ п+1/2
п+1/2
*п+1/2 / \ п+1/2
ъ2 )1+1/2 =(ик - Ск )г-+1/2 , ("3 );+1/2 1/2 1/2
1/2 к )1+1/2 '
1/2 1/2 1/2 1/2 (^1,к )1 +1/2 =(^2,к )г+1/2 4(^/рА )],1/2 , (й,к )г-+1/2 = 0,
(IX -+1/г, (¡2 -+1/2, (¡з X -+1/2 - локальные римановы инварианты в ячейке i + 1/2
слоя k на интервале времени [?я, ?я+ 1 ]:
(71 )к,+1/2 = ик +СМ/2Л +Вм2,крк , (¡2 )к
/к ,1+1/2
_ _гп+1/2 7 _ р.п+1/2 /т \ _
= ик °г+1/2,кПк 1/2,крк , (¡3 )к ,+1/2 = рк
о
1+1/2,к
д
1+1/2,к
ёК 2ркСк
Отметим, что при переменной плотности собственные числа матрицы А^ остаются такими же, как и при постоянной плотности [24], тогда как собственные левые векторы модифицируются, что приводит к появлению в локальных инвариантах дополнительного члена, пропорционального плотности.
Первое действие фазы 2 состоит в вычислении предварительных значений локальных инвариантов на новом временном слое методом линейной экстраполяции с учетом знака характеристической скорости (Хт ) , т = 1, 2, 3 :
ka:,;,-Cm»;., f (а:;;, > о and (к с, * «■
kj+i
\я+1/2 _( Т )n
~m'k,¿+1/2 V m)k,i
2 _(Тт)П,г+2 if < 0 and * 0,
/ \ n+1/2 /r \ n+1/2 .r /. \ n+1/2 \ n+1/2 _
(Im)n,i+3/2 +(In)n,+1/2 /2 f (K)n++3/2 x(Mn++1/2 < 0
s.n+1/2 ) k,¿+1/2 \ n+1/2 ) k,¿+3/2
n+1/2 / k ,¿+3/2
n+1/2 "m) k ,¿+3/2 n+1/2 ) k,¿+1/2
\ n+1/2
-J k ,¿+1/2
Следующее действие состоит в коррекции полученных значений в соответствии с принципом максимума для уравнений (3). Вычислим для каждой ячейки минимальные и максимальные значения локальных инвариантов на текущем временном слое:
min
max
min
max
in (Im,k),+1/2 = min {(I:,k)г+l, (С ,k)г+l/2, (Ilk\^ (I:,k)+1/2 = max {(I:,k)г+l, (Г: kL^ (In,k)
in (I:,k )¿+1 = min {min (I:,k),+1/2 , min (I:,k )¿+3/2 }, (I:,k )¿+1 = maX {maX (I:,k )¿+1/2 , maX (I~k)i+m }
и правые части уравнений переноса (3) по известным левым частям:
(In+1/2) _(r) (In) _(In)
/ \n+1/2 V : 'k,¿+1/2 V :>k,¿+1/2 , /Л yn+1/2 V :)k,¿+1 V 4k,i , r>( A
(Q:,k)¿+1/2 =-T/2-+ (^1+1/2 -^-+ 0(T ' '
At
(e3,k)n::=o.
Из уравнений переноса по характеристикам (3) следует, что если величины
/ \n+1/2
(Q: k) равны нулю и число Куранта - Фридрихса - Леви (CFL) меньше единицы, то должны выполняться условия:
min (I , ) , max (I m,\ if (A, )"+1/2 > 0 and (A )"+12 > 0,
V 'n'k Ji+1/2 V )i+1/2 ^ : 'fc,>+1/2 V Wt,¿+3/2 '
™n (I:0,+3/2 ] f (АЛТъ12 < 0 and (^ 0,
™n (L, )j+! , max (,+! ] f (AЛ]X ( A: )k ■ +1/2 < 0
/ \ n+1/2
При отличных от нуля (Qm А ) интервалы допустимых значения смещаются:
(I Y
\ rn/k
е <
max (lmk )I+1/2 = max (lm,k )I+1/2 + T (Qmk )"Z ,
П1т (1 т.к ),i2 = П1т (1 т.к ),|2
\ n+\
Im)k заключается в следующем: если предварительное значение инварианта {Jm)k ] удовлетворяет этим неравенствам, то
инвариант не меняется; если нет, то ему присваивается значение, соответствующее ближайшей границе допустимого интервала. По откорректированным значениям локальных инвариантов на новом временном слое из системы уравнений
(pk ГМ 4 ,
\Я+1 s-,* /т \n+i __.* / \n+i / , \n+i
(uk )"++i + g* ( к )n+; + D* (Pk )n+; =( ii ^,
(uk )n++ii - G2,k( hk )n++ii - D*k (pk )n++ii =( 12 )n++i вычисляются новые потоковые переменные:
G* и ( I* ) "+i + G*t ( I* ) "+i
f \"+i -i г \n+i i v+i _ ik У 2 /k,i+i 2'M i /k,i+i (Pk) i+i -( 13)k,,+i , (uk) i+i - r* j-r*
Gi,k + G2,k
« \ n+i / * \ n+i
IJ
( к )П+ -
(*\n+i / *\n
i) -( I* )
Wk,i+i V 4k,i+i
k i+i /-1* /-1* Gi,k + G2,k
где
(Ji* )n++ -(iiC+- D*k ып, ( i2*) n+i-(i 2C++Dik (px;,
Gj G2 k принимают значения:
te ,D-£ ) if (^ )n+;+/i2/2 > 0 and (^ ^ * 0, K,k, Dl,k Hte ,DS ) if (^m C2 < 0 and (^ C2 < 0, (G^, D-f ) if (К^ x(^C2 < 0,
где m -1,2, G^ - (g££ + G^k ) / 2, D^2 - (D^ + D^k ) /2.
Фильтрация потоковых переменных. В областях потери гиперболичности происходит генерация коротковолновых возмущений. Для регуляризации решения вычисленные потоковые переменные подвергаются фильтрации:
\и+1 / \и+1 /1 \ /_ \и+1 /_ \и+1 / \и+1 / \и+1
МД = = (и* )<+1 +("*),_!
(рД = МрД дрд = [(рД+1 +(р*)м ]/2'
(А^Г = (\)Г -М, (мХ1 =[(АЛк):!1 ]/2, аи,о4,о^[0,1].
По регуляризованным толщинам слоев определяются новые значения их уровней:
^г1), =0с0г +(4)ги+1,(z;;11)г = я,я/" =(^+1)г ,к = \,...,к,1 = \,...,кх.
Следует отметить, что нелинейная коррекция локальных инвариантов происходит не всегда, а только в областях резких градиентов. Фильтры (4) не приводят к понижению порядка аппроксимации.
Фаза 3. По найденным значениям потоковых переменных находятся величины консервативных переменных на новом временном слое:
Н:+1
С ,к
- сг,(чп:к -(ч
):+1
Ь,к
т/2
Ах
= 0,
(рА)С1-(р^СГ ,(рки)Тк -(рки)
:+1 'Ь,к
т/2
Ах
(р^с;: -мст | (^+1/2 )Г-(»л+1/2 )Г
= 0,
т/2
+-
:+1 :+1 :+1 :+1 РК,к+1 + РЬ,к+1 +1 - ¿Ь,к+1
Р
Ах
:+1
Ах
Я,к
:+1 :+1 :+1
2
Ах
Ах
= 0,
где /о* = 2о/:+1 + (1 - 2о)/:, о е[1/2, 3].
Разностные уравнения фазы 3 аппроксимируют исходные дифференциальные со вторым порядком по пространственной переменной и с первым -по времени. Если сложить эти уравнения с уравнениями фазы 1, то результирующая система при о = 0,5 будет иметь второй порядок аппроксимации и при отсутствии фильтрации (ои = о = ор = 0) обладать свойством временной обратимости. Изменяя параметр о, можно управлять диссипативными свойствами схемы. При о > 1 схема становится сверхнеявной [23], хотя вычислительный алгоритм остается явным. Данный алгоритм при постоянной плотности является сбалансированным: состояние покоя над произвольным донным рельефом не нарушается [24].
Выбор величины шага по времени. Условие устойчивости схемы КАБАРЕ для каждого слоя имеет вид [22]
( с0:+1/2 +(| иЛ1
1+1/2
ах,
< 1,
+1/2
откуда
т = СЕЬ х тт
/ чя+1/2 /| |\я+1/2
( ^ ,+1/2 +(1 "Л ,+1/2
Перестройка границ слоев. Фильтрация потоковых переменных, не затрагивающая консервативных величин, не делает полную задачу корректной на достаточно продолжительных промежутках времени. Приходится корректировать и консервативные переменные, что приводит к обмену массой и импульсом между слоями.
Перестройка сетки в предлагаемом алгоритме проводится два раза в одном временном цикле. Первый раз - после вычисления консервативных переменных на промежуточном слое в фазе 1, второй раз - после фазы 3.
На фазе 1 вычисляются промежуточные консервативные скорости и^1/^
и толщины слоев , по которым находятся новые высоты границ раздела
гуП+1/2 _ гуп 7 П+1/2 гуП+1/2 _ гуп _ ту гуП+1/2 _ Т_Т"И+1/2 7 _ 7 ДТ
А'+1/2,к = А'+1/2,к+1 + "г+1/2,к, ^¡+1/2,^+1 = А'+1/2, N+1 = в;+1/2 , А'+1/2,1 _ Н;+1/2 , к = 1,.", В общем случае эти координаты не являются окончательными и подлежат коррекции.
Рассмотрим три способа задания , что не исчерпывает всех
возможностей.
1. Высоты не перестраиваются и ¿П+1 = . В этом случае обмен
массой и импульсом между слоями отсутствует и задача остается некорректной.
2. Новые высоты границ раздела задаются таким образом, чтобы отношения толщин слоев в каждом вертикальном столбце оставались одинаковыми (сигма-сетка):
4п+1/2 _ уп , гп+1/2 ¿п+1/2 / гтп+1/2 _ ту ~\/дТ
^¡+1/2,к = ^г+1/2,к+1 + °к"+1/2 , = (Н¡+1/2 Вг+1/2 )/1У ,
N
к = 1,...,N, Ок>0, X°к= N.
к=1
3. Нижняя граница слоя к = 8, N > 8 > 1 остается неподвижной (эйлеровой), все границы слоев при к > 8 также неподвижны, а толщины слоев при к < 8 имеют заданные пропорции (сигма-сетка):
4»+1/2 _ 70 /. V . 7»+1/2 _ 7»+1/2 Г„+1/2 /. |
51 к=\
Во втором и третьем случаях законы сохранения приводят к обмену массой и импульсом между слоями. Потоки между слоями можно аппроксимировать как с первым, так и со вторым порядком точности.
Потоки первого порядка точности вычисляются методом донорной ячейки, который заключается в следующем. Рассмотрим самый нижний слой (к = N).
Если ЬП+гы — К+гпы, то толщина нижнего слоя уменьшится на величину
ДЯ& = - > 0 и часть ее массы Лтмп = Р^А/^Лхмп
и импульса А(ти).+1/2 = Рn++l1//22лUГ+ш2лА/т^Ах,.+1/2 перейдет в ячейку лежащего выше слоя. В результате получим:
"п+1/2 _ п+1/2 77п+1/2 — )/П+1/2 7п+1/2 _ п , £п+1/2
Р,+= р1+1/2^, г +1/2,N = и,+1/2, N, Z¡+1/2,N = В,+1/2 + 1/2^,
п 1/2 п 1/2 п 1/2 п 1/2
Ч+1/2 Р, +1/2, N-1й,+1/2, N-1 +Р1+1/2,N Ай,+1/2, N
Р,+1/2,N-1 , п+1/2 д 7 п+1/2 ,
"¡+1/2, N-1 + АП,+1/2, N
п+1/2 п+1/2 7 п+1/2 п+1/2 п+1/2 д 7 п+1/2
Лп+1/2 _ Р,+1/2,N-1г +1/2,N-1+1/2,N-1 + Р,+1/2,NUi+1/2,N Ай,+1/2,У и+1/2,N-1
п+1/2 7 п+1/2 п+1/2 д7
Р,+1/2, N-1",+Ш^Ч + Р,+1/2, N ЛП1
¡+
В противоположном случае, когда /£+ 1 /22ЛГ > 1 ^, верхняя ячейка отдает часть своей массы и импульса нижней:
"п+1/2 _ п+1/2 Лп+1/2 _ п+1/2 уп+1/2 _ п , £п+1/2
Р, +1/2,N-1 = Р,+1/2,N-1 , и,+1/2,N-1 = и,+1/2,N-1 , 1/2,N = +1/2 + "¿+1/2,N,
п 1/2 Р, +1/2,N :
п+1/2 1,п+1/2 _ п+1/2 к! п+1/2 Р, +1/2, N"1+1/2, N Р1+1/2, N-1^,+1/2, Л
/п+1/2 _ А/п+1/2 '4+1/2,N А/ ,+1/2,N
п+1/2 п+1/2 7 п+1/2 _ п+1/2 7/п+1/2 Дг,п+1/2 ^п+1/2 Рг+1/2, NU¡+1/2,Nh¡+1/2,N Р,+1/2, N-1^+1/2, N-1А^,1+1/2, N
и
п+1/2 ^п+1/2 _ п+1/2 Д7 п+1/2 Р¡+1/2,Nh¡+1/2,N Рг+1/2, N-1Аh¡+1/2,N
Отметим, что в методе донорных ячеек плотности и скорости считаются кусочно-постоянными функциями по вертикали.
Потоки второго порядка точности находятся из условия, согласно которому плотность и горизонтальная скорость являются линейными функциями на отрезке, соединяющем середины соседних слоев. Рассмотрим ячейки двух смежных слоев (с индексами В - нижний слой, Т- верхний слой). Обозначим через Аг = /в - /в гечшгей изменение границы между слоями, через Ах -ширину этих ячеек, через Ат, Ар - изменение массы и импульса слоев. Тогда получим следующие выражения:
ат = Рагах +1 ( ' + аг )—!——^— агах,
Нв 2^ в % +')/2
ар = оки„ агах +1 (' + аг) —агах. р Рв в г в ' (' + ')/2
Новые значения плотностей и скоростей слоев вычисляются по формулам:
Рг ' Ах + Ат рв ' Ах -Ат
Р^ , Рв ,
К ^шгеё Ах "В ^шгеё Ах
Рг Ах + Ар рв "¿ив Ах -Ар
и^ — Л , и £ — Л .
РThT,гeqшгedЛx РВЙВ,геди1геёАх
Пересчет консервативных переменных на других слоях осуществляется аналогично, от нижних слоев к верхним. Затем откорректированные консер-
/Л к+1/2
, принимаются за окончательные значения /к+112 и используются в последующей фазе 2. Аналогичная коррекция консервативных переменных осуществляется и после фазы 3.
Искусственная (турбулентная) вязкость. В задачах газовой динамики для регуляризации численного решения при наличии ударных волн используется метод искусственной вязкости [25]. В численный алгоритм добавляется разностный оператор, аппроксимирующий так называемую вторую - объемную вязкость, коэффициент при которой пропорционален пространственному шагу сетки и находится из условия наибыстрейшего затухания самой высокочастотной гармоники [26]. В областях волн разрежения, где дивергенция скорости положительна, искусственная вязкость обнуляется. Это обеспечивает монотонность численного решения на волнах сжатия и не портит волн разрежения.
В случае многослойной гидродинамики включение линейной искусственной вязкости сводится к модификации в фазе 1 и фазе 3 потоковой сеточной функции давления на вертикальных гранях расчетной сетки. Для фазы 1:
(Рк+1/2 ) ^ (Р+1/2 )П = (Рк+1/2 ) " 0ку (РС)] (<у+1/2 " <у-1/2 ) ,
д* _|0 К,/ +1/2 - К,у—1/2 < 0, У ""[О г/ < у +1/2 - «к, у—1/2 ^ 0,
где с. = ф) | /рк - локальная скорость звука, 9 ~ 1 - безразмерный параметр.
Для фазы 3:
(Р+1/2Г ^ (^+1/2) Г" (Р+1/2Г — б"к, у (рС2 (иут — <+—у,
* _ I 0 <+/+1/2 — ик+'-1/2 < 0
к' [О г/ — ик+--2/2 * 0.
Обобщение на случай двух пространственных измерений. Переход от описанного выше одномерного алгоритма к многослойному двумерному осуществляется аналогично тому, как это описано в работе [27] для однослойной мелкой воды на сфере. Учет силы Кориолиса в схеме КАБАРЕ, оперирующей двумя типами переменных (консервативные и потоковые), не вызывает затруднений и не накладывает дополнительных ограничений на шаг по времени.
Верификация алгоритма
Целью верификации является проверка алгоритма на отсутствие ошибок программирования, вычислительную устойчивость и сходимость на задачах, допускающих сравнение с аналитическим решением либо с численным решением, полученным другими авторами. Для одного слоя и постоянной плот-
ности предлагаемый алгоритм совпадает с описанным в [24]. При наличии одного слоя и переменной плотности новый многослойный алгоритм совпадает с описанным в [28].
Особенности многослойных моделей проявляются уже на задачах с двумя или тремя слоями [29-31]. В этих работах приведены примеры расчетов модельных задач, на которых была протестирована и наша методика. Следует отметить, что в указанных работах приводятся результаты решения на определенный момент времени на фиксированных расчетных сетках и не изучается поведение численного решения на больших временных интервалах и на сетках различной густоты.
Проблемы с устойчивостью численного решения многослойных уравнений мелкой воды при отсутствии обмена массой и импульсом между слоями проиллюстрируем на двухслойной задаче А. Курганова из [30]. Начальные данные (рис. 3) следующие:
/ ч / ч С1 / 2 >\А > 1
щ (Х,0) = 0,4, р = 0,98, к (х,0) = <! . . ,,
и ' 1 и ' [1 - 0,25 8т (2лх) 1/ |х| <1,
/ ч / ч С1 1/ 2 > IX > 1,
щ (х, 0) = -0,4, р, = 1, к (х,0) = < , ч ,,
2К ' П ' [1 + 0,258т(2лх) 1/ IX <1.
Рельеф дна В(х) = -2, х е [-2, 2], g = 10. Расчетная сетка - 801 узел по оси x.
Р и с. 3. Начальные данные Fig. 3. Initial conditions
Результаты расчетов по схеме КАБАРЕ показаны на рис. 4. Расчеты проводились при следующих параметрах: аи = ай = а = 2/3, а* = 3, 0 = 0,
СЕЬ = 0,3. Графики границ раздела слоев и их скоростей на момент времени t = 0,5 представлены на рис. 4, a, Ь. Они практически совпадают с приве-
денными в работе [30]. При продолжении расчета происходит вырождение толщин слоев и при ^ = 0,65 расчет прекращается (рис. 4, с,
Р и с. 4. Расчет по схеме КАБАРЕ на сетке в 800 ячеек: границы и скорости слоев на момент t = 0,5 (a, b) и на момент t = 0,65 (c, d)
Fig. 4. Calculation on the 800 cell grid by the CABARET scheme: boundaries and velocity of the layers at t = 0.5 (a, b) and t = 0.65 (c, d)
Если использовать более густую сетку с числом узлов 1601, то потеря устойчивости происходит значительно раньше. На более грубой сетке (при N= = 201) неустойчивость развивается медленнее (рис. 5), но плохая обусловленность задачи неизбежно приводит к аварийному прекращению расчета.
Перестройка расчетной сетки и возникающий при этом обмен массой и импульсом между слоями в значительной степени регуляризируют задачу [15, 32-34], однако и в этом случае надлежащий выбор параметров фильтрации и искусственной вязкости остается весьма важным.
Р и с. 5. Расчет по схеме КАБАРЕ на сетке в 200 ячеек: границы и скорости слоев на момент t = 0,5(a, b) и на момент t = 1(c, d)
Fig. 5. Calculation on the 200 cell grid by the CABARET scheme: boundaries and velocity of the layers at t = 0.5 (a, b) and t = 1 (c, d)
Расчет баротропного течения по многослойной модели с учетом обмена массой и импульсом между слоями
Рассмотрим задачу о колебаниях жидкости в бассейне с ровным дном и возмущенной в начальный момент формой свободной границы. Параметры задачи, начальные и граничные условия определены следующим образом:
х е[-5,5], B(х) = -2, u(x, t0) = 0, u(-5, t) = u(5, t) = 0, p = 1, g = 10,
H (x, t0) = ^
1 ( (2nx^ 5 5 —I 1 + cos I--if — < x <—,
2 ^ ^ 5 )) 2 2
0 else.
Расчет проводился по схеме КАБАРЕ по однослойной модели и при числе слоев N = 10 на сигма-сетке с одинаковыми толщинами слоев до момента времени I = 6, соответствующего более чем одному периоду колебания по-
верхности. На рис. 6, а приведена форма свободной поверхности в начальный момент, на рис. 6, Ь - свободная поверхность при I = 6 и числе расчетных ячеек N = 128.
Р и с. 6. Форма поверхности жидкости: a - начальные данные; b - при t = 6. Сплошная линия - однослойная вода, маркеры - многослойная вода с числом слоев N = 10 Fig. 6. Form of the liquid surface: a - initial conditions, b - at t = 6. Solid line denotes one-layer water, the markers - multilayer water (number of the layers N = 10)
Расчеты проводились при параметрах ои = oh = о = 2/3, о* = 2, 0 = 0,
CFL = 0,3 . При увеличении времени расчета и сгущении сетки алгоритм не теряет устойчивости и демонстрирует сходимость по сетке.
Заключение
Многослойные гидростатические модели течений со свободной поверхностью без учета обмена массой и импульсом между слоями в большинстве случаев являются некорректными и не могут использоваться в практических расчетах. Перестройка расчетной сетки на каждом временном шаге способна регуляризировать задачу, однако для обеспечения гарантированной устойчивости необходимо использовать фильтрацию потоковых переменных и искусственную вязкость, моделирующую турбулентное перемешивание. Дальнейшее развитие предложенной модели будет направлено на разработку правил применения описанных методов регуляризации решения, минимизирующих численную диссипацию.
СПИСОК ЛИТЕРАТУРЫ
1. Марчук Г. И., Дымников В. П., Залесный В. Б. Математические модели в геофизической гидродинамике и численные методы их реализации. Л. : Гидрометеоиздат, 1987. 296 с.
2. Numerical simulation of large-scale ocean circulation based on the multicomponent splitting method / V. B. Zalesny [et al.] // Russian Journal of Numerical Analysis and Mathematical Modelling. 2010. Vol. 25, iss. 6. Р. 581-609. https://doi.org/10.1515/RJNAMM.2010.036
3. Информационно-вычислительные технологии - новый этап развития оперативной океанографии / Г. И. Марчук [и др.] // Известия Российской академии наук. Физика атмосферы и океана. 2013. Т. 49, № 6. С. 629-642. https://doi.org/10.7868/S0002351513060114
4. Доценко С. Ф., Залесный В. Б., Санникова Н. К. В. Модульный подход к расчету циркуляции и приливов в Черном море // Морской гидрофизический журнал. 2016. № 1. С. 3-19. https://doi.org/10.22449/0233-7584-2016-1-3-19
5. Залесный В. Б., Гусев А. В., Фомин В. В. Численная модель негидростатической морской динамики, основанная на методах искусственной сжимаемости и многокомпонентного расщепления // Океанология. 2016. Т. 56, № 6. С. 951-971. https://doi.org/10.7868/S0030157416050178
6. Новые алгоритмы вычислительной гидродинамики для многопроцессорных вычислительных комплексов / В. М. Головизнин [и др.]. М. : Изд-во МГУ, 2013. 467 с.
7. Audusse E. A multilayer Saint-Venant model: Derivation and numerical validation // Discrete & Continuous Dynamical Systems - B. 2005. Vol. 5, iss. 2. P. 189-214. https://doi.org/10.3934/dcdsb.2005.5.189
8. Audusse E., Bristeau M.-O. Finite-Volume Solvers for a Multilayer Saint-Venant System // International Journal of Applied Mathematics and Computer Science. 2007. Vol. 17, no. 3. Р. 311-320. https://doi.org/10.2478/v10006-007-0025-0
9. Овсянников Л. В. Модели двухслойной "мелкой воды" // Прикладная механика и техническая физика. 1979. Т. 20, № 2. С. 3-14.
10. Duchene V. A note on the well-posedness of the one-dimensional multilayer shallow water model. 2013. URL: https://hal.archives-ouvertes.fr/hal-00922045/document (date of access: 12.07.2019).
11. Chandrasekhar S. Hydrodynamic and Hydromagnetic Stability. London : Oxford University Press, 1961. 652 p.
12. Чухарев А. М., Руновский К. В., Кульша О. Е. Моделирование статистического распределения турбулентных пятен в стратифицированных слоях океана // Морской гидрофизический журнал. 2017. № 5. С. 35-46. https://doi.org/10.22449/0233-7584-2017-5-35-46
13. Методы расчета турбулентных течений / Дж. Ламли [и др.]. М. : Мир, 1984. 463 с.
14. A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation / E. Audusse [et al.] // ESAIM: Mathematical Modelling and Numerical Analysis. 2011. Vol. 45, no. 1. P. 169-200. https://doi.org/10.1051/m2an/2010036
15. A fast finite volume solver for multi-layered shallow water flows with mass exchange / E. Audusse [et al.] // Journal of Computational Physics. 2014. Vol. 272. P. 23-45. https://doi.org/10.1016/jjcp.2014.04.026
16. Hirt C. W., Amsden A. A., Cook J. L. An Arbitrary Lagrangian Eulerian computing method for all flow speeds // Journal of Computational Physics. 1974. Vol. 14, iss. 3. P. 227-253. https://doi.org/10.1016/0021-9991(74)90051-5
17. Ringler T. D, Randall D. A. The ZM Grid: An Alternative to the Z Grid // Monthly Weather Review. 2002. Vol. 130, no. 5. P. 1411-1422. https://doi.org/10.1175/1520-0493(2002)130<1411:TZGAAT>2.0.m;2
18. Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model: Kinetic interpretation and numerical solution / E. Audusse [et al.] // Journal of Computational Physics. 2011. Vol. 230, iss. 9. P. 3453-3478. https://doi.org/10.1016/jjcp.2011.01.042
19. Stewart A. L., Dellar P. J. Multilayer shallow water equations with complete Coriolis force. Part 1. Derivation on a non-traditional beta-plane // Journal of Fluid Mechanics. 2010. Vol. 651. P. 387-413. https://doi.org/10.1017/S0022112009993922
20. Bermudez A., Vazquez Ma. E. Upwind methods for hyperbolic conservation laws with source terms // Computers & Fluids. 1994. Vol. 23, iss. 8. P. 1049-1071. https://doi.org/10.1016/0045-7930(94)90004-3
21. Toro E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Berlin : SpringerVerlag, 2009. 724 p. doi:10.1007/b79761
22. Karabasov S. A., Goloviznin V. M. Compact accurately boundary-adjusting high-resolution technique for fluid dynamics // Journal of Computational Physics. 2009. Vol. 228, iss. 19. P. 7426-7451. https://doi.org/10.1016/jjcp.2009.06.037
23. Головизнин В. М., Канюкова В. Д., Самарская Е. А. Сверхнеявные разностные схемы газовой динамики // Дифференциальные уравнения. 1983. Т. 19, № 7. С. 1186-1197. URL: http://www.mathnet.ru/links/23381b18a4caaff660036527679f7476/de4901.pdf (дата обращения: 12.07.2019).
24. Головизнин В. М., Исаков В. А. Применение балансно-характеристической схемы для решения уравнений мелкой воды над неровным дном // Журнал вычислительной математики и математической физики. 2017. Т. 57, № 7. С. 1142-1160. https://doi.org/10.7868/S0044466917070092
25. Рихтмайер Р. Д., Мортон К. Разностные методы решения краевых задач. М. : Мир, 1972. С. 300-304.
26. Головизнин В. М. Об одном способе введения искусственной диссипации в вариационно-разностные схемы магнитной гидродинамики // Журнал вычислительной математики и математической физики. 1982. Т. 22, № 1. С. 144-150. URL: http://mi.mathnet.ru/zvmmf5796 (дата обращения: 11.07.2019).
27. Goloviznin V. M., Solovjov A. V., Zalesny V. B. A new algorithm for solving the shallow water equations on the sphere based on the cabaret scheme // Journal of Physics: Conference Series. 2018. Vol. 1128. 012091. https://doi.org/10.1088/1742-6596/1128/1/012091
28. Evtushenko Y. G., Gorchakov A. Y., Goloviznin V. M. Fast automatic differentiation in problems variations four-dimensional data assimilation (4Dvar) // Journal of Physics: Conference Series. 2018. Vol. 1128. 012001. https://doi.org/10.1088/1742-6596/1128/1/012001
29. Ведерников А. Б., Холодов А. С. Численное моделирование течений двух- и трехслойной жидкости в рамках модели мелкой воды // Математическое моделирование. 1990. Т. 2, № 6. С. 9-18. URL: http://mi.mathnet.ru/mm2375 (дата обращения: 11.07.2019).
30. Kurganov A., Petrova G. Central-Upwind Schemes for Two-Layer Shallow Water Equations // SIAM Journal on Scientific Computing. 2009. Vol. 31, iss. 3. Р. 1742-1773. https://doi.org/10.1137/080719091
31. Елизарова Т. Г., Иванов А. В. Квазигазодинамический алгоритм численного решения двухслойных уравнений мелкой воды // Препринты ИПМ им. М. В. Келдыша. 2016. № 69. 27 с. doi:10.20948/prepr-2016-69
32. Shankar N. J., Cheong H. F., Sankaranarayanan S. Multilevel finite-difference model for three-dimensional hydrodynamic circulation // Ocean Engineering. 1997. Vol. 24, iss. 9. P. 785-816. https://doi.org/10.1016/S0029-8018(96)00036-4
33. Bouchut F., Zeitlin V. A robust well-balanced scheme for multi-layer shallow water equations // Discrete and Continuous Dynamical Systems-Series B. 2010. Vol. 13, iss. 4. P. 739758. http://doi.org/10.3934/dcdsb.2010.13.739
34. Couderc F., Duran A., Vila J.-P. An explicit asymptotic preserving low Froude scheme for the multilayer shallow water model with density stratification // Journal of Computational Physics. 2017. Vol. 343. P. 235-270. https://doi.org/10.1016/j.jcp.2017.04.018
Об авторах:
Головизнин Василий Михайлович, заведующий отделом, Институт проблем безопасного развития атомной энергетики РАН (115191, Россия, г. Москва, ул. Большая Тульская, д. 52), доктор физико-математических наук, профессор, ResearcherID: Y-2025-2018, gol@ibrae.ac.ru
Майоров Павел Александрович, инженер, Институт проблем безопасного развития атомной энергетики РАН (115191, Россия, г. Москва, ул. Большая Тульская, д. 52), Researcher-ID:AAD-6169-2019, pavel.a.mayorov@gmail.com
Майоров Петр Александрович, инженер, Институт проблем безопасного развития атомной энергетики РАН (115191, Россия, г. Москва, ул. Большая Тульская, д. 52), Researcher-ID: AAD-6173-2019, maiorov.peter@gmail.com
Соловьев Андрей Валерьевич, ведущий научный сотрудник, Институт проблем безопасного развития атомной энергетики РАН (115191, Россия, г. Москва, ул. Большая Тульская, д. 52), кандидат физико-математических наук, ResearcherID: X-4858-2019, solovjev@ibrae.ac.ru
