CC BY
19Original Russian Text O The Authors, 2019 published in MORSKOY GIDROFIZICHESKIY ZHURNAL, Vol. 35, Iss. 6 (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. Solovjev
Nuclear Safety Institute, Russian Academy of Sciences, Moscow, Russian Federation
* 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.
Conclusions. 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
DOI: 10.22449/1573-160X-2019-6-528-546
© 2019, V. M. Goloviznin, Pavel A. Maiorov, Petr A. Maiorov, A. V. Solovjev © 2019, Physical Oceanography
Introduction
The achievements of recent decades in the field of modeling and forecasting of sea currents are primarily associated with the rapid development of highperformance computing systems and new mathematical algorithms [ 1-6]. Modern problems of ocean hydrodynamics are described by nonlinear partial differential equations and are solved in complex regions with high spatial resolution [3, 4].
528 ISSN 1573-160X PHYSICAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019)
(J) ® The content is available under Creative Commons Attribution-NonCommercial 4.0 nw^ International (CC BY-NC 4.0) License
In many cases, these problems are not classical and require the development of special methods for their spatial approximation, as well as methods for solving by time, including explicit, implicit and splitting schemes [1, 2, 5].
One of the main advantages of using explicit schemes is their efficient parallelization on multi-core computing platforms. In order to solve prognostic problems, the gain due to the efficiency of parallelization often exceeds the loss due to a decrease in the time step.
For modeling the hydrodynamics of the seas and oceans, the multilayer hydrostatic approximation is widely used. Despite the rich history of development and a large number of numerical implementations, the task of further improving computational algorithms for the equations of multilayer hydrodynamics with a free surface remains relevant [7, 8].
The equations of multilayer shallow water are not unconditionally hyperbolic [9, 10]. This leads to the fact that the initial-boundary value problem in the process of its solving becomes incorrect and the algorithms for solving unconditionally hyperbolic equations lose stability *. It is considered that the loss of hyperbolicity occurs with the development of Kelvin - Helmholtz instability at the layer interfaces [11]. In these cases, an intense exchange of mass and momentum between the layers should occur, which is forbidden in the classical multilayer approximation [8]. This disadvantage can be eliminated by the inclusion of the so-called turbulent viscosity in multilayer equations. It often depends on the parameters of the calculated current and contains empirical parameters setting to various types of currents [12, 13]. The disadvantage of this approach is the critical dependence of the calculation results on the calculator experience.
Another approach to the solution of this problem is to modify the classical equations of multilayer shallow water - to remove the prohibition of mass and momentum exchange between layers [14, 15]. In fact, this means a transition from the Lagrangian description of the vertical motion of the layer boundaries to the mixed Euler-Lagrangian description [16]. In this case, the fluxes of mass and momentum arise between the layers, which regularize the problem.
When implementing this approach, we use the method of splitting into the physical processes. Firstly, the values of physical variables on a new time layer are calculated according to the classical equations of multilayer shallow water with a density variable in each layer. Secondly, the predetermined vertical coordinates of the boundaries of layers are specified and the exchanges of mass and momentum fluxes that arise between them are found. When constructing the prescribed vertical boundaries of the layers, it is assumed that the upper layer with a free boundary remains Lagrangian and the lower one remains motionless. If the coordinates of the intermediate boundaries for each vertical are proportional to each other, then such computational grids can be called sigma-coordinate [1, 17]. Other oftenly used in oceanology is z-coordinate, when the boundaries of all layers, except for the upper one, remain motionless [1, 17].
For the numerical solution of the classical equations of multilayer shallow water with different layer densities, many algorithms that have their own advantages and disadvantages [18, 19] were developed. Usually, computational algorithms are based on the finite volume method, and some variant of downstream transferring of the values is used to calculate convective flows [20].
* Kabanikhin, S.I., 2009. Inverse and Ill-Posed Problems. Novosibirsk: Siberian Scientific Publishing House, 457 p.
It should be noted that classical methods based on solving the Riemann Problem for one layer shallow water [21] are not applicable in the multilayer case. This is due to the fact that the Riemann Problem in the multilayer case does not have a simple solution.
In the proposed new numerical technique, at the first stage of splitting by extrapolation of local Riemannian invariants horizontal fluxes are calculated, as is done in CABARET scheme [6]. For this, the multilayer model is represented as the sum of single-layer (hyperbolic decomposition) ones interacting by means of reaction forces applied to the interface. For each of the single-layer models, CABARET scheme with nonlinear flux correction based on the maximum principle is recorded [22]. In this case, the incompressibility condition and the laws of mass and momentum conservation are satisfied. As a result (in a new time layer) there are flux and conservative variables without taking into account the exchange of mass and momentum between the layers. Since this task is incorrect, new flux variables undergo filtering that does not violate conservation laws. The resulting algorithm has a second order of approximation and has the well-balance property -to maintain resting state over any bottom relief at a constant density [23].
At the second stage of splitting into physical processes the vertical boundaries of the layers are set and the fluxes between the layers are calculated.
In order to expand the stability domain of the algorithm, two procedures are used that lower the order of approximation in time to the first. This is the inclusion of linear viscosity, which leads to the fastest attenuation of the highest frequency harmonic itself [24], and the super implicit approximation of the pressure gradient [25].
Equations of single-layer shallow water with variable density, with regard to external pressure and bottom relief
In the case of one spatial measurement, we write out the integral conservation
(balances) laws for an arbitrary subdomain x e[ x1, x2 ] (Fig. 1). We introduce the following notation:
V = J hdx, M = J phdx, n = j puhdx, h = H - B,
where Vis an area occupied by liquid; M is a mass; n is momentum H(x, t) is a free surface; B(x) is a bottom relief; h(x, t) is a layer thickness; p(x, t) is a density; u(x, t) is horizontal velocity of liquid.
The equation for variation of the area occupied by the liquid has the following form
dV = —[hdx = -uhl + uh\ =- [2Ldx. dt dt x x x dx
Fig. 1. One-layer model
The variation of mass and momentum is described respectively by the equations
dM d
= — J phdx = - p + p uh|_^ =-|
dphu
dt dt
11 J dx
dx,
-= — Jpuhdx = -(pu2h + pgh2 /2 + PTh) +(pu2h + pgh2 /2 + PTh)
dt dt
x2 dB x2
J p -dx + J P
dH
dx
dx.
Here g is a free fall acceleration; Pt is a pressure on a free surface; Pb = Pt + pgh is a pressure on the bottom.
The corresponding conservation laws in differential form follow from the given integral equations:
5ph dphu
dh dhu
— +-= 0,
dt dx
dt dx
■ = 0,
dpuh dphu2 g dph2 dPTh ^ „ dH „ dB —— + —-+ -—-— + —— = F = PT--PB —.
dt dx 2 dx dx dx dx
After some transformations these equations can be represented as ^ + A x^d^ = d, ( = (h,p, u )T, d = (0,0, F/ph),
where
A =
A
u 0
0 h u 0 gh dp
Pb (ph )-1 ^
bk ' 2p dx
u
f=p dH - p dB-„dPT
dx
dx dx
All matrix A eigenvalues are valid:
\ = u + c, A,2 = u - c, A,3 = u, c = iJPB/p,
whence the so-called characteristic form of the initial equations follows:
„v
du c dh — +--+
dt h dt 2pc dt
gh dp , J du + (u + c
du dt
c dh h dt
gh dp 2pc dt
/ v
A / Jcu + (u - c)
c dh gh dp
— +--+ —---
dx h dx 2pc dx
v
dx
c dh h dx
gh dp 2pc dx
= ( F/ph ), = ( F/ ph),
dp dp
— + u — = 0. dt dx
Equations of multilayer shallow water with variable density with regard to external pressure and bottom relief
The equations of multilayer shallow water with variable density in the absence of mass and momentum exchange between the layers can be represented as the sum of single-layer equations connected by pressure and reaction forces affecting the interface (Fig. 2):
d(phu\
d t k = 1,
, d(huI
dt dx
d(p hu 2 I
= 0,
d(ph2)
5pA , d(phu )k d t
, d(hkPk )
dx
= 0,
, N,
dx Z = H,
2
dx = B,
dx
P
dZ,
dx
■- P
dZt
hk = Zk Zk+1,
P = P
1 k+1 1 k
dx
+ gp A.
= 0, (1)
Here k is a layer number counted from the free surface; Zk is the coordinates of the upper boundary of the layer; P\ = Pt is an external pressure on a free surface Z1. The so-called simple form of the system of equations (1) will have the form:
- _ , st
— + G = D, ¥ = (h1,■■■, hN ,p1,-,pN , U1,■■■, uN ) ,
where G is a matrix of N * N dimensions and D is some right part. It is known well that already at N = 2 the matrix G can have complex roots and the system will not be unconditionally hyperbolic [9], which gives rise to known computational difficulties. For a larger number of layers the situation is exacerbated. Already in the presence of two layers the calculation of the eigenvalues of the matrix (4 * 4) is a rather difficult task. With a large number of layers it becomes practically insoluble and the direct use of balance-characteristic methods [23] is impossible. For overcoming this difficulty we use a technique that we will call the hyperbolic decomposition of the problem.
Fig. 2. Multilayer model 532
If we consider the forces affecting the interface to be external, then (1) can be represented as:
d(pk dt
"Ak x
d(pk dx
= dk, (k = ( hk ,pk, uk), dk = ( 0,0 Fkl pkhk), k = 1,..., N,
where
Ak =
0
Pk+1 ( ph )
0 K
u,, 0
2pk
F = P — - P
dx
dZk+1 - dPL k+1 ~ 'lk~
dx dx
which leads to a system of independent characteristic equations:
duk , ck dhk , ghk %
dt hk dt 2pkck dt
f duk - dk - ghk % dt hk dt 2pkck dt
dpk , ., = 0
+ (uk + ck)
+ (uk - ck)
duk , ck dh , ghk dpk
dx hk dx 2pkck dx
du
ck dh ghk dpk ^
______
dx hk dx 2pkck dx
= ( Fkl pkhk), = (FjpA ), (2)
dt
u
ex
This makes it possible to use algorithms that have proven themselves in the single-layer case for the numerical solution of a complete system of multilayer equations. The incorrectness of the complete system will not disappear and will manifest itself in high-frequency distortion of the boundaries. Using special filters the distortion can be significantly reduced while maintaining the approximation order and without violating the conservation laws.
CABARET scheme for the layered solution of the equations of multilayer shallow water
We cover the domain of the problem solution with an uneven computational grid with xt, i = 1, Nx coordinates of the nodes. In the nodes at the initial moment of time we set the vertical coordinates of nki layers:
(Z0)i= nk,, $)= H0, (Z04= B,, i = 1,..., N.
The constant coordinates of the nodes along the x axis and the time-dependent vertical coordinates of Znki layers determine the computational grid on the (x, z) plane.
Two types of variables are used in the CABARET scheme: flux and conservative [22] ones. Flux variables refer to the midpoints of the vertical faces of computational cells, conservative ones - to their centers. Flux variables related to the layer with the number k will be denoted as p"ki, u"ki, hnki = Znki - Zl+li,
conservative - as pnki+1/2, unki+1/2, hlt+1/2. At the initial moment flux variables are set. Conservative variables are determined by flux ones as follows:
P0,i+1/2 = 0 5 (Pi + P0,i+1 ) , Uk,i+1/2 = 0 5 (i + i+1 ) , Ki,+1/2 = 0,5 ( hi, + h0 i + ).
The computational algorithm consists of the following elements: calculation of conservative variables on an intermediate time layer according to conservative finite-volume schemes (phase 7); calculation of flux variables on a new time layer by extrapolation of local Riemann invariants (phase 2); calculation of conservative variables on a new layer using finite-volume schemes (phase 3).
Phase 1. In order to calculate the intermediate conservative variables, we use:
klÏ2 - hekk.(hu L-(hu I
= 0,
t/2 at
(phfi/2 - (ph)c,k (phu )"R,k- (phu )
L,k
= 0,
T/2 AX
(phur/2 -(phu)n _ (Phu2)" -(Phu21
(hkPk+1/2 )R (hkPk+1/2 )L
t/2
Ax
Ax
T>n I T>n r7n r7n T)n I T>n ryn ryn
PR,k+1 + PL,k+1 ZR,k+1 — ZL,k+1 PR,k + PL,k ZR,k — ZL,k
AX
Ax
= 0,
where Pk+1/2 = (P + 0,5pgh. These equations have the first order of
approximation in time and the second in spatial variable.
Phase 2. In order to find new flux variables, a linearized system of characteristic equations (2), which can be written in the following form, is used:
u h ck L+1/2 dhk + ghk L+1/2
dt \ . hk _ .+1/2 dt _ 2PkCk _ i+1/2
dPk dt
\
+ [(uk +ck )]
L + 1/2 i +1/2
duk + ck L+1/2 dK
dx \ . hk _ i+1/2 &
2PkCk
L+1/2 ^
dPk
,+1/2 &
= [( FJ P khk )]
L+1/2 i+1/2 '
duk - ck L+1/2 dhk ghk
dt . K. +1/2 dt _ 2PkCk _
L+1/2
i+1/2
dPk dt
\
+[( uk -ck )]
L+1/2 i+1/2
duk - ck L+1/2 % __ ghk l+1/2 A dPk
dx \ . hk _ i+1/2 & _ 2PkCk _ i+1/2 & ,
= [( FJ P A )]
L+1/2 i+1/2 '
^ + \uk ]L+1/2 ^ = 0, dt L kJl+1/2 cX
which can be written down as
where
k,i+1/2 , \n+H2 A ) k,i+1/2 I r\ \"+1// +( x 1 )i+i/2 a, =( ßu).+1/2 ,
d(l2 )k,i+1/2 y+1/2 d( 12 )k,i+1/2 \"+1/2 ■ + ( 12 )i+1/2 -^-= (ß2k )i+1/2 '
dt v /i+1/2 dx
d( h ) k,i+1/2 , s.n+1/2
k+1/2
+( ^ c2^12=(ß3,k )
a v J/i+1/2 ex
i+1/2
( L ) , (/, ) , (L ) are local Riemann invariants in a cell i + 1/2 of k
V I Jk,i+1/2 ' 2 >k,i+1/2 ' V3 )k,i+1/2
layer at time interval [t:, t:+1 ] :
( I1 )
= 11 + G:+1/2 k ,i+1/2 Uk^ i+1/2,k
:+1/2 ^ _ :+1/2
:+1/2
_ "k Gi+1/2,k
hk + D^Pk , (12 )
hk _ Di+112,kpk , (/3 )k i+1/2 = Pk
k ,i+1/2
G
:+1/2
D
:+1/2
i+1/2
2PkCk
:+1/2
i+1/2
Note that at variable density the eigenvalues of Ak matrix remain the same as with constant density [24], while left eigenvetctors are modified which leads to the appearance of an additional term proportional to the density in local invariants.
The first action of phase 2 consists in calculating the preliminary values of local invariants on a new time layer using the linear extrapolation method taking
into account the sign of the characteristic velocity (Xm Y^l, m = 1, 2, 3 :
( îm )
:+1 k ,i+1
' k ,i
_ / \:+1/2 / j \î V m)k,i+1/2 V m)i •\:+1/2 / j \:
k,i+3/2 V m)k,i +2 \:+1/2 / j \:+1/2 'k ,i +3/2 +( m ) k,i+1/2
2 ( Im I ( Im I
f ( 1 m ):+++12/2 > 0 a:d ( 1m \
\:+1/2
if ( 1 m )
/2 if ( K c
k, i+1/2 V m/k ,i +3/2
\:+1/2 „ 7 /a \:+1/2
, < 0 and ( 1 m ),
k,i+ 3/2 v m / k,i+1/2
> 0, < 0,
'A 1m )
^:+1/2 k, i+1/2
< 0.
The next step consists in the correction of the obtained values in accordance with the maximum principle for equations (3). We calculate the minimum and maximum values of local invariants on the current time layer for each cell:
min (^ )+1/2 = mm {(II, )M, (II, );+i/2, (II, )}, max (^ )+i/2 = max {(C,k )M,(C,k (C,k )},
min (Im,k )+i = min {min (lm,k )i+i/2 , ^ (lm,k ),+3/2}' maX (Im,k ),+i = maX {maX (Im,k ),+i/2 , ™X (lm,k ),+3/2 }
and the right sides of the transport equations (3) from the known left sides:
(r+i/2) -(I") (I") -(r)
/ \«+l/2 V m )k,i+i/2 V m)k,i+i/2 , /» N«+i/2 V m)k,i+i V m>k,i 2\
(Q )i+i/2 =-T/2-+ (^ii+i/2 -^-+ 0(T )'
/ \n+i/2 m = l2, Q ),,/2=0.
From the transport equations by characteristics (3) it follows that if
/ \«+i/2
(Qm k) values are equal to zero and the Courant - Friedrichs - Lewy number (CFL) is less than unity, then the following conditions must be satisfied:
( 4 )
n +1
m Jh ,i+1
in (Im,k )I+1/2 , max (Im,k )I+1/2 J if (^m )k1/2 > 0 aLd (^m )k^+3/2 ^ 0,
in (Im, )I+ 3/2 , max (Im,k ), + 3/2 ] f (^m )kkI1+/32/2 < 0 ^ (^m )k^2 ^ 0,
min (Im,k ),+1 , max (Im,k ),+1 ] if (^m )kkI1+/32/2 X (^m )k^2 < 0.
/ \ l+1/2
For non-zero (Qmk ) intervals of permissible values are shifted:
max (Im,k )i +1/2 = max (Im,k )I+1/2 + T (Qm,k ^ ,
min (Im,k )i+1/2 = min (Im,k ^ + T (Qm,k ^ . \ l+1
Im ) invariants consists in the following: if the preliminary value
/ k ,i+1
\ l+1
I ) invariant satisfies these inequalities, then the invariant does not change;
/ k ,i+1
if not, it is assigned a value corresponding to the closest boundary of the allowable interval. According to the adjusted values of local invariants on a new time layer from the system of equations
/ \L+1 / T \L+1
(Pk )i+1 =( I3 )ki+1,
(uk )L++11 + G*k ( hk )L++11 + D* k ( Pk )L++11 = ( I1 )L+i+1, (uk )L++11 - G*k (hk )L++i - D2,k (Pk )L++i = (I2 )
n+1 h ,i+1
new flux variables are calculated:
/ \;+1 / j \ñ+1/ \;+1 (pk ),■+! = ( h ),,,+! , (U ),+! =■
(*\;+1 / *\;
G k(^ c+(/;)
* \ ;+1 k ,1+1
G1,k + G2,k
(hk ^ =
\ ;+1
'k ,i+1
G1, k + G2, k
where
(7* );+1 = (7i CI 1 " D* (Pk;, (* = (^);++1 + D* (Pk;, G1*k, G2* k take the values:
(Gm,k , Dm,k ) :
fe, D+Z ) if (^m C2 > 0 a;d (Xm ^ > 0, (G^, D^k ) if (X„ C2 < 0 a;d ^ < 0,
(G^, D-f ) if (Xm )k Z X(^m )k+TL < 0, where m=1,2, G-f=(g^+G^ ) /2, D-f=(D-- + D^ ) /2.
Filtering of flux variables. In the areas of hyperbolicity loss, short-wave disturbances are generated. To regularize the solution, the calculated flux variables are filtered:
/2,
(4)
\ n+1 / \ n+1 /1 \ /_ \ n +1 /_ \ n+1 / \n+1 / \n+ J
ukI = Ou (Uk\ +(1 - o„)(uk\ , (MkI =[(uk)i+1 +(uk)i_i
¡~ \n+1 / \n+1 W— \n+W— \n+1 f/ \n+1 / \n+1~|
(PkI = MPkI +11 -op)(PkI ,(Pk),. =L(pkL + (PkL J/2, (hk )n+J =(hk )n+1 +(Ahk f, (Ahk f = Oh (Ahk )n+ +(1 - Oh )(Ahk r , (Ahk)n+1 =(hkr-(hk)n, (Ahk)n+1 =[(Ahk)n: +(Ahk)n+1 ]/2, Ou,Oh,o^[0,1].
Regularized thicknesses of the layers determine new values of their levels:
(zr1) =( z£ J +(hk )n+J, ( =b,, h"+1 =( ZJ ), k = 1,..., N,t = 1,..., n .
It should be noted that nonlinear correction of local invariants does not always occur, but only in the areas of sharp gradients. Filters (4) do not lead to a decrease in the approximation order.
Phase 3. According to the found values of flux variables, the values of conservative variables are found on the new time layer:
hff - K7 , (hu -(hu= 0
t/2 ar '
(ph JCl -(phET. (phu C-(phu }
= 0,
t/2 Ar
\n+1 / 0\n+1
(phu )nc+;-(phu c _ (Phu 2):;-(Phu 2);; _ (Kpk+V2 %_(hkpk+y21
t/2 Ar Ar
pn + 1 . pn + 1 ryn + 1 ryn + 1 pn + 1 , pn + 1 ry n + 1 ryn + 1
+ P:,k+1 + P;,k+1 Z:,k+1 ~ Z;,k+1 _ P:,k + P;,k Z:,k ~ Z;,k = 0
2 Ar 2 Ar '
where f = 2o/n+1 + (1 _ 2a)fn, o e [1/2,3].
Phase 3 difference equations approximate the initial differential equations with the second order in the spatial variable and with the first in time. If we add these equations to the ones of phase 1, then the resulting system at o = 0.5 will have a second order of approximation and, in the absence of filtering (ou = oh = ap = 0),
have the property of temporary reversibility. By changing the parameter o, one can control the dissipative properties of the scheme. At a > 1 the scheme becomes super implicit [23], although the computational algorithm remains explicit. This algorithm at constant density is balanced: the resting state above an arbitrary bottom relief is not violated [24].
The selection of time step value. The stability condition for CABARET scheme for each layer has the form [22]
(^ )n+1/2 +(lM)i
+1/2
AXi+1/2
t
-< 1,
whence
t = CF; x min
/ \n+1/2 /1 i\n+1/2 (Ck )i+1/2 +(luk|)i+1/2
Rearrangement of layer boundaries. Filtering of flux variables without affecting conservative values does not make the complete task correct for sufficiently long periods of time. We have to correct conservative variables, which leads to the mass and momentum exchange between the layers.
The grid reconstruction in the proposed algorithm is carried out twice in one time cycle. The first time is after calculating the conservative variables on the intermediate layer in phase 1, the second time - after phase 3.
In phase 1 intermediate conservative velocities uin+1//2:2k and layer thicknesses
K+mk are calculated. Along them new heights of the interfaces
ryn+1/2 _ ryn i^n+1/2 ryn+1/2 _ ryn _ T) ryn+1/2 _ Tjn+1/2 i _ 1 \T
Zi+1/2,k ~ Zi+1/2,k+1 + hi+1/2,k , Zi+1/2,^+1 = Zi+1/2,^+1 ~ Bi+1/2, Zi+1/2,1 H i+1/2 , k = 1,...,N .
are found. In general case these coordinates are not final and are subject to
correction.
538 PHYSICAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019)
1
We consider three ways to set Z^/A that does not exhaust all the possibilities.
1. The heights are not rearranged and Z^1/^ = Z^^ . In this case, there is no
mass and momentum exchange between the layers and the problem remains incorrect.
2. New heights of the interfaces are set so that the ratio of the layer thicknesses in each vertical column remains the same (sigma-coordinate):
Zn+1/2 = Zn + L"+ 1/2 L"+1/2 = ZJn+1/2 _ B \ In
7i+1/2,k = 7 i+1/2,k+1 + 0khi+1/2 , hi+1/2N i+1/2 Bi+1/2 j/ N ,
k = 1,...,N, ok >0, £Ok=N
k=1
3. The lower boundary of k = s, N > s > 1 layer remains motionless (Euler), all the boundaries of the layers at k > 5 are also motionless and the thicknesses of the layers at k < 5 have the specified proportions (sigma-coordinate):
Z n+1/2 = Z0 k > S 7 n+1/2 - 7 n+1/2 + o hn+1/2 k - 1 S
i+1/2,k+1 i+1/2,k+1' ' i+1/2,k i+1/2,k+1 Uk"i+1/2,k' _
-(h^2 _ z°1/2,k+1 )/s, k -1,...,5, Ok > o, £ Ok - 5.
k-1
The conservation laws lead to the mass and momentum exchange between the layers in the second and third cases. The fluxes between the layers can be approximated with both the first and second order of accuracy.
The first order accuracy fluxs are calculated by the donor cell method, which
is as follows. We consider the lowest layer (k - N). If h"^1/^ < Ld/22N, then the lower layer thickness decreases by the value AL"4//^ - L"4//^ _ > 0 and
a part of its mass Ami+l/2 - p"++11//2NAhi"+11//2:iNAr,.+1/2 and momentum A(mu),+1/2 - p"++l1//22,NuГ++l1/2:/NALmHnAxt+1/2 enters the cell of the layer lying above. As a result, we obtain:
pn+1/2 - pn+1/2 un+1/2 - u n+1/2 7n+1/2 - b . Ln+1/2 pi+1/2,N - pi+1/2,N , ui+1/2,N - ui+1/2,N , Zi+1/2,N - Bi+1/2 + hi+1/2,N ,
pn+1/2 Ln+1/2 . pn+1/2 A Ln+1/2 " n+1/2 - pi+1/2,N_1hi+1/2,N_1 + pi+1/2,N ALi+1/2,N pi+1/2, N _1 - Ln+1/2 . ALn+1/2 ,
hi+1/2, N _1 + Ahi+1/2, N
pn+1/2 u n+1/2 Ln+1/2 i pn+1/2 u n+1/2 A Ln+1/2
"n+1/2 pi+1/2, N _1ui+1/2, N _1hi+1/2, N _1 + p,+1/2, Nui+1/2, N Ahi+1/2, N
u -_-_-_-_-_-_-_
"i+1/2,N_1 pn+1/2 Ln+1/2 i pn+1/2 A Ln+1/2
pi+1/2, N _1hi+1/2, N _1 + pi+1/2, N Ahi+1/2, N
In the opposite case whenh"+y2N >h"+y2N, the upper cell gives up part of its mass and momentum to the lower one:
pn+1/2 _ pn+l/2 "n+1/2 _ . n+1/2 zb+1/2 _ n , ;_n+1/2
P/+1/2,iV-1 _ Pi+1/2,N-1, "i+1/2,N-1 _ "i+1/2,N-1, Zi+1/2,N _ Bi+
-¡+1/2' "i+1/2,N '
pp b+1/2 _ Pi+1/2,N _
pn+1/2 hb+1/2 _ pn+1/2 A^B+1/2 Pi+1/2, Nhi+1/2, N pi+1/2, N+1/2, N
hn+1/2 - Ahn+1/2 i+1/2, N Z-S";+1/2,N
pn+1/2 "B+1/2 hn+1/2 - pn+1/2 "B+1/2 Ahn+1/2
"n+1/2 Pi+1/2, N i+1/2, Nhi+1/2,N Pi+1/2, N-l"i+1/2, N-l^i+1/2, N
Pn+1/2 hn+1/2 - Pn+1/2 Ahn+1/2 Pi+1/2, Nhi+1/2,N Pi+1/2, N-^Zihi+1/2, N
It should be pointed out that in the method of donor cells the densities and velocities are considered piecewise constant functions along the vertical.
The second-order accuracy fluxes are found from the condition according to which the density and horizontal velocity are linear functions on the segment connecting the midpoints of neighboring layers. We consider the cells of two adjacent layers (with the indices B is the lower layer, T is the upper layer). By Az _ hB - hB required we denote the variation of the boundary between the layers, by Ax - the width of these cells, by Am, Ap - the mass and momentum variation of the layers. Then we get the following expressions:
Am _ pBAzAx + — (hB + Az-Pt—PB— AzAx, FB 2 \hT + hB)/2
Ap _ pr"rAzAx + — (hB + Az)-p"T—PB-BAzAx.
F WB B 2 \hT + hB)/2
New values of densities and velocities of layers are calculated by the following formulas:
PThT Ax + Am PBhB Ax -Am
Pt , Pb ,
hT ,required Ax hB,required Ax
PThTuT Ax + Ap P BhBuB Ax -Ap
"t — ^ , "b — ^ .
P t hT, required Ax PBhB,required Ax
Recalculation of conservative variables on other layers is carried out similarly, from the lower layers to the upper ones. Then, the corrected conservative variables
and layer thicknesses marked with fB+1/2 symbol are taken as final values fB+112 and are used in the next phase 2. A similar correction of conservative variables is also carried out after phase 3.
Artificial (turbulent) viscosity. In gas dynamics problems the artificial viscosity method is used to regularize the numerical solution in the presence of shock waves [25]. A difference operator is added to the numerical algorithm, which approximates the so-called second - bulk viscosity, the coefficient at which is proportional to the spatial step of the grid and is found from the condition of the fastest attenuation of the highest frequency harmonic [26]. In the areas of rarefaction waves where the velocity divergence is positive, the artificial viscosity is zeroed. This provides the monotonicity of the numerical solution on compression waves and does not spoil rarefaction waves.
In the case of multilayer hydrodynamics, the inclusion of linear artificial viscosity is reduced to the modification in phase 1 and phase 3 of the flux grid pressure function on the vertical faces of the computational grid. For phase 1 :
(Pk+1/2 )k )k=(p+1* )k - ei (p^rkJ (
j \Uk, j+1/2 Uk, j-1/2
),
e* = I e if uk, j+1/2 uk, j-1/2 < 0, k,j = [0 if ulj+1/2 - uk j-1/2 > 0,
where cj jPk+1/pk is a local speed of sound, e ~ 1 is a dimensionless parameter.
For phase 3:
( Pk+1/2 )k+1 H P+1/2 )k+1 =( Pk+1/2 )k+1 - ek,j (pC2 (
k+1/2 / k+1/2
U, . - U
k+1/2 k, j+1/2 Uk,j-1/2
) ■
e*,j =
P if uk!,+V2 - u
.k+1/2
k, j+1/2 k, j-1/2
< 0,
n „.k+1/2 „,k+1/2 r\
0 if uk,j+1/2 - uk,j-1/2 > 0.
k+1/2
A generalization for the case of two spatial dimensions. A transition from the one-dimensional algorithm described above to the two-dimensional multilayer one is carried out similarly to that described in [27] for single-layer shallow water on a sphere. The consideration of Coriolis force in the CABARET scheme, which operates with two types of variables (conservative and flux), does not cause difficulties and does not impose additional restrictions on the time step.
Algorithm Verification
The purpose of verification is to check the algorithm for the absence of programming errors, computational stability and convergence on problems that can be compared with an analytical solution or with a numerical solution obtained by other authors. The proposed algorithm coincides with that described in [24] for one layer and constant density. In the presence of one layer and variable density the new multilayer algorithm coincides with that described in [28].
The features of multilayer models are already evident in problems with two or three layers [29-31]. In these works, the examples of calculations of model problems, on which our methodology was tested, are given. It should be noted that in the mentioned papers the results of the solution for a certain moment of time on fixed computational grids are presented, and the behavior of the numerical solution on large time intervals and on the grids of various densities is not studied.
We will illustrate the stability problems of the numerical solution of multilayer shallow water equations in the absence of mass and momentum exchange between layers on the two-layer problem of A. Kurganov from [30]. The initial data (Fig. 3) are as follows:
/ v / v fl if 2 > Ixl > 1,
ui(x,0) = 0,4,pi = 0,98,hi(x,0) = \i 025.(2 ) f ,, ''
11 -0,25sin (2nx) if x <1,
/ ^ / ^ f1 if 2 > Ixl > 1,
U2 (x,0) = -0,4, p 2 =1, h2 (x,0) = j 1 + 025.(2 ) f ,, <1
11 + 0,25sin (2nx) if x <1.
Bottom reliefB(x) = -2, x e [-2,2], g = 10. Computational grid - 801 node on x-axis.
Fig. 3. Initial conditions
The results of calculations according to the CABARET scheme are shown in Fig. 4. The calculations were carried out with the following parameters: ou = °h = °P = 2/3, o* = 3, 0 = 0, CFL = 0,3. Graphs of the layer interfaces and
their velocities at time t = 0.5 are given in Fig. 4, a, b. They practically coincide with those presented in [30]. When the calculation is continued, the thickness of the layers degenerates and at t = 0.65 the calculation ceases (Fig. 4, c, d).
b
a
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)
If we use a denser grid with the number of nodes 1601, then the loss of stability occurs much earlier. On a coarser grid (at N = 201) the instability develops more slowly (Fig. 5) but poor problem conditioning inevitably leads to an abortion of calculation.
The computational grid rearrangement and the resulting exchange of mass and momentum between the layers largely regularize the problem [15, 32-34], however, in this case the proper choice of filtrering parameters and artificial viscosity remains very important.
H
■1.6 -0.8 0 0.6 1.6 X
------ a
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)
Calculation of the barotropic flux according to the multilayer model with regard to the mass and momentum exchange between the layers We consider the problem of fluid oscillations in a basin with a flat bottom and a free boundary form perturbed at the initial moment. The problem parameters, initial and boundary conditions are defined as follows:
x e [-5, 5], B(x) = -2, u(x, t0) = 0, u(-5, t) = u(5, t) = 0, p = 1, g = 10,
H ( x, t0) = ^
f
1 + cos
0
2nx
., 5 5
if--< x <—,
2 2 '
else.
The calculation was carried out according to the CABARET scheme according to a single-layer model and at N _ 10 number of layers on the sigma-coordinate with the same layer thicknesses up to t _ 6 point in time corresponding to more than one period of surface oscillation. In Fig. 6, a the shape of the free surface at
b
the initial moment is represented, in Fig. 6, b - the free surface at t = 6 and Nx = 128 number of computational cells.
a b
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)
The calculations were carried out at au = ah = op = 2/3, a* = 2, 0 = 0,
CFL = 0,3 parameters. As the calculation time increases and the grid thickens, the algorithm does not lose stability and demonstrates the grid convergence.
Conclusion
Multilayer hydrostatic model with a free surface excluding the mass and momentum exchange between the layers are incorrect in most cases and cannot be used in practical calculations. Rearrangement of the computational grid at each time step can regularize the problem. However, in order to provide guaranteed stability it is necessary to use filtering of flux variables and artificial viscosity simulating the turbulent mixing. Further development of the proposed model will be aimed at developing rules for the application of the described methods for regularizing the solution that minimizes the numerical dissipation.
REFERENCES
1. Marchuk, G.I., Dymnikov, V.P. and Zalesny, V.B., 1987. [Mathematical Models in Geophysical Hydrodynamics and Methods of their Numerical Implementation]. Leningrad: Gidrometeoizdat, 296 p. (in Russian).
2. Zalesny, V.B., Marchuk, G.I., Agoshkov, V.I., Bagno, A.V., Gusev, A.V., Diansky, N.A., Moshonkin, S.N., Tamsalu, R. and Volodin, E.M., 2010. Numerical Simulation of Large-Scale Ocean Circulation Based on the Multicomponent Splitting Method. Russian Journal of Numerical Analysis and Mathematical Modelling, [e-journal] 25(6), pp. 581-609. https://doi.org/10.1515/RJNAMM.2010.036
3. Marchuk, G.I., Paton, B.E., Korotaev, G.K. and Zalesny, V.B., 2013. Data-Computing Technologies: A New Stage in the Development of Operational Oceanography. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 49(6), pp. 579-591. doi:10.1134/S000143381306011X
4. Dotsenko, S.F., Zalesny, V.B. and Sannikova, N.K.V., 2016. Block Approach to the Simulation of Circulation and Tides in the Black Sea. Physical Oceanography, [e-journal] (1), pp. 3-19. doi:10.22449/1573-160X-2016-1-3-19
5. Zalesny, V.B., Gusev, A.V. and Fomin, V.V., 2016. Numerical Model of Nonhydrostatic Ocean Dynamics Based on Methods of Artificial Compressibility and Multicomponent Splitting. Oceanology, [e-journal] 56(6), pp. 876-887. doi:10.1134/S0001437016050167
6. Goloviznin, V.M., Zaytsev, M.A., Karabasov, S.A. and Korotkin, I.A., 2013. [New Algorythms of Numerical Hydrodynamic for Multiprocessor Computational System]. Moscow: MSU, 467 p. (in Russian).
7. Audusse, E., 2005. A Multilayer Saint-Venant Model: Derivation and Numerical Validation. Discrete & Continuous Dynamical Systems - B, [e-journal] 5(2), pp. 189-214. https://doi.org/10.3934/dcdsb.2005.5.189
8. Audusse, E. and Bristeau, M.-O., 2007. Finite-Volume Solvers for a Multilayer Saint-Venant System. International Journal of Applied Mathematics and Computer Science, [e-journal] 17(3), pp. 311-320. https://doi.org/10.2478/v10006-007-0025-0
9. Ovsyannikov, L.V., 1979. Two-Layer "Shallow Water" Model. Journal of Applied Mechanics and Technical Physics, [e-journal] 20(2), pp. 127-135. https://doi.org/10.1007/BF00910010
10. Duchêne, V., 2013. A Note on the Well-Posedness of the One-Dimensional Multilayer Shallow Water Model. Available at: https://hal.archives-ouvertes.fr/hal-00922045/document [Accessed: 12 July 2019].
11. Chandrasekhar, S., 1961. Hydrodynamic and Hydromagnetic Stability. London: Oxford University Press, 652 p.
12. Chukharev, A.M., Runovsky, K.V. and Kulsha, O.E., 2017. Modeling of Turbulent Patches Statistical Distribution in the Stratified Ocean Layers. Physical Oceanography, [e-journal] (5), pp. 31-41. doi: 10.22449/1573- 160X-2017-5-31 -41
13. Kollman, V., ed., 1980. Prediction Methods for Turbulent Flows. Washington, D.C.: Hemisphere Pub., 468 p.
14. Audusse, E., Bristeau, M.-O., Perthame, B. and Sainte-Marie, J., 2011. A Multilayer Saint-Venant System with Mass Exchanges for Shallow Water Flows. Derivation and Numerical Validation. ESAIM: M2AN, [e-journal] 45(1), pp. 169-200. doi:10.1051/m2an/2010036
15. Audusse, E., Benkhaldoun, F., Sari, S., Seaid, M. and Tassi, P., 2014. A Fast Finite Volume Solver for Multi-Layered Shallow Water Flows with Mass Exchange. Journal of Computational Physics, [e-journal] 272, pp. 23-45. https://doi.org/10.1016/jjcp.2014.04.026
16. Hirt, C.W., Amsden, A.A. and Cook, J.L., 1974. An Arbitrary Lagrangian Eulerian Computing Method for All Flow Speeds. Journal of Computational Physics, [e-journal] 14(3), pp. 227-253. https://doi.org/10.1016/0021-9991(74)90051-5
17. Ringler, T.D. and Randall, D.A., 2002. The ZM Grid: An Alternative to the Z Grid. Monthly Weather Review, [e-journal] 130(5), pp. 1411-1422. https://doi.org/10.1175/1520-0493(2002)130<1411:TZGAAT>2.0.œ;2
18. Audusse, E., Bristeau, M.-O., Pelanti, M. and Sainte-Marie, J., 2011. Approximation of the Hydrostatic Navier-Stokes System for Density Stratified Flows by a Multilayer Model: Kinetic Interpretation and Numerical Solution. Journal of Computational Physics, [e-journal] 230(9), pp. 3453-3478. https://doi.org/10.1016/jjcp.2011.01.042
19. Stewart, A.L. and Dellar, P.J., 2010. Multilayer Shallow Water Equations with Complete Coriolis Force. Part 1. Derivation on a Non-Traditional Beta-Plane. Journal of Fluid Mechanics, [e-journal] 651, pp. 387-413. https://doi.org/10.1017/S0022112009993922
20. Bermudez, A. and Vazquez, M.E., 1994. Upwind Methods for Hyperbolic Conservation Laws with Source Terms. Computers & Fluids, [e-journal] 23(8), pp. 1049-1071. https://doi.org/10.1016/0045-7930(94)90004-3
21. Toro, E.F., 2009. Riemann Solvers and Numerical Methods for Fluid Dynamics. Berlin: Springer, 724 p. doi:10.1007/b79761
22. Karabasov, S.A. and Goloviznin, V.M., 2009. Compact Accurately Boundary-Adjusting High-Resolution Technique for Fluid Dynamics. Journal of Computational Physics, [e-journal] 228(19), pp. 7426-7451. https://doi.org/10.1016/jjcp.2009.06.037
23. Goloviznin, V.M., Kanyukova, V.D. and Samarskaya, E.A., 1983. Super-Implicit Difference Schemes of Gas Dynamics. Differential Equations, 19(7), pp. 1186-1197.
24. Goloviznin, V.M. and Isakov, V.A., 2017. Balance-Characteristic Scheme as Applied to the Shallow Water Equations over a Rough Bottom. Computational Mathematics and Mathematical Physics, [e-journal] 57(7), pp. 1140-1157. https://doi.org/10.1134/S0965542517070089
25. Richtmyer, R.D. and Morton, K.W., 1967. Difference Methods for Initial-Value Problems. New York: Interscience Publishers, 405 p.
26. Goloviznin, V.M., 1982. On a Method of Introducing Artificial Dissipation into the Variational-Difference Schemes of Magnetohydrodynamics. USSR Computational Mathematics and Mathematical Physics, [e-journal] 22(1), pp. 150-156. https://doi.org/10.1016/0041-5553(82)90172-0
27. Goloviznin, V.M., Solovjov, A.V. and Zalesny, V.B., 2018. A New Algorithm for Solving the Shallow Water Equations on the Sphere Based on the Cabaret Scheme. Journal of Physics: Conference Series, [e-journal] 1128(1), 012091. https://doi.org/10.1088/1742-6596/1128/1/012091
28. Evtushenko, Y.G., Gorchakov, A.Y. and Goloviznin, V.M., 2018. Fast Automatic Differentiation in Problems Variations Four-Dimensional Data Assimilation (4Dvar). Journal of Physics: Conference Series, [e-journal] 1128(1), 012001. https://doi.org/10.1088/1742-6596/1128/1/012001
29. Vedernikov, A.B. and Kholodov, A.S., 1990. Numerical Investigations of Flows of the Double-and Threelayer Fluid in Framework of "Shallow Water" Model. Matematicheskoe Modelirovanie, 2(6), pp. 9-18 (in Russian).
30. Kurganov, A. and Petrova, G., 2009. Central-Upwind Schemes for Two-Layer Shallow Water Equations. SIAM Journal on Scientific Computing, [e-journal] 31(3), pp. 1742-1773. https://doi.org/10.1137/080719091
31. Elizarova, T.G. and Ivanov, A.V., 2016. Quasi-Gasdynamic Algorithm for Numerical Solution of Two-Layer Shallow Water Equations. KIAM Preprint No. 69. 27 p. (in Russian). doi:10.20948/prepr-2016-69
32. Shankar, N.J., Cheong, H.F. and Sankaranarayanan, S., 1997. Multilevel Finite-Difference Model for Three-Dimensional Hydrodynamic Circulation. Ocean Engineering, [e-journal] 24(9), pp. 785-816. https://doi.org/10.1016/S0029-8018(96)00036-4
33. Bouchut, F. and Zeitlin, V., 2010. A Robust Well-Balanced Scheme for Multi-Layer Shallow Water Equations. Discrete & Continuous Dynamical Systems - B, [e-journal] 3(4), pp. 739758. http://doi.org/10.3934/dcdsb.2010.13.739
34. Couderc, F., Duran, A. and Vila, J.-P., 2017. An Explicit Asymptotic Preserving low Froude Scheme for the Multilayer Shallow Water Model with Density Stratification. Journal of Computational Physics, [e-journal] 343, pp. 235-270. https://doi.org/10.1016/jjcp.2017.04.018
About the authors:
Vasiliy M. Goloviznin - Head of Department, Nuclear Safety Institute of the Russian Academy of Sciences, (52, Bolshaya Tulskaya str., Moscow, 115191, Russian Federation), Dr. Sci. (Phys.-Math.), Professor, ResearcherID: Y-2025-2018, gol@ibrae.ac.ru
Pavel A. Maiorov - Engineer, Nuclear Safety Institute of the Russian Academy of Sciences, (52, Bolshaya Tulskaya str., Moscow, 115191, Russian Federation), ResearcherID: AAD-6169-2019, pavel.a.mayorov@gmail.com
Petr A. Maiorov - Engineer, Nuclear Safety Institute of the Russian Academy of Sciences, (52, Bolshaya Tulskaya str., Moscow, 115191, Russian Federation), ResearcherID: AAD-6173-2019, maiorov.peter@gmail.com
Andrey V. Solovjev - Leading Research Associate, Nuclear Safety Institute of the Russian Academy of Sciences, (52, Bolshaya Tulskaya str., Moscow, 115191, Russian Federation), Ph. D. (Phys.-Math.), ResearcherID: X-4858-2019, solovjev@ibrae.ac.ru
Contribution of the co-authors:
Vasiliy M. Goloviznin - general research guidance, study initiation, scientific leadership, formulation of research goals and objectives, development and scientific substantiation of the concept, preparation of the initial version of the paper text
Pavel A. Maiorov - correction of the mathematical model and calculations, analysis of the applied research methodology, processing and description of research results, computer implementation of algorithms, refinement of the text, visualization / representation of the data in the text, representation of the data in the text, their analysis
Petr A. Maiorov - selection and justification of numerical methods for solving equations, development of methodology, analysis and generalization of research results, collection of available materials on the research topic, development and debugging of a computer program for solving a problem, building tables, graphs, diagrams, computer work
Andrey V. Solovjev - study of the concept, development of a mathematical model, qualitative and quantitative analysis of the results, analysis of materials on the research topic, analysis of materials in domestic and foreign sources according to research methods, analysis of materials on the research topic, literature review on the research problem, literature review on foreign sources, the search for analytical materials in domestic and foreign sources, analysis of literature data, critical analysis, writing the annotations, critical analysis and revision of the text, editing and addition of the paper text, the study of the concept
All the authors have read and approved the final manuscript.
The authors declare that they have no conflict of interest.
