Novel entropy stable schemes for 1D and 2D fluid equations
Novel entropy stable schemes for 1D and 2D fluid equations Eitan Tadmor1 and Weigang Zhong2 1
2
Department of Mathematics, Center for Scientific Computation and Mathematical Modeling (CSCAMM), Institute of Physical Science and Technology, University of Maryland, College Park, MD 20742 [email protected] Department of Mathematics, Center for Scientific Computation and Mathematical Modeling (CSCAMM), University of Maryland, College Park, MD 20742 [email protected]
Summary. We present a systematic study of the novel entropy stable approximations of a variety of nonlinear conservation laws, from the scalar Burgers equation to 1D Navier-Stoke and 2D shallow water equations. This new family of second-order difference schemes avoid using artificial numerical viscosity in the sense that their entropy dissipation is dictated solely by physical dissipation terms. The numerical results of 1D compressible Navier-Stokes equations equations provide us a remarkable evidence for different roles of viscosity and heat conduction in forming sharp monotone profiles in the immediate neighborhoods of shocks and contacts. Further implementation in 2D shallow water equations is realized dimension by dimension.
1 The Inviscid Burgers’ Equation 1.1 Entropy conservative schemes We begin with the inviscid Burgers equation, ∂u ∂ + f(u) = 0, ∂t ∂x
f(u) =
1 2 u . 2
(1)
It is equipped with a family of entropy functions, Up (u) = u2p, p = 1, 2, · · ·, such that solutions of (1) satisfy, at the formal level, ∂ ∂ Up (u) + Fp (u) = 0. ∂t ∂x
(2)
These are additional conservation laws balanced by the corresponding entropy flux functions Fp (u) = 2pu2p+1/(2p + 1). Spatial integration then yields (ignoring boundary contributions)
2
Eitan Tadmor and Weigang Zhong
Z
u2p (x, t) dx = x
Z
u2p(x, 0) dx.
(3)
x
We turn to the discrete framework. Discretization in space yields the semidiscrete scheme, d 1 (4) uν (t) + fν+ 12 − fν− 12 = 0. dt ∆x Here, uν (t) denotes the discrete solution along the gridline (xν , t) with xν := ν∆x, ∆x being the uniform meshsize, and fν+ 1 := f(uν−r+1 , · · · , uν+r ) is a 2 consistent numerical flux based on a stencil of 2rP+ 1 neighboringPgrid values, that makes (4) conservative in the sense that ν uν (t)∆x = ν uν (0)∆x. Fix p. We seek a semi-discrete scheme that conserves the entropy Up (u) = u2p in the sense of satisfying the discrete analogue of (2)p , d 1 Up (uν (t)) + (F 1 − Fν− 1 ) = 0. 2 dt ∆x ν+ 2 Here, Fν+ 12 is a consistent numerical entropy flux. According to [Ta1987, Theorem 5.2], such 3-point scalar entropy conservative schemes are uniquely ∗ determined by the entropy conservative numerical flux fν+ 12 = fν+ 1 given by, 2
2p+1
∗ fν+ 12 = fν+ 1 = 2
−1 2p − 1 (uν+1/uν ) · u2ν · . 2p−1 2(2p + 1) (uν+1 /uν ) −1
(5)
The resulting scheme (4), (5) is entropy conservative in the sense that the discrete analogue of total entropy conservation (3) is satisfied, X X u2p u2p ν (t) ∆x = ν (0) ∆x ν
ν
Of course, all the above manipulations are at the formal level. To recover the physical relevant entropy inequality, ∂t Up (u) + ∂x Fp (u) ≤ 0, one can add numerical dissipation, d 1 ∗ ∗ = d(u ) − 2d(u ) + d(u ) , > 0. uν (t) + fν+ 1 − fν− 1 ν+1 ν ν−1 2 2 dt ∆x (∆x)2 This serves as an approximation to the vanishing viscosity regularization ut + f(u)x = d(u)xx, d0(u) > 0. Sum this scheme against vν := Up0 (uν ) = 2pu2p−1: the resulting entropy balance that follows reads, d X X ∆dν+ 12 Up (uν (t))∆x = − (∆vν+ 12 )2 ≤ 0, dt ν ∆x ν ∆vν+ 1
(6)
2
since
∆dν+ 1
2
∆vν+ 1
2
:=
d(uν+1 )−d(uν ) vν+1 −vν
> 0 for d0(u) > 0. Observe that the amount of
entropy dissipation on the right is completely determined by the dissipation
Novel entropy stable schemes for 1D and 2D fluid equations
3
term d(u). No artificial viscosity is introduced by the convective term. If we exclude any dissipative mechanism ( = 0), the entropy conservative solutions admit dispersive oscillations interesting for their own sake, consult [La1986, LL1996].
Scalar Burgers equation,entropy−consertive scheme w/ U(u)=u2,∆ t=0.005,∆ x=0.005
Scalar Burgers equation,entropy−consertive scheme w/ U(u)=u8,∆ t=0.005,∆ x=0.005
2
1.5 t=0.0 t=0.125 t=0.25
1.5
t=0.0 t=0.125 t=0.25 1
1 0.5
u(x)
u(x)
0.5
0
0
−0.5 −0.5 −1 −1 −1.5
−2
0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
−1.5
1
0
0.1
0.2
0.3
(a) p=1
0.4
0.5 x
0.6
0.7
0.8
Scalar Burgers equation,entropy−consertive scheme w/ U(u)=u32,∆ t=0.005,∆ x=0.005
Scalar Burgers equation,entropy−consertive scheme w/ U(u)=u64,∆ t=0.005,∆ x=0.005 1.5
t=0.0 t=0.125 t=0.25
t=0.0 t=0.125 t=0.25
1
1
0.5
0.5
u(x)
u(x)
1
(b) p=4
1.5
0
0
−0.5
−0.5
−1
−1
−1.5
0.9
0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
(c) p=16
0.9
1
−1.5
0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
(d) p=32
Fig. 1. Scalar Burger’s equation, sine initial condition, entropy-conservative schemes, 200 spatial grids, U (u) = u2p , ∆t = 5 × 10−3 , ∆x = 5 × 10−3
1.2 Numerical experiments The semi-discrete entropy conservative scheme (4),(5) is integrated with the following third-order Runge-Kutta (RK3), consult [GST2001].
1
4
Eitan Tadmor and Weigang Zhong
u(1) = un + ∆tL(un ), 3 1 u(2) = un + u(1) + 4 4 1 n 2 (2) n+1 u = u + u + 3 3
1 ∆tL(u(1) ), 4 2 ∆tL(u(2) ), 3
(7)
1 ∗ ∗ where [L(u)]ν = − ∆x (fν+ ). We note that this explicit RK3 time dis1 −f ν− 12 2 cretization produces a negligible amount of entropy dissipation. For a general framework of entropy conservative fully discrete schemes, consult [LMR2002]. We solve the inviscid Burgers equation with the sine initial condition, u(0, x) = sin(2πx) and periodic boundary. In Fig.1, we display the numerical solutions for (7) with the numerical flux (5) for different choices of p. Observe that the amplitude of the spurious dispersive oscillations decreases as p increases. Indeed, as we increase p, the control of a constant Up entropy, 1 P 2p 2p approaches the control of L∞ -norm. ν uν (t) ∆x
2 The 1D Navier-Stokes Equations 2.1 Entropy dissipation We consider the Navier-Stokes equations (NSE) governing the density ρ = ρ(x, t), momentum m = m(x, t), and energy E = E(x, t), ∂ ∂ ∂2 u+ f (u) = 2 d(u), ∂t ∂x ∂x
u = [ρ, m, E]>.
(8)
> They are driven by the convective flux f (u) = m, qm + p, q(E + p) , to > > gether with the dissipative flux d(u) = (λ + 2µ) 0, q, q2/2 + κ 0, 0, θ which stands for the combined viscous and heat fluxes. Here, represents the amplitudes of the viscosity and conductivity. These fluxes involve the velocity 2 q := m/ρ, the pressure p = p(x, t) = (γ − 1)e where e := E − m 2ρ , and the absolute temperature, θ = θ(x, t) > 0, such that Cv ρθ = e. Here γ > 1, λ, µ are fixed and κ 7→ κ/Cv with Cv = 1. The viscous and heat fluxes are dissipative in the sense that they are responsible for the dissipation of total entropy, U (u) = −ρS with S = ln(pρ−γ ) being the specific entropy, 2 ∂ qx2 ∂ θx ≤ 0. (9) (−ρS) + (−mS + κ(ln θ)x ) = −(λ + 2µ) − κ ∂t ∂x θ θ Spatial integration of (9) then yields the second law of thermodynamics, Z 2 Z Z 2 d θx qx (−ρS) dx = −(λ + 2µ) dx ≤ 0. (10) dx − κ dt x θ θ x x
Novel entropy stable schemes for 1D and 2D fluid equations
5
In fact, (10) specifies the precise entropy decay rate. In the case of the Euler equations, λ = µ = κ = 0, total entropy is precisely conserved, Z Z −ρS(x, t) dx = −ρS(x, 0) dx. x
x
This corresponds to the scalar entropy conservation (3). 2.2 Entropy stable schemes for the Navier-Stokes equations We now turn our attention to the construction of entropy stable schemes for NSE (8). We use the conservative differences for convective flux and standard centered differences for the dissipative terms on the RHS. d 1 ∗ ∗ = d , (11) − 2d + d uν (t) + fν+ 1 − fν− 1 ν+1 ν ν−1 2 2 dt ∆x (∆x)2 Here, dν := d(uν ). As in the scalar case, we seek entropy conservative flux ∗ fν+ 1 , so that the entropy decay will be dictated solely by the dissipation terms 2 ∗ on the right of (11). The construction of the entropy conservative flux fν+ 1 2 follows [Ta2004] and [TZ2006], and is summarized in the following Algorithm 1 If
∗ fν+ 1 = f (vν ); else
uν = uν+1 then
2
• Set ∆uν+ 12 := uν+1 − uν . Starting with u1ν+ 1 := uν , compute recussively 2 the intermediate states, E D j j b b uj+1 rjν+ 1 , j = 1, 2, 3, . ` 1 , ∆u 1 (12) 1 = u 1 + ν+ ν+ ν+ ν+ 2
2
2
2
2
j {b `ν+ 12 }
Here, and {b rjν+ 1 } are the left and right eigensystems of the Roe 2 matrix A(uν , uν+1)(see [Roe1981]). • Set rjν+ 1 := v(uj+1 ) − v(ujν+ 1 ) and let {`j }3j=1 be the corresponding ν+ 12 2 2 orthogonal system. Compute the entropy-conservative numerical flux, ∗ fν+ 1 2
3 mj+1 − mjν+ 1 X ν+ 12 2 D E `jν+ 1 , = (γ − 1) j 2 j=1 `ν+ 1 , ∆vν+ 12
D
j+1 j `jν+ 1 , vν+ 1 − v ν+ 1 2
2
E
= δjk
2
2
(13) Here, v(u) := Uu (u) = [−E/e − S + γ + 1, q/θ, −1/θ]> are the entropy variables, and {mj } are intermediate momentum values along the path. We now arrive at our main result of NSE corresponding to the Burgers statement (6). Theorem 1 ([TZ2006], Theorem 3.6). The semi-discrete scheme (11) ∗ with the entropy conservative numerical flux fν+ 1 in (12)-(13) and d(uν ) 2
being the dissipative NS flux, is entropy stable in the sense that 3
√ We let zbν+ 1 = zν + zν+1 /2 and zeν+ 1 = zν zν+1 . 2
2
3
6
Eitan Tadmor and Weigang Zhong
X d X [−ρν (t)Sν (t)] ∆x = − dt ν ∆x ν
*
= − (λ + 2µ)
∆vν+ 1 , 2
∆vν+ 12 2
X ∆qν+ 1 2
∆x
ν
−κ
∆dν+ 12
X ∆θν+ 1 2 2
ν
∆x
g 1/θ
∆vν+ 1
+
2
d 1/θ
2 ν+ 12
ν+ 12
(14)
∆x
∆x ≤ 0.
This statement is a discrete analog of the entropy balance statement (10). Here is our main point: we introduce no excessive entropy dissipation due to spurious, artificial numerical viscosity. According to (14), the semi-discrete scheme contains the precise amount of numerical viscosity to enforce the correct entropy dissipation dictated by NSE. More can be found in [Ta2004, TZ2006]. 2.3 Numerical experiments We consider ideal polytropic gas equations as an approximation of air with γ = 1.4,
Cv = 716,
κ = 0.03,
λ + 2µ = 2.28 × 10−5 .
We simulate the Sod’s shocktube problem, where the Euler and NSE are solved over the interval [0, 1] subject to Riemann initial conditions (1.0, 0.0, 2.5) 0 < x ≤ 0.5 (ρ, m, E)t=0 = (0.125, 0.0, 0.25) 0.5 < x < 1. In Fig.2, we display the numerical solutions for the fully discrete scheme (11) with RK3 method (7) and the numerical flux (13). The density fields of four different cases are recorded. Density field of the Euler equations 2(a) demonstrates the purely dispersive character of the entropy conservative schemes. Dispersive oscillations on the mesh scale are observed in shocks and contact regions due to the absence of any dissipative mechanism. These oscillations approach a modulated wave envelop, consult [La1986, LL1996] for more discussions on dispersive oscillations. For the results of NSE in 2(b)-2(d), the presence of heat flux causes the oscillations to be dramatically reduced around the contact discontinuity and the shock in 2(b). The viscous flux in NSE, on the other hand, is doing a better job than heat flux in damping oscillations around the shock in 2(c), but we still can observe an oscillatory behavior around the contact discontinuity. In 2(d), not only the oscillations around the shock are damped out by viscosity, but the oscillations around the contact discontinuity are significantly reduced due to the heat flux.
Novel entropy stable schemes for 1D and 2D fluid equations
1−D Euler,density,γ=1.4,Cv=716,κ=0,entropy=−ρ ln(pρ−γ),∆t=2.5e−005,∆x=0.00025,Tmax=0.1 t=0.0 t=0.05 t=0.1
0.9
1−D Euler,density,γ=1.4,Cv=716,κ=0.03,entropy=−ρ ln(pρ−γ),∆t=2.5e−005,∆x=0.00025,Tmax=0.1
1
1
7
t=0.0 t=0.05 t=0.1
0.9
0.8
0.7
0.7
0.6
0.6
ρ
ρ
0.8
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
1
0
(a) Euler
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
1
(b) Navier-Stokes with heat only
1−D N−S,density,λ+2µ=2.28e−005,γ=1.4,Cv=716,κ=0,entropy=−ρ ln(pρ−γ),∆t=2.5e−005,∆x=0.00025
1
0.1
1−D N−S,density,λ+2µ=2.28e−005,γ=1.4,Cv=716,κ=0.03,entropy=−ρ ln(pρ−γ),∆t=2.5e−005,∆x=0.00025
1 t=0.0 t=0.05 t=0.1
0.9
t=0.0 t=0.05 t=0.1
0.9
0.8
0.7
0.7
0.6
0.6 ρ
ρ
0.8
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
(c) Navier-Stokes with viscosity only
1
0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
(d) Navier-Stokes with viscosity and heat
Fig. 2. Density field Sod problem with 4000 spatial gridpoints, U (u) = −ρ ln pρ−γ , ∆t = 2.5 × 10−5 , ∆x = 2.5 × 10−4
3 The 2D Shallow Water Equations We turn to 2D shallow water equations, ∂ ∂ ∂ ∂ ∂ ∂ ∂ u+ f (u) + g(u) = η h d(u) + η h d(u) , ∂t ∂x ∂y ∂x ∂x ∂y ∂y
(15)
with u = [h, uh, vh]> being the vector of conserved variables balanced by the flux vectors f = [uh, u2h + gh2 /2, uvh]> ,
g = [vh, uvh, v2 h + gh2 /2]>,
1
8
Eitan Tadmor and Weigang Zhong
and the viscous flux vector d = [0, u, v]> . Here, h = h(x, t) is the total water height, (u(x, t), v(x, t)) are the depth-averaged velocities along x and y direction. Finally, g is the constant acceleration due to gravity, and η > 0 is the constant eddy viscosity which models the turbulence stress in the flow. The total energy U (u) = (gh2 +u2 h+v2 h)/2 serves as an entropy function, ∂ ∂ ∂ U (u) + F (u) + G(u) = −ηh(u2x + u2y + vx2 + vy2 ), ∂t ∂x ∂y 3
2
2
(16) 3
h where F (u) = guh2 + u h+uv − huux and G(u) = gvh2 + u vh+v 2 2 are the entropy fluxes. Spatial integration of (16) yields Z Z Z Z d U (u) dxdy = −η h(u2x + u2y + vx2 + vy2 ) dxdy. dt y x y x
h
− hvvy
(17)
For the inviscid case (η = 0), the global entropy conservation is satisfied, Z Z Z Z U (u(x, t)) dxdy = U (u(x, 0)) dxdy y
x
y
x
Arguing along the same line as the above NSE dimension by dimension, we obtain the entropy-stable semi-discrete schemes (recall z[ ν+ 1 := (zν+1 +zν )/2) 2
d 1 ∗ 1 ∗ ∗ uν, µ (t) + (fν+ 1 , µ − fν− (g∗ 1 1 − g ν, µ− 12 ) , µ) + 2 2 dt ∆x ∆y ν, µ+ 2 η \ dν+1, µ − dν, µ dν, µ − dν−1, µ = (h 1 − h\ ) ν− 12 , µ ∆x ν+ 2 , µ ∆x ∆x η \ dν, µ+1 − dν, µ dν, µ − dν, µ−1 1 + (h − h\ ), ν, µ− 12 ∆y ν, µ+ 2 ∆x ∆x
(18a)
∗ ∗ with the entropy-conservative fluxes fν+ and gν, outlined in Algorithm 1 µ+ 12 2,µ 1 along x and y direction, respectively,
∗ fν+ 1 2,µ
j+1 j 2 j+1 2 j 3 g X (hν+ 12 , µ ) uν+ 12 , µ − (hν+ 12 , µ ) uν+ 12 , µ xj D j E `ν+ 1 , µ , = 2 2 j=1 `x 1 , ∆v 1 ν+ 2 , µ
∗ gν, = µ+ 1 2
g 2
3 X j=1
(hj+1 )2uj+1 ν, µ+ 12 ν, µ+ 12 D
(18b)
ν+ 2 , µ
− (hjν, µ+ 1 )2 ujν, µ+ 1 j 2 2 E `yν, µ+ 1 , yj 2 `ν, µ+ 1 , ∆vν, µ+ 12
(18c)
2
Here, uν, µ (t) denotes the discrete solution at the grid point (xν , yν , t), dν, µ := d(uν, µ ), and v := Uu = [gh − 12 (u2 + v2 ), u, v]> is the entropy variable. Numerical flux f ∗ and g∗ are constructed separately along two different phase j j pathes dictated by two sets of vectors {`x } and {`y }. {hj } and {uj } are intermediate values of height and velocity along the path. The above difference scheme (18a)-(18c) is an entropy stable scheme with no artificial viscosity in the sense that the following discrete entropy balance is satisfied,
Novel entropy stable schemes for 1D and 2D fluid equations
(
"
9
2 #
2
X ∆uν+ 1 , µ ∆vν+ 1 , µ d X 2 2 h\ U (uν, µ (t))∆x∆y = −η + ν+ 12 , µ dt ν, µ ∆x ∆x ν, µ " #) ∆uν, µ+ 12 2 ∆vν, µ+ 12 2 ∆x∆y. (19) +h\ + ν, µ+ 12 ∆y ∆y
(19) is a discrete analogue of the entropy balance statement (17). Equipped by RK3, we test the entropy stable scheme (18a)-(18c) by the 2D partial Dam-Break problem with free-slip boundary described in [FC1990]. Both the inviscid and viscous case are tested. The water surface profiles at time t = 7.2s are recorded in Fig.2.
Height at t=7.2 10 9 8 7 6 5 4 3 11 10 9 8 7 6 5 4 3 2
0
5
10
40 35 30 25 20 15 10 15
20
25
30
(a) η = 0
5 35
40 0
11 10 9 8 7 6 5 4 3 2
Height at t=7.2 10 9 8 7 6 5 4
10 9 8
10
7
9
6
8
5
7 6
4 3
5 4
40 35 30 25 20 15 10
3
0
5
10
15
20
25
30
5 35
40 0
(b) η = 10m2 · s−1
Fig. 3. Shallow water equations, Dam-Break, 200×200 m2 basin, free-slip boundary, ∆x = ∆y = 5 m, ∆t = 5 × 10−3 s
Comparing 3(b) to 3(a), we observe the improvements in smoothness of the numerical solutions. There is no analytical reference solution for this test case, but other numerical results are available in [FC1990].
References [FC1990] Fennema, R.J., Chaudhry, M.H.: Explicit methods for 2D transient freesurface flows. J. Hydraul. Eng. ASCE, 116, 1013–1034 (1990) [GST2001] Gottlieb, S., Shu, C.-W., and Tadmor, E.: Strong stability-preserving high-order time discretization methods. SIAM Review, 43, 1, 89–112 (2001) [La1986] Lax, P.D.: On dispersive difference schemes. Phys. D, 18, 250–254 (1986) [LMR2002] LeFloch, P.G., Mercier, J.M., and Rohde, C.: Full discrete, entropy conservative schemes of arbitrary order. SIAM J. Numer. Anal., 40, 5, 1968– 1992 (2002)
10
Eitan Tadmor and Weigang Zhong
[LL1996] Levermore, C.D., Liu, J.-G.: Oscillations arising in numerical experiments. Phys. D, 99, 191–216 (1996) [Roe1981] Roe, P.L.: Approximate Riemann solvers, parameter vectors and difference schemes. J. Comp. Phys., 43, 357–372 (1981) [Ta1987] Tadmor, E.: The numerical viscosity of entropy stable schemes for systems of conservation laws, I. Math. Comp., 49, 91–103 (1987) [Ta2004] Tadmor, E.: Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. Acta Numerica, 451–512 (2004) [TZ2006] Tadmor, E. and Zhong, W.: Entropy stable approximations of NavierStokes equations with no artificial numerical viscosity. J. of Hype. Diff. Equa., 3, 3, 529–559 (2006)
2
Department of Mathematics, Center for Scientific Computation and Mathematical Modeling (CSCAMM), Institute of Physical Science and Technology, University of Maryland, College Park, MD 20742 [email protected] Department of Mathematics, Center for Scientific Computation and Mathematical Modeling (CSCAMM), University of Maryland, College Park, MD 20742 [email protected]
Summary. We present a systematic study of the novel entropy stable approximations of a variety of nonlinear conservation laws, from the scalar Burgers equation to 1D Navier-Stoke and 2D shallow water equations. This new family of second-order difference schemes avoid using artificial numerical viscosity in the sense that their entropy dissipation is dictated solely by physical dissipation terms. The numerical results of 1D compressible Navier-Stokes equations equations provide us a remarkable evidence for different roles of viscosity and heat conduction in forming sharp monotone profiles in the immediate neighborhoods of shocks and contacts. Further implementation in 2D shallow water equations is realized dimension by dimension.
1 The Inviscid Burgers’ Equation 1.1 Entropy conservative schemes We begin with the inviscid Burgers equation, ∂u ∂ + f(u) = 0, ∂t ∂x
f(u) =
1 2 u . 2
(1)
It is equipped with a family of entropy functions, Up (u) = u2p, p = 1, 2, · · ·, such that solutions of (1) satisfy, at the formal level, ∂ ∂ Up (u) + Fp (u) = 0. ∂t ∂x
(2)
These are additional conservation laws balanced by the corresponding entropy flux functions Fp (u) = 2pu2p+1/(2p + 1). Spatial integration then yields (ignoring boundary contributions)
2
Eitan Tadmor and Weigang Zhong
Z
u2p (x, t) dx = x
Z
u2p(x, 0) dx.
(3)
x
We turn to the discrete framework. Discretization in space yields the semidiscrete scheme, d 1 (4) uν (t) + fν+ 12 − fν− 12 = 0. dt ∆x Here, uν (t) denotes the discrete solution along the gridline (xν , t) with xν := ν∆x, ∆x being the uniform meshsize, and fν+ 1 := f(uν−r+1 , · · · , uν+r ) is a 2 consistent numerical flux based on a stencil of 2rP+ 1 neighboringPgrid values, that makes (4) conservative in the sense that ν uν (t)∆x = ν uν (0)∆x. Fix p. We seek a semi-discrete scheme that conserves the entropy Up (u) = u2p in the sense of satisfying the discrete analogue of (2)p , d 1 Up (uν (t)) + (F 1 − Fν− 1 ) = 0. 2 dt ∆x ν+ 2 Here, Fν+ 12 is a consistent numerical entropy flux. According to [Ta1987, Theorem 5.2], such 3-point scalar entropy conservative schemes are uniquely ∗ determined by the entropy conservative numerical flux fν+ 12 = fν+ 1 given by, 2
2p+1
∗ fν+ 12 = fν+ 1 = 2
−1 2p − 1 (uν+1/uν ) · u2ν · . 2p−1 2(2p + 1) (uν+1 /uν ) −1
(5)
The resulting scheme (4), (5) is entropy conservative in the sense that the discrete analogue of total entropy conservation (3) is satisfied, X X u2p u2p ν (t) ∆x = ν (0) ∆x ν
ν
Of course, all the above manipulations are at the formal level. To recover the physical relevant entropy inequality, ∂t Up (u) + ∂x Fp (u) ≤ 0, one can add numerical dissipation, d 1 ∗ ∗ = d(u ) − 2d(u ) + d(u ) , > 0. uν (t) + fν+ 1 − fν− 1 ν+1 ν ν−1 2 2 dt ∆x (∆x)2 This serves as an approximation to the vanishing viscosity regularization ut + f(u)x = d(u)xx, d0(u) > 0. Sum this scheme against vν := Up0 (uν ) = 2pu2p−1: the resulting entropy balance that follows reads, d X X ∆dν+ 12 Up (uν (t))∆x = − (∆vν+ 12 )2 ≤ 0, dt ν ∆x ν ∆vν+ 1
(6)
2
since
∆dν+ 1
2
∆vν+ 1
2
:=
d(uν+1 )−d(uν ) vν+1 −vν
> 0 for d0(u) > 0. Observe that the amount of
entropy dissipation on the right is completely determined by the dissipation
Novel entropy stable schemes for 1D and 2D fluid equations
3
term d(u). No artificial viscosity is introduced by the convective term. If we exclude any dissipative mechanism ( = 0), the entropy conservative solutions admit dispersive oscillations interesting for their own sake, consult [La1986, LL1996].
Scalar Burgers equation,entropy−consertive scheme w/ U(u)=u2,∆ t=0.005,∆ x=0.005
Scalar Burgers equation,entropy−consertive scheme w/ U(u)=u8,∆ t=0.005,∆ x=0.005
2
1.5 t=0.0 t=0.125 t=0.25
1.5
t=0.0 t=0.125 t=0.25 1
1 0.5
u(x)
u(x)
0.5
0
0
−0.5 −0.5 −1 −1 −1.5
−2
0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
−1.5
1
0
0.1
0.2
0.3
(a) p=1
0.4
0.5 x
0.6
0.7
0.8
Scalar Burgers equation,entropy−consertive scheme w/ U(u)=u32,∆ t=0.005,∆ x=0.005
Scalar Burgers equation,entropy−consertive scheme w/ U(u)=u64,∆ t=0.005,∆ x=0.005 1.5
t=0.0 t=0.125 t=0.25
t=0.0 t=0.125 t=0.25
1
1
0.5
0.5
u(x)
u(x)
1
(b) p=4
1.5
0
0
−0.5
−0.5
−1
−1
−1.5
0.9
0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
(c) p=16
0.9
1
−1.5
0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
(d) p=32
Fig. 1. Scalar Burger’s equation, sine initial condition, entropy-conservative schemes, 200 spatial grids, U (u) = u2p , ∆t = 5 × 10−3 , ∆x = 5 × 10−3
1.2 Numerical experiments The semi-discrete entropy conservative scheme (4),(5) is integrated with the following third-order Runge-Kutta (RK3), consult [GST2001].
1
4
Eitan Tadmor and Weigang Zhong
u(1) = un + ∆tL(un ), 3 1 u(2) = un + u(1) + 4 4 1 n 2 (2) n+1 u = u + u + 3 3
1 ∆tL(u(1) ), 4 2 ∆tL(u(2) ), 3
(7)
1 ∗ ∗ where [L(u)]ν = − ∆x (fν+ ). We note that this explicit RK3 time dis1 −f ν− 12 2 cretization produces a negligible amount of entropy dissipation. For a general framework of entropy conservative fully discrete schemes, consult [LMR2002]. We solve the inviscid Burgers equation with the sine initial condition, u(0, x) = sin(2πx) and periodic boundary. In Fig.1, we display the numerical solutions for (7) with the numerical flux (5) for different choices of p. Observe that the amplitude of the spurious dispersive oscillations decreases as p increases. Indeed, as we increase p, the control of a constant Up entropy, 1 P 2p 2p approaches the control of L∞ -norm. ν uν (t) ∆x
2 The 1D Navier-Stokes Equations 2.1 Entropy dissipation We consider the Navier-Stokes equations (NSE) governing the density ρ = ρ(x, t), momentum m = m(x, t), and energy E = E(x, t), ∂ ∂ ∂2 u+ f (u) = 2 d(u), ∂t ∂x ∂x
u = [ρ, m, E]>.
(8)
> They are driven by the convective flux f (u) = m, qm + p, q(E + p) , to > > gether with the dissipative flux d(u) = (λ + 2µ) 0, q, q2/2 + κ 0, 0, θ which stands for the combined viscous and heat fluxes. Here, represents the amplitudes of the viscosity and conductivity. These fluxes involve the velocity 2 q := m/ρ, the pressure p = p(x, t) = (γ − 1)e where e := E − m 2ρ , and the absolute temperature, θ = θ(x, t) > 0, such that Cv ρθ = e. Here γ > 1, λ, µ are fixed and κ 7→ κ/Cv with Cv = 1. The viscous and heat fluxes are dissipative in the sense that they are responsible for the dissipation of total entropy, U (u) = −ρS with S = ln(pρ−γ ) being the specific entropy, 2 ∂ qx2 ∂ θx ≤ 0. (9) (−ρS) + (−mS + κ(ln θ)x ) = −(λ + 2µ) − κ ∂t ∂x θ θ Spatial integration of (9) then yields the second law of thermodynamics, Z 2 Z Z 2 d θx qx (−ρS) dx = −(λ + 2µ) dx ≤ 0. (10) dx − κ dt x θ θ x x
Novel entropy stable schemes for 1D and 2D fluid equations
5
In fact, (10) specifies the precise entropy decay rate. In the case of the Euler equations, λ = µ = κ = 0, total entropy is precisely conserved, Z Z −ρS(x, t) dx = −ρS(x, 0) dx. x
x
This corresponds to the scalar entropy conservation (3). 2.2 Entropy stable schemes for the Navier-Stokes equations We now turn our attention to the construction of entropy stable schemes for NSE (8). We use the conservative differences for convective flux and standard centered differences for the dissipative terms on the RHS. d 1 ∗ ∗ = d , (11) − 2d + d uν (t) + fν+ 1 − fν− 1 ν+1 ν ν−1 2 2 dt ∆x (∆x)2 Here, dν := d(uν ). As in the scalar case, we seek entropy conservative flux ∗ fν+ 1 , so that the entropy decay will be dictated solely by the dissipation terms 2 ∗ on the right of (11). The construction of the entropy conservative flux fν+ 1 2 follows [Ta2004] and [TZ2006], and is summarized in the following Algorithm 1 If
∗ fν+ 1 = f (vν ); else
uν = uν+1 then
2
• Set ∆uν+ 12 := uν+1 − uν . Starting with u1ν+ 1 := uν , compute recussively 2 the intermediate states, E D j j b b uj+1 rjν+ 1 , j = 1, 2, 3, . ` 1 , ∆u 1 (12) 1 = u 1 + ν+ ν+ ν+ ν+ 2
2
2
2
2
j {b `ν+ 12 }
Here, and {b rjν+ 1 } are the left and right eigensystems of the Roe 2 matrix A(uν , uν+1)(see [Roe1981]). • Set rjν+ 1 := v(uj+1 ) − v(ujν+ 1 ) and let {`j }3j=1 be the corresponding ν+ 12 2 2 orthogonal system. Compute the entropy-conservative numerical flux, ∗ fν+ 1 2
3 mj+1 − mjν+ 1 X ν+ 12 2 D E `jν+ 1 , = (γ − 1) j 2 j=1 `ν+ 1 , ∆vν+ 12
D
j+1 j `jν+ 1 , vν+ 1 − v ν+ 1 2
2
E
= δjk
2
2
(13) Here, v(u) := Uu (u) = [−E/e − S + γ + 1, q/θ, −1/θ]> are the entropy variables, and {mj } are intermediate momentum values along the path. We now arrive at our main result of NSE corresponding to the Burgers statement (6). Theorem 1 ([TZ2006], Theorem 3.6). The semi-discrete scheme (11) ∗ with the entropy conservative numerical flux fν+ 1 in (12)-(13) and d(uν ) 2
being the dissipative NS flux, is entropy stable in the sense that 3
√ We let zbν+ 1 = zν + zν+1 /2 and zeν+ 1 = zν zν+1 . 2
2
3
6
Eitan Tadmor and Weigang Zhong
X d X [−ρν (t)Sν (t)] ∆x = − dt ν ∆x ν
*
= − (λ + 2µ)
∆vν+ 1 , 2
∆vν+ 12 2
X ∆qν+ 1 2
∆x
ν
−κ
∆dν+ 12
X ∆θν+ 1 2 2
ν
∆x
g 1/θ
∆vν+ 1
+
2
d 1/θ
2 ν+ 12
ν+ 12
(14)
∆x
∆x ≤ 0.
This statement is a discrete analog of the entropy balance statement (10). Here is our main point: we introduce no excessive entropy dissipation due to spurious, artificial numerical viscosity. According to (14), the semi-discrete scheme contains the precise amount of numerical viscosity to enforce the correct entropy dissipation dictated by NSE. More can be found in [Ta2004, TZ2006]. 2.3 Numerical experiments We consider ideal polytropic gas equations as an approximation of air with γ = 1.4,
Cv = 716,
κ = 0.03,
λ + 2µ = 2.28 × 10−5 .
We simulate the Sod’s shocktube problem, where the Euler and NSE are solved over the interval [0, 1] subject to Riemann initial conditions (1.0, 0.0, 2.5) 0 < x ≤ 0.5 (ρ, m, E)t=0 = (0.125, 0.0, 0.25) 0.5 < x < 1. In Fig.2, we display the numerical solutions for the fully discrete scheme (11) with RK3 method (7) and the numerical flux (13). The density fields of four different cases are recorded. Density field of the Euler equations 2(a) demonstrates the purely dispersive character of the entropy conservative schemes. Dispersive oscillations on the mesh scale are observed in shocks and contact regions due to the absence of any dissipative mechanism. These oscillations approach a modulated wave envelop, consult [La1986, LL1996] for more discussions on dispersive oscillations. For the results of NSE in 2(b)-2(d), the presence of heat flux causes the oscillations to be dramatically reduced around the contact discontinuity and the shock in 2(b). The viscous flux in NSE, on the other hand, is doing a better job than heat flux in damping oscillations around the shock in 2(c), but we still can observe an oscillatory behavior around the contact discontinuity. In 2(d), not only the oscillations around the shock are damped out by viscosity, but the oscillations around the contact discontinuity are significantly reduced due to the heat flux.
Novel entropy stable schemes for 1D and 2D fluid equations
1−D Euler,density,γ=1.4,Cv=716,κ=0,entropy=−ρ ln(pρ−γ),∆t=2.5e−005,∆x=0.00025,Tmax=0.1 t=0.0 t=0.05 t=0.1
0.9
1−D Euler,density,γ=1.4,Cv=716,κ=0.03,entropy=−ρ ln(pρ−γ),∆t=2.5e−005,∆x=0.00025,Tmax=0.1
1
1
7
t=0.0 t=0.05 t=0.1
0.9
0.8
0.7
0.7
0.6
0.6
ρ
ρ
0.8
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
1
0
(a) Euler
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
1
(b) Navier-Stokes with heat only
1−D N−S,density,λ+2µ=2.28e−005,γ=1.4,Cv=716,κ=0,entropy=−ρ ln(pρ−γ),∆t=2.5e−005,∆x=0.00025
1
0.1
1−D N−S,density,λ+2µ=2.28e−005,γ=1.4,Cv=716,κ=0.03,entropy=−ρ ln(pρ−γ),∆t=2.5e−005,∆x=0.00025
1 t=0.0 t=0.05 t=0.1
0.9
t=0.0 t=0.05 t=0.1
0.9
0.8
0.7
0.7
0.6
0.6 ρ
ρ
0.8
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
(c) Navier-Stokes with viscosity only
1
0
0.1
0.2
0.3
0.4
0.5 x
0.6
0.7
0.8
0.9
(d) Navier-Stokes with viscosity and heat
Fig. 2. Density field Sod problem with 4000 spatial gridpoints, U (u) = −ρ ln pρ−γ , ∆t = 2.5 × 10−5 , ∆x = 2.5 × 10−4
3 The 2D Shallow Water Equations We turn to 2D shallow water equations, ∂ ∂ ∂ ∂ ∂ ∂ ∂ u+ f (u) + g(u) = η h d(u) + η h d(u) , ∂t ∂x ∂y ∂x ∂x ∂y ∂y
(15)
with u = [h, uh, vh]> being the vector of conserved variables balanced by the flux vectors f = [uh, u2h + gh2 /2, uvh]> ,
g = [vh, uvh, v2 h + gh2 /2]>,
1
8
Eitan Tadmor and Weigang Zhong
and the viscous flux vector d = [0, u, v]> . Here, h = h(x, t) is the total water height, (u(x, t), v(x, t)) are the depth-averaged velocities along x and y direction. Finally, g is the constant acceleration due to gravity, and η > 0 is the constant eddy viscosity which models the turbulence stress in the flow. The total energy U (u) = (gh2 +u2 h+v2 h)/2 serves as an entropy function, ∂ ∂ ∂ U (u) + F (u) + G(u) = −ηh(u2x + u2y + vx2 + vy2 ), ∂t ∂x ∂y 3
2
2
(16) 3
h where F (u) = guh2 + u h+uv − huux and G(u) = gvh2 + u vh+v 2 2 are the entropy fluxes. Spatial integration of (16) yields Z Z Z Z d U (u) dxdy = −η h(u2x + u2y + vx2 + vy2 ) dxdy. dt y x y x
h
− hvvy
(17)
For the inviscid case (η = 0), the global entropy conservation is satisfied, Z Z Z Z U (u(x, t)) dxdy = U (u(x, 0)) dxdy y
x
y
x
Arguing along the same line as the above NSE dimension by dimension, we obtain the entropy-stable semi-discrete schemes (recall z[ ν+ 1 := (zν+1 +zν )/2) 2
d 1 ∗ 1 ∗ ∗ uν, µ (t) + (fν+ 1 , µ − fν− (g∗ 1 1 − g ν, µ− 12 ) , µ) + 2 2 dt ∆x ∆y ν, µ+ 2 η \ dν+1, µ − dν, µ dν, µ − dν−1, µ = (h 1 − h\ ) ν− 12 , µ ∆x ν+ 2 , µ ∆x ∆x η \ dν, µ+1 − dν, µ dν, µ − dν, µ−1 1 + (h − h\ ), ν, µ− 12 ∆y ν, µ+ 2 ∆x ∆x
(18a)
∗ ∗ with the entropy-conservative fluxes fν+ and gν, outlined in Algorithm 1 µ+ 12 2,µ 1 along x and y direction, respectively,
∗ fν+ 1 2,µ
j+1 j 2 j+1 2 j 3 g X (hν+ 12 , µ ) uν+ 12 , µ − (hν+ 12 , µ ) uν+ 12 , µ xj D j E `ν+ 1 , µ , = 2 2 j=1 `x 1 , ∆v 1 ν+ 2 , µ
∗ gν, = µ+ 1 2
g 2
3 X j=1
(hj+1 )2uj+1 ν, µ+ 12 ν, µ+ 12 D
(18b)
ν+ 2 , µ
− (hjν, µ+ 1 )2 ujν, µ+ 1 j 2 2 E `yν, µ+ 1 , yj 2 `ν, µ+ 1 , ∆vν, µ+ 12
(18c)
2
Here, uν, µ (t) denotes the discrete solution at the grid point (xν , yν , t), dν, µ := d(uν, µ ), and v := Uu = [gh − 12 (u2 + v2 ), u, v]> is the entropy variable. Numerical flux f ∗ and g∗ are constructed separately along two different phase j j pathes dictated by two sets of vectors {`x } and {`y }. {hj } and {uj } are intermediate values of height and velocity along the path. The above difference scheme (18a)-(18c) is an entropy stable scheme with no artificial viscosity in the sense that the following discrete entropy balance is satisfied,
Novel entropy stable schemes for 1D and 2D fluid equations
(
"
9
2 #
2
X ∆uν+ 1 , µ ∆vν+ 1 , µ d X 2 2 h\ U (uν, µ (t))∆x∆y = −η + ν+ 12 , µ dt ν, µ ∆x ∆x ν, µ " #) ∆uν, µ+ 12 2 ∆vν, µ+ 12 2 ∆x∆y. (19) +h\ + ν, µ+ 12 ∆y ∆y
(19) is a discrete analogue of the entropy balance statement (17). Equipped by RK3, we test the entropy stable scheme (18a)-(18c) by the 2D partial Dam-Break problem with free-slip boundary described in [FC1990]. Both the inviscid and viscous case are tested. The water surface profiles at time t = 7.2s are recorded in Fig.2.
Height at t=7.2 10 9 8 7 6 5 4 3 11 10 9 8 7 6 5 4 3 2
0
5
10
40 35 30 25 20 15 10 15
20
25
30
(a) η = 0
5 35
40 0
11 10 9 8 7 6 5 4 3 2
Height at t=7.2 10 9 8 7 6 5 4
10 9 8
10
7
9
6
8
5
7 6
4 3
5 4
40 35 30 25 20 15 10
3
0
5
10
15
20
25
30
5 35
40 0
(b) η = 10m2 · s−1
Fig. 3. Shallow water equations, Dam-Break, 200×200 m2 basin, free-slip boundary, ∆x = ∆y = 5 m, ∆t = 5 × 10−3 s
Comparing 3(b) to 3(a), we observe the improvements in smoothness of the numerical solutions. There is no analytical reference solution for this test case, but other numerical results are available in [FC1990].
References [FC1990] Fennema, R.J., Chaudhry, M.H.: Explicit methods for 2D transient freesurface flows. J. Hydraul. Eng. ASCE, 116, 1013–1034 (1990) [GST2001] Gottlieb, S., Shu, C.-W., and Tadmor, E.: Strong stability-preserving high-order time discretization methods. SIAM Review, 43, 1, 89–112 (2001) [La1986] Lax, P.D.: On dispersive difference schemes. Phys. D, 18, 250–254 (1986) [LMR2002] LeFloch, P.G., Mercier, J.M., and Rohde, C.: Full discrete, entropy conservative schemes of arbitrary order. SIAM J. Numer. Anal., 40, 5, 1968– 1992 (2002)
10
Eitan Tadmor and Weigang Zhong
[LL1996] Levermore, C.D., Liu, J.-G.: Oscillations arising in numerical experiments. Phys. D, 99, 191–216 (1996) [Roe1981] Roe, P.L.: Approximate Riemann solvers, parameter vectors and difference schemes. J. Comp. Phys., 43, 357–372 (1981) [Ta1987] Tadmor, E.: The numerical viscosity of entropy stable schemes for systems of conservation laws, I. Math. Comp., 49, 91–103 (1987) [Ta2004] Tadmor, E.: Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. Acta Numerica, 451–512 (2004) [TZ2006] Tadmor, E. and Zhong, W.: Entropy stable approximations of NavierStokes equations with no artificial numerical viscosity. J. of Hype. Diff. Equa., 3, 3, 529–559 (2006)
