Developments of entropy-stable residual distribution methods for conservation laws I: Scalar problems

Accepted Manuscript Developments of Entropy-Stable Residual Distribution Methods for Conservation Laws I: Scalar Problems

Farzad Ismail, Hossain Chizari

PII: DOI: Reference:

S0021-9991(16)30573-3 http://dx.doi.org/10.1016/j.jcp.2016.10.065 YJCPH 6943

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

29 January 2016 7 October 2016 29 October 2016

Please cite this article in press as: F. Ismail, H. Chizari, Developments of Entropy-Stable Residual Distribution Methods for Conservation Laws I: Scalar Problems, J. Comput. Phys. (2016), http://dx.doi.org/10.1016/j.jcp.2016.10.065

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Developments of Entropy-Stable Residual Distribution Methods for Conservation Laws I: Scalar Problems Farzad Ismaila , Hossain Chizari1 a

School of Aerospace Engineering, Engineering Campus, Universiti Sains Malaysia, 14300 Pulau Pinang, Malaysia

Abstract This paper presents preliminary developments of entropy-stable residual distribution methods for scalar problems. Controlling entropy generation is achieved by formulating an entropy conserved signals distribution coupled with an entropy-stable signals distribution. Numerical results of the entropystable residual distribution methods are accurate and comparable with the classic residual distribution methods for steady-state problems. High order accurate extensions for the new method on steady-state problems are also demonstrated. Moreover, the new method preserves second order accuracy on unsteady problems using an explicit time integration scheme. The idea of the multi-dimensional entropy-stable residual distribution method is generic enough to be extended to the system of hyperbolic equations, which will be presented in the sequel of this paper. Keywords: entropy conservation, entropy-stable, residual distribution, semi-discrete, scalar advection

Email addresses: [email protected], Tel: 604-599 5902, Fax: 604-594 1025 (Farzad Ismail), hc12 [email protected], [email protected] (Hossain Chizari)

Preprint submitted to Journal of Computational Physics

October 31, 2016

1. Introduction In the CFD folklore, it has been widely believed that the residual distribution (RD) methods are more accurate on irregular triangular grids relative to a finite volume (FV) methods. The work of [9, 16] presented a formal justification to this claim. The upwind RD method was first developed by Roe [30] and continued to be developed by other researchers [1, 10, 11, 35, 43]. More recent work which includes unsteady problems [5, 28], high order methods [6, 8, 21, 23, 33] and advection-diffusion problems [25, 26, 29] had also been done but will not be discussed here. The details of these RD developments can be obtained in [2, 3] and the references therein. In hind sight, the RD method is designed based on accurately capturing multi-dimensional wave-characteristics without having to rely on the grid-orientation. On the other hand, the FV approach would generally would seek waves ”normal” to a grid interface to solve a local one-dimensional Riemann-problem even for multi-dimensional settings. As a result, there is the unwanted diffusion particularly when the flow is diagonal to the grids unlike the RD methods. In addition, RD methods have a much compact stencil compared to the FV methods and this presents a nice platform for parallel computation. Compact stencils would also allow the construction of a powerful and efficient implicit solver based on Newton’s method as recently reported by [23]. Overall, there is very strong appeal to develop numerical methods based on the RD approach. Most recent work in RD includes the extension to compressible Navier-Stokes equations [12] and a high order RD methods on arbitrary grids [23]. The concept of entropy-control was first studied by Tadmor [37, 38, 39] 2

defining the concept of entropy-conserved and entropy-stable finite volume methods. This was followed by his later work [40]. Tadmor and Zhong [41] continued the developments of an entropy-stable methods using physical diffusion coupled with the entropy-conserved discrete method. Ismail and Roe [19] presented an explicit and well-formed entropy-conserved fluxes for the system of Euler equations coupled with Barth’s [7] entropy-stable diffusion part. Further work of entropy-stable FV methods have been done, particularly the extension to high order accuracy [13, 14], to other systems of equations [15, 20, 44] and for different boundary conditions [27, 36]. An extension of the entropy-stable approach for the compressible Navier-Stokes equations using fully first-order-systems has also been done [24]. Almost all of these work are based on the FV approach although the foundation of [19] is based on a one-dimensional residual-distribution method. The work of [4] provided an entropy-stable RD method where the signals are based on the entropy variables. However, their proposed method does not precisely address entropy conservation for smooth flows and requires a specific averaging within each element to ensure conservation of the primary variables. [34] introduced an ”entropy-fix” to avoid non-physical shocks, similar to the approach done to the original Roe-flux in the FV approach. This paper aims to develop a multi-dimensional entropy-stable residual distribution method to solve hyperbolic-type equations. There would be developments to include discrete entropy-conservation when entropy should be physically conserved and discrete entropy production where there is physical entropy being generated, both enforced locally within an element. Since this is a preliminary attempt to study a multi-dimensional entropy-stable RD ap-

3

proach, the focus would be on scalar equations on the interior computational domain. Note that entropy-stability can be lost at the boundaries [27] if inappropriate boundary conditions are used but will not be done herein. The paper begins with a concept review of entropy conservation and entropy-stability in Section 2 followed by a brief discussion of entropy-control using a one-dimensional RD approach (Section 3). An alternative residual distribution approach based on ”flux-difference” will be introduced in Section 4 to bridge the gap between the one dimensional and multidimensional entropy-stable RD methods. The developments of an entropy-conserved and entropy-stable multi-dimensional RD approach will be done on the Burgers’ equation in Section 5. Numerical results will be included to verify the analytical work in Section 6 followed by some concluding remarks in Section 7. 2. Entropy-Conservation and Entropy-Stability We define a hyperbolic conservation law equation in the form of ut + ∇ · f = 0.

(1)

The hyperbolic conservation law has an entropy function (U ) and its fluxes (F ) satisfying a mathematical entropy condition given by Ut + ∇ · F ≤ 0.

(2)

The inequality is to cater for discontinuous flow where entropy is produced. The mathematical convention is that entropy is a decreasing function which is opposite to the nature of physical entropy. Based on the work of [18], there 4

exists a valid entropy pair (U, F ) which must satisfy ∂U ∂ f ∂ F = . ∂u ∂u ∂u

(3)

The entropy production in a domain Ω with boundary ∂Ω is defined as        ˙ Ut + ∇ · F dAdt = U dA − F dt , (4) U= Ω

∂Ω

for which U˙ = 0 implies entropy conservation and that U˙ < 0 denotes entropy is decreasing (or entropy-stable). For smooth flows, the equality of Eq. (2) must be satisfied which implies zero entropy production. A discrete approach satisfying this proposition is deemed as an entropy-conserved method. On the other hand, a discrete method that satisfies the inequality of Eq. (2) is defined as entropy-stable. Entropy-stability however, only guarantees that the discrete approach produces the correct sign of entropy production but not necessarily the correct amount. An entropy-consistent method produces the ”correct” amount of entropy generation. Achieving entropy-consistency is less straightforward to achieving entropyconservation and entropy-stability since mainly, we do not know the correct amount of entropy that has to be generated at a shock. Furthermore, unlike entropy-conservation and entropy-stability, the problem is no longer local since a shock may be captured over several cells. Overall, it is quite hard in general to design a numerical method that has the precise amount of entropy generation. Too much of entropy generation would cause a smearing of the discontinuous profile while too little entropy production would give rise to oscillatory behavior around the discontinuity. Thus, we will focus mainly in achieving entropy-stability in this paper. 5

hL

hR

hL

f∗

f∗

Δt

L

hR

R

L

(a) Finite volume

Δt R

(b) Residual distribution

Figure 1: Control volume for establishing entropy production. Cell domain is showing by thick lines.

3. Residual Distribution Approach in One Dimension Consider Eq. (1) in one dimension. Following [19], the two neighboring discrete states (uL,R ) with dual cell area (hL,R ) are discretized semi-discretely as h L ∂ t uL = f L − f ∗ = φ L , h R ∂ t uR = f ∗ − f R = φ R .

(5)

This has two kinds of interpretations. One, which is a finite volume interpretation where the left and right states are cell-averaged values separated by a flux-interface

∗

(Fig. 1a). The other represents point values at vertices

that surround a linear element centered at

∗

(Fig. 1b). This is a residual

distribution scheme where the total element residual (fL − fR ) is split as (fL − f ∗ ) + (f ∗ − fR ) and distributed to the left (φL ) and right (φR ) signals respectively. Define the entropy variables as v =

∂U ∂u

which gives the one-to-one map-

ping from the primary variable (u) space to the entropy function (U ) space. 6

We intend to locally impose a semi-discrete update of u in Eq. (5) which will simultaneously update U (semi-discretely) in Eq. (2). In other words, the update for u on each node (L, R) will give the local entropy update U for the respective nodes through the use of the mapping function v given by following equations. vL hL ∂t uL = hL ∂t UL = vL φL = vL (fL − f ∗ ), vR hR ∂t uR = hR ∂t UR = vR φR = vR (f ∗ − fR ).

(6)

For a closed system, entropy is only generated through the fluxes and it is our intention to correctly model the primary fluxes f such that it would mimic the entropy-fluxes F . To ensure discrete entropy conservation, the total element update which is computed as the sum of two nodal entropy values within the element has to be equal to the discrete changes in the entropy-fluxes such that ∂t (hL UL + hR UR ) = −[vf ] + [v]f ∗ = −[F ].

(7)

Note that [·] = (·)R − (·)L . If f ∗ satisfies the previous equation, then we say that the flux is an entropy-conserved flux (f C ). Entropy may be spuriously produced if f ∗ does not satisfy the previous equation. [19] defined an entropy production U˙ to be the difference between the entropy produced by an arbitrary f ∗ and the entropy-conserved flux f C . In other words, the entropy production U˙ is defined as U˙ = f ∗ (uL , uR )[v] − [vf ] + [F ],

(8)

which is numerically speaking, an alternative to Eq. (4). For the one dimensional Burgers’ equation, selecting an entropy pair of (U = u2 , F = 7

2u3 ) 3

as in [24] would yield an entropy variable v = 1 (u2L 6

∂U ∂u

= 2u. Choosing f ∗ =

+ uL uR + u2R ) = f C [19] yields zero entropy production hence f C is an

entropy-conserved method. There is a need to generate physical entropy around discontinuous flows. The entropy-conserved flux f C alone cannot do this, it needs another numerical term. Coupling the entropy-conserved flux f C with a dissipation term which depends on [v] [19] provides this avenue such that f∗ = fC − where a ¯ =

uL +uR 2

|¯ a| du [v], 2 dv

(9)

is the average local wave speed. Inserting the previous

equation into Eq. (8) would yield an entropy production of |a| du [v], U˙ = ([v]f C − [vf ] − [F ]) − [v] 2 dv |a| du [v] ≤ 0, = −[v] 2 dv and is decreasing since

du dv

(10)

is always positive. Hence we have an entropy-stable

flux function. Further details of this entropy-stable method will not be discussed herein but can be readily found in [19]. Instead, we would like to focus on extending this one dimensional entropy-conserved and consequently the entropy-stable residual distribution idea in two dimensions with the inviscid Burgers’ equation. But before that, we intend to introduce an alternative residual distribution method to bridge the gap between the one dimensional RD to the multi-dimensional RD approach. 4. Residual Distribution Method Using Flux-Difference The idea is to borrow the ”flux-difference” approach from finite volume and implement it in a residual distribution setting. This section will also 8

nj

k ni φT

i j

nk Figure 2: Residual Distributed Over Triangular Mesh.

serve to establish first and second order accuracy of the new RD method for steady and unsteady problems, with potential to achieve high-order accuracy for at least steady-state problems. 4.1. Isotropic Signal Distribution From the ideas of one-dimensional residual distribution approach as done in [19], we first present an alternative two-dimensional residual distribution method. Consider Eq. (1) in two dimensions. Assume a piecewise linear elements with data located at the nodes (p = i, j, k) of a triangle (T ) as shown in Fig. 2. The residuals over element (T ) satisfying Eq. (1) can be

9

computed using the trapezoidal rule such that   (∇ · f)dxdy = − (f dy − gdx) φT = − T

∂T

1 1 1 ≈ (fi + fj )(yj − yi ) + (fj + fk )(yk − yj ) + (fk + fi )(yi − yk ) 2 2 2 1 1 1 − (gi + gj )(xj − xi ) − (gj + gk )(xk − xj ) − (gk + gi )(xi − xk ) 2 2 2 1 = (fi (yj − yk ) + fj (yk − yi ) + fk (yi − yj ) 2 −gi (xj − xk ) − gj (xk − xj ) − gk (xi − xj ))  1 1  ∗ 1   = fp − f∗ · np , fp · np − f · np = 2 p 2 p 2 p  

(11)

=0

where

∗

is some averaged quantity which is arbitrary and denotes the degree

of freedom that is influencing the signals to each node p within element (T ). The np is the inward normal to the side of the opposite node p as given in Fig. 2. We propose a signal distribution to each node (p = i, j, k) similar to the one dimensional approach in Eq. (5). 1 φp = ((fp , gp ) − (f ∗ , g ∗ )) · np = φiso p . 2

(12)

This isotropic signals (φiso ) distribution is a central-type flux difference scheme of which the total residual of Eq. (11) is equally distributed to each of the three nodes within element T . The φiso also depends on f(u) rather than u, which is one of the key features of this alternative RD method that is quite different from the ones currently available in the literature. Similar to a FV approach, conservation of the primary variables (u) is automatic for the isotropic signals integrated over each local element since the summation 10

k (

2

−

u

β (u k 2

γ

ui

k

) k

− uj

φT

)

i

i

α 2 (uj

j (a) Isotropic signal going to i.

− ui )

j

(b) Artificial diffusion terms.

Figure 3: Isotropic and artificial signals going to point i.

of any f∗ would be zero within each element. For the time being we choose an arithmetic average of three nodal values for f∗ within the element. 4.2. Artificial Signal Distribution The isotropic signal distribution is central-type method, equally distributing the total residual to each node (Fig. 3a). We propose adding artificial terms as a means for numerical stability. The idea is to add the artificial terms such that the primary variable u is discretely conserved over a local element (cell) but at the same time, will discretely augment u on each node. Let us focus on a local element T which has nodes i, j, k. Define [·]ji = (·)j − (·)i so for node i, the newly proposed signals are distributed as 1 φi = (fi − f ∗ , gi − g ∗ ) · ni −α[u]ji − β[u]kj − γ[u]ik .  

2

(13)

φart i

The φart denotes the artificial terms on node i to offset the isotropic signals i (φiso i ) as it is illustrated in Fig. 3b. This can be be a platform to introduce 11

upwinding by cancelling some of the portion of the isotropic signals to the respective nodes to ensure the characteristics waves propagate in the correct manner. However, we shall save this exploration for future work and just focus on a central approach to ensure proper entropy-control. We assume that α, β, γ are the same for each node within element T . Similarly, the signals for nodes (j, k) are distributed as 1 φj = (fj − f ∗ , gj − g ∗ ) · nj −α[u]kj − β[u]ik − γ[u]ji ,  

2

(14)

1 φk = (fk − f ∗ , gk − g ∗ ) · nk −α[u]ik − β[u]ji − γ[u]kj .  

2

(15)

φart j

φart k

For the time being, we have not properly defined the precise form for (α, β, γ). We shall do so in the next subsection but for now, we just introduce a general framework of the coefficients. Based on Eq. (15), (α, β, γ) should contain a characteristic speed and a length vector to conform with the correct units. Note that the summation of all the artificial terms cancel out within an element since the summation of the artificial signals in Eqs. (13-15) is zero. Thus, conservation of u is still achieved in the element since the summation of isotropic signals is the total residual as defined in Eq. (11). The artificial terms will nevertheless influence u on each node. Although the artificial terms do not influence the u-space within an element (cell), we shall demonstrate in Section 5.2 that they will affect the entropy(U)-space to ensure physically correct entropy behavior within each cell. This is one of the unique features of this ”flux-difference” RD approach. Overall, the scheme is a family of RD methods of which (α, β, γ) and f ∗ will be the free parameters. In appendix A, we shall demonstrate the 12

specific values of (α, β, γ) which will recover the classic RD methods such as the Narrow (N) scheme, Low Diffusion Advection (LDA), Lax-Friedrich and Lax-Wendroff for the two dimensional linear advection equation. If we consider other than arithmetic averaging for f ∗ , we could enforce another physical constraint as will be shown in Section 5. Before we do that, we shall investigate the accuracy of the ”flux-difference” RD method. 4.3. Truncation Error Analysis for the Flux-Difference RD Method To establish the order-of-accuracy for the ”flux-difference” approach, first we need to determine its truncation error (TE) in general form with arbitrary (α, β, γ). The TE analysis will be done on a two dimensional linear advection equation with a constant characteristic speed λ = (a, b). In terms of Eq. (1),     √     we have f = λu and for this case, λ = λcell  = a2 + b2 . The TE analysis will be based on Taylor’s expansion on uniform triangular grids[9]. We are aware of the limitations of Taylor’s expansion on irregular unstructured grids. However, we intend to get an approximation of the orderof-accuracy of the ”flux-difference” RD approach on the simplest uniform triangular grids. Based on Fig. 4, the structured triangular grids would have uniform grid length and height with dimensions (h1 , h2 ). For conciseness, we rewrite the grid length in terms of h and the height in terms of a stretching factor (s) of the grid length. h1 = h,

h2 = hs.

(16)

Recall that the ”flux-difference” signal distribution has an isotropic and ar-

13

h

4

3

2 I

h2 = hs

II

II

5

1

0 I

I II

6

7

8

Figure 4: Right-running grid topology. The characteristic line is y = ab x

tificial terms given by 1 φi = (fi − f ∗ , gi − g ∗ ) · ni −α[u]ji − β[u]kj − γ[u]ik . 2 

 

(17)

φart i

φiso i

We chose (f ∗ , g ∗ ) = (f¯, g¯) as the arithmetic mean of three nodal values within an element. As it is mentioned before, (α, β, γ) contain a length scale from the cell to ensure consistent units of Eq. (17). Thus, we define the ˜ γ˜ ) as the following. dimensionless parameters (˜ α, β,     α=α ˜ λcell  h,

    β = β˜ λcell  h,

    γ = γ˜ λcell  h.

(18)

The TE analysis would depend on the isotropic and artificial signals about node 0 which depends on the neighboring nodes, TE = TEiso + TEart .

(19)

In a structured (or right-running) triangular grids, there are two types of triangle (I, II) which are shown in Fig. 4. The tangential to streamline (t) and normal to streamline (n) coordinate system will be utilized. 14

For conciseness, we have included the detailed TE analysis in appendix B, where the overall truncation error of the isotropic signals are given as TEiso =



a2 uttt + suttn + s2 utnn





+ab suttt + 2s2 − 2 uttn − sutnn   2

h2 2 √ + O h4 . +b s uttt − suttn + utnn 6 a2 + b 2

(20)

Note that uttt is the third derivative tangential to the streamline. This implies that the isotropic signals are second order accurate in general, which is to be expected since the isotropic part is based on trapezoidal integration. A highorder accurate isotropic signals can be achieved using high order numerical quadrature for integration but will not be done herein since we intend to keep a compact stencil. Nevertheless, the current isotropic signal is a good foundation for unsteady problems where preserving second order is an issue using classic RD methods. For steady-state conditions, we expect no changes in the derivatives tangential to the streamline, thus the following is recovered.   4 4 3 3 3 4 s − 5a bs + 5ab s − 2b ) ab (2a TEiso unnnnn h4 + O h5 . (21) − RR = 5/2 360 (a2 + b2 ) The isotropic signals is fourth order accurate for steady-state conditions. From appendix B.2, the truncation error for the artificial signals is



TEart = − a2 s2 unn + sutn + utt + ab 2s2 utn + s (utt − unn ) − 2utn  

α ˜ II − γ˜I − γ˜II   ˜I + α 2 +b (s (sutt − utn ) + unn ) λcell  h s (a2 + b2 )

(22)

+O (h2 ) .

From the previous equation, note that the artificial signals are second order 15

accurate when we select α = γ. For the steady state case, the truncation error reduces to      λcell  unn h      γII −2β˜I +2β˜II )(b−as)  I −˜ 2 λ + ab(α˜ I −α˜ II +˜γ2(a   cell  unnn h 2 +b2 )3/2      γII )(b2 −abs+a2 s2 )2  I −˜ λ − (α˜ I +α˜ II −˜γ12(a  unnnn h3   2 cell 2 +b2 ) s     ab(a3 s3 −2a2 bs2 +2ab2 s−b3 )(α ˜ I −α ˜ II −2β˜I +2β˜II +˜ γI −˜ γII )   − λ    unnnnn h4 cell 5/2 24(a2 +b2 )

TEart = −



(α ˜ I +α ˜ II −˜ γI −˜ γII )(b2 −abs+a2 s2 ) (a2 +b2 )s

(23)

+O (h5 ) .

which recovers second order accuracy if α = γ. The truncation error of the artificial signals limits the overall ”flux-difference” order of accuracy in steady-state conditions since the isotropic signals are fourth order accurate. We shall demonstrate in section 5.3.2 a neat way to achieve third and fourth order accurate ”flux-difference” approach for steady-state problems. Nevertheless, the overall isotropic and artificial signals can be at least second order accurate for both steady and unsteady problems. 4.4. Time Integration We can use any consistent ODE-type solvers to update node p such that (

∂u Δt  e )p = − φ, ∂t 2Ap e p

(24)

where φep is the signals at node p coming from cell (element) e and Ap is the median dual area of the point p as shown in Fig. 5. Note that the ”flux-difference” RD approach is still compact. It must be emphasized that there is also more to be explored in the ”fluxdifference” residual distribution approach in terms of u, particularly in terms of upwind-developments which we hope to investigate in subsequent papers. 16

i

m p

j

l

k Figure 5: Dual median area of a point (Ap is the shaded area).

Herein, we intend to further develop the ”flux-difference” approach using the the entropy variable v as part of entropy-control. 5. Entropy-Control in a Residual Distribution Method for the 2D Burgers’ Equation The two-dimensional Burgers’ equation satisfies Eq. (1), with f = (f, g) = 2

( u2 , u) as the respective fluxes. As done in [24], we shall define the entropy function U for the Burgers’ equation as the ”total-energy” of the system. If we select U = u2 as in [24] , we could determine the entropy-fluxes using 3 Eq. (3) to be F = (F, G) = ( 2u3 , u2 ) respectively. The entropy function U

and its fluxes (F, G) would satisfy Eq. (4). For this choice of entropy function U = u2 , the entropy variable is given exactly like the one dimensional Burgers’ equation, where v=

∂U = 2u. ∂u 17

(25)

Using the ”flux-difference” RD method presented in the previous section, we shall develop an RD method which will have two distinctive parts. One, would be the part where the newly developed RD method would conserve entropy. The other part would be designed to produce entropy with the correct sign (entropy-stability). 5.1. An Entropy-Conserved Signal Distribution Assume that each node has signals as in Eq. (12) where the isotropic signals are given by φp = 12 (fp − f∗ ). To create an entropy conserving scheme, the change of entropy due to the residuals in element (T ) must be identical to the change of entropy governed by the entropy conservation law similar to the 1D approach in Eq. (7). From the principles of [19], we have the following constraints to ensure entropy conservation.  p

(vp fp · np ) − f∗



(vp np ) =

p



Fp · np .

(26)

p

This is an equation with two unknowns f∗ = (f ∗ , g ∗ ). Note that the difference between summation of RHS and LHS represents the entropy production U˙ within the element T . The dot product of the fluxes and their normal quantities can be rewritten as follows. 

fp · np =fi (yk − yj ) + fj (yi − yk ) + fk (yj − yi )

p

+ gi (xj − xk ) + gj (xk − xi ) + gk (xi − xj ) =(fi − fj )yk + (fk − fi )yj + (fj − fk )yi + (gj − gi )xk + (gi − gk )xj + (gk − gj )xi .

18

(27)

Focusing on the f - flux, the entropy conservation constraint of Eq. (26) can be written in terms of discrete differentials as the following. (vi fi − vj fj )yk + (vk fk − vi fi )yj + (vj fj − vk fk )yi −f ∗ ((vi − vj )yk + (vk − vi )yj + (vj − vk )yi ) = (Fi − Fj )yk + (Fk − Fi )yj + (Fj − Fk )yi .

(28)

Recall that for the 2D Burgers’ equation, we have (f, g) = ( 12 u2 , u) with the entropy variable v = 2u and its respective entropy-fluxes (F, G) = ( 23 u3 , u2 ). To achieve entropy conservation within element T , we need to solve for f ∗ satisfying the previous equation. [vf − F ]ij yk + [vf − F ]ki yj + [vf − F ]jk yi [v]ij yk + [v]ki yj + [v]jk yi 3 [u]ij yk − [u]3ki yj + [u]3jk yi = 6 ([u]ij yk + [u]ki yj + [u]jk yi ) 



f∗ =

[u]ij yk = Df

 

Df

u2i

   + ui uj + u2j [u]ki yj u2k + ui uk + u2i + 6 Df 6  



∗ fij

[u]jk yi + Df

 u2j + uj uk + u2k . 6  



∗ fki

(29)

∗ fki

Note that the values inside the parentheses in the last line of the previous equation represent the exact one dimensional entropy-conserved flux on each edge as given by [19]. In other words, we can view the two-dimensional entropy-conservation constraint as an integrated entropy-conserved fluxes from each element edge. Fig. 6 depicts the entropy-conservation constraints on each edge of the element. From this perspective, the multi-dimensional 19

k

−

= [F ] jk

i

∗

f∗

f jk [v ] j k

i

k

] ki

] jk −

]k [v f

[v ]

=

[v f

∗

ki if

[F

[v f ]ij − [v ]ij f ∗ ij

= [F ]ij

j

Figure 6: Entropy conservation constraints in an element.

entropy-conservation constraint reduces to a one-dimensional constraint for which many entropy-conserved fluxes have been derived from various systems of equations. The challenge is to determine the exact weights for each edge. From the previous equation, we can see that the weights from each element edge is f wij =

[u]ij yk f [u]jk yi f [u]ki yj , wjk = , wki = . Df Df Df

(30)

We could add the value of = 10−6 in the denominator of the weights to ensure well-posedness when the velocity difference is small. However, there is a more elegant way to overcome this which is presented in appendix C. The exact integration using the entropy-conserved fluxes on each element edge to determine f ∗ within element T is, f ∗ f ∗ f ∗ f ∗ = wij fij + wjk fjk + wki fki = f C .

(31)

We can use a similar process for g. g ∗ g ∗ g ∗ g ∗ = wij gij + wjk gjk + wki gki = g C .

20

(32)

Define Dg = [u]ij xk + [u]ki xj + [u]jk xi , the weights are g = wij

[u]ij xk g [u]jk xi g [u]ki xj , wjk = , wki = , Dg Dg Dg

1 1 1 ∗ ∗ = (uj + uk ), gki = (uk + ui ). gij∗ = (ui + uj ), gjk 2 2 2

(33) (34)

We have derived an entropy-conserved flux f ∗ , g ∗ for the 2D Burgers’ equation and consecutively an entropy-conserved signal distribution as given in Eq. (12). The next step is to develop an entropy-stable method. 5.2. Entropy-Stability We propose adding artificial terms as means for entropy-stability as shown in Fig. 3b. The idea is to add the artificial terms such that the primary variable u is discretely conserved over a local element (cell) but at the same time, there will be an entropy production within the same element to decrease the mathematical entropy function U . For node i, the newly proposed signals which include the isotropic and artificial signals are distributed as 1 φi = (fi − f C , gi − g C ) · ni − α[u]ji − β[u]kj − γ[u]ik 2 1 du du du = (fi − f C , gi − g C ) · ni − α[v]ji ( )ji − β[v]kj ( )kj − γ[v]ik ( )ik 2 dv dv dv 1 α β γ = (fi − f C , gi − g C ) · ni − [v]ji − [v]kj − [v]ik , (35) 2 2 2

 2 φart i

since

du dv

=

1 . 2

The choice of [v] instead of [u] is to ensure that entropy-

stability can be achieved [7], as explained in Section 3. Similarly, the signals

21

for nodes (j, k) are distributed as 1 α β γ φj = (fj − f C , gj − g C ) · nj − [v]kj − [v]ik − [v]ji , 2 2 2

 2

(36)

1 α β γ φk = (fk − f C , gk − g C ) · nk − [v]ik − [v]ji − [v]kj . 2 2 2

 2

(37)

φart j

φart k

The artificial terms will influence u on each node but not within the element. Although the artificial terms do not affect u within a local element, we shall show that the artificial terms will influence the entropy production U˙ within the element. 5.2.1. Conditions for Entropy-Stability Entropy-stability is achieved if the discretization of u would fulfill the discrete entropy inequality as given in Eq. (2). We now intend to demonstrate the effects of (α, β, γ) on the entropy production. Similar to the one dimensional case, we define the multi-dimensional entropy production within element T as U˙ =

 p

(Fp , Gp ) · np −



vp (fp , gp ) · np − (f ∗ , g ∗ )

p



(vpnp ).

(38)

p

If we select an entropy-conserved flux (f ∗ , g ∗ ) = (f C , g C ) derived in Section 5.1, we recover the semi-discrete equality of Eq. (2) since there is zero entropy production. For any other choices of (f ∗ , g ∗ ), there will be an entropy

22

production within the element which can be rewritten in terms of the signals. U˙ =



(Fp , Gp ) · np −

p

=





(vp φp )

p

(Fp , Gp ) · np −

p



vp (fp , gp ) · np + (f C , g C )

p



(vpnp ) +

p

vj vi = − (α[v]ij + β[v]jk + γ[v]ki ) − (α[v]jk + β[v]ki + γ[v]ij ) 2 2 vk − (α[v]ki + β[v]ij + γ[v]jk )  2 

α−γ =− (vi − vj )2 + (vj − vk )2 + (vk − vi )2 4 β − (vi vj − vj vi + vj vk − vk vj + vk vi − vi vk ) . 2



vp φart p

p

(39)

If we choose α − γ ≥ 0, then U˙ ≤ 0 thus entropy will always be decreasing. Note that there will be more than one combination for α, γ to ensure entropystability. We choose, γ = −α,

α ≥ 0.

(40)

There is no condition for β since the terms associated with it cancels out at the entropy production level. The point is that although the artificial diffusion terms do not affect the element in terms of the total integration of f(u), but (excluding the β terms) they affect the entropy production U˙ within the element. The entropy-stable condition differs from the choice of α = γ in Sec 4.3.2 to achieve second order accuracy for the artificial terms. We shall present another way to achieve second order accuracy in the next subsection.

23

5.3. An Entropy-Stable Residual Distribution Approach In this section, an entropy-stable methods will be developed using the artificial signals based on the conditions obtained from the previous subsection. We also intend to establish the first, second, limited and high order methods. The idea is to propose an entropy-stable approach which has the capability of extension to the system of Euler equations. Refer to Eqs. (35) to (37). The isotropic signals will be kept fixed since the conservation of the primary variable u may be lost if the isotropic terms are altered. Furthermore, the isotropic signals are at least second order accurate for transient calculations, and fourth-order accurate for steady problems as shown in Section 4.3. However, the artificial signals are at best first order accurate unless the α = γ condition is met, for which contradicts the α = −γ entropy-stability condition as shown in the previous subsection. We intend to augment the artificial signals specifically by varying α to achieve second order accuracy for unsteady problems and perhaps at least third order accurate for steady-state cases while maintaining the entropy-stability condition. Conservation of u will be retain if we insist (α, β, γ) to be locally invariant parameters within each element. The philosophy of achieving low order and high order terms can be viewed from the perspective of entropy-control. For smooth-type flows where entropy is mostly conserved, we expect minimal effects from the artificial signals thus we impose φart → 0. By this notion, we recover almost only the isotropic signals which conserve entropy and are at least second order accurate for unsteady problems and fourth order for steady-state cases. In theory, we could achieve higher accuracy if we use a high order numerical quadrature 24

for the integration of the isotropic signals as mentioned in Section 4.3, but it will be part of future work. For shock-dominated flows, we expect entropy to be produced and therefore φart to have more influence. In this case, we want first order accuracy to be implemented to overcome Godunov’s theorem, and maybe this where the current form of the artificial signals is appropriate. Overall, the variation between first and second (or higher) order methods depends only on changing the artificial signals φart without altering the entropy-conserved (isotropic) signals. The limiting procedure will be an adaptive approach between the entropy-stable (first) and entropy-conserved (high order) methods, but also specifically by controlling φart . This idea of limiting was first presented by [22] for FV methods. 5.3.1. Designing the Artificial Signals for First-Order, Second-Order and Limited Methods From the entropy stable condition in Eq. (40), we have γ = −α and that β is a free parameter which does not influence entropy generation within an element hence we set β = 0. Thus, we could simplify the artificial signals in Eq. (35)-(37) as

⎧ ⎪ ⎪ φart = − 12 α ([v]ji − [v]ik ) , ⎪ ⎨ i

1 φart j = − 2 α ([v]kj − [v]ji ) , ⎪ ⎪ ⎪ ⎩ φart = − 1 α ([v] − [v] ) . ik kj k 2

(41)

The artificial terms alter the u-space for each node but not alter the u-space overall within each element. We therefore have flexibility to alter its form, particularly the (α) to achieve the appropriate order-of-accuracy which is measured at the nodes. 25

Recall that from Eq. (B.12) the dominant error terms originate from the artificial signals as given by



TEart = − a2 s2 unn + sutn + utt + ab 2s2 utn + s (utt − unn ) − 2utn  

α ˜ II − γ˜I − γ˜II   ˜I + α +b2 (s (sutt − utn ) + unn ) λcell  h s (a2 + b2 )

(42)

+O (h2 ) .

The artificial signals are fictitious terms which are trivial within each element provided the same α, β, γ are used for each node within that element. We intend to exploit this fact. If we multiply the first order error terms by some factor of O(h), we could increase the order of accuracy to second order. To ensure proper units and length scales for α we propose,  q   h   α= λ  cell  h, Lr  

(43)

α ˜

where Lr is the reference length for numerical case. This is yet another degree of freedom, that is the average value of λcell , for which the arithmetic mean within the element is used for simplicity. For fully unstructured grid, the h is calculated from the local grid size, h=

1 li , 3

(44)

where li is the length of ith edge of the cell. From Eq. (43), it shows that α ˜  O (hq ). Choosing α ˜ with any q and γ˜ = −˜ α, will yield an overall order of accuracy for the artificial signals to be O (hq+1 ). Therefore, we could produce first and second order entropy-stable methods accordingly based on the choice of q. Note that we can only attain at best second order accuracy in unsteady cases even if q > 1 is selected since the current isotropic signals (trapezoidal) 26

are at best second order accurate. The magnitude of the artificial signals which produce entropy is inversely proportional to q since

h Lr

<< 1. This is

consistent with our assumption of high order accurate methods will produce less entropy production. First Order: To construct a first order entropy-stable method, we can just use the original φart with α = −γ and that β = 0. In a more structured form, this can viewed by choosing α such that  α

1st

=

h Lr

q1st     λcell  h,

q 1 → 0+ . st

(45)

The overall ”flux-difference” signals are art,1 φ1i = φiso , i + φi st

st

st 1 st φart,1 = − α1 ([v]ji − [v]ik ) . i 2

(46)

Second Order: To produce a second order method, we use  α

2nd

=

h Lr

q2nd     λcell  h,

q2

→ 1− .

(47)

1 nd = − α2 ([v]ji − [v]ik ) . 2

(48)

nd

Consequently, φ2i

nd

art,2 = φiso , i + φi nd

φart,2 i

nd

Limited Approach: The current limiting procedure in RD [31] is a postprocessing approach to alter the signals in which entropy-stability may be lost. We shall construct a limited approach that blends the low and high order methods which will still retain entropy-stability. To do this, we shall design a locally adaptive selection of different values of q to generate different 27

order of accuracy. For a start, we shall focus on only limiting second order accurate solutions.  q1st   h st   α1 = λcell  h, Lr

 α

2nd

=

h Lr

q2nd     λcell  h.

(49)

Analytically, the q → 0+ will generate the first order accuracy while the q → 1− will produce second order solutions. Controlling the order of accuracy is implemented on the artificial signals only since we might lose conservation of the primary variable within each element if we alter the isotropic signals. Numerical experiments in this paper use the lower-bound values q 1 and and st

upper bound values q 2

nd

to denote limited methods.

The limited entropy-stable approach uses local variations of entropy as a sensor to indicate decide whether to use first order, second order or anywhere in between. This is where we differ from the current limiting approach for which the sensor is based on the changes of the primary variables u, either in terms of signals [31] or basic TVD methods [17]. We propose a limited approach by blending low and high order artificial terms. Thus, art,lim φlim = φiso , i i + φi

(50)

where, φart,lim = ψ(S)φart,1 + (1 − ψ(S))φart,2 . i i i st

nd

(51)

and φart,high are calculated by q low and q high , respectively Note that, φart,low i i (Eqs. (46) and (48)). Therefore, they produce low order and high order results. In the scalar equation the slope of the (negative) entropy function (S = −U ) is chosen as the local parameter to decide between low and high order approach. We do this to be consistent with the Euler approach which 28

k mki mmax i

mij

mjk

j

Figure 7: The entropy sensor topology.

we intend to present in the sequel, where we used the physical entropy S (negative of mathematical entropy function) as the indicator. In principle, for relatively smooth flow we expect very little change in entropy, thus the limiter would produce closer toward second order accuracy. However, when there is discontinuity, we expect entropy being generated hence the sensor would detect these entropy generations to limit the scheme to first order. The sensor is built on the spatial slope of entropy calculated from each neighboring cells. According to Fig. 7, for a compact limiting the absolute slope of three neighbors within each element is determined at first. Of course, we could extend the stencil size to perform the slope limiting for shock detection but will not do so herein to focus on a compact approach. We are only interested in the magnitude of the slope difference in developing the sensor. For instance, the calculation on nodes (j, k) is given by the following.    S j − Sk  , mjk =  li 

(52)

where S = −u2 is the increment from the inflow which will be always positive. The maximum of three slopes is taken as the input to the sensor such that mmax = max(mij , mjk , mki ). 29

(53)

1

iter

-1

0.8 0.6

rLi

m

ite

0.4

2

Lim

ψ(S)

y 2 = Cx y = Cx2

Low Order

0.2 1 mcritical, max

0

2 mcritical, max

High Order

mmax 1, 2 Figure 8: The function ψ(S) for the sensor with mcritical, = 200, 400 to construct max

limiter-1 and limiter-2, respectively. Note that ψ(S) = 1 is for low order (dotted line); and, ψ(S) = 0 is for high order (long dashed line). The two functions shown are the non-linear possibilities for the limiter.

The ψ(S) depicts the portion of low and high order methods as demonstrated in Fig. 8, which is defined as ⎧ ⎪ ⎨ 1 ψ(S) = mmax ⎪ ⎩ critical mmax

mmax ≥ mcritical max , mmax < mcritical max .

(54)

The choice of mcritical defines the type of limiter. If the slope (mmax ) becomes max very small ψ(S) will be close to 0 hence, the solution will be high order art,2 (φart ). As the mmax is increasing ψ(S) would move from 0 through i → φi nd

1. Therefore, the method will be more toward first order (φart → φart,1 ). i i st

The ψ(S) could also be constructed with a parabolic or cubic function as shown in Fig. 8. Note that this the limited-entropy-diagram is similar to the classic TVD/TVB diagram for which other type of limiters can be developed. 30

5.3.2. High Order Methods for Steady-State Conditions From Section 4.3, it is clear that the obstacle to achieve high order spatial accuracy in steady-state advection problems is due to the artificial signals. From Eq. (34), selecting β = 0 and imposing α = −γ will yield the truncation error of the artificial signals in steady-state condition to be   −˜ γII )(b2 −abs+a2 s2 ) TE = − (α˜ I +α˜ II −˜γI(a unn h 2 +b2 )s   γII )(b2 −abs+a2 s2 )2 I −˜ unnnn h3 + O (h4 ) . − (α˜ I +α˜ II −˜γ12(a 2 2 +b2 ) s

(55)

To develop a compact high order (beyond second order) method, the main concept still remains the same. By controlling entropy generation produced by the artificial signals, one might be able to construct a high order approach. Thus, we could achieve third order by choosing  qhigh   h   high α = λcell  h, Lr

(56)

where q high → 2− . Fourth order accuracy can be theoretically achieved if q high → 3− . Observe that by increasing the q high , the value of α decreases because the h/Lr is mostly less than one. This will lead to the artificial signal becomes very small which may yield numerical instability. To resolve this issue and remain the order of accuracy we propose to use a constant coefficient (κ) as follows.  α

high

=κ

h Lr

qhigh     λcell  h,

κ > 0.

(57)

Overall, the order of accuracy will be third and fourth order. It should be noted that the κ will be determined numerically. In theory, we could achieve beyond fourth order accuracy for the artificial signals if we choose q > 3 but now the restricting factor will be the isotropic 31

signals. At best, the trapezoidal integration would give fourth order accuracy for steady-state problems. We could attempt for even higher order accuracy (beyond fourth order) by using a high order numerical quadrature to integrate the isotropic signals but will not do so for this paper since we will lose compactness of the RD method. However, the extension to very high order isotropic signals will be presented in the subsequent papers. 6. Results and Discussion We shall investigate the numerical performance of the newly proposed RD methods in this section. For the linear advection case we assume the speed to be (a, b) = (1, 1) which is computed in [0, 1] × [0, 1] domain. For the Burgers’ equation, consider the characteristic speed of (a, b) = (u, 1) with the computational domain a square area of [−1, 0]×[0, 1]. Various irregular grids were used in all of the test cases (Fig. 9) with three different grid sizes which are 8989, 35175 and 78862 elements for the order-of-accuracy study, similar to the sizes done in [2]. Excluding the unsteady linear advection case, the boundary condition for the other test cases is set to be the exact solution to eliminate the boundary condition effects. We shall use forward Euler and second order Runge-Kutta time integration methods for steady-state and unsteady computations respectively, with CFL ≤ 0.1 where its definition is obtained from [42]. For higher order methods, the choice of κ is selected based on the optimum values determined empirically as shown in appendix D.

32

(a) 0%

(b) 20%

(c) 40%

Figure 9: Irregular grids with different randomized percent.

6.1. 2D Linear Advection 6.1.1. Steady This is sine-cosine problem where the exact solution will be, u(x, y) = − cos(2π(bx − ay)),

a = b = 1.

(58)

In Fig. 12 a selective number of q is analyzed and the L2 error is shown in a logarithmic scale verses the log(h). Additionally, the order-of-accuracy for different values of q is also presented in Fig. 10. The numerical order of accuracy is tested on the equilateral (0%) grids to ensure consistency with the analytical TE analysis. Overall, there is a good agreement between the analytical order-of -accuracy and the numerical order-of-accuracy as q increases. The second order accuracy is established at q = 1, which is identical to the analytical prediction. Although not shown in Fig. 10, numerical order-of-accuracy matches the analytical order-of-accuracy up to fourth-order accurate. First order accuracy is numerically obtained at q ∼ 0.3 as opposed to 33

q = 0 from the TE analysis. For the q = 0, the difference in order-of-accuracy is larger than q = 1. This might be due two reasons. First, the analytical results (TE analysis) is based on structured (right-running) triangular grids but the numerical results are based on equilateral grids. Second, the analytical derivation is based on the size of one edge of the element whereas the numerical results are based on the average length of the element. For q ≥ 1, the error is in O(h)q+1 hence the effects of the two reasons above are perhaps negligible. However, the difference is larger when q ≤ 0.8 since the difference corresponds to roughly O(h). The order-of-accuracy of the entropy approach for various q on the irregular grids are shown in Fig. 12. Excluding the fourth-order ES method, the other ES methods seem to preserve their respective order-of-accuracies even on highly randomized grids. The order-of-accuracy behavior of classic RD methods are well understood thus omitted herein, but interested readers can learn more about it in [9]. Additionally, the cross section of the linear advection is also shown in Fig. 11 for different values of q compared to some classic RD methods, which in general depicts good agreement. The results of first-order ES method (q = 0.3) is very similar to the Lax-Friedrich, thus omitted for brevity. This is not surprising since both are first order central methods. Perhaps a first order upwind entropy-stable method would be more accurate than a central entropy-stable approach, maybe as good as the N-scheme.

34

Exact Equilateral Grid

Order of Accuracy

2

1.5

1

0.5 0

0.1

0.2

0.3

0.4

0.5 q

0.6

0.7

0.8

0.9

1

Figure 10: The linear advection order of accuracy for the ES approach with different values for q in equilateral (0%) grid types.

1

u

0.5 0 −0.5 −1 0

Exact q = 1.0 q = 2.0 q = 3.0 N LDA LxF LxW 0.1 0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Figure 11: Cross section at y = 0.5 with different methods for 2D linear advection case based on the ES method.

35

q = 0.3, 1.03 q = 2.0, 2.90

100

q = 1.0, 1.96 q = 3.0, 3.95

q = 0.3, 1.06 q = 2.0, 2.68

q = 1.0, 1.81 q = 3.0, 3.00

L2 error

10−1 10−2 10−3 10−4

10−2.4 10−2.2 10−2 Grid Distance log(h)

10−2.4 10−2.2 10−2 Grid Distance log(h)

Figure 12: Numerical L2 error verses the grid distance in logarithmic scale for 2D linear advection case for equilateral (left side) and 40% randomized element (right side). The order of accuracy is shown in the legend.

6.1.2. Unsteady This is a rotating cosine problem as defined by [2]. The advection velocity will be, (a, b) = (−2πy, 2πx),

(59)

which constructs a rotation around (0, 0) in which after t = 1 the exact solution will be in the initial condition. The cosine wave for the initial condition is, ⎧ ⎨ u=

1 2

⎩ 0

(1 + cos (4πr)) r <

1 4

r≥

1 4

 ,

r=

1 x− 2

2 + y2.

(60)

Note that, the domain will be [−1, 1] × [−1, 1]. Fig. 13-14 demonstrates that the entropy-stable RD methods preserve second order-accuracy even on unsteady problems, unlike the LDA. The compact-limited entropy-stable RD 36

L2 error

10−1

q = 1.0, 1.84 q = 3.0, 2.55 Lim-1, 1.37 N, 0.48 LDA, 0.47 LxF, 0.34 LxW, 1.72

10−2

10−3

10−2.4

10−2.3 10−2.2 10−2.1 Grid Distance log(h)

10−2

Figure 13: Numerical L2 error verses the grid size in logarithmic scale for unsteady 2D rotating cosine case for equilateral element. The order of accuracy is shown in the legend.

method has an order of accuracy of about 1.4, but a less stringent limiter can perhaps increase its order of accuracy. This can be done by increasing mcritical max , or exploring other functions in the limited-entropy-diagram. The Lax-Wendroff method here is based on the work of [32], which also preserves second order accuracy. 6.2. The 2D Burgers’ Equations 6.2.1. Steady: Expansion Problem In the expansion case, the exact solution is u(x, y) =

Ur + x(Ur − Ul ) , 1 + y(Ur − Ul )

Ur > Ul ,

(61)

where Ur and Ul are the right and left bottom corners values. The condition Ur > Ul with the inlet and outlet boundary conditions in Fig. 15 ensures that the solution will expand smoothly in the interior domain. In this study, 37

Exact q = 1.0 q = 2.0 q = 3.0 Lim-1 LDA LxW

1 0.8

u

0.6 0.4 0.2 0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Figure 14: Cross section of the second order methods for unsteady rotating cosine case for 20% randomized grid.

Ur = 1.5 and Ul = −0.5 are chosen. The bottom boundary condition as is calculated from the above equation for y = 0. Moreover, the exact solution is also shown in Fig. 15. In order to compare the numerical results, the cross section of y = 0.6 is shown for different methods. The logarithmic L2 error convergence study is also done for this test case using different q values. Fig. 16 demonstrates that the entropy-stable RD methods preserve their respective order-of-accuracy even for the nonlinear Burgers’ equation. The error convergence for the classic RD methods are omitted to avoid cluttering, but their results are similar with the entropystable methods which is confirmed by the cross-sectional results shown in Fig. 17. The grid irregularities also have very little influence on the quality of the results. Note that we do not include the limited ES method for brevity

38

u=0

Inlet Boundary

Inlet Boundary

(a) Initial Condition

(b) Exact Solution

Figure 15: The Burgers’ expansion case.

but it produces similar results trend as before. 6.2.2. Steady: Shock-Tree Problem Assume the Burgers’ equation with the inflow boundary at the bottom, u(x, 0) = 1.5 − 2x,

x ∈ [−1, 0].

The steady state exact solution in the [−1, 0] × [0, 1] domain is, ⎧ x− 3 ⎪ −0.5, y ≤ 12 & y− 14 < ⎪ ⎪ 2 ⎪ ⎨ 3 x− u(x, y) = 1.5, y ≤ 12 & y− 14 > ⎪ 2 ⎪    ⎪ ⎪ ⎩ max −0.5, min 1.5, x− 34 (elsewhere). y− 1

(62)

1 2 1 2

(63)

2

This is an opportunity to test the various limited approach in section 5.3.1. The choice of mcritical which is the critical slope of the entropy for limiting max approach is defined as Lim-1 (mcritical = 200) and Lim-2 (mcritical = 400). max max Note that, the high order of the limited entropy approach is based on q high = 1.0, 2.0. Figs. 18-21 show the results in two cross sections at y = 0.3 and 39

L2 error

10−2

q = 0.3, 1.16 q = 1.0, 1.93 q = 2.0, 2.90

10−3

10−4

10−5

10−2.4

10−2.3 10−2.2 10−2.1 Grid Distance log(h)

10−2

Figure 16: Numerical L2 error verses the grid distance in logarithmic scale for 2D expansion case for equilateral elements. The order of accuracy is shown in the legend.

Exact q = 0.3 q = 1.0 q = 2.0 q = 3.0 N LDA LxF LxW

u

0.5

0

−1 −0.8−0.6−0.4−0.2 0 −1 −0.8−0.6−0.4−0.2 0 x x Figure 17: Cross section at y = 0.6 with different methods for nonlinear expansion case in the entropy stable methods. Equilateral on the left and 40% grid disturbance on the right.

40

1.5

u

1 0.5 0

Exact q = 0.3 N −0.5 LxF −1 −0.8 −0.6 −0.4 −0.2 x

Exact q = 1.0 q = 2.0 LDA LxW 0 −1 −0.8 −0.6 −0.4 −0.2 x

0

Figure 18: Cross section at y = 0.3 for first order methods with 20% randomized grid.

y = 0.6 for the smooth expansion and shock regions. Clearly, the first order methods are diffusive but the second order approaches have spurious oscillations around the shock. The limiters remove these oscillations yet produce sharp profiles for both shock and expansion regions. 6.2.3. Unsteady: Shock Expansion Case The initial condition has two stationary shocks defined as, ⎧ ⎪ ⎨ 1, − 23 < x < − 13 u= ⎪ ⎩ −1, −1 ≤ x ≤ − 2 & − 1 < x < 0. 3 3

(64)

It should be mentioned that this case is tested in a 2D domain. This is a problem where there will be strong expansion wave which may interact with a shockwave if we let the time run long enough. We consider the unsteady solution at a time where expansion wave does not meet with the shockwave. 41

1.5

u

1 0.5 0 Exact Lim-1 −0.5 Lim-2 −1 −0.8 −0.6 −0.4 −0.2 x

Exact Lim-1 Lim-2 0 −1 −0.8 −0.6 −0.4 −0.2 x

0

Figure 19: Cross section at y = 0.3 for limited methods with 20% randomized grid. Note that in the left figure q low = 0.3 and q high = 1.0; and, in the right side q low = 0.3 and q high = 2.0.

1.5

u

1 0.5 0

Exact q = 0.3 N −0.5 LxF −1 −0.8 −0.6 −0.4 −0.2 x

Exact q = 1.0 q = 2.0 LDA LxW 0 −1 −0.8 −0.6 −0.4 −0.2 x

0

Figure 20: Cross section at y = 0.6 for the first and second order methods with 20% randomized grid.

42

1.5

u

1 0.5 0

Exact Lim-1 −0.5 Lim-2 −1 −0.8 −0.6 −0.4 −0.2 x

Exact Lim-1 Lim-2 0 −1 −0.8 −0.6 −0.4 −0.2 x

0

Figure 21: Cross section at y = 0.6 for limited methods with 20% randomized grid. Note that in the left figure q low = 0.3 and q high = 1.0; and, in the right side q low = 0.3 and q high = 2.0.

The particular time is determined as follows.  1  1 1 3 = . tfinal = 2 umax 6

(65)

The cross-sectional results for different methods are demonstrated in Figs. 22 and 23 for the equilateral elements and 40% randomized grids, respectively. Notice that, all the entropy-stable RD methods produce a smooth expansion region which demonstrates the effect of having a good entropy-control. However, almost all of the classic RD methods show a jump (”dog-leg”) in the middle of the expansion region, especially the Lax-Wendroff method. This gets worse as on the irregular grids. The spurious oscillations also grow around the shock for the Lax-Wendroff method on irregular grids. On the contrary, LDA method does not produce spurious oscillations around the shock since it becomes at best first order accurate for unsteady problems. 43

1

u

0.5 0 −0.5

Exact q = 0.3 q = 1.0 q = 2.0 q = 3.0 N LDA LxF LxW

−1 −1

−0.9

−0.8

−0.7

−0.6 x

−0.5

−0.4

−0.3

−0.2

Figure 22: Cross section of different methods for unsteady burgers case for equilateral element.

The second order ES methods produce very small oscillation around the shock, which is easily removed by the limited ES approach. 7. Concluding Remarks We have developed a multi dimensional entropy-stable RD method to solve the Burgers’ equation. The newly proposed signals distribution approach has five key features: 1. A family of isotropic and artificial signals distribution that ensures conservation within each element by default, which can be designed for upwind or central type RD methods. 2. An isotropic entropy-conserved signal distribution that is based on entropy-conservation on each element edge (leading to entropy conser44

1

u

0.5

Exact q = 1.0 q = 2.0 q = 3.0 Lim-2 LxW

0 −0.5 −1 −1

−0.9

−0.8

−0.7

−0.6 x

−0.5

−0.4

−0.3

−0.2

Figure 23: Cross section of different methods for unsteady burgers case for 40% randomization.

vation within the element), which is very similar to the one dimensional entropy-conserved finite volume methods. 3. The entropy-stable part is based on the artificial terms, to offset the isotropic (entropy-conserved) signals to ensure proper entropy generation. The entropy-stable RD method is up to fourth-order accurate in space for steady problems, and its high order versions are at least second-order accurate for unsteady problems while retaining a compact stencil. 4. The newly developed limiter is based on entropy-control, for which low order accuracies commensurate to artificial (entropy-stable) signals being dominant whereas high order accuracies reflect minimal entropygeneration since the entropy-conserved signal is dominant.

45

Method

ES

N

LDA

LxF

LxW

Time (s)

8.68 × 10−3

4.35 × 10−3

5.76 × 10−3

5.83 × 10−3

5.85 × 10−3

Table 1: The computational cost for one iteration in second.

In addition, the entropy-stable (ES) RD methods are about as accurate as classic RD methods for steady-state problems with about an additional 45 percent computational cost as shown in Tab. 1. However, the entropy-stable RD methods are far superior on unsteady problems due to their second-order spatially accurate preserving feature unlike most of the classic high-order RD methods. The entropy-stable methods do not need an ”entropy-fix” to breakup the expansion shock and also have a smoother profile prediction of the expansion region. Most of the classic RD methods tested herein suffer from having the ”dog-leg” feature when computing the expansion shock problem. Although the work herein is based on scalar problems, the fundamental ideas of the these key features above can be straightforwardly extended to the systems of hyperbolic equations due to the generality of the approach. The extension to the system of Euler equations is almost done and will be presented in the sequel of this paper. Acknowledgements We would like to thank Universiti Sains Malaysia for financially supporting this research work under the University Research grant (No: 1001 / PAERO / 814152) and to Malaysian Ministry of Higher Education Fundamental Research Grant (No: 203 / PAERO / 6071316).

46

[1] R. Abgrall. Toward the ultimate conservative scheme: Following the quest. Journal of Computational Physics, 167(2):277–315, 2001. [2] R. Abgrall. Residual distribution schemes: Current status and future trends. Computers & Fluids, 35(7):641–669, 2006. [3] R. Abgrall. A Review of Residual Distribution Schemes for Hyperbolic and Parabolic Problems: The July 2010 State of the Art. Communications in Computational Physics, 11(4):1043–1080, 2012. [4] R. Abgrall and T. Barth. Resdidual Distribution Scheme for Conservation Laws via Adaptive Quadrature. SIAM Journal on Scientific Computing, 24(3):732–769, 2002. [5] R. Abgrall and M. Mezine. Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow. Journal of Computational Physics, 188(1):16–55, 2003. [6] R. Abgrall and P.L. Roe. High-order fluctuation schemes on triangular meshes. Journal of Scientific Computing, 19(1):3–36, 2003. [7] T.J. Barth. An Introduction to Recent Developments in Theory and Numerics of Conservation Laws. In Numerical Methods for Gas-dynamic Systems On Unstructured Meshes. 1999. [8] D. Caraeni and L. Fuchs. Compact third-order multidimensional upwind scheme for Navier Stokes simulations. Theoretical and Computational Fluid Dynamics, 15:373–401, 2002.

47

[9] H. Chizari and F. Ismail. Accuracy Variations in Residual Distribution and Finite Volume Methods on Triangular Grids. Bulletin of Malaysian Mathematical Sciences Society, 22:1–34, 2015. [10] H. Deconinck, P.L. Roe, and R. Struijs. A Multidimensional generalization of Roe’s difference splitter for the Euler equations. Computers & Fluids, 22:215–222, 1993. [11] H. Deconinck, R. Struijs, G. Bourgeois, and P.L. Roe. Compact advection schemes on unstructured mehses. VKI Lecture Series: Computational Fluid Dynamics, 1993. [12] J. Dobes, M. Ricchiuto, R. Abgrall, and H. Deconinck. On hybrid residual distribution Galerkin discretizations for steady and time dependent viscous laminar flows. Computer Methods in Applied Mechanics and Engineering, 283:1336–1356, 2015. [13] T. Fischer and M. Carpenter. High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains. Journal of Computational Physics, 252(1):518–557, 2013. [14] U. Fjordholm, S. Mishra, and E. Tadmor. Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws. SIAM Journal on Numerical Analysis, 50(2):544–573, 2012. [15] G. Gassner, A. , Winters, and D. Kopriva. A well balanced and entropy conservative discontinuous Galerkin spectral element method for

48

the shallow water equations. Applied Mathematics and Computation, 272:291–308, 2016. [16] S. Guzik and C. Groth. Comparison of solution accuracy of multidimensional residual distribution and Godunov-type finite-volume methods. International Journal of Computational Fluid Dynamics, 22:61–83, 2008. [17] A. Harten. High resolution schemes for hyperbolic conservation laws. Journal of Computational Physics, 49:357–393, 1983. [18] T. Hughes, L. Franca, and M. Mallet. A new finite element formulation for compressible fluid dynamics: I. Symmetric forms of the compressible Euler and Navier Stokes equations and the second law of thermodynamics. Computer Methods in Applied Mechanics and Engineering, 54:223–234, 1986. [19] F. Ismail and P.L. Roe.

Affordable Entropy Consistent Euler flux

functions II: Entropy production at shocks. Journal of Computational Physics, 228(15):5410–5436, 2009. [20] S. Kumar and S. Mishra. Entropy stable numerical schemes for twofluid plasma equations. Journal of Scientific Computing, 52(2):401–425, 2012. [21] A. Lerat, R. Abgrall, and M. Richiuto. Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes. Journal of Computational Physics, 230(1):4103– 4136, 2011. 49

[22] Youqiong Liu, Jianhu Feng, and Jiong Ren. High resolution, entropyconsistent scheme using flux limiter for hyperbolic systems of conservation laws. Journal of Scientific Computing, 64(3):914–937, 2015. [23] A. Mazaheri and H. Nishikawa.

Improved second-order hyperbolic

residual-distribution scheme and its extension to third-order on arbitrary triangular grids. Journal of Computational Physics, 300(1):455– 491, 2015. [24] A. Mohamed and F. Ismail. Entropy Consistent Methods for the NavierStokes Equations: A First Order System Approach. Journal of Scientific Computing, 63(2):612–631, 2015. [25] H. Nishikawa. A first-order system approach for diffusion equation. I: Second-order residual-distribution schemes. Journal of Computational Physics, 227(1):315–352, 2007. [26] H. Nishikawa. A first-order system approach for diffusion equation. II: Unification of Advection-diffusion. Journal of Computational Physics, 229(11):3989–4016, 2010. [27] M. Parsani, M. Carpenter, and E. Nielsen. Entropy stable wall boundary conditions for the three-dimensional compressible Navier-Stokes equations. Journal of Computational Physics, 292(1):88–113, 2015. [28] M. Ricchiuto, R. Abgrall, and H. Deconinck. Very high order fluctuation splitting schemes for unsteady scalar advection. VKI Lecture Series, 22:215–222, 2003.

50

[29] M. Ricchiuto, N. Villedieu, R. Abgrall, and H. Deconinck. On uniformly high-order accurate residual distribution schemes for advectiondiffusion.

Journal of Computational and Applied Mathematics,

215(2):547–556, 2008. [30] P. Roe. Characteristic-based schemes for the Euler equations. Ann. Rev. Fluid Mechanics, 19:337–365, 1986. [31] P. L. Roe and D. Sidilkover. Optimum positive linear schemes for advection in two and three dimensions. SIAM Journal on Numerical Analysis, 29(6):1542–1568, 1992. [32] G. Rossiello. High order accurate Residual Distribution methods for unsteady compressible flows. PhD thesis, 2004. [33] G. Rossiello, P. De Palma, G. Pascazio, and M. Napolitano. Third-orderaccurate fluctuation splitting schemes for unsteady hyperbolic problems. Journal of Computational Physics, 222(1):332–352, 2007. [34] K. Sermeus and H. Deconinck. An entropy fix for multi-dimensional upwind residual distribution schemes. Computers & Fluids, 34:617—640, 2005. [35] R. Struijs, H. Deconinck, and P.L. Roe. Fluctuating Splitting Schemes for the 2D Euler equations. VKI Lecture Series: Computational Fluid Dynamics, 22, 1991. [36] M. Svard and H. zcan. Entropy-Stable Schemes for the Euler Equations with Far-Field and Wall Boundary Conditions. Journal of Scientific Computing, 58(1):61–89, 2014. 51

[37] E. Tadmor.

Skew-Self Adjoint Form For Systems Of Conservation

Laws. Journal of Mathematical Analysis and Applications, 103(2):428– 442, 1984. [38] E. Tadmor. Entropy Functions For Symmetric Systems Of Conservation Laws. Journal of Mathematical Analysis and Applications, 122(2):355– 359, 1987. [39] E. Tadmor. The numerical viscosity of entropy stable schemes for systems of conservation laws. Mathematics of Computation, 49(179):91– 103, 1987. [40] E. Tadmor. Entropy stability theory for difference approximations of nonlinear Conservation Laws and related time dependent problems. In Acta Numerica. 2002. [41] E. Tadmor and W. Zhong. Entropy Stable Approximations of Navier Stokes Equations with No Artificial Numerical Viscosity. Journal of Hyperbolic Differential Equations, 3:529–559, 2006. [42] E. van der Weide. Compressible Flow Simulation on Unstructured Grids using Multi-dimensional Upwind Schemes. PhD thesis, 1998. [43] E. van der Weide. Compressible flow simulation on unstructured grids using multidimensional upwind schemes. PhD thesis, Delft University of Technology, 1998. [44] A Winters and G. Gassner. An Entropy Stable Finite Volume Scheme for the Equations of Shallow Water Magnetohydrodynamics. Journal of Scientific Computing, pages 1–26, 2015. 52

Appendix A. Classic RD Recovery We could exactly recover the classic RD methods for the two dimesional linear advection equation using various α and β values. We assume f ∗ to be the arithmetic mean within three points within an element. Appendix A.1. Lax-Friedrich Method The signal distribution for Lax-Friedrich using α and β, αLxF = β

LxF

=

2|k|max −kj −kk 6 kk −kj 6

+

+

([u]ij +[u]ik )(kj [u]ik +kk [u]ij ) 3([u]2ij +[u]2jk +[u]2ki )

[u]jk (kj [u]ik +kk [u]ij ) [u]2ij +[u]2jk +[u]2ki

(A.1)

Appendix A.2. Lax-Wendroff Recovery

αLxW = −

kj +kk 6

+

β LxW = −

kj −kk 6

+

6Δtφ2T +2A([u]ij +[u]ik )(2[k]jk [u]jk −2φT ) 12A([u]2ij +[u]2jk +[u]2ki ) ΔtφT (ki [u]jk +kj [u]ki +kk [u]ij ) 2A([u]2ij +[u]2jk +[u]2ki )

+

[u]jk ([k]jk [u]jk −φT )

(A.2)

[u]2ij +[u]2jk +[u]2ki

where A is the cell area. Appendix A.3. N-scheme Recovery For two-target, αN = N

β =

([u]ji +[u]ki )φT 6([u]2ij +[u]2jk +[u]2ki ) kj −kk 6

+

(A.3)

kj [u]2ik −kk [u]2ij [u]2ij +[u]2jk +[u]2ki

For one-target cell, αOneTarget = β OneTarget =

φT (4ui +uj +uk )−3(ki u2i +kj u2j +kk u2k )



kj −kk 6

3([u]2ij +[u]2jk +[u]2ki )  6[u]2

jk [u]2ij +[u]2jk +[u]2ki

53

−1



(A.4)

Appendix A.4. LDA Recovery For two-target, αLDA = − β

LDA

=

kj +kk 6

k −k − j6 k

+ +

3φ2T +2(kj +kk )[k]jk [u]jk ([u]ij +[u]ik ) 6(kj +kk )([u]2ij +[u]2jk +[u]2ki )

φT (ki [u]jk +kj [u]ki +kk [u]ij ) (kj +kk )([u]2ij +[u]2jk +[u]2ki )

+

(kj +kk )[k]jk [u]2jk

(A.5)

(kj +kk )([u]2ij +[u]2jk +[u]2ki )

Note that the one-target LDA is identical with one-target N scheme. Appendix B. Truncation Error Analysis Following the work of [9], the first step would be to determine general spatial update equation prior to the TE analysis. The equation could be discretely written in the following form.  n  Δt un+1 w i ui + = uni − w j uj . i Ai j

(B.1)

Note that Ai is the median-dual cell area. Assume j denotes the neighboring nodes and wj is the coefficient distribution of the residual to that node. It is also assumed for this analysis, that (a, b) > 0 and b a

>

h2 h

b a

<

h2 . h

The analysis for

follows the same steps and will not be shown here for conciseness.

In the limit of steady state, un+1 → uni . Thus, the terms inside the i parentheses in Eq. B.1 will be the truncation error (TE). TE = wi ui +



w j uj

(B.2)

j

Using the Taylor series expansion of the neighboring points about the main point of interest (node 0), the TE will be determined. For a right

54

running (RR) triangular grid, the Taylor series expansion of the neighboring points about the main node 0 is given as the following [9].   ∞ d  d! ∂du 1  j d−k j k uj = (l ) (ln ) d! k=0 (d − k)!k! ∂td−k ∂nk t d=0

(B.3)

where, ltj and lnj are the tangential to streamline and normal to streamline distances respectively for node j from the main node of interest i. The ∂du ∂td−k ∂nk

is the dth order partial derivative with respect to tangential direction

along the streamline, t and normal, n. For an arbitrary ”flux-difference” RD method, node 0 will be dependent on all the surrounding element signals (φ012 , φ023 , ..., φ071 ). Appendix B.1. Truncation Error Analysis for Isotropic Signals The isotropic signal truncation error is given as the following. TEiso =

1  iso φ Ai e e

(B.4)

We could expand the signals coming from each element in terms of u. For instance, the signals from elements 012, 023 and 032 are written as follows.  ahs u0 − =− 2  bh iso φ023 = − u0 − 2  ahs − bh u0 − = 2

φiso 012

φiso 032

55

 1 (u0 + u1 + u2 ) 3  1 (u0 + u2 + u3 ) 3  1 (u0 + u3 + u5 ) 3

(B.5) (B.6) (B.7)

We could do a similar procedure for other signals coming from other elements thus, the TE reduces to the following. TEiso =

1 (as (2u1 + u2 − u3 − 2u5 − u6 + u7 ) 6hs

(B.8)

+b (−u1 + u2 + 2u3 + u5 − u6 − 2u7 )) By performing a Taylor series expansion about node 0 which includes both the derivatives tangential to streamline and normal to streamline, the overall truncation error of the isotropic signals can be written as below. TEiso =



a2 uttt + suttn + s2 utnn





+ab suttt + 2s2 − 2 uttn − sutnn   2

h2 2 √ + O h4 +b s uttt − suttn + utnn 6 a2 + b 2

(B.9)

Note that un is the first derivative of the solution normal to the streamline. Appendix B.2. Truncation Error Analysis for Artificial Signals Based on Eq. (18), the truncation error for the artificial terms of the signals could also be determined. TEart =

1  art φ Ai e e

(B.10)

By expanding each signal coming from the respective element, TEart =

||λcell || hs

(˜ α1 (3u0 − u1 − u3 − u6 ) + α ˜ 2 (3u0 − u2 − u5 − u7 )

+β˜1 (u1 − u2 + u3 − u5 + u6 − u7 ) −β˜2 (u1 + u2 − u3 + u5 − u6 + u7 ) γ2 (u0 + u1 + u3 + u6 )) −3˜ γ1 (u0 + u2 + u5 + u7 ) − 3˜ (B.11) 56

The overall truncation error for the artificial signals is given as the following.



TEart = − a2 s2 unn + sutn + utt + ab 2s2 utn + s (utt − unn ) − 2utn   

α ˜I + α ˜ II − γ˜I − γ˜II   +b2 (s (sutt − utn ) + unn ) λcell  h s (a2 + b2 ) + a3 (sutnn + uttn ) + a2 b (−sunnn + 2suttn − 2utnn + uttt )

(B.12)

2

+ab (−2sutnn + suttt + unnn − 2uttn )    α ˜ II − 2β˜I + 2β˜II + γ˜I − γ˜II   2 ˜I − α +b3 (utnn − suttn ) λcell  h + O h3 3/2 2 2 2 (a + b )

From the previous equation, note that the artificial signals are second order accurate in general when we select α = γ. Appendix C. The Denominator Issue of Entropy Conservative Flux The exact flux is, 

      vp fp · np − vp f∗ · np = Fp · np

p

p

(C.1)

p

In the scalar equation vp is a scalar value hence,  p





vp fp nxp − f ∗

y p vp g p n p

−g

∗

p



vp nxp =

y p vp np

=

 p



Fp nxp

(C.2)

y p Gp np

Therefore,  ∗

f =

p



(vp fp − Fp ) nxp  , x p vp n p

∗

g =

p

(vp gp − Gp ) nyp  y p vp n p

(C.3)

Note that, based on differential form of the denominator in Eq. (27)  p

vp nxp  0 or



vp nyp  0

p

57

⇒

vi  vj  vk

(C.4)

Since (v, f, F ) all depends on the same u values, this implies the following. vi fi − Fi  vj fj − Fj  vk fk − Fk

(C.5)

Thus in the limit form, Eq. (C.3) can be reduced to, f∗ =

vi fi − Fi , vi

g∗ =

v i g i − Gi vi

which for the Burgers’ equation one could find,  2   3    (2u) u2 − 2u3  u2  (2u)(u) − (u2 )  u  ∗ ∗  f =  = 6 , g =  = 2 i 2u 2u  i i

(C.6)

(C.7)

i

If all [v]ij , [v]jk , [v]ki ≤ , we use the above equation for (f ∗ , g ∗ ). Otherwise, we use Eq. (31)-(32). Note that = 10−6 . Moreover, for the general case, if vi , vj , vk → 0 then, based on L’Hopital’s rule,   ∂F  ∂G  ∗ ∗ , g = gi − f = fi − ∂v i ∂v i

(C.8)

, ∂G ) can be related back to u using the chain rule. ( ∂F ∂v ∂v Appendix D. The Entropy Generation Control for High Order Based on Eq. 57, the minimum value of κ might be find according to the numerical test in order to stabilize the solution. The numerical investigations are shown in Fig. D.24.

58

104

κ

103 102 101 100 0.8

1

1.2

1.4

1.6

1.8

2 q

2.2

2.4

2.6

Figure D.24: The κ vs the q for high order methods.

59

2.8

3

3.2