A high-order monotonicity-preserving scheme for hyperbolic conservation laws

Computers and Fluids 144 (2017) 86–116

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

A high-order monotonicity-preserving scheme for hyperbolic conservation laws G. Capdeville Ecole Centrale de Nantes, LHEEA Lab. (ECN/CNRS), Nantes, France

a r t i c l e

i n f o

Article history: Received 14 December 2015 Accepted 29 October 2016 Available online 8 December 2016 Keywords: Finite-volume Hyperbolicity Monotonicity MUSCL High-order accuracy Diagonalization Gasdynamics Maximum principle Least-square Coupling term Characteristic variable

a b s t r a c t This work aims at designing a numerical model for discretizing non-linear hyperbolic systems of conservation laws, with high-order accuracy, while verifying a positivity property of the discretization. From a finite-volume discretization of the PDEs, we identify three positivity constraints for the resulting numerical scheme to be a convex combination of the stencil supporting the discretization. These constraints are then introduced into the MUSCL-like least-square interpolation procedure by using limiters conveniently sized. Lastly, a high-order stability-preserving time integration method (SSPRK) enables to generate the final algorithm. Extensive numerical tests on scalar problems help to assess the benefits and potentialities of the resulting algorithm. Then, upon using a local diagonalization, we extend the scope of the algorithm to non linear hyperbolic systems. To this aim, we propose two strategies for identifying local directions that promotes this diagonalization. Numerical tests on standard problems of gas dynamics, help to appreciate the advantages and drawbacks of each strategy. The resulting numerical scheme is a fourth-order least-square positive scheme (LSQP4) and it is designed and formulated for discretizing unsteady multi-dimensional hyperbolic conservation laws, over non-Cartesian meshes. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction High-order accuracy along with monotonicity-preservation are mathematical properties difficult to simultaneously assess when we wish to design a numerical model for the discretization of conservation laws. However, these properties are desirable as long as we want to compute highly varying solutions with embedded discontinuities that often arise in gas dynamics, aero-acoustics or magneto-hydrodynamics, for example. Since numerous years, many attempts have been made to conciliate these constraints in order to generate robust numerical methods (see [14], for a good overview on this matter). With this paper, we propose a numerical model that explores a new way of introducing positivity constraints into a high-order accuracy discretization method. To begin, we set the foundations for the numerical method by modeling non linear scalar conservation laws.

E-mail address: [email protected] http://dx.doi.org/10.1016/j.compfluid.2016.10.029 0045-7930/© 2016 Elsevier Ltd. All rights reserved.

In the framework of a finite-volume formulation combined with a MUSCL-like least-square procedure, we can withdraw three sufficient conditions for the resulting discretization to be a convex combination of the numerical stencil typifying the high-order discretization. Based upon the original work of S. Spekreijse [10], on positive schemes, T. J. Barth, early in the 90s (see [1] for practical details), derived conditions for generating a monotonicitypreserving numerical scheme for scalar conservation laws. The solution we promote in this paper closely follows this line of thought – at least for the scalar case. However, instead of generating a only one limiter from these conditions, as T. J. Barth did in [1], we split the problem of monotonicity preservation in three sub-problems: each sub-problem is then handled with a specific limiter. By doing so, we intend to alleviate the henceforth classical deficiency of the family of limiters inspired from Barth’s work: the loss of accuracy at local and regular extrema. To reach this goal, we generate a “positivity monitor”, which is designed to check the positivity-preservation in time of the resulting discretization. When the numerical scheme is positive, there is no need to activate the conditions that would be detrimental for the accuracy at extrema. So, those conditions are not activated.

G. Capdeville / Computers and Fluids 144 (2017) 86–116

On the other side, when there is a local loss of positivity in the numerical scheme, the specific condition is activated. In addition, the splitting in three sub-problems makes it possible to identify two supplementary sufficient conditions: a condition that prevents the possible oscillatory behavior of the MUSCL procedure at discontinuities and a condition that ensures that this latter procedure preserves a “physical“ meaning into the interpolated quantities. Limiters are specifically formulated so as to deal with those conditions without destroying the formal accuracy of the scheme. Lastly, by employing a stability-preserving high-order Runge–Kutta algorithm (SSPRK, [5]), the monotonicity preservation in time follows. Numerical tests on scalar problems enable to check the correct behavior of the resulting numerical model. Usually, the extension of a scalar numerical procedure to nonlinear systems, is far from being straightforward. Indeed, unless the non-linear system is diagonalizable, there is no rigorous mathematical theory for such an extension. In the last two decades, among the huge amount of solutions displayed to circumvent this fact, H. Deconinck et al., provided an interesting and original approach by designing an approximated diagonalization of the Euler equations (see [3,11], for a good and practical overview on this technique). Unfortunately, the method they promoted was only derived for steady computations and experimented spurious non-linear oscillations that prevented the convergence of the algorithm towards the steady-state (see [11], for example). Since we aim at computing unsteady solutions, our point of view is slightly different. Moreover, by capitalizing upon the experience of the pre-cited authors, we could derive two numerical solutions that permit to generate an approximate diagonalization of the Euler equations, this without destroying the monotonicity properties of the initial scalar procedure, for the second solution. The first solution is inspired from [3] and adapted to our framework. The second solution is a simplified form of the previous one; we use a least-square approach for generating a local characteristic direction that minimizes the source term coupling the equations. Both solutions are computed so as to ensure a smooth transition between the local direction of convection given by the flow, and directions given by the gradients of the solution. In that way, we are able to prove that, if the hyperbolic system of equations is not fully diagonalized, then this results in a decrease in range of the stability condition (CFL condition). The source term that couples the equations is incriminated in this phenomenon. Consequently, provided that the time-step is conveniently reduced according to the magnitude of the source term, we can prove – for the linearized form of the equations – the monotonicity properties of the scheme. Numerical tests on classical problems of gas dynamics help to point out the advantages but also some drawbacks of our approach. Upon using such an algorithm, the numerical scheme so generated is a fourth-order accurate scheme (in time and space), which is formulated in a finite-volume/MUSCL least-square approach to deal with multi-dimensional problems on non-Cartesian meshes; this scheme is referred to as the “LSQP4” scheme (Least-SquarePositive-fourth-order). The logical structure of this paper is the following one: Section 2 is devoted to the derivation of the algorithm for discretizing scalar conservation laws over non Cartesian meshes. Then, we do numerical tests to check the main features of the resulting method and to understand the behavior of the corrections we introduced. Section 3 starts from theoretical considerations for diagonalizing a non linear hyperbolic system of PDEs. From these considerations, extension of the scalar work is then

87

investigated and the practical solutions (approximate diagonalizations) that promote this extension, are derived. We end this section by doing numerical checks and comparisons between the alternatives we proposed. This work ends with a summary of the salient features of the LSQP4 scheme. Finally, perspectives about the possible improvements and extensions of the numerical model, are drawn. Some theoretical results concerning the diagonalization approach of the Euler equations of gas-dynamics, are given in Appendix A, whereas Appendix B gives details about the error made in linearizing a non-linear scalar conservation law. 2. Discretization of the scalar problem We consider the following non-linear scalar conservation law, defined over the computation domain, D:

− →→ → ut + ∇ . F (u ) = 0,∀ x ∈ D

(1)

or, alternatively, in the convective form: →

→

ut + λ (u ).∇ u = 0

→  → λ (u ) ≡ ∂ F (u )/∂ u

(2)

To model some features of the Euler equations, the generalized →

t

flux, F (u ) ≡ ( f (u ), g(u ) ) , is supposed to have convexity properties, (f  (u) ≥ 0, g  (u) ≥ 0). If we integrate the conservative form, (1), over the discrete cell, Ii , we get the following semi-discrete algebraic form:



|Ii | ×

du¯ dt



=−

 

  → ωq × f˜ uLq , uRq ; n × Sij

(3)

j∈Ni q∈Qij

i

with the following rth order approximation (r ≥ p), upon using a Gauss quadrature over the cell face, Sij :



 →     → → → ωq × f˜ uLq , uRq ; n × Sij = F u x .nij dS + O (hr )

(4)

Sij

q∈Qij

where Qij represents the set of quadrature points over Sij and Ni is the set of points directly linked to i. In a 2D context, we use a 2-point Gaussian integration formula for computing (4). This formula is exact for a polynomial up to degree 3 (r = 4), if we take: ωq = 12 , ∀q, then q1 = w × A + (1 − w ) × √

B and q2 = w × B + (1 − w ) × A, with w ≡ 12 (1 + 33 ) (see Fig. 1). Therefore, the resulting discretization (3) is at best fourth-order in space.  L R uq , uq represent the local interpolated values at quadrature points over the cell face Sij and at both sides of this face (see Fig. 1). In all what follows, we suppose that (uLq , uRq ) are interpolated →

→

values at x = xq , from, respectively, the average values (u¯ ni , u¯ nj ), this with an error O (h p ) (p>1). (uLq , uRq ) are obtained from a high-order least-square interpolation procedure that will be explained in a following section. Now, upon using an approximate Riemann solver, we must pre→ cise the way of computing the approximate flux, f˜(uL , uR ; n ). q

q

To this aim, we select the classical HLLE Riemann solver [4], of which we briefly recall the main monotonicity properties in what follows. 2.1. Monotonicity properties of the HLLE scheme. For the considerations that will follow in the next sub-section, we re-interpret the HLLE scheme as a flux-difference splitting

88

G. Capdeville / Computers and Fluids 144 (2017) 86–116

The quantities (λnL , λnR ) are now defined by extending the definitions (7):

⎧    →  → ⎪ L ¯ ⎪ ⎨λnL ≡ min F u .nij , λn , 0    →  → ⎪ R ⎪ ⎩λnR ≡ max F u .nij , λn , 0

(12)

Summarizing, if we model problem (1) by Eqs. (3), (9), (10)–(12), then, provided that an explicit CFL-like condition is L ensuredR for the time-step, the resulting first-order scheme u ≡ u¯ i , u ≡ u¯ j is a monotonicity-preserving scheme, with all the stability properties derived therefrom [14]. In other terms, the discrete solution u¯ n+1 is a convex combination of all the quantities, u¯ nj , j ∈ Ni : this result was early quoted as a “positive scheme” by Spekreijse [10]. In the following sub-section, we derive the conditions that permit to extend this result to a high-order accuracy. 2.2. Monotonicity-preservation conditions 



Fig. 1. : Metrics and interpolated values uLq , uRq at quadrature points (q1 , q2 ) for a finite-volume algorithm.

From Eq. (3) and definition (9), we can write, now, the following result that discretizes (1):

 |Ii | ×

scheme, according to the one-dimensional formula:









 R



f˜ uL , uR = f uL + aL × u − uL = f u

R

− aR × u − u

where the mean quantity, condition” (see [4]):

u ,

 R





(5)

is calculated from the “conservativity

  L

aL ≡ min a u , a¯ , 0

  R

where we defined the new quantity, λqL , by using the Rankine– Hugoniot relationship, linking the states typified by uLq and u¯ i :

 →

→

(6)

a¯ ≡





(7)

( ) ( ) , ifuR −uL   ε ε ≡ O10−12   R

f uR − f uL uR −uL

f u , else



(8)



then, f˜ uL , uR , computed according to Eqs. (5), (7) and (8) is a Lipschitz continuous function of is arguments [4]. These definitions and properties constitute the main features of the so-called “HLLE Riemann solver” derived in [4] and, for these reasons, we select this solver as a building block in our model. Now, it is left to extend these properties to multi-dimensional computations. This is simply done by using the characteristic di→

rections of the mesh, nij , in the following way:



 →

f˜ uL , uR ; n



→



→

R L ≡ F uL .nij + λ− n × u −u



(9)

with the following definition obtained by inserting Eq. (6) into formula (5):



λnR − λ¯ n λ ≡ λnL × λnR − λnL



− n

and:

λn ≡

→ R F u

→  L

−F u

uR − uL



→

(10)

→

.nij (11)



(14)

To derive the monotonicity conditions, Eq. (13) is re-written as follows:

|Ii | ×

du¯ dt



=−

 





−λ

− qL

uLq − u¯ i

−

λ



u¯ j − u¯ i

 − nq



ωq × λ+qL

j∈Ni q∈Qij

i

The set of conditions (7) for the characteristic velocities is the one determined by Einfeldt [4], in the context of Euler equations. Moreover, if we practically define the mean characteristic velocity, a¯ , as:



j∈Ni q∈Qij

i

F uLq . n = F (u¯ i ). n + λqL uLq − u¯ i



aR ≡ max a u , a¯ , 0

    ωq × λqL uLq − u¯ i + λ−nq uRq − uLq × Sij

 

=−

→

 L

f u −f u aR u R − aL u L u ≡ − aR − aL aR − aL





(13)





with:

du¯ dt

uRq − uLq u¯ j − u¯ i



uLq − u¯ i u¯ i − u¯ k

× u¯ i − u¯ j







× (u¯ i − u¯ k )



× u¯ i − u¯ j





× Sij

(15)

where k is an index such that: k ∈ Ni , k = j. This latter result can be written, more synthetically, as:



|Ii | ×

du¯ dt



=− i





Cij+ × u¯ i − u¯ j



(16)

j∈Ni

Therefore, if the following three sufficient conditions are verified:

⎧ R  uq − uLq ⎪ ⎪ ≥ 0, ∀q ∈ Qij , j ∈ Ni : condition 1 ⎪ ⎪ u j − ui ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ uL − u  i q ≥ 0, ∀q ∈ Qij , j ∈ Ni : condition 2 u − u ⎪ j i ⎪ ⎪ ⎪ ⎪ ⎪ ⎪  L  ⎪ ⎪ uq − ui ⎪ ⎪ ≥ 0, ∀q ∈ Qik , k ∈ Ni , k = j : condition 3 ⎩

(17)

ui − uk

then, we have the following property: Cij+ ≥ 0, ∀ j ∈ Ni . Those conditions were already obtained by Barth in [1]. In our paper, the difference with his work appears with the way of introducing (17) into Eq. (16).

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Now, by integrating in time Eq. (13) over the range t ≡ tn+1 − tn , we get the forward Euler discretization of (1):

 u¯ ni +1 = u¯ ni ×

1−





t t × C+ + × C + × u¯ nj |Ii | j∈N ij |Ii | j∈N ij i



t  + 1− × C |Ii | j∈N ij

(18)

i

From this latter result, we can  t verifies the CFL− like condition:

deduce

that,

provided



 0, ∀i

(19)

then, the scheme typified by (16) can be termed as a “positive” scheme with the following monotonicity property:



u¯ ni +1 ∈ [mi , Mi ], ∀i with : Mi ≡ max u¯ ni , u j j∈Ni

n





and mi ≡ min u¯ ni , u j

n



j∈Ni

(20) with as ultimate consequences:

i

L2



 0,∀ j ∈ Ni

u¯ j − u¯ i

q∈Qij

 (1 + O (h ) ) × u¯ n L . 2

(21)

This is also a sufficient condition that has the advantage to be simpler and cheaper to compute than the original condition; for this main reason, we retain this assumption in what follows. Now, introducing the limiter, ϕ ij , which is designed to enforce condition (21), we can define the following corrected interpolated value, u˜Lq , such that: →

u˜Lq ≡ u¯ i + ϕij × ui (xq ; u¯ ) ≡ u¯ i + u1i (xq ; u¯ )

→ 



(22)



with: ui xq ; u¯ ≡ uLq − u¯ i . Therefore, we want to devise ϕ ij such that the following condition is locally ensured, inside the cell Ii :



uRq − u˜Lq



 0, ∀q ∈ Qij , j ∈ Ni

u¯ j − u¯ i

(23)

→

of ui (xq ; u¯ ), we get the following inequality:

1. u¯ ni +1 ∈ [m, M], ∀i with : M ≡ max(u¯ ni ) and m ≡ min u¯ i



uRq − uLq

Then, by using this condition with Eq. (22) and the definition

 n





ωq ×

→

i

2. u¯ n+1 

If we suppose that: λnq ≡ Cte over Sij , then, the first condition may be weakened by the following one:





89

  ϕij × uLq −u¯ i  uRq −u¯ i

i

(24)

This condition is verified if we define ϕ ij , over Sij , as follows:

The first property tells us that a global maximum principle is ensured, at each time-step, for un+1 ; the second property is more accurate since it provides an asymptotic stability in L2 -norm of the numerical solution. Note that this latter property holds only if Cij+ is a Lipschitz continuous function of u¯ n (see [7], for a detailed demonstration): this is, of course, no more true in regions where discontinuities appear.

   R  uq −u¯ i   , 1 j ∈ Ni ϕij ≡ Min  L q∈Qij uq −u¯ i 

(25)

or, by favoring condition (21):

⎧   ⎫ R ⎨ q∈Q ωq uq −u¯ i  ⎬ ij , 1 j ∈ Ni ϕij ≡ Min   ωq uLq −u¯ i  ⎭ ⎩ q∈Q

(26)

ij

Remark 1. Conditions (17) do not act on the same level. Indeed, when the solution remains smooth, the first condition imply an O (h p−1 ) quantity whereas condition 2 and 3 use O (1 ) quantities: clearly, in that case (smooth solution and p > 1), it becomes useless to activate condition 1. Remark 2. The only necessary condition to ensure positivity of the scheme, is: Cij+ ≥ 0, ∀ j ∈ Ni . Conditions (17) are only sufficient conditions; other different conditions might be derived: this point will be exploited in the following subsection. Remark 3. The maximum-principle comes from the positivity condition: Cij+ ≥ 0, ∀ j ∈ Ni . However, this latter condition does not tell us about the bounds of the locally interpolated quantities  Lanything uq , uRq . This means that those interpolated quantities may be locally out of the range [m, M] if no specific correction is made to remedy that. This point is particularly crucial for Riemann solvers, such as HLLE or equivalent, that use intermediate states to compute the numerical flux. These states should maintain a maximumprinciple so that the resulting numerical flux is a convex combination of “physical” quantities (see [4], for example) Now, the main difficulty is to practically derive a numerical procedure from conditions (17) in order to ensure the monotonicity preservation of (18).

Remark 4. The local assumption λnq ≡ Cte over Sij that enables to derive condition (21), is not so crucial. Indeed, in deriving formula (13), we only need to compute the variations of the solution in →

the local normal direction to Sij given by n ij . Therefore, we can logically assume that variations in transverse directions, along Sij , do not play a significant role into the positivity of Eq. (18). Now, let us have a glimpse at the numerical behavior of ϕ ij . Let us suppose that, for a given point jo ∈ Ni , linked to i, we have the following behavior for uLq :

 q∈Qijo



ωq ×

uRq − uLq u¯ jo − u¯ i

<0

and suppose that this latter negative quantity is less than 1, in absolute value. Then, definition (26) gives, for ϕijo :

   R    q∈Q ωq uq −u¯ i   ijo  ϕijo ≡    L −u ¯ ω u  q q i  q∈Qijo 

From this latter result, two possibilities may arise: 1. If the solution remains smooth in the region centered around i, then, according to the order of accuracy of the interpolation   polynomial, we have: q ωq uLq − q ωq uRq = O (h p ). p −1 Therefore ϕijo = 1 + O (h ), hence we have: →

→

ϕijo × ui (xq ; u¯ ) = ui (xq ; u¯ ) + O (h p ).

2.3. Derivation of practical conditions for positivity

→

Consequently, we can deduce from the definitions of ui (xq ; u¯ ), uRq −uLq

Let us begin with the first condition from (17): ( u¯

0, ∀q ∈ Qij , j ∈ Ni (see Fig. 1, for the notations).



j −u¯ i

)

→ u ( xq , tn )

uLq and u˜Lq , that: u˜Lq = + O (h p ). In other words, in such a case, there is no local loss of accuracy when ϕ ij is activated.

90

G. Capdeville / Computers and Fluids 144 (2017) 86–116

2. Now, if there exists a discontinuity between points i and jo , we   get henceforth the following result: q ωq uLq − q ωq uRq = O (1 ) →

We obtain, consequently, the property: u˜Lq = u(xq , tn ) + O (h ) and the scheme becomes locally first-order accurate. Therefore, the first part of the monotonicity-preserving procedure gives us a corrected interpolated value, u˜Lq , such that:

→ 

→ 

u˜Lq = u¯ i + u1i xq ; u¯ ≡ u¯ i + ϕij × ui xq ; u¯

∀q ∈ Qij , j ∈ Ni

Now, let us adapt the second condition of (17) to the interpouL −u¯

lation procedure. We recall that this condition writes as: ( u¯q −u¯ i )  j

i

0, ∀q ∈ Qij , j ∈ Ni . In words, this latter condition simply stipulates that the interpolation procedure that produces uLq must follow the monotonicity of the mean solution, characterized by u¯ i . With the same assumption as before concerning the behavior of λqL , the simplest way for introducing this condition is to impose that, ∀q ∈ Qij :



i f sign



ωq 

u1i

→ 





→ 



= sign u j − ui ⇒ u1i xq ; u = 0,

xq ; u

q∈Qij

(27) Of course, if we directly apply such a condition, the local accuracy of the numerical solution is lost, whatever the smoothness of this solution. To attenuate this detrimental behavior, as much as possible, while maintaining the monotonicity properties of the scheme, we introduce a “positivity monitor”, which is designed to control the local positivity of the solution. To this end, this latter quantity, named α i , is formally defined over the discrete cell, Ii , and for each time-step, as follows:

 

αin+1 ≡

u¯ ni +1 −min u¯ nj

 

j∈Ni

 

(28)

max u¯ nj − min u¯ nj

If this quantity is out of the range [0, 1], then there is a local loss of positivity for the discretization (18) and a corrective procedure must be employed; otherwise, the scheme remains locally positive and no corrections are needed. Definition (28) is only formal. Indeed, the interpolation procedure is computed at t = tn ; so, at this discrete time, there is no way to know the quantity u¯ ni +1 that depends on this latter procedure. For this reason, we practically compute α i , at t = tn+1 , by the following third-order time extrapolation procedure:



∂α ∂t

n

+

t

i

2

2



×

∂ α ∂t2 2

n



+ O t 3



(29)

i

Upon using a backward discretization of the time derivatives, we get the following final approximation for α i , at t = tn+1 :

αin+1

of the scheme, is not entirely evacuated and there may still exist some possible loss of accuracy in regions of smoothness of the solution, when this latter correction is activated. Nevertheless, the advantage in using α i is that, even though the quantity that appears in condition 2 is negative, the related correction (27) may not be activated, provided that the scheme remains locally positive: in this sense, we can consider that the global accuracy of the scheme is improved by this method. The proper functioning of this procedure will be particularly checked in the numerical tests that follow.

j∈Ni

j∈Ni

αin+1 = αin + t ×

Fig. 2. One-dimensional interpolation procedure at a regular extremum.

1 n 5αi −4αin−1 + αin−2 2

(30)

Practically, α i is used to introduce the correction displayed by (27) into (18), by modifying the interpolation procedure: if α i lies outside the range [0, 1], then the modification given by (27) is activated and so the second positivity condition from (17) is verified; otherwise, the scheme remains locally positive and, consequently, there is no need to modify the interpolation procedure. Remark 5. Introducing α i into the interpolation procedure is a way to ensure that some positivity conditions are applied, only when it becomes necessary. However, this modification is done independently of the smoothness of the solution. This means that the harmful behavior of correction (27) for the formal accuracy

Now, it is left to introduce the condition 3 of positivity into the interpolation procedure. In order to better understand the meaning of this condition, let us examine the one-dimensional problem. In this simplified framework, conditions 2 and 3, respectively write as follows, at the interpolated point x = xi+1/2 :

uLi+1/2 − u¯ i u¯ i+1 − u¯ i

 0 and

uLi+1/2 − u¯ i u¯ i − u¯ i−1

0

(31)

To illustrate the role of these conditions in that context, let us consider the usual problematic case of a regular extremum, pictured in Fig. 2, which follows. Clearly, condition 2 in (17), is verified at x = xi+1/2 . However, if we look at condition 3, it appears immediately that this latter condition is no more verified: a correction may become necessary in order to avoid a possible local loss of positivity. To this end, we use exactly the same corrective procedure as the one devised for applying condition 2, ∀q ∈ Qij :



i f sign



ωq 

u1i

→  xq ; u



→ 

= sign (ui − uk ) ⇒ u1i xq ; u = 0,

q∈Qij

(32) Once again, this correction is activated only if αin+1 lies outside the range [0, 1]. As one can see it from Fig. 2, condition 3 of positivity occurs at regular extrema: it is this condition that is mainly responsible for the classical “clipping” phenomenon that appears at regular extrema of the numerical solution when applying positivity

G. Capdeville / Computers and Fluids 144 (2017) 86–116

conditions. By using α i , our intent is to reduce the scope of this local attenuation, over the global accuracy of the scheme. Lastly, taking into consideration Remark 3, we introduce an additional limiter, ψ i , into the already modified interpolation procedure. This limiter is aimed at imposing a maximum-principle on the interpolation procedure so that the resulting interpolated values, uLq , q ∈ Qi , do not lie outside a specified range that owns a  physical meaning (Qi ≡ {quadrature points for Ii } ≡ j∈N Qij ).

91

2.4. Numerical algorithm for monotonicity preservation In this subsection, we describe the numerical algorithm we implemented for computing the corrected pointwise values (u˜Lq , u˜Rq ) that appear in the scalar Eq. (3). Those values are designed so as to generate a monotonicity preserving scheme for discretizing the scalar conservation law (1). Algorithm for a scalar problem.

i

Specifically, we want that:

• Step 1: Computation of condition 1.

→

m u¯ i + ψi × u1i (xq ; u¯ )  M, ∀q ∈ Qi

(33)

Then, taking into account all the possibilities, according to the →

sign of u1i (xq ; u¯ ), it is sufficient to define ψ i , as follows:

 ⎧  →  →  ⎪   ¯ M− u ⎪ i→  if maxu1 xq ; u¯ > 0 and minu1 xq ; u¯ > 0 ⎪  ⎪ i i ⎪ q∈Qi q∈Qi u1i xq ;u¯  ⎪ max q∈Qi ⎪ ⎪ ⎪ ⎪ ⎪   ⎪  ⎨ →  →   m−u¯ i   ψi ≡  if minu1i xq ; u¯ < 0 and maxu1i xq ; u¯ < 0  → q∈Qi u1i xq ;u¯  q∈Qi  min ⎪ ⎪ q∈Qi ⎪ ⎪ ⎪ ⎪    

 ⎪ ⎪     ⎪ ⎪     ¯ M− u¯ i m − u ⎪ i  ,   →  , 1 , otherwise ⎪  → ⎩  maxu1 xq ;u¯   minu1 xq ;u¯  i i q∈Q q∈Q i

i

(34) Thus, by this definition, it is easy to check that the quantity, →

ψi × u1i (xq ; u¯ ), lies inside the range [m − u¯ i , M − u¯ i ], ∀q ∈ Qi . The bounds (m, M) may be defined, for example, from the →initial solu→ → tion, u0 ( x ), of (1): M ≡ max u0 ( x ) and m ≡ min u0 x . → → x ∈D

lution is supposed to be regular, we can assert that: u˜Li+1/2 = M +

O (h p ). Therefore, using the fact that u1i (xi+1/2 ; u¯ ) ≡ u˜Li+1/2 −u¯ i , we have the following result:

= 1 + O (h p−1 )



u˜Li+1/2 = u xi+1/2 + O (h p ) In words, this means that the limiter, ψ i , does not destroy the accuracy of the polynomial interpolation in regions of smoothness of the numerical solution. Remark 6. It is important to note that the limiter, ψ i , acts upon → the already limited increment, u1i (xq ; u¯ ), → limiting the initial increment, ui (xq ; u¯ ) ≡

• Step 2 : Computation of the “positivity monitor”, αin+1 For all i, do: ! 1. αin =

u¯ ni −min u¯ nj −1 j∈N

!i

max u¯ nj −1 −min u¯ nj −1

!

j∈Ni

j∈Ni

2. if n > 2 then  compute: αin+1 = 12 5αin −4αin−1 + αin−2

Else αin+1 ≡ 1 • Step 3: Computations of conditions 2 and 3, according to the value of αin+1 For all i, do: If αin+1 ∈ / [0, 1]then 

1. s1 j ≡ sign 

s2 j ≡ sign

 q∈Qij

 q∈Qij

 →  ωq u1i xq ;u¯ ×sign(u¯ j −u¯ i ), j ∈ Ni

 →  ωq u1i xq ;u¯ ×sign(u¯ i −u¯ k ), k ∈ Ni , k = j

2. if s1j < 0 or s2 j < 0 ⇒ u1i

→ 

xq ;u¯ =0,

j ∈ Ni , q ∈ Qij

αni +1

Endif • Step 4: Computation of ψ i from equation (34) • Step 5: Computation of u˜Lq : For all i, do: →  u˜Lq ≡ u¯ i + ψi × u1i xq ;u¯ , j ∈ Ni , q ∈ Qij



Consequently,  we can write that:  ψi × u1i xi+1/2 ; u¯ = u1i xi+1/2 ; u¯ + O (h p ) and we can conclude that:



⎧   ⎫ ⎨ q∈Q ωq uRq −u¯ i  ⎬ ij , 1 j ∈ Ni 1. ϕij ≡ Min   ⎩ q∈Q ωq uLq −u¯ i  ⎭  →  ij →   1 2. ui xq ; u¯ ≡ ϕij × ui xq ; u¯ ≡ ϕij × uLq −u¯ i , j ∈ Ni , q ∈ Qij

x ∈D

To understand the behavior of ψ i , in regions of smoothness, let us return to the one-dimensional framework. For example, let us suppose that, at some cell interface, x = xi+1/2 , the interpolated value is such that: u˜Li+1/2 > M. Since the so-

M − u¯ i L u˜i+1/2 −u¯ i

→

For all i, compute u1i (xq ; u¯ ), q ∈ Qi , as follows:

which was obtained by uLq -u¯ i , by the local lim-

iter, ϕ ij . Moreover, the order in which these both limiters are used, is not practically important: numerical tests we managed proved that the introduction of ψ i into the interpolation procedure does not destroy the modifications imposed by conditions (17). Conversely, if we apply conditions (17) after having imposed the maximumprinciple condition based upon ψ i , the positivity of the resulting scheme is unaltered. A this point, it is important to summarize the whole procedure that ensure the monotonicity of discretization (18) and to describe the numerical algorithm we employed to enforce conditions (17).



Note. For the computation of the numerical flux, f˜ u˜Lq , u˜Rq , at the cell interface, Sij , the quantity, u˜Rq , is computed from the jth adjacent cell, j ∈ Ni : this is simply the quantity, u˜Lq , which was computed in the discrete cell Ij . Now, in the subsection that follows, we briefly present the  →way  of computing the unlimited increment over the cell Ii : ui xq ;u¯ ≡ uLq −u¯ i ,

j ∈ Ni , q ∈ Qij .

2.5. Least-square interpolation procedure To set the base interpolation procedure, all we need is to com→ → pute the local increment ui (xq ;u¯ )≡ui (xq ;u )−u¯ i ≡uLq −u¯ i , j ∈ Ni , q ∈ Qij ; later, this increment will be limited according to the algorithm presented just above. For the sake of simplicity, we just present in what follows the bi-dimensional formulation of the polynomial interpolation. We begin by defining the following polynomial of degree ( p−1 ), over the discrete cell, Ii :

→ 

ui x ; u ≡ u¯ i +

 p−1 k k=1

1 l=0 hk k!

  ×

k l

"



× (x−xi )l × (y−yi )k−l − xl yk−l →

× Dl,k−l , ∀ x ∈ Ii

# i

(35)

92

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Fig. 4. Quadrature points, qj , and points, j, linked to i, across the cell interfaces, Sij , j ∈ Ni ≡ {1, 3, 5, 7}. Fig. 3. 2D discretization stencil, Ni ,for fourth-order accuracy.

Practically, the unknown coefficients of this interpolation, namely: Dl,k−l ≡ hk × Dl,k−l = (

∂ k u ) + O h p , are computed ac( ) ∂ xl ∂ yk−l i

cording to the following conservative conditions:

u¯ j ≡

1

 × I j 

 

→ 



 

→

u j x ; u¯ × d x , j ∈ Ni

Ij



(36)



The term xl yk−l is a metric term introduced so that Eq. i (35) verify condition (36), for each discrete cell, Ii .





In other words, we have the following definition for xl yk−l :



xl yk−l

i

≡

1 × |Ii |

 Ii

 →   (x−xi )l × (y−yi )k−l × d x 

i

(37)

This being said, if we impose conditions (36), for all j ∈ Ni , Ni ≡ {points of the discretization stencil}, (see Fig. 3), to formula (35), we get the following system, for the unknowns, Dl,k−l :

⎧ k  p−1 ⎪  ⎪ ⎪ al,k−l × Dl,k−l , j ∈ Ni ⎨ u¯ j −u¯ i = k=1 l=0   " ⎪   # ⎪ k 1 ⎪ × xl yk−l − xl yk−l ⎩al,k−l ≡ hk k! × ij i

(38)

l

xl yk−l



ij

≡

1

 × I j 



Ij

xl yk−l



"

ij

=

k−1 l  p=0

q=0



k−l p

 

× (x−xi ) × (y−yi ) − q

p

l q

Ni ≡ {1, 3, 5, 7} and Qi ≡ {q1 , q2 , .., q8 } For example, along the cell-face, Si1 , we have: Qi1 ≡ {q1 , q2 }. Therefore, the system (38) can be written, equivalently, as:

 ⎧ AX = B; with : dimA = m, n ≡ ⎪ ⎪ ⎪ ⎨  A ≡ al,k−l , k ∈ {1, p}, l ∈ {0, k} ⎪ ⎪ ⎪  t  t ⎩

xl −q yk−l−p

# j

, j ∈ Ni

p( p + 1 ) 2

u¯ j −u¯ i

p( p+1 ) 2



(41)

, j ∈ Ni

Practically, this is accomplished by using the “matrix normal method”, which amounts to solving the following symmetric definite positive system:

 (39)

Consequently, in order to get a determined system for the Dl,k−l ’s, we need to impose the following sufficient condition:

card Ni ≡ m 

In addition, according to Fig. 4, below, we can precise the following notations:

J (X ) ≡ | AX −B|2

 

×

(40)

This is an over-determined system (m > n) of which the solution may be computed by minimizing the function:

 →   (x−xi )l × (y−yi )k−l × d x , j ∈ Ni

Practically, this latter quantity is computed according to the following formula:



Ni ≡ N0 = {1, 2, .., 12} ⇒ card Ni ≡ m = 12

X ≡ Dl,k−l ; B ≡

with the metric term:



For example, if we select p = 4 (fourth-order accuracy), then we must have m ≥ 10. For all the bi-dimensional computations that follow, we made this particular choice. For that purpose, we selected the stencil pictured by Fig. 3, below, which is centered around the discrete point “0”, labelled “i”, throughout this paper. In a cell-center finite-volume approach, each discrete point, j, represents the barycenter of the jth cell, Ij . With such a choice, we have, obviously, the following identity:

SX = V ; S ∈ Rn,n

S ≡ t A.A;V ≡ t A.B

(42)

The numerical result, X, represents, for each discrete cell, Ii , the minimum in L2 -norm of J(X). To end with the features of this model, it is left to precise the time-integration method we selected in order to integrate the semi-discrete form, (3).

G. Capdeville / Computers and Fluids 144 (2017) 86–116 Table 1 Linear advection  sin π2 (x + y ) .

2.6. Preserving positivity time-integration method In order to preserve the positivity properties of the scheme in its high-order time-integration version, we made the choice of the fourth-order, five stages, SSPRK(5,4) scheme of Gottlieb [5]. Indeed, this scheme is designed as a convex combination of forward Euler discretizations such as Eq. (18). Therefore, if the discretization (18) is positive, so is the fourthorder time integration using the SSPRK(5,4) scheme. By writing, formally, Eq. (13) as:



du¯ dt

 ≡ F (u¯ (t ) )

93

ut + ux + uy = 0,

test:

u0 (x, y ) =

N

L1 -error

Order

L∞ -error

Order

10 20 40 80 160 320

1.71 × 10−3 1.08 × 10−4 6.78 × 10−6 4.25 × 10−7 2.65 × 10−8 1.66 × 10−9

– 4 4 4 4 4

2.65 × 10−3 1.70 × 10−4 1.08 × 10−5 6.65 × 10−7 4.17 × 10−8 2.61 × 10−9

– 4 4 4 4 4

Table 2 Linear advection test: ut + ux + uy = 0, u0 (x, y ) =  sin π2 (x + y ) conditions 2 and 3 activated, whatever α i .

(43)

i

integration of this equation by the SSPRK(5,4) scheme, gives the following algorithm: Algorithm SSPRK(5,4), [5]. 1 ui( ) ≡ u¯ ni + 0.391752226571890 × t × F (u¯ (tn ) )

N

L1 -error

Order

L∞ -error

Order

10 20 40 80 160 320

8.90 × 10−3 2.33 × 10−3 4.64 × 10−4 1.45 × 10−4 4.18 × 10−5 9.94 × 10−6

– 1 2.2 1.7 2.1 2.1

1.47 × 10−2 1.25 × 10−2 4.23 × 10−3 2.20 × 10−3 8.94 × 10−4 3.42 × 10−4

– <1 1.5 1 1.3 1.3

2 ui( ) ≡ 0.444370493651235 × u¯ ni + 0.555629506348765

 1 ×ui( ) + 0.368410593050371 × F u(1 )

(3 )

ui

≡ 0.620101851488403 ×

u¯ ni

+ 0.379898148511597 × ui

+0.251891774271694 × F (4 )

ui

≡ 0.178079954393132 ×

u¯ ni



u (2 )





u (3 )

i

(3 )



2 3 ≡ 0.517231671970585 × ui( ) + 0.096059710526147 × ui( )

+0.063692468666290 × F



u (3 )



+0.386708617503269 × u(4 ) + 0.226007483236906



×F u(4 )

A Cartesian grid (h ≡ x = y ) is used for all the computations of this section. The CFL number is defined as:

CFL ≡ max

+ 0.821920045606868 × ui

+0.544974750228521 × F u¯ ni +1

(2 )



The advantage in using this algorithm is twofold: firstly, the resulting scheme is henceforth fourth-order accurate in time and space and, moreover, it remains monotonicitypreserving; secondly, this time-integration method is more efficient than the popular SSPRK(3,3) method (see [5]), with an effective CFL number (≡ CFL/m, m:number of computations of F (u¯ (t ) ) required per time step) of 0.377 instead of 1/3 for SSPRK(3,3). Now, we have designed a numerical model of which the essential features are the following ones: 1. A fourth-order least-square and monotonicity-preserving polynomial reconstruction of the discrete variable. 2. A fourth-order monotonicity-preserving Runge–Kutta time integration method. The result is the “Least-Square Positive Fourth-order accurate scheme” (LSQP4). In the subsection that follows, we study the essential properties of this model by some numerical experiments.

2.7. Numerical tests We use for all the computations that follow, the LSQP4 scheme with positivity conditions (17) and the SSPRK(5,4) time-integration algorithm. The positivity constraints 2 and 3 from (17) are conditionally activated, according to the magnitude of the “positivity monitor”, α i .

   n     n   f u¯ i , g u¯ i  × t/h

Test 1: Bi-dimensional linear advection. We begin this series of numerical tests with an accuracy study. To this aim, we compute the linear form of Eq. (1) :

ut + ux + uy = 0∀(x, y ) ∈ [−2, 2]2 with the following initial solution:

u0 (x, y ) ≡ u(x, y, t = 0 ) = sin

π 2

 (x + y )

Computations are run up to the dimensionless time, T = 0.20, with periodic boundary conditions. Table 1, displays the numerical results we obtained by using the LSQP4 scheme with a CFL number of 0.25. The Cartesian mesh is gradually refined from N = 10 to N = 320 grid points, in each space direction. As can be seen, the accuracy expected is reached, both in L1 and L∞ norms, even on the coarsest meshes. The LSQP4 scheme is uniformly accurate and α i remains everywhere and every time in the “positivity range” [0, 1]. It would be interesting to see what would happen if conditions 2 and 3, from (17), were activated, whatever the value taken by α i . Results obtained by such a choice, are displayed in Table 2. As predicted by the theoretical analysis, there is a local loss of accuracy that comes from errors generated at extrema of the solution. This error (L∞ error) is accumulated during the time iterations and modifies the overall accuracy of the scheme by generating a global error (L1 error). From results of Table 2, two points should be underlined: firstly, when the solution remains smooth, the conditions 2 and 3 of positivity are not necessary to get a positive scheme; secondly, even though α i is approximately computed, this quantity permits to check, efficiently, the positivity of the LSQP4 scheme. However, these first findings should be confirmed before drawing teachings about the overall behavior of positivity conditions (17). Hence, we pursue this accuracy analysis, with the more irregular initial solution: 4

u0 (x, y ) = sin

π 2

 (x + y )

94

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Fig. 5. Two-dimensional linear advection of a discontinuity at T = 2. a) h = 1/10. b) h = 1/40.

Table 3 Linear advection  4 sin π2 (x + y ) .

test:

ut + ux + uy = 0,

u0 (x, y ) =

N

L1 -error

Order

L∞ -error

Order

20 40 80 160 320

8.0 × 10−3 8.24 × 10−4 5.38 × 10−5 3.39 × 10−6 2.13 × 10−7

– 3 4 4 4

2.10 × 10−2 2.26 × 10−3 9.25 × 10−5 5.81 × 10−6 3.63 × 10−7

– 3.10 4.6 4 4

All the remaining parameters are left unchanged. Then, we get the following results, typified by Table 3: Once again, the numerical scheme reaches the theoretical accuracy with, however, a higher error level than before; if we discard the use of α i , then we get the same tendencies as previously: an order of accuracy of 2, in L1 error, and an order slightly greater than 1, in L∞ error. Now, we study the linear advection of an initial discontinuous solution. This latter solution is as follows:



u0 (x, y ) =

1 if (x, y ) ∈ S 0 otherwise

∀(x, y ) ∈ [−1, 1]2

 √ √  where S ≡ (x, y )/|x − y| < 1/ 2, |x + y| < 1/ 2 represents an unit square centered at the origin and rotated by an angle of π /4. For this problem, we select two Cartesian mesh sizes: h = 1/10 and h = 1/40. Computations are run up to T = 2 (one period in time) and the optimum CFL number is obtained for the value 0.25. To be as accurate as possible, both the maximum and minimum values of the numerical solution and the error in L1 -norm, are computed, every time. The first results are displayed by Fig. 5 for an horizontal cut along the line y = 0. In addition, we plot, on the same graph, the profiles for α i (dashedline), at the final time, T = 2, of the simulation. As one can see it, the monotonicity of the numerical solution is ensured; this is true, even on the coarse mesh. The maximum principle is verified by the numerical solution, since its numerical values lie inside the range [umin = 0, umax =1]. Lastly, the “positivity monitor”, α i , is everywhere inside the range [0,1]: this means that the positivity conditions 2 and 3 are not activated during the computations. Now, it is interesting to study the behavior of the numerical solution when α i is not employed in the numerical procedure, the

remaining parameters of the computations being, moreover, unchanged. Doing so, we get the results pictured by Fig. 6. The damping behavior associated to the use of conditions 2 and 3 is obvious: the L1 -norm of the discretization error makes it possible to quantify this fact. Once again, these latter conditions do not appear to be necessary to preserve the monotonicity of the scheme. If we look, now, at the significance of positivity condition 1 in the monotonicity of the numerical solution, we get, then, the results displayed by Fig. 7. Clearly, the monotonicity of the numerical solution, is lost. In addition, the maximum principle is no more ensured because the maximum and minimum values of the solution fall outside the range [0,1]. Lastly, the “positivity monitor”, α i , is locally (just before and just after the discontinuity) outside the interval [0,1]: the loss of positivity of the scheme is confirmed. Therefore, condition 1 is necessary, for this test-case, so that the numerical scheme preserves the monotonicity of the exact solution. Pursuing this series of tests, we consider, henceforth, the following linear form for the scalar Eq. (1):

ut + yux − xuy = 0 with the Gaussian initial data: u0 (x, y ) = exp [ − 100 × [(x − 1/2 )2 + y2 ]] over the square domain [−1, 1]2 . This test-case is significant to evaluate the dissipative and dispersive properties of the LSQP4 scheme. Indeed, usually slight oscillations appear behind the convected hill, when using a highorder scheme. These oscillations come from phase errors generated by the scheme and illustrate a local loss of positivity of that scheme. The mesh size selected is h = 1/40, computations are run up to T = 2π (one period) and the optimized CFL number is set to 1.25. Firstly, we begin these computations with results displayed by Fig. 8, in order to illustrate the dispersive properties of the scheme, without any positivity correction. In other words, conditions (17) are deactivated for obtaining these first results. Although slight, the numerical oscillations behind the hill, clearly appear (Fig. 8b). In addition, α i is lower than zero (dashed line), at this location: there is a loss of positivity at the place where oscillations appear. Therefore, to remedy this problem, the solution is to enforce positivity conditions, (17).

G. Capdeville / Computers and Fluids 144 (2017) 86–116

95

Fig. 6. Two-dimensional linear advection of a discontinuity at T = 2. Activation of positivity conditions 2 and 3.

Fig. 7. Two-dimensional linear advection of a discontinuity at T = 2. Deactivation of positivity condition 1.

Fig. 8. Rotation of a Gaussian hill over one period (h = 1/40) Conditions (17) deactivated.

By doing that, we get the numerical results illustrated by the Fig. 9. A first remarkable result must be emphasized: the maximum value of the solution is left unchanged when passing from an uncorrected interpolation procedure (Fig. 8) to a corrected one (Fig. 9).

Indeed, we have umax = 0.904, for both results. Besides, the discretization error of the corrected interpolation (err−L1 = 5.74 × 10−4 ) is lower than the uncorrected one (err−L1 = 6.18 × 10−4 ). There are no more visible oscillations and α i remains everywhere into the “positivity range”: this means that conditions 2 and 3 remain deactivated, throughout the time evolution since umax is the same as the one obtained without any correction.

96

G. Capdeville / Computers and Fluids 144 (2017) 86–116 Table 4 2D Burgers equation. u0 (x, y ) = 0.30.

1 4

+

1 2





sin π2 (x + y ) . T =

N

L1 -error

Order

L∞ -error

Order

10 20 40 80 160 320

8.28 × 10−4 9.26 × 10−5 1.03 × 10−5 8.90 × 10−7 6.70 × 10−8 4.49 × 10−9

– 3.1 3.2 3.5 3.7 3.9

2.28 × 10−3 5.29 × 10−4 7.39 × 10−5 9.31 × 10−6 7.44 × 10−7 5.14 × 10−8

– 2.1 2.9 3 3.6 3.9

Table 5  2D Burgers equation. u0 (x, y ) = 14 + 12 sin π2 (x + y ) . T = 0.30 conditions 2 and 3 activated, whatever α i .

Fig. 9. Rotation of a Gaussian hill over one period (h = 1/40) Conditions (17) activated.

N

L1 -error

Order

L∞ -error

Order

10 20 40 80 160 320

4.95 × 10−3 1.76 × 10−3 3.87 × 10−4 7.31 × 10−5 1.43 × 10−5 2.86 × 10−6

– 1.5 2.0 2.4 2.4 2.3

1.95 × 10−2 9.38 × 10−3 3.25 × 10−3 1.35 × 10−3 4.56 × 10−4 1.58 × 10−4

– 1 1.5 1.6 1.6 1.5

case, the non-linear form of Eq. (1) writes as:



ut +

u2 2





+ x

u2 2



with the initial data:

u0 (x, y ) =

Fig. 10. Rotation of a Gaussian hill over one period (h = 1/40) Positivity conditions 2 and 3 activated.

Therefore, upon considering this unique benchmark, it appears that condition 1, from (17), is necessary, even when the solution remains smooth; clearly, this condition is a non-oscillatory condition that acts similarly to a high-frequency filter, without any damage for accuracy of the smooth part of the solution. Lastly, if we activate conditions 2 and 3 by removing the influence of α i , we get the result illustrated by Fig. 10 Once again, with a significant damping of the extremum (umax = 0.744), the dissipative nature of conditions 2 and 3, is obvious; as a consequence, the error in L1 -norm, is increased (err−L1 = 1.01 × 10−3 ). However, one can notice that α i remains, everywhere, in the “positivity range”. Test 2: Bi-dimensional non-linear advection. Now, we want to check the previous conclusions in a non-linear context. To this aim, we study the 2D Burgers equations. In such a

1 1 + sin 4 2

=0

∀(x, y ) ∈ [−2, 2]2

y

π 2

 (x + y )

The optimized CFL number is set to 0.50, every time, and computations are completed over Cartesian grids. The exact solution for this problem is obtained from the “method of characteristics”. Firstly, the L1 - and L∞ -norms of the error, are collected for the dimensionless time, T = 0.30, when the exact solution is still smooth. Table 4, below, synthesizes those first findings. The error convergence is slower than in the linear case; however, the theoretical order of accuracy is asymptotically reached over the finest meshes, both in L1 and L∞ norms. During the computations, the “positivity monitor”, α i , remains, everywhere, in the “positivity range”. Now, by discarding the influence of this latter quantity, we get the results typified by Table 5, below: As foreseen, those results show the dissipative influence of the positivity conditions 2 and 3 and are similar to results obtained in Table 2. Now, if we compute the numerical solution at T = 0.80, a moving shock develops and interacts with a rarefaction wave. The numerical solution is plotted in Fig. 11, for h = 1/20 and CFL = 0.50. As we can see it, the monotonicity of the solution is preserved between the shocks and these shocks are accurately captured without any oscillation. In addition, α i , remains, everywhere and every time, in the “positivity range”. For this result, the L1 -norm of the error is: 1.44 × 10−2 , and umax = 0.749, umin = −0.248. By activating the positivity conditions 2 and 3, we then get the results visualized by Fig. 12, just below: This time, conditions 2 and 3 do not deteriorate the solution; if we look at the L1 -norm of the error, we find: errorL = 1.41 × 1

10−2 . This slight improvement comes from the discontinuity that is slightly better described. By enforcing those two conditions, we get, then, umax = 0.747, umin = −0.248. Lastly, if we remove the positivity condition 1 from the set of conditions, (17), we get absolutely the same results as that obtained in Fig. 11: this means that, for this problem, the positivity

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Fig. 11. 2D Burgers equation. Solution at T = 0.80.

Fig. 12. 2D Burgers equation. Solution at T = 0.80. Positivity conditions 2 and 3 activated

condition 1 is not necessary to ensure the positivity of the numerical scheme. Up to now, we did not say any word about the influence of limiter, ψ i , defined by Eq. (34). Indeed, according to Eq. (15), this limiter is not mandatory to get a positive scheme; we introduced it so as to impose a maximum-principle into the interpolation procedure. For all the benchmarks we computed where the numerical solution remains smooth, and everywhere in the computation domain, this limiter has no influence. This is different when the solution becomes discontinuous. To illustrate this fact, let us compute the 2D Burgers equation without using limiter (34) but in using all the positivity conditions from (17) weighted by the use of α i . The result is given by Fig. 13, below:

97

Fig. 13. 2D Burgers equation. Solution at T = 0.80. Limiter (34) deactivated.

Fig. 14. 2D Burgers equation. Solution at T = 0.80. Limiter (34) deactivated and positivity conditions (17) activated.

As can be seen, there appears an overshoot just behind the shock (umax = 0.876) and, besides, the parameter, α i , becomes greater than one, at this place (dashed line). Considering those results, one might conclude to the inefficiency of conditions (17)! Fortunately, this is not true and this can be proved by considering the results displayed by Fig. 14, below: To obtain those results, we enforced positivity conditions (17), by discarding the influence of α i , and we kept the limiter (34) deactivated. As one can see, we find again almost the results of Fig. 11: errorL = 1.41 × 10−2 and umax = 0.758, umin = −0.248. 1 Those latter results tell us that the failure illustrated by Fig. 13 is only due to the parameter, α i , which does not detect always, with sufficient accuracy, the positivity range of the numerical solution. Indeed, the computed values of α i are prone to computation errors; moreover, α i computed by formula (30) is only an approximation of the exact quantity, αin+1 : this brings about a

98

G. Capdeville / Computers and Fluids 144 (2017) 86–116

supplementary error. This shortcoming is especially sensitive in regions of discontinuity of the solution. Therefore, in that case, the limiter (34) becomes useful by correcting the deficiencies of the procedure that activates conditions 2 and 3, according to the computed value of α i .

2.8. Synthesis At this point, we think necessary to pause and summarize our findings. Those numerical results demonstrate that it is possible to design a high-order numerical model while ensuring monotonicity preservation properties. The positivity conditions (17) were designed for that purpose. However, those conditions do not act on the same level and at the same time.  RThe L positivity condition 1 is based 3 upon the ratio u − u /(u¯ i − u¯ j ), and, therefore, it is an O h condition (for a fourth-order interpolation procedure) when the solution is smooth, and an O (1 ) condition, otherwise. As proved by the benchmarks, this is a non-oscillatory condition that has a damping behavior similar to high-order dissipative terms, when the numerical solution is smooth, and that prevents unbounded oscillations when the numerical solution is discontinuous. However, this condition is not always active and is not always necessary to get a positive scheme. Conditions 2 and 3 are always O (1 ) conditions, whatever the regularity of the numerical solution. These conditions introduce second-order dissipative terms into the numerical procedure →

(ui (xq ; u¯ ) = 0) and can alter the numerical accuracy unless they are selectively activated according to a “positivity monitor”. The limiter (34), which introduces a maximum-principle into the interpolation procedure, acts to thwart the deficiencies due to the use of the “positivity monitor”. Therefore, upon considering those results one might conclude that the best strategy would be to impose, only, the positivity condition 1 along with the limiter (34). We are not in favor of such a solution. Firstly, because it does not strictly follow the theoretical analysis we made and so, the positivity of the resulting scheme is not strictly ensured; secondly, because the benchmarks we selected are not exhaustive and, maybe, there can exist situations where condition 2 and 3 become substantial; lastly, the numerical results obtained for the 2D Burgers equation, proved the efficiency of conditions 2 and 3, when activated. For all those reasons, the remainder of this article uses all the conditions from (17) combined with limiter (34). Finally, one should note a significant advantage in using the algorithm we have just detailed: the high-order interpolation polynomial is absolutely free of the positivity conditions, (17). In other words, we may use any interpolation polynomial of any degree, obtained with any other method than a least-square method, the way of introducing the positivity conditions will always be unchanged. Things are different if we use classical schemes such as ENO or CWENO, because conditions of monotonicity are strongly linked with the shape of the polynomial interpolation. So, having designed and validated the scalar procedure that defines the LSQP4 scheme, the next step, now, is to extend it to nonlinear hyperbolic systems of conservation laws.

3. Discretization of a non-linear system of conservation laws As an example of an hyperbolic non-linear system, we make the choice of the Euler equations of gas-dynamics.

However, the reasoning and methods presented in this section are not limited to those equations and are directly transferable to any non-linear hyperbolic system. Therefore, we want to generate a monotonicity-preserving numerical algorithm for non-linear systems, which uses the numerical procedures designed in the preceding section, for a scalar problem. To begin, we give a brief description of our framework and we put the question of the choice of the formulation for generating the equivalent of Eq. (15). 3.1. Finite-volume discretization We consider the following hyperbolic system of conservation laws, for the vector of conservative variables: U ≡ [ρ , ρv, E ]t : →→

→

Ut + ∇ . F (U ) = 0, ∀ x ∈ D

(44)

With the same accuracy level and the same hypothesis as in the scalar problem, this latter system can be approximated by the following semi-discrete form:



|Ii | ×

dU¯ dt



=− i

 

  → ωq × f˜ VqL , VqR ; n × Sij

(45)

j∈Ni q∈Qij

Formally, this equation is identical to Eq. (3), with, now, the following approximation:



 →     → → → ωq × f˜ VqL , VqR ; n × Sij F U xq .nij dS, ∀ j ∈ Ni

(46)

Sij

q∈Qij

→

→

→

f˜(VqL , VqR ; n ) represents the numerical approximation of F (U ). n at

→

→

x = xq , which is computed by using the HLLE flux of [4] and the vector of interpolated values, VqL , VqR , computed from the vector →

of primitive variables, V ≡ [ρ , v , p]t . Note 1. As can be seen, we arbitrarily selected the vector of the interpolated primitive variables, VqL,R , for the definition and the com→ putation of the numerical flux, f˜(V L , V R ; n ). This choice is mainly q

q

dictated by practical worries concerning the positivity of the interpolated quantities such as the pressure and the density. However, the consistency of the numerical flux is not affected → → → by that choice and, of course, we have always: f˜(V¯ , V¯ ; n ) ≡ F (U¯ ). n Practically, we formulate the numerical flux as follows:

⎧→ → ⎪ F VqL . n if λLn > 0 ⎪ ⎪ ⎪ ⎪ → → → ⎪ ⎪ ⎪ λRn F VqR .→ n − λLn F VqL . n ⎨ λR λL   → − Rn nL R L f˜ VqL , VqR ; n ≡ λn − λn λn − λn ⎪  R L ⎪ L R ⎪ × V −V if λ  0  λ ⎪ q q n n ⎪ ⎪ ⎪ ⎪ ⎩→ R → F Vq . n if λRn < 0

(47)

For  L the Euler equations, if we define the characteristic velocities, λn , λRn , as:

→L →  ⎧ → → ⎨λLn ≡ min vq . n − aLq , v i+1/2 . n − ai+1/2   → → ⎩λR ≡ max v→ R .→ R n + a , v . n + a q i+1/2 i+1/2 n q

(48)

where the notation (. )i+1/2 , stands for the Roe average, then, the resulting flux is a monotone flux that computes accurately the speed of a shock wave (see [4]), provided that this latter discontinuity is aligned with the mesh. For any hyperbolic system, these properties should be re-demonstrated.

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Therefore, to ensure that discretization (45) is positive, it is left to produce the positivity conditions that act upon the interpolation procedure generating the components of VqL , ∀q ∈ Qi , in the cell Ii . Unfortunately, the form (45) based upon the conservative formulation (44), is difficult to handle for generating positivity conditions. Indeed, the more natural way to extend our scalar procedure to a system is to diagonalize this system. For that purpose, the convective form is preferable. Thus, we present in the subsection that follows, the procedure we designed for extending the results of the scalar analysis; this latter procedure mainly draws on the earlier work of H. Deconinck et al. [3].

get the following semi-discrete form, for each component, w¯ np , of ¯ n (average of Wn , over Ii ): W



|Ii | ×

dwnp dt



i

→

(49)



|Ii | ×

dwnp dt

 

ωq ×

$

→ np

λ



→

U i . nij



(55)

Upon following the same logic as in the scalar case, we can rewrite those equations in the following way:

dwnp dt



=−

 

$

ωq ×

j∈Ni q∈Qij

i

 ×

wnpi − wnpk

 −

wLnpq − wnpi

→ np

λ



→

U i . nij



→ np

λ

 

Ui .

→ nij

wnpi − wnpk

− 

+



wRnpq − wLnpq

 

wnpj − wnpi

wnpi − wnpj



    %  → − wLnpq − wnpi  → − λ np U i . nij wnpi − wnpj w −w npj

Having introduced this approximation, we discretize (52) by closely following the scalar procedure detailed in Section 2. Thus, upon integrating (52) over the discrete cell, Ii , and up → winding the spatial term according to the sign of λnp U¯ i .nij , we

wLnpq − wnpi



npi



(51)

(53)

+ 

j∈Ni q∈Qij

i

(50)

(52)

(54)

× Sij − Qwp U i × |Ii |, p ∈ {1, 2, 3, 4}

|Ii | ×

With that choice, we can define the vector of characteristic vari ¯ U ables, Wn , by the following relationship: Wn ≡ R−1 .U, or in a i n more usable form for what follows:

  −1 ¯ ¯ ∂ Wn ≡ R−1 n Ui .∂ U ≡ Rn Vi .∂ V

% →  → − λ np U i . nij × wRnpq  × Sij − Qwp U i × |Ii |, p ∈ {1, 2, 3, 4} =−

or, equivalently, for each component, wnp , of Wn and, Qwp , of QW :

 ∂ wnp →  ¯ → + λ np Ui .∇ wnp + Qwp U¯ i = 0, p ∈ {1, 2, 3, 4} ∂t

× wLnpq

→  → −  λ np U i . nij wRnpq − wLnpq % →  → −  L + λ np U i . nij wnpq − wnpi



Except for the source term, Qwp , this latter form is very close to the scalar form, (2). However, such a form is not of practical use. Indeed, because of the non-linearity of system (49) or its characteristic form, (50), the explicit computation of the wnp ’s is feasible only in special circumstances (isentropy, for example, for the case of Euler equations); in the general case, we are reduced at employing an approximation based upon a local linearization so as to compute the wnp ’s. Of course, by doing so, we introduce an extra error into the discretization procedure. However, one can prove (see Appendix B, for the proof) that this error does not alter significantly the positive properties of the scheme we are going to derive. Therefore, we replace Eq. (51) by its locally linearized form that follows:

+

j∈Ni q∈Qij



Now, if we define characteristic directions, say nα , along which we diagonalize system (49) (see Appendix A, for technical details concerning the derivation), we get the following equivalent system, for Wn , vector of characteristic variables:

∂ wnp → → + λ np .∇ wnp + Qwp = 0, p ∈ {1, 2, 3, 4} ∂t

→

U i . nij

+

→

→ ∂ Wn → + n (Wn ).∇ Wn + QW = 0 ∂t

λ

This result can be formulated in the more attractive form:

→

A (V ): generalized Jacobian matrix for the convective variables. In the smooth case, this system is equivalent to (44).



→ np

ωq ×

+

Now, instead of considering the conservative form, (44), we start from the following convective form, for the primitive variables: →

$

 

=−

3.2. Positive discretization of the characteristic form of an hyperbolic system

Vt + A (V ).∇ V = 0

99

× Sij − Qwp U i × |Ii |, p ∈ {1, 2, 3, 4}

(56)

This result should be brought closer to that of Eq. (15), for the scalar case. Therefore, on the basis of this writing, one can produce the fol¯ n: lowing condensed form, for each component, w¯ np , of W



|Ii | ×

dw¯ np dt



=− i

 j∈Ni



D+ × w¯ npi −w¯ npj pij





− Qwp U¯ i × |Ii |, p ∈ {1, .., 4}

(57)

D+ pij

where is directly identified from (56). Thus, if the following set of sufficient conditions is ensured, for each component, w¯ np , ∀ p ∈ {1, 2, 3, 4}:

⎧ R  wnpq −wLnpq ⎪ ⎪  0, ∀q ∈ Qij , j ∈ Ni : condition 1 ⎪ ⎪ w¯ npj −w¯ npi ⎪ ⎪ ⎨ wL −w¯  npi npq  0, ∀q ∈ Qij , j ∈ Ni : condition 2 ¯ ¯ − w w ⎪ ⎪  Lnpj npi  ⎪ ⎪ wnpq −w¯ npi ⎪ ⎪  0, ∀q ∈ Qik , k ∈ Ni , k = j: condition 3 ⎩

(58)

w¯ npi −w¯ npk

then, we have the result: D+  0, ∀ j ∈ Ni , p ∈ {1, 2, 3, 4}. pij In kind, this set of sufficient conditions is identical to conditions (17) obtained for the scalar problem. Therefore, if we use the scalar procedure ensuring monotonicity of the interpolating polynomial, for each component, wnp , of the characteristic vector, Wn , then we can conclude that conditions (58), are verified, whatever p ∈ {1, 2, 3, 4}.

100

G. Capdeville / Computers and Fluids 144 (2017) 86–116

In that case, we may integrate (62) in time to produce the following forward Euler discretization, ∀p ∈ {1, 2, 3, 4}:



+1 = w¯ nnpi × w¯ nnpi

1−

+ t ×





t × D+ −t × |Ii | j∈N pij i

D+ pij

×

 

Qwp U¯ i w¯ nnpi

w¯ nnpj

(59)

j∈Ni

As can be seen from this result, if the time-step, t, verifies the following CFL-like condition:



t ×



     1, ∀i, ∀ p ∈ {1, 2, 3, 4} (60) 

 Qwp U¯  1 i  × D+ + pij |Ii |  w¯ nnpi j∈Ni

then the resulting scheme is a positive scheme, with the following associated properties:



+1 1. w¯ nnpi ∈ [m p , M p ], ∀i with : M p ≡ max w¯ nnpi



min w¯ nnpi

i



2. w¯ nnp+1 

L2





and m p ≡

i





 (1 + O (h ) ) × w¯ nnp  . L





D+ , p ∈ {1, .., 4} pij

(61)

Doing so, Eq. (63) can be written in a more compact form:



|Ii | ×

¯n dW dt



=− i









¯ nj −QW U¯ i × |Ii | ¯ ni −W D+ × W ij

j∈Ni

(62)



Multiplying this equation by Rn U¯ i , in order to return to the conservative variables, U¯ i , and using the linearized relationship, (53), we get the following result:

      ¯ ¯ ¯ |Ii | × ddtU¯ = − Rn U¯ i .D+ij .R−1 n Ui × Ui −U j i   j∈Ni −Rn U¯ i .QW U¯ i × |Ii | Lastly, if we introduce the notations:

Cij+









=− i









Cij+ × U¯ i −U¯ j −QU U¯ i × |Ii |





¯ ¯ ¯ ¯ ≡ Rn U¯ i .D+ .R−1 n Ui and QU Ui ≡ Rn Ui .QW Ui ij

(63)

Now, Cij+ , is a diagonal matrix, similar to D+ , of which all ij the eigenvalues are positive quantities, provided that conditions (58) are verified. By comparing this latter equation with the finite-volume discretization, (45), we see that the right-hand term of this latter equation represents the consistent approximation of the flux in& → → tegral: − S F (U ).nij dS. ij

Note, however, that discretization, (65), is not a conservative discretization: it is simply used for deriving the proofs of the positivity of the scheme. Therefore, by integrating in time, we obtain the discrete conservative form, equivalent to the characteristic form, (59):



  t  + ¯ I− × C −t Qi × I U¯ in + t × Cij+ × U¯ jn (66) |Ii | j∈N ij j∈N i

where we introduced the following useful definition, for the source term:



Q¯ i ≡

t

QU U¯ i .U i

 2 U¯ i 

(67)

Consequently, if the time-step verifies the global CFL-like condition:



I  | i | t  max i maxC +  + |Ii | × Q¯ i p



(68)

pij

then one can state that the linearized conservative discretization of system of equations, (44), is positive; according to the proofs given in [7], we have the following stability property:

 n+1    U¯   (1 + O (h ) )U¯ n  L L 2

(69) 2

for the forward Euler discretization, (66). Remark 7. To obtain this latter result, we implicitly used the theoretical result saying that the only functional to be bounded for symmetrizable linear systems of hyperbolic equations, is energy, represented by the L2 − norm of the solution (see [7]). Unfortunately, the system of equations studied in this paper, is non-linear, and not always symmetrizable. Strictly speaking, this means that the hypothesis and the conclusion for the boundedness of the L2 norm of the solution are no longer valid. However, we adopt the welcome heuristic reasoning of P. D. Lax [7], who explained why it might exist “plausible” reasons – at least for the discretization of Euler equations – for a successful extension of the linear theoretical proof of stability to non-linear problems. In other words, the only result we can acknowledge with Eq. (66), is that the discretization of the locally linearized form of the system, (44), is presumably stable in L2 − norm, when we enforce conditions (58) and (68). Although this remark gives rise to a minimalist finding, we think that those first results constitute a reliable basis for the design of effective monotonicity-preserving numerical schemes. Now, we pursue our development of a positive scheme, for nonlinear systems, by explaining the way of computing the charac→

(64)

(65)

j∈Ni

i

To begin, we introduce the following definition:



dU¯ dt

2

3.3. Positive discretization of the conservative form of an hyperbolic system

≡ diag



|Ii | ×

U¯ in+1 =

These results are not surprising since they come directly from the definition of a positive scheme. On the contrary, the CFL condition, (60), is more interesting. Indeed, the coupling term, QW ≡ [{Qwp }, p ∈ {1, . . . , 4}]t , acts so as to limit the variation range of the time-step: the greater this term is, the smaller must be t in order to ensure the positivity condition. Or, in other words, if the characteristic system, (50), is strongly coupled, then the generation of a high-order positive scheme discretizing that system, becomes problematic, simply because we have to use small time-steps to balance the significance of QW . Therefore, there is a great interest, for the efficiency of the numerical method, in minimizing, in some norm, the quantity, QW ; this point will be emphasized in a later subsection. Up to this point, one can state that, if conditions (58) and (60) are warranted, then the high-order and upwind discretization of system (52) is monotonicity preserving, according to the above properties. However, we are interested by the discretization of the conservative form, (44). Now, the question that arises is the following one: how the positivity result we have just announced can be introduced into the discretization of the conservative form so as to generate a positive scheme? This point is discussed in the following subsection.

D+ ij

we can simply write:

teristic directions, nα . These latter directions are needed for the discretization of the characteristic form and the implementation of positivity conditions, (58).

G. Capdeville / Computers and Fluids 144 (2017) 86–116

3.4. Computation of the characteristic directions 3.4.1. Computation of the characteristic directions: the solution of Deconinck et al. [3] →

The problem is to find the characteristic directions, nα , which minimize, in some norm, the source term, QW , which appears into the system of Euler equations. Looking at the components of this quantity, Eq. (97) of Appendix A, we can see that this minimization can be achieved in a simple way, by enforcing the following conditions:

⎧→ → ⎨n2 such that →s .∇ p = 0   → → →  ⎩n3 such that→s →s .∇ v  ≡ Min

(70)

The physical meaning of the first condition is obvious since it →

indicates that n2 must be aligned with the pressure gradient. So, if →

we define this first characteristic direction by n2 ≡ (cos θ , sin θ ) , then we have for the definition of θ :

θ ≡ Atan

t

p  y

(71)

px

The second condition in (70), is more difficult to assess and to translate in a physical meaning. To clarify this condition, let us de→ → →→

velop s. ( s .∇ v ), in the 2D Cartesian frame. We get, then: →

→ → →

s . s .∇ v

= s2x ux + (vx + uy )sx sy + s2y vy

(72)

→

or equivalently, by using the fact that s is a unit vector:





→ → →→

s. s .∇ v

= s2x ux ± (vx + uy )sx

'



1−s2x + 1−s2x

vy

1 2

1 2

ux + vy ≡ Min 2 (74)

If we derive this condition with respect to θ , weget the new condition for the existence of an extremum:



−(ux −vy ) sin 2θ ± (vx + uy ) cos 2θ = 0 ⇔ tg2θ = ±

vx + uy ux −vy



 v¯

u¯ i

(78)

with the quantities (A, B), being defined according to the condition  considered in (70), and θi ≡ Atan AB . Thus, if we consider the first condition of (70), we define:

A ≡ ( px )i , B ≡ ( py )i

(79)

whereas for the second condition of (70), we set:

A ≡ (ux − vy )i , B ≡ (vx + uy )i

(80)

That way, assuming that all the gradients are approximated from the coefficients of the polynomial interpolation in each discrete cell, Ii , we have the following numerical properties:





1. in regions of discontinuity: |A| = |B| = O h1 ⇒ θ(i = Atan AB + O (h ) 2. in regions of smoothness: |A| = |B| = O (1 ) ⇒ θ(i = β ×  Atan AB + (1−β ) × Atan uv¯¯ , β ∈ [0, 1] i  3. in regions of uniformity: |A| = |B| ≡ 0 ⇒ θ( = Atan v¯ u¯ i

Doing this, the characteristic directions are computed with less numerical noise; in addition, in regions of uniformity of the solution, the characteristic directions coincide with the convection direction. Apart from the spurious oscillations, the procedure we have just detailed has also the drawback to be somewhat cumbersome and, above all, restricted to the Euler equations of gasdynamics. For this reason, we detail in the next subsection, an alternative procedure that is more general and easier to implement into a numerical algorithm.

(75)

(76)

direction, given by n = (cos θ , sin θ ) , which minimizes, in some norm, the tangential derivatives of the numerical solution. Mathematically, formulated for a 2D problem, we search θ , in each discrete cell, Ii , which is such that:

For this solution, the function in (74) is an extremum. To guarantee that this extremum is a minimum, we select the value of θ that gives the absolute minimum value among the solutions proposed by (74). →

This being done, the second characteristic direction, n3 , is simply given by: n3 ≡ (−sinθ , cos θ ) . The physical meaning of this latter solution, is analyzed in [3], where it is demonstrated that this computation of θ amounts to t

→

choosing s along the principal axis corresponding to the minimum normal strain of the flow. → → Consequently, if we compute (n2 , n3 ), according to formula (71) and (76), then we have the following result for QW : t

(|A| + |B| ) × θi + Atan 1 + (|A| + |B| )

→

⎧   vx + uy 1 ⎪ ⎪ + kπ , ∀k ∈ Z ⎪ Atan ⎨ 2 ux −vy θ=   ⎪ v x + uy ⎪ π 1 ⎪ 2 − 2 Atan + kπ , ∀k ∈ Z ⎩ ux −vy

QW ≡ [0, 0, |Min|, |Min|]

θ(i ≡

3.4.2. Computation of the characteristic directions: least-square formulation For the sake of simplicity, we look for a single characteristic

The analytic solution for this problem is given by:

→

where |Min| is drawn from Eq. (74). Thus, if QW is minimized, so is Q¯ i , which is defined by (67). However, as already mentioned, this method of computing the characteristic direction may introduce spurious oscillations into the numerical solution. These oscillations may come from computations of the solution gradients in regions of uniformity, for example. To attenuate this likely flaw, we practically compute the characteristic directions, symbolized by the generic quantity, θ i , by the following heuristic formula:

i

(73)

and, finally, if we set sx ≡ cos θ , the condition 2 from (70) is tantamount to writing the following minimization problem:

∃θ such that (ux −vy ) cos 2θ ± (vx + uy ) sin 2θ +

101

(77)



∂V → ∂s





t

∂V ∂V ≡ − sin θ × + cos θ × ∂ x ∂y i



=0

(81)

i

Practically, this condition upon the parameter θ is ensured in a least-square sense, by solving the following equation:

tgθ =

 ∂V

∂V ∂y , ∂x

 ∂ V 2  



(82)

∂x

where ( ∂∂Vy , ∂∂Vx ) ≡ ( ∂∂Vy )t . ∂∂Vx , represents the Euclidean scalar product. The quantities ( ∂∂Vx )i and ( ∂∂Vy )i are computed from the knowledge of the coefficients, D0,1 and D1,0 , of the interpolation polynomial, for each component of V. This choice of computation for θ not only simplifies the numerical procedure but also, in the case of the Euler equations, it enables a minimization of the source term, Q¯ i .

102

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Indeed, if we look at the meaning of the components that appear in QW , in such a case (see Appendix A), we can see that all the gradients, into these components, are expressed along the tangen→

tial direction, s ; therefore, by computing θ according to condition (82), we minimize the tangential terms of QW . In the more general case of any system, this property should be demonstrated. Again, the characteristic direction is practically computed from the weighting procedure, (78), to produce the useful direction given by θ(i . So far, we have detailed all the salient features of the numerical model for discretizing, with monotonicity constraints, any hyperbolic system of equations. Now, we summarize and organize these latter features in the algorithm presented in the section that follows.

3.5. Numerical algorithm for monotonicity-preservation Algorithm for a non-linear hyperbolic system →

• Step 1: Computation of Vi (xq ; V¯ ) ≡ VqL − V¯i , q ∈ Qi , for all i. • Step 2 : Computation of the characteristic directions, θ˜i , from (78) according to (76) (Euler equations) or (82), for all i. • Step 3: Computation of condition 1 from (58). → For all i, compute W 1 (xq ; V¯ ), q ∈ Qi , as follows:

⎧  ni   →   ⎫ →   ⎪ ⎨ q∈Q ωq lnp (V¯i ). V¯ j −V¯i +V j xq ;V¯  ⎪ ⎬  ij  →  1. ϕijp ≡ Min  , 1 j ∈ Ni , p ∈ →  ⎪  ⎪ ωq lnp (V¯i ).Vi xq ;V¯ ⎩  ⎭ q∈Qij → → → {1,.., 4}, lnp V¯i : pth left eigenvector of : A V¯i . n →  "  → ! #t 1 x ; V¯ ≡ 2. Wni w1npi xq , p ∈ {1, .., 4} q → →  → with : w1npi xq ≡ ϕijp × lnp V¯i .Vi xq ; V¯ j ∈ Ni , q ∈ Qij

• Step 4 : Computation of the “positivity monitor”, αin+1 . For all i, do:  ! 1. αin =

 ¯ n−1  U j   ! i  !  ¯ n−1    max U j  −min U¯ nj −1 

U¯ in −min j∈N

j∈Ni

j∈Ni

2. if n > 2 then  compute: αin+1 = 12 5αin −4αin−1 + αin−2 Else αin+1 ≡ 1

• Step 5: Computations of conditions 2 and 3 from (58), according to the value of αin+1 For all i, do: If αin+1 ∈ / [0, 1]then



1.

s1jp ≡ sign

 q∈Qij

 →

ωq w1npi xq

j ∈ Ni , p ∈ {1, .., 4}.



s2jp ≡ sign

 q∈Qij

 →

ωq w1npi xq

k = j, k ∈ Ni , p ∈ {1, .., 4}.

→  



→  



× sign lnp V¯i . V¯ j −V¯i

× sign lnp V¯i . V¯i −V¯k

,

→

,

2. if s1jp < 0 or s2jp < 0 ⇒ w1npi xq = 0, j ∈ Ni , q ∈ Qij , p ∈ {1, .., 4}. Endif αni +1 • Step 6: Computation, for all i, of :  → → →   Vi1 xq ; V¯ ≡ w1npi xq × rnp V¯i , j ∈ Ni , q ∈ Qij →  rnp V¯i

p∈{1,..,4}

→

→

: pth right eigenvector of : A V¯i . n

• Step 7: Computation of ψ pi from equation (34), for each component, vp , of V, p ∈ {1, .., 4}. • Step 8: Computation of V˜qL : For all i, do: →  V˜qL ≡ V¯i + ψi × V 1 xq ; V¯ , j ∈ Ni , q ∈ Qij i

with ψ i ≡ diag{ψ pi , p ∈ {1, .., 4}}

Remark 8. →

a) Step1: The computation of Vi (xq ; V¯ ), suppose that the fourth-order least-square interpolation polynomial designed in section 2.5 has been used to compute each interpolated component, vLpq , of the vector VqL . With the computation of this polynomial, we get the coefficients, Dl,k−l , which are also used into the computations of the characteristic directions, θ˜i . b) Step2: The computation of the characteristic directions, θ˜i , make it possible to compute, then, the matrices of left and →→



¯ ¯ right eigenvectors of A . n , R−1 n Vi and Rn (Vi ), respectively. →

c) Step3: The vector quantity, (V¯ j − V¯i + V j (xq ; V¯ )), that appears into the computation of ϕijp , could be written, also, as: VqR −V¯i , and represents the extra information that comes from the cell Ij , j ∈ Ni . d) Step 3: The condition q ∈ Qi is tantamount to the condition: q ∈ Qij , j ∈ Ni . We choose a formulation or another depending on our needs. e) Step 4: The quantity, αin , is computed from the vector of conservative variables, U¯ , because we want to test the positivity of the finite-volume discretization, (45). The Euclidean   scalar product is used to compute this quantity: U¯ in  ≡



1/2

. Logically, one should have computed αin , for each component of U¯ . However, the need to simplify the overall procedure and, above all, the numerical tests we did, which have proved that this simplification is not crucial for the monotonicity of the resulting scheme, led us to this solution. f) Step 6–8: Instead of keeping the use of characteristic variables into those steps, we have preferred to return to the primitive variables, V¯ . Indeed, the limitation procedure using the limiter, ψ i , from (34), was initially designed so as to introduce a “physical” maximum principle into the interpolation polynomial. In this spirit, it is therefore more logical to enforce this principle to the density and the pressure that must remain positive quantities. This constraint is enforced into step 8 of the algorithm. g) In comparison with a numerical scheme that would have used the primitive variables throughout the interpolation algorithm, the extra cost is around 25% in using the characteristic variables. U¯ in , U¯ in

Therefore, this algorithm leads to the computation of the corrected pointwise quantity, V˜qL , in each discrete cell, Ii , and for all the quadrature points, q, of this cell; or, said in other words, this procedure enables the computation of the quantities V˜qL , V˜qR , at each quadrature point of the cell interface, Sij . This being done, those quantities can be introduced in the finite-volume discretization, (45), in order to generate a conservative and positive discretization numerical method. The HLLE algorithm  L R from Section 3.1 enables to select between the quantities V˜q , V˜q in computing the numerical flux. The semi-discrete form, (45), is then integrated in time by the SSPRK(5,4) algorithm precised in Section 2.6.

G. Capdeville / Computers and Fluids 144 (2017) 86–116 Table 7 2D Euler equations: isentropic vortex advection. T = 0.10.

Table 6 → 2D Euler equations. ρ0 ( x ) = 1 + 0.20 × sin [π (x + y )]. T = 2. N

L1 -error

Order

L∞ -error

Order

10 20 40 80 160

3.45 × 10−3 2.19 × 10−4 1.38 × 10−5 8.62 × 10−7 5.40 × 10−8

4 4 4 4

5.92 × 10−3 3.70 × 10−4 2.25 × 10−5 1.39 × 10−6 8.60 × 10−8

– 4 4 4 4

Consequently, the resulting numerical model is designed to be a fourth-order accurate and monotonicity-preserving conservative scheme: similarly to the scalar case, we call it the LSQP4 scheme. Now, it is time to make some numerical tests in order, not only to validate the theoretical results, but also to compare the alternatives we proposed. 3.6. Numerical experiments In all that follows, these tests are applied to the 2D formulation of the Euler equations for gas-dynamics. Although the least-square interpolation procedure is formulated to be used over curvilinear meshes, we restrict ourselves to Cartesian grids, with h ≡ x = y. According to this latter choice, the CFL condition, (68), becomes:



t  max i

 

h2

 





(83)

max u¯ ni  + a¯ ni , v¯ ni  + a¯ ni + h2 × Q¯ i

We begin this series of tests with an accuracy analysis. Test 1: Time evolution of a sinusoidal solution. We compute the Euler equations of gas-dynamics, (44), with the following initial data:

 ⎧ → ⎪ ρ x , t = 0 = 1 + 0.20 × sin [π (x + y )] ⎪ ⎪ ⎨ →  → u x , t = 0 = 0.70∀ x ∈ [0, 2]2 ⎪   ⎪ ⎪ ⎩ →

N

L1 -error

Order

L∞ -error

Order

10 20 40 80 160

2.72 × 10−3 1.73 × 10−4 1.10 × 10−5 7.52 × 10−7 5.61 × 10−8

– 4 4 3.9 3.8

4.42 × 10−2 3.03 × 10−3 2.61 × 10−4 4.11 × 10−5 5.14 × 10−6

– 3.9 3.5 2.9 3.0

An isotropic vortex perturbation is propagated in a mean flow, given by: (ρ , u, v, p) = (1, 1, 1, 1 ). The perturbation is characterized by:

⎧ (1−r2 ) ε ⎪ × e 2 × (−(y − y0 ), (x − x0 ) ) ⎨(δ u, δv ) ≡ 2π γ − 1 )ε 2 ⎪ 2 ( ⎩δ T ≡ − × e(1−r ) 8γ π 2

with: r 2 ≡ (x−x0 )2 + (y−y0 )2 , T ≡ p/ρ , γ = 1.40. The computation domain is taken as [0, 10]2 . We set (x0 , y0 ) = (5, 5 ) and the vortex parameter value is ε ≡ 5. Computations are run up to the dimensionless time, T = 0.10, and the CFL number is 1. Periodic boundary conditions are used. The exact solution, at any time, is the initial solution translated according to the advection velocity (1, 1)t . Computations are run with the least-square solution, (82), for the characteristic direction. Table 7, below shows the convergence rate of the LSQP4 scheme. As can be seen, the theoretical accuracy of the numerical scheme is slightly deteriorated on the finest grids: this problem is mainly caused by small spurious oscillations that appear at the extremum of the solution; these oscillations are caused, either by the extra non-linearities introduced by the computation of the characteristic direction, or by cumulative errors generated by the multistep SSPRK method. Those results are unchanged if we use formulas, (71) and (76), for the computation of the characteristic directions.   On the finest grids, we get: maxQ¯ i  0.50. This result means i

p x,t = 0 = 1

The exact solution is simply the passive advection of the →

density perturbation in the mean flow: ρ ( x , t ) = 1 + 0.20 × sin [π (x + y−t )]. Computations are run up to the dimensionless time, T = 2, and the optimized CFL number that enables to ensure condition, (83), is set to 1. N discretization points are used in each space direction and periodic boundary conditions are enforced. The L1 - and L∞ -norms of the discretization error are reported in Table 6, below: Those results have been obtained by using a least-square computation of the characteristic directions, formula (82); however, the choice of characteristic directions computed according to formulas (71) and (76), does not change anything to those results. As expected, we get an uniformly accurate scheme, even at smooth extrema of the solution. The designed fourth-order accuracy is reached, even on the coarsest meshes. In addition, the following result, for the coupling term:   we have maxQ¯ i  = O 10−14 . This means that the decoupling of the Euler i

equations is effective when using the least-square solution, (82); moreover, the “monotonicity monitor”, α i , is everywhere and every time, in the range [0,1]. Consequently, the LSQP4 algorithm is in agreement with the theoretical results, for this test-case. Test 2: Advection of an isentropic vortex [8].

103

that the equations are not entirely decoupled; however, since it is multiplied by h2 , Q¯ i has a minor influence over the computation of the time-step, (83). For all the computations, we have: α i ∈ [0, 1]. Again, the monotonicity properties of the numerical scheme, are confirmed, for a smooth solution. Test 3: Rotated Riemann problems. We consider a series of classical Riemann problems that we decided to rotate in order to emphasize the multi-dimensional features of the LSQP4 scheme. The principle of this rotation (angle θ ) is illustrated by Fig. (15): now, the one-dimensional direction in which the Riemann problem is defined, is supposed to be the ξ -direction; the orthogonal direction (η-direction), simply becomes a periodic solution. For each solution investigated, the CFL selected is the optimum quantity in terms of the L1 -norm of the discretization error. Two values of the parameter, θ , are selected for this study: θ = 0 (one-dimensional problem) and θ = Atan(2 ) (bi-dimensional problem). Example 1. Lax’s problem This classical problem is solved over the bi-dimensional domain [0, 1]2 , with N=100 grid points in the x-direction and only 5 points in the y-direction. For this test, the Riemann problem, in the ξ -direction, is given by the following initial data:



(ρ , u, v, p)t=0 =

(0.445, 0.698, 0, 3.528 ) if ξ < 0.50 (0.5, 0, 0, 0.571 )

if

ξ > 0.50

104

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Fig. 15. Bi-dimensional rotated Riemann problem.

Fig. 16. Lax’s problem at T = 0.130. N = 100 points, CFL = 0.25. Computation of characteristic directions from (82). a) θ = 0 ; b) θ = Atan(2 )

The CFL number is 0.25 and computations are run up to T = 0.130. Fig. (16), below, display the density profiles obtained with both θ = 0 (1D computations) and θ = Atan(2 ) (2D computations) and for a characteristic direction computed with formula (82). As one can see, the solution is free of spurious oscillations and the maximum principle is locally verified, except for the last point behind the shock wave: ρmax = 1.313 instead of 1.304, for the exact solution on the same grid. However, when the mesh is regularly refined, ρ max converges towards the exact value: this means that this error is not liable to a Gibbs-like oscillation but is likely linked to the linearization error we made for deriving the positivity conditions (see Appendix B). As expected, when the Riemann problem is rotated, the discretization error increases, mainly because the contact discontinuity  is strongly damped; for this case (Fig. 16b), we have: maxQ¯ i  = 58. This latter result means that the characteristic sysi

tem of equations is not entirely decoupled; in addition, the size of the admissible is reduced (CFL condition (83)) because  time-step h2 × Q¯ i = O 10−1 . Of course, we check that Q¯ i ≡ 0, when we compute a one-dimensional problem. The “monotonicity monitor”, α i , is everywhere in the range [0,1], even just behind the shock: this might explain, also, the sin-

gular point behind the shock since at this location, the conditions 2 and 3 from (58), are not activated. Fortunately, according to its design principle, the limiter (34) acts so as to introduce a maximumprinciple into the interpolation procedure: the result is, especially, a limitation in the maximum value of the solution in the region close to the shock-wave. Now, if we compute the same problem with the characteristic directions from (71) and (76), we get the result typified by Fig. 17, below: There is no noticeable difference between this result and the one of Fig. 16b: in spite of a different computation of the characteristic directions, the contact discontinuity remains excessively damped. The problem comes, rather, from the numerical flux that is computed along directions which have no physical meaning. However,   if we look at the magnitude of the coupling term, we get maxQ¯ i =36: this time, the decoupling procedure based upon i

solutions (71) and (76) is more efficient than the least-square solution, (82). It is interesting to look at the regions where Q¯ i is maximum and, also, to check if the characteristic directions are correctly calculated.

G. Capdeville / Computers and Fluids 144 (2017) 86–116

105

Now, if we only change the way of computing the characteristic direction (formulas (71) and (76) instead of (82)), we get the results illustrated in Fig. 19, below: The first major difference comes from the region of the rarefaction wave, Fig. 19 (top). Indeed, this region is more irregular than in the previous case; in addition, in spite of an overall Q¯ i lower than in Fig. 18 (top), this latter quantity is higher in the region of the rarefaction wave. The second major difference comes from the characteristic direction, Fig. 19 (bottom), which shows some irregularities in the uniform region between the shock-wave and the rarefaction wave; those irregularities might explain the behavior of Q¯ i in the rarefaction wave region.

Example 2. Sod’s problem modified by Toro [12]. For this test, the initial data become:

 (ρ , u, v, p)t=0 = Fig. 17. Lax’s problem at T = 0.130. N = 100 points, CFL = 0.25. Computation of characteristic directions from (71) and (76), with θ = Atan(2 ).

This is done in Fig. 18, below, for the characteristic direction computed according to (82): Clearly, the maximum value for Q¯ i is in the region close to the shock-wave; however the magnitude of Q¯ i is also important in the smooth region defined by the rarefaction wave. In what concerns the characteristic direction computed according to the leastsquare procedure, (82), one can see that the physical direction, θ = Atan(2 ), is correctly identified, Fig. 18 (bottom), without any oscillation in regions of uniformity of the solution.

(1, 0.75, 0, 1 ) if ξ < 0.50 (0.125, 0, 0, 0.10 ) if ξ > 0.50

This test is a modified version of Sod’s test in order to generate a left sonic rarefaction wave, instead. The computations are run up to T = 0.20 and all the remaining parameters are unchanged. Then, we get the following results typified by Fig. 20, below, with the characteristic direction computed according to formula (82) : This time, the monotonicity of the numerical solution is preserved, everywhere and whatever the value of the rotation angle, θ. However, “the monotonicity monitor”, α i , is locally out of the range of [0,1]; this happens in the intermediate region between the shock-wave and the contact discontinuity: in that case, conditions

Fig. 18. Lax’s problem for N = 100 points and T = 0.130. Contour levels of Q¯ i (top). Contour levels of density with characteristic direction computed according to (82), (bottom).

106

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Fig. 19. Lax’s problem for N = 100 points and T = 0.130. Contour levels of Q¯ i (top). Contour levels of density with characteristic direction computed according to (71) and (76), (bottom).

Fig. 20. Sod’s problem at T = 0.20. N = 100 points, CFL = 0.25. Computation of characteristic directions from (82). (a) θ = 0 ; (b) θ = Atan(2 )

2 and 3 from (58), are activated for the characteristic variables concerned.   For θ = Atan(2 ), Fig. 20b, we have maxQ¯ i =83, with always a i

maximum value located in the region of the shock-wave. If we change the way of computing the characteristic direction by selecting, now, the formulas (71) and (76), we get exactly the same results as before except for the maximum value of Q¯ i , since we have, max Q¯ i =25: once again, the method of H. Deconinck i

et al., is more effective in minimizing the coupling terms of the characteristic system. Example 3. To check the robustness of the model in presence of highly varying smooth solutions, with density and pressure near

vacuum, we study the following Riemann problem:



(ρ , u, v, p)t=0 =

(1, −2, 0, 3 ) if ξ < 0.50 (1, 2, 0, 3 ) if ξ > 0.50

This is the “1-2-0-3” problem, initially defined in [4]. The solution of this problem consists in two strong rarefaction waves, travelling to the left and to the right. The initial conditions are selected to produce vacuum, at x = 0. Computations are run up to T = 0.150 with h = 1/100 and the CFL number is set to 0.125. For a characteristic direction computed according to (82), we get the results that follow, for the density profile. Fig. 21 show the resulting distribution of density for both cases, θ = 0 (Fig. 21a) and θ = Atan(2 ) (Fig. 21b). Whatever the case con-

G. Capdeville / Computers and Fluids 144 (2017) 86–116

107

Fig. 21. “1-2-0-3” problem at T = 0.150. Density profiles for N = 100 points. Computation of characteristic directions from (82). (a) θ = 0 ; (b) θ = Atan(2 )

sidered, the density is nowhere below zero and everywhere below 1. The same result is obtained for the pressure. One should note that the solution is not very sensitive to the grid orientation, although a slight increase in the L1 -norm of the error can be noticed when passing from θ =  0 to θ = Atan(2 ). For this choice of parameters, we have maxQ¯ i  = 6, whereas α i remains, everywhere, into i

the positivity range [0,1]. Now, if we use the solution of H. Deconinck for the characteristic direction, formula (71) and (76), there are no changes in the new results, including the maximum value of Q¯ i . Test 4: Two-dimensional Riemann problem [9]. In the square domain, [0, 1]2 , we consider the following bidimensional Riemann problem (“configuration 3”):

(ρ , u, v, p)t=0

⎧ (1.5, 0, 0, 1.5 ) if x > 0.50, y > 0.50 ⎪ ⎪ ⎪ ⎨(0.5323, 1.206, 0, 0.30 ) if x < 0.50, y > 0.50 = ⎪ (0.138, 1.206, 1.206, 0.029 ) if x < 0.50, y < 0.50 ⎪ ⎪ ⎩ (0.5323, 0, 1.206, 0.30 ) if x > 0.50, y < 0.50

The grid resolution is h = 1/1120, and computations are performed up to the dimensionless time, T = 0.30, with CFL = 1.20 in every cases. We begin the computations with the formula (82), for the computation of the characteristic direction. Fig. 22, shows the density contour lines. As can be seen, the oblique shock-waves are well resolved and the solution remains symmetric. The contact discontinuities which terminate in a jet, with its “mushroom” shape, are also clearly identified. There are no noticeable oscillations behind the discontinuities and the regions of uniformity of the solution are not altered by the slowly moving shock waves. If look at the characteristic directions in the computation domain as well as the contour levels of Q¯ i , we get the following figures (computed with h = 1/560, for the sake of readability): Firstly, one can see that the shock-wave directions are clearly identified by the least-square procedure (82), Fig. 23a, (one should notice that all directions are not plotted: this is why the horizontal shock wave for x ∈ [0.15, 0.25] seems to be not detected). The same tendency occurs for the contact discontinuities and the jet, although the identification of those phenomena is less marked. Looking at the contour lines for Q¯ i , Fig. 23b, we see that the maximum values are obtained in the regions defined by the four

Fig. 22. Bi-dimensional Riemann problem: “configuration 3”. Density contours at T = 0.30, for h = 1/1120. Computation of the characteristic direction with (82).

 

normal shock waves: we have maxQ¯ i  = 824, for h = 1/1120, at i

this location. However this latter result does not destroy the monotonicity properties of the numerical scheme. This point is illustrated by Fig. 24, below, which displays the contour lines of the “positivity monitor”, α i . This quantity is, everywhere, in the range [0,1], including the discontinuities and the jet: according to the theoretical meaning of α i , one can conclude that there is nowhere a loss of positivity of the numerical scheme. Note that Fig. 24 is plotted at a given instant, T = 0.30; at the prior moments, α i may have values that fall outside the range [0,1]. Therefore, monotonicity conditions 2 and 3 from (58) may be activated during the computation of the intermediate states of the numerical solution. This is, for example, what happens in the region close to the point (0.5,0.5), at the beginning of the numerical procedure. Now, if we change the way of computing the characteristic directions, all the rest being unchanged, we get results of Fig. 25.

108

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Fig. 23. Bi-dimensional Riemann problem: “configuration 3”. (a) Characteristic directions according to formula (82). (b) Levels for the coupling term, Q¯ i .

Fig. 24. Bi-dimensional Riemann problem: “configuration 3”. Levels of the “positivity monitor”, α i , for h = 1/1120. Computation of the characteristic direction with (82).

Fig. 25. Bi-dimensional Riemann problem: “configuration 3”. Density contours at T = 0.30, for h = 1/1120. Computation of the characteristic direction with (71) and (76).

Even though the overall structure of the solution is unchanged, there are noticeable modifications. Indeed, some Kelvin–Helmholtz instabilities begin to appear at the contact discontinuities; the consequence of these instabilities is to destabilize the jet that ends those discontinuities (Fig. 25). Fig. 26, enable to look at the characteristic directions computed according to (71) and (76) as well as the levels of Q¯ i . We note that there are some strong irregularities, in regions of uniformity, for the computation of the characteristic directions; this is especially true in what concerns the direction computed along the normal strain of the flow, Fig. 26a. The direction of pressure gradient is less concerned by those oscillations, Fig. 26b. In spite of these oscillations, the discontinuities are, overall, well identified by procedures (71) and (76): the pressure gradient direction does not detect the contact discontinuity, Fig. 26b, but the direction along the normal strain does, Fig. 26a. Note that if the smoothing process, (78), were removed, the oscillations would be more pronounced.

In addition, if we look at the contour levels of Q¯ i , we note that the regions of maximum value for this latter quantity, Fig. 26c, are complementary of those identified by Fig. 23b.   The maximum value is also lower since we have maxQ¯ i  = 439, i

for h = 1/1120. Lastly, Fig. 27 produces the contour lines for the “positivity monitor”, α i . In comparison with results of Fig. 24, the range of variation of α i is very similar. Thisquantity remains, everywhere, in the range [0,1]: the monotonicity preservation is verified, at least for the simulation time illustrated by Fig. 27. Test 5: Forward facing step problem [2,13]. This benchmark constitutes a good test for high-order methods that aim at computing complex flows with discontinuities. A right-going Mach3 flow enters a wind tunnel of 1 unit wide and 3 units long. The step is 0.2 unit high and is located 0.6 unit from the left-hand end of the tunnel.

G. Capdeville / Computers and Fluids 144 (2017) 86–116

109

→

Fig. 26. Bi-dimensional Riemann problem: “configuration 3”. (a) Characteristic directions according to Atan( (vx + uy )/(ux −vy ) ). (b) Characteristic directions according to ∇ p. (c) Contour levels for the coupling term, Q¯ i .

Reflective boundary conditions are applied along the walls of the tunnel; at the left is a flow-in boundary condition whereas, at the right, all gradients are assumed to vanish. Initially, the wind tunnel is filled with a gas which has, everywhere, the following dimensionless quantities as physical parameters:

(ρ , u, v, p)t=0 ≡ (1.40, 3, 0, 1 ) The corner of the step is the center of a rarefaction fan and, hence, is a singular point of the flow; unlike [13], we have not modified the schemes near the corner in any way; similarly, we did not use any mesh refinement at this location, as it is done in [2]. The mesh spacing selected is h = 1/160 and we used a CFL number of 0.90 for all the computations. Simulations are run up to the dimensionless time T = 4. We start the computations with formula (82) for defining the characteristic directions. Fig. 28 displays the density distribution in the wind tunnel.

The general shape and the position of the shock-waves are correctly computed. However, one can note some oscillations behind the stationary part of the bow shock. The contact discontinuity behind the upper Mach stem, although somewhat damped, clearly appears. There is no spurious Mach stem at the bottom wall as it may happen when there is no special treatment at the corner of the step (see [13] or [2]). If we look at the characteristic directions and the contour levels of Q¯ i , computed by using formula (82), we get Fig. 29 (obtained for h = 1/80). As one can see, Fig. 29a, the least-square procedure (82) correctly detects all the gradients, including the centered rarefaction   fan. In what concerns the coupling term, Q¯ i , we have maxQ¯ i  = i

120456 (for h = 1/160) and this maximum value is obtained in the region of the upper Mach stem, Fig. 29b. Lastly, Fig. 30 shows the contour levels for the “positivity monitor”, α i .

110

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Fig. 27. Bi-dimensional Riemann problem: “configuration 3”. Levels of the “positivity monitor”, α i , for h = 1/1120. Computation of the characteristic direction with (71) and (76).

1, 0  y  6; Neumann at 1  x  13, y = 0 and 0  x  13, y = 11. Again, unlike B. Cockburn [2], no special numerical treatment was enforced at the corner. The intensity of the rarefaction wave that appears at the corner, is such that it is easy to get negative density and/or pressure values for a high-order numerical scheme that would not preserve the monotonicity of the solution. For this particular benchmark, we selected the following numerical parameters: CFL = 1.10, h = 1/40 and T = 2.0, for the simulation time. To begin, computations are run with the formula (82) for the computation of the characteristic directions. Fig. 31 displays the contour lines of density. The contact surface that separates the part of the flow field processed by the diffracting shock wave from that processed by the incident rarefaction wave, is correctly predicted. The numerical solution is free of any instability and density and pressure remain, everywhere and at every time, positive quantities. Upon looking at the characteristic directions and the distribution of Q¯ i , in the computation domain, we get Fig. 32. As can be seen from Fig. 32a, the gradients are correctly identified   by the characteristic directions. In addition, we have maxQ¯ i  = 111 for h = 1/40, and this maximum is situated in the i

One can note some oscillations of this quantity behind the upper Mach stem; however α i remains, everywhere, in the positivity range [0,1]. The monotonicity conditions 2 and 3 from (58) are not activated, accordingly: this likely explains the permanency of numerical oscillations behind the Mach stem, Fig. 28. However, these conditions are sometimes activated, during the calculations, for some discrete points in the computation domain. Now, if we change the way of computing the characteristic directions by using formulas (71) and (76), all the remaining parameters being unchanged, the computations result in a failure. Indeed, at the very beginning of the numerical simulation, some strong numerical oscillations appear in computing the gradients in the zone close to the corner of the step; these oscillations cannot be entirely damped by the smoothing procedure (78) and lead to a rapid divergence of the numerical algorithm. A lowering of the CFL number value, has no effect on this trouble. Test 6: Shock diffraction problem, [2]. A right-moving shock of Mach = 5.09, initially located at x = 0.50 and 6 ≤ y ≤ 11, moves into undisturbed air ahead of the shock, with a dimensionless density of 1.40 and a dimensionless pressure of 1. The boundary conditions are: inflow at x = 0, 6  y  11; outflow at x = 13, 0  y  11; reflective at 0  x  1, y = 6 and at x =

region of the normal part of the incident shock, Fig. 32b. Lastly, Fig. 33, displays the contour levels for α i . Now, there are some points where α i is out of the range [0,1]: for example, αi = 1.675, just behind the upper normal shock wave and αi = −1.045, for a point near the corner of the step. This means that the monotonicity conditions 2 and 3 from (58), are activated at those points. However, as one can see from Fig. 33, the main part of the computation domain is not concerned by those positivity corrections. By using formulas (71) and (76) to compute the characteristic directions, all the remaining parameters being unchanged, we get numerical results typified by Fig. 34, below. The results are almost identical with those of Fig. 31. The characteristic directions and the distribution of Q¯ i , in the computation domain, are displayed by Fig. 35 that follow. Both characteristic directions identify, correctly, the gradients of the solution, Fig. 35a and b. There are no significant oscillations for these directions and the transition with regions of uniformity is smoothly ensured. The region of significance   of Q¯ i , Fig. 35c, is the same as for Fig. 32b and we have maxQ¯ i  = 100, for h = 1/40. i

In what concerns the contour lines for α i , we get exactly the same graph as the one of Fig. 33, with the same properties.

Fig. 28. Forward-facing step problem (h = 1/160). Density contours at T = 4. Characteristic directions according to formula (82).

G. Capdeville / Computers and Fluids 144 (2017) 86–116

111

Fig. 29. Forward-facing step problem at T = 4. (a) Characteristic directions according to formula (82). (b) Contour Levels for the coupling term, Q¯ i .

Fig. 30. Forward-facing step problem at T = 4. Contour Levels of the “positivity monitor”, α i , for h = 1/160.

4. Summary and perspectives

Fig. 31. Shock wave diffraction (h = 1/40). Contour lines of density at T = 2.0. Computation of characteristic directions according to (82)

By initiating this work, our aim was to design a numerical model that enables the computation of complex compressible flows, with a high accuracy and without the numerical artifacts classically associated. For that purpose, we mainly have used the notion of “positive scheme”, early defined by S. Spekreijse [10], and the results of P. D. Lax [7], on positive hyperbolic systems. Upon extending the method proposed by T. J. Barth in [1], we have been able, thus, to derive three positivity conditions for discretizing a scalar hyperbolic conservation law. Through the use of numerical tests, we have validated the theoretical basis of our work and we have been able to analyze and to classify each condition. Having validated the scalar theory, we have used it in order to generate a numerical model for discretizing any non-linear hyperbolic system of conservation laws. To manage this generalization, we have had to find a method to overcome the problems induced by the coupling between equations constituting the system. The seminal work of H. Deconinck [3], has given us a key for reducing the significance of this coupling. Thus, we have been able to derive a new method based upon a least-square process that

112

G. Capdeville / Computers and Fluids 144 (2017) 86–116

Fig. 32. Shock wave diffraction (h = 1/40). (a) Characteristic directions according to formula (82). (b) Levels for the coupling term, Q¯ i .

Fig. 33. Shock wave diffraction (h = 1/40). Contour levels of the “positivity monitor”, α i .

minimizes this coupling. Then, numerical tests have helped us to highlight the benefits and shortcomings of our generalization. Upon examining this work, one can do some comments. In what concerns our initial aim, we may consider that this one is reached. Indeed, the numerical results we have obtained are in close agreement with our theoretical analysis: the LSQP4 scheme we have derived from the theory, is monotonicity-preserving while being fourth-order accurate. In addition, we think that the resulting algorithm is easier to improve than more classical schemes such as WENO or ENO, since the procedure to get high-order of accuracy and that to get a positive scheme, are totally decoupled. The point of view is more delicately shaded in what concerns the modeling of non-linear systems. Firstly, for practical considerations, we have had to limit our theoretical work to a locally linearized version of the initial system.Therefore, the theoretical positivity conditions we have derived from this simplification, do not guarantee the positivity of

Fig. 34. Shock wave diffraction (h = 1/40). Contour lines of density at T = 2.0. Computation of characteristic directions according to (71) and (76)

the non-linear algebraic form. This discrepancy is even more emphasized by the extra linearities introduced by the gradients based upon characteristic directions. On the other side, the way of minimizing the coupling term, is partially successful. Indeed, the first solution proposed by H. Deconinck [3], is not sufficiently robust to be systematically employed; however, the least-square alternative, although often generating coupling terms of greater magnitude, is more promising. The numerical results we have obtained with this latter approach are free of spurious oscillations and the gradient directions are clearly identified, whatever the complexity of the numerical solution. By relying upon this positive observation, we think that our model should be improved in, at least, two directions.

G. Capdeville / Computers and Fluids 144 (2017) 86–116

113

→

Fig. 35. Shock wave diffraction (h = 1/40). (a) Characteristic directions according to Atan( (vx + uy )/(ux −vy ) ). (b) Characteristic directions according to ∇ p. (c) Contour levels for the coupling term, Q¯ i .

Firstly, the characteristic directions given by the least-square procedure, might be used to improve the computation of the numerical flux. This way, it might be possible to introduce some multi-dimensional features into a genuinely uni-dimensional algorithm. Secondly, upon using the flexibility of the least-square approach generating the high-order interpolation polynomial, one might introduce a multi-dimensional bias into this interpolation, by the way of weighting coefficients depending on the gradient directions, for example. By increasing the influence of normal derivatives to gradients of the solution, we think that these two latter solutions might enable a minimization of the numerical effects associated to the coupling terms; this way, we could be nearer of an effective decoupling of the equations, and, consequently, we could approach, more closely, the performances of the scalar procedure.

Acknowledgment The author wishes to thank warmly, Adrien Grellier (ECNFrance), who provided a technical support essential to this work. Appendix A Approximate diagonalization of the Euler equations in two space dimensions [3]. In this appendix, we explain the way of deriving Eq. (50), starting from the convective form (49). These forms are linked to the Euler equation of gasdynamics. This derivation is only made for two dimensional problems. For an extension to three-dimensions, see, for example, C. Hirsch [6]. We start from the convective form: →

→

Vt + A (V ).∇ V = 0

(84)

114

G. Capdeville / Computers and Fluids 144 (2017) 86–116

This system is an hyperbolic system. Therefore, we have the mathematical property that the roots (characteristic velocities) of:

→ →



det A . n −λn I = 0

(85) →

are all real, whatever n , real. For the 2D Euler equations, these characteristic velocities are the following ones:

⎧ →→ ⎨λn1,2 ≡ v . n →→ λ ≡ v.n + a ⎩ n3 → → λn4 ≡ v . n −a

(86)

→→

saying that the matrix, A . n , is similar to the real diagonal matrix, →

n , whatever n , real. Mathematically, one can write, then, that it exists a real matrix, Rn (V), that can be inverted, to give:

A . n = Rn (V ).

−1 n Rn

(V )

(87) →→

Rn (V), is the matrix of right eigenvectors of A . n and its inverse, →→

R−1 n (V ), the matrix of left eigenvectors of A . n . We write those quantities as:

0

0

−1/a2

n2y

−n2x

0

n3x

n3y

−n3x

−n3y

1

⎜ ⎜0 ⎜ R−1 n (V ) = ⎜ ⎜0 ⎝ 0

→ 2

⎧→ → ⎪ λ n1,2 ≡ v ⎨→ → → λ n3 ≡ v + an3 ⎪ → ⎩→ → λ n4 ≡ v − an3

(94)

along with the definition:

∂ Wn → → + n .∇ Wn + QW = 0 ∂t

(88)

∂ (.) may represent a spatial increase as well as a time variation. Consequently, upon using definition (88) in the convective form (84), the differential system, using Wn , can be written as follows:

∂ Wn ∂ Wn ∂ Wn + R−1 + R−1 =0 n ARn n BRn ∂t ∂x ∂y

(89)

Unfortunately, the Jacobian matrices A and B, do not share the same set of eigenvectors: the consequence is that system (84) cannot be strictly diagonalized. However, if we follow the presentation of [3], system (89) can be formulated to be as close as possible of a diagonalized form. Indeed, for each component, wnp , p ∈ {1, 2, 3, 4}, of Wn , we can expand (89) in the following way:

⎧ → → ∂ wn1 ⎪ ⎪ ∂ t + v .∇ wn1 = 0 ⎪ ⎪ ⎪   ⎪ ⎪ → → → → ⎪ → ∂ wn2 ⎪ 1 2 ⎪ ⎨ ∂ t + v .∇ w n2 − ρ ıy ∧ n .∇ p   →  → → → →→ → ⎪ ∂ wn3 3 3 3 ⎪ .∇ wn3 + a × s . s .∇ v = 0 ⎪ ∂ t + v + an ⎪ ⎪ ⎪   →  ⎪ → → ⎪ → →→ → ⎪ ⎪ ⎩ ∂ w∂ tn4 + v − an3 .∇ wn4 + a × s3 . s3 .∇ v = 0

⎜0 ⎜ Rn (V ) = ⎜ ⎝0 0

n2y /P

n3x /2P

−n2x /P

n3y /2P

0

ρ a/2

ρ



→ ıy

→ 2



∧n

→ 3

→

→



→→

→ 3

→

%t

→→

.∇ p, a × s . s .∇ v , a × s . s .∇ v 3

3

(97) QW is the source term of the characteristic form of the system of Euler equations. This source term prevents an entire decoupling of the initial system, (84). To appreciate which form is taken by the coupling term, QW , in primitive variables, one can transform Eq. (96) to primitive variables by multiplying by Rn (V)and rearranging the spatial terms. Doing so, we get the following result [11]:

∂V ∂V ∂V + Rn (V ) x R−1 + Rn (V ) y R−1 + QV = 0 n (V ) n (V ) ∂t ∂x ∂y

(90)



→→

⎛

ρ

(98)

⎞

⎜ 0 ⎟ 0 ⎠ ρ a2

×⎝

(99)

→

→ →→

Therefore, QV is proportional to the term: s3 .(s3 .∇ v ). As early demonstrated in [3], this term represents the normal →

⎞

⎟ → → ⎟ , P ≡ n2 .n3 ⎟ −n3y /2P⎠ ρ a/2 −n3x /2P

QW ≡ 0, −

1

QV = s3 . s3 .∇ v

→

ρ /2a

$

→ →

In this formulation, we used the possibility for the shear wave →

(96)

with the following components of QW , identified from (90):

with for QV :

direction, n2 , to be different from the acoustic wave direction, n3 . With such a choice, the matrices Rn (V) and R−1 n (V ), linked to the primitive variables, can be precised. We get, then, the following results [11]:

ρ /2a

(95)

one can write, the final, semi-diagonalized form, which is condensed from system (90):

∂ Wn ≡ R−1 n (V ).∂ V

0

(93)

(see [3], for a physical meaning of each component) Lastly, upon using the following notations:

From those definitions, we can, now, define the vector of char→ acteristic variables, in direction n , Wn , by the following relationship:

1

(92)

→

"→ → → → # →

n ≡ diag λ n1 , λ n2 , λ n3 , λ n4

⎛

⎟ ⎟ ⎟ ⎟ 1/ρ a ⎟ ⎠ 1/ρ a

%t ∂ p →3 → ∂ p →3 → ∂ Wn = ∂ − 2 ,s .∂ v , + n .∂ v , −n .∂ v a a a ∂p

#t R−1 lnp , p ∈ {1, 2, 3, 4} and n (V ) ≡ " →  # Rn (V ) ≡ rnp , p ∈ {1, 2, 3, 4}

" → 

⎞

Hence, the above derivations make it possible to formulate the components of ∂ Wn , in two space dimensions:

$

The hyperbolicity of (84), can be formulated in another way by

→→

⎛

strain rate of the flow field in the direction s . For any hyperbolic system, the same procedure must be employed so as to formulate the equivalent coupling source term, QW . Appendix B Errors made in locally linearizing the scalar conservation law: →

(91)

→

ut + λ (u ).∇ u = 0. Our intent in this section is to evaluate the impact of a local linearization over the procedure designed for the monotonicity preservation of the LSQP4 scheme.

G. Capdeville / Computers and Fluids 144 (2017) 86–116

We start from the following assertion: in the smooth case, since the quasi-linear and the conservative forms are equivalent (they describe the same physical properties), the sufficient conditions for positivity we can draw from the discretizations of these both forms, must be the same, provided that both discretizations have the same numerical features. Therefore, let us begin our study with the quasi-linear form of the scalar conservation law: →

→

ut + λ (u ).∇ u = 0

(100)

By locally linearizing around the point, i, we get the modified equation: →

→

ut + λ (u¯ i ).∇ u = 0

(101)

Now, if we integrate this latter results over the discrete cell, Ii , then we get the following result:



|Ii | ×

du¯ dt



→

+ λ (u¯ i ).

→  →   ∇ ud x  = 0

 Ii

i

(102)

then, we can re-write the quasi-linear form, (100), as follows: →

|Ii | ×

du¯ dt





→

+ λ (u¯ i ). i

→ unij dS

Si

=0

(103)

Using the same notations and approximations as previously, the integral quantity can be approximated in the following way: →

λ (u¯ i ).



Si

→ unij dS



=





R

ωq λ+nij uLq + λ−nij uq Sij

j∈Ni q∈Qij

   ωq λ+nij uLq −u¯ i + λ−nij uRq −u¯ i Sij

(104)

j∈Ni q∈Qij

→ λ (u¯ i ).nij .

with the notation: λnij ≡ This latter result can also be re-written in a form similar to Eq. (15):

λ (u¯ i ).



→

Si

unij dS

$  L  u −u¯ i ωq λ+nij q (u¯ i −u¯ k )

 j∈Ni q∈Qij



− λ− nij

uLq −u¯ i

 −

 

u¯ j −u¯ i

λ−nij

uRq −uLq



u¯ j −u¯ i

u¯ i −u¯ k

u¯ i −u¯ j



%

u¯ i −u¯ j



Sij

(105)

Consequently, we can write the following semi-discrete approximation:



→

→

ut + λ (u ).∇ u



i





du¯ dt



+ i

  1 × D+ u¯ i −u¯ j ij |Ii |

(106)

j∈Ni

D+ , ij

with the coefficients, directly identified from Eq. (105). Thus, upon considering Eq. (105), the sufficient conditions to get D+  0, ∀ j ∈ Ni , are exactly the same as conditions (17), genij erated to warrant the positivity of the conservative discretization. Unfortunately, since we used a local linearization to get this result, this means that we introduced an error into the discretization procedure; the consequence appears with the fact that D+ and Cij+ , ij from Eq. (16), are not formally identical. So, let us evaluate the influence of this linearization error over the positivity condition of the discretization. For that purpose, let us start with the following local expansion, over the discrete cell, Ii : →

→

→





→

.∇ u + O h2 = 0

(108)

i

Upon integrating this result  over the discrete cell, Ii , and by neglecting the integrated O h2 residual term, compared to the two first leading terms, we get the following semi-discrete form:

 →   → ∂λ   + λ (u¯ i ). ∇ ud x  + . (u − u¯ i ) |Ii | × ∂u Ii Ii i i →  →   × ∇ ud x  = 0 

du¯ dt





→

→

(109)

The third term of Eq. (109) permits to size the leading error made by linearizing the quasi-linear form, (100). More precisely, we can re-write Eq. (109) as:

|Ii | ×

du¯ dt



 +

i

→ × nij dS

Si

→



→ u¯ i .nij dS

uλ ( )

1 + 2

→

∂λ ∂u



 . i

Si

(u − u¯ i )2

=0

(110)

And, by using the same quadrature rules for approximating the surface integrals, we get:

 du¯   ⎧ + ⎨|Ii | × dt i + j∈Ni Dij u¯ i −u¯ j  → →  ∂ λ .n × ω uL − u¯ 2 × S = 0 ⎩+ 1   q i ij ij q j∈Ni q∈Qij 2 ∂u

(111)

i

→

→

→

∂λ × ∂u





→

→

ut + λ (u ).∇ u = ut + λ (u¯ i ).∇ u + (u − u¯ i )

or, written equivalently:



115

λ (u ) = λ (u¯ i ) + (u − u¯ i ) ×



→

∂λ ∂u



 2

+O h i

(107)

More synthetically, this latter result can be written as

⎧      2 ⎪ |Ii | ddtu¯ i + D+ij u¯ i −u¯ j + Eij u¯ i −u¯ j = 0 ⎪ ⎪ j∈Ni j∈Ni ⎪ ⎨  →  →   L 2 u −u¯ ⎪ ⎪ Eij ≡ 12 Sij ∂∂ λu .nij × ωq u¯qj −u¯ ii ⎪ ⎪ i ⎩ q∈Qij and, finally, by [tn , tn+1 ] we get:

integrating,

explicitly,

over

(112)

the



u¯ ni +1

=

u¯ ni ×

 t  + t  + t × 1− × Dij + × D × u¯ nj + I | i | j∈N |Ii | j∈N ij |Ii | i i   n n 2

range

Eij u¯ i −u¯ j

(113)

j∈Ni

For bi-dimensional computations, the quantity, Eij , is an O (h ) quantity, since Sij = O (h ). Therefore, provided that the CFL condi tion is verified (t = O (h )), the last term of Eq. (113) is an O h2 quantity. From this latter result, we can conclude that if we apply the sufficient positivity conditions, (17), to the locally linearized form of Eq. (1) or, equivalently in the smooth case, to Eq. (100), then we do an O h2 error; in other words, we have henceforth the following result:

 2  2  ⎧ n+1  ⎨u¯ ∈ mi + O h , Mi + O h     ⎩mi ≡ min u¯ n , Mi ≡ max u¯ n j∈Ni

j

j∈Ni

(114)

j

Consequently, in the limit h → 0, the maximum principle is asymptotically verified.

116

G. Capdeville / Computers and Fluids 144 (2017) 86–116

References [1] Barth TJ. Aspects of unstructured grids and finite-volume solvers for the Euler and Navier–Stokes equations. VKI-Lect Ser 1994;1994-05. [2] Cockburn B, Shu C-W. The Runge–Kutta discontinuous Galerkin method for conservation laws. J Comp Phys 1988;141. [3] Deconinck H. A survey of upwind principles for the multi-dimensional euler equations. VKI-Lect Ser 1987;1987-02. [4] Einfeldt B, et al. On Godunov-type methods near low densities. J Comp Phys 1991;92. [5] Gottlieb S. On high-order strong stability preserving Runge–Kutta and multi step time discretizations. J Sci Comp 2005;25. [6] Hirsch C. Numerical computation of internal and external flows, Vol. 2. Wiley and Sons; 1992. Chapter 16:177-191 [7] Liu X-D, Lax PD. Positive schemes for solving multi-dimensional hyperbolic systems of conservation laws, II. J Comp Phys 2003;187.

[8] Shu CW. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. ICASE 1997;97-65. [9] Shulz-Rinne CW. Numerical solutions of the riemann problem for two-dimensional gas dynamics. SIAM 1993;14. [10] Spekreijse S. Multigrid solution of monotone second-order discretizations of hyperbolic conservation laws. Math Comp 1987;49. [11] Struijs R. A multi-dimensional upwind discretization method for the euler equations on unstructured grids, Ph.D. thesis. Delft University; 1994. [12] Toro EF. Riemann solvers and numerical methods for fluid dynamics. Springer-Verlag; 1997. [13] Woodward P, Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. J Comp Phys 1984;54. [14] Zhang X, Shu C-W. Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments. Proc R Soc 2011;467.