# Staggered Residual Distribution scheme for compressible flow

This paper is focused on the approximation of the Euler equation of compressible fluid dynamics on a staggered mesh. To this aim, the flow parameter are described by the velocity, the density and the internal energy. The thermodynamic quantities are described on the elements of the mesh, and this the approximation is only L^2, while the kinematic quantities are globally continuous. The method is general in the sense that the thermodynamical and kinetic parameters are described by arbitrary degree polynomials, in practice the difference between the degrees of the kinematic parameters and the thermodynamical ones is equal to 1. The integration in time is done using a defect correction method. As such, there is no hope that the limit solution, if it exists, will be a weak solution of the problem. In order to guaranty this property, we introduce a general correction method in the spirit of the Lagrangian stagered method described in <cit.>, and we prove a Lax Wendroff theorem. The proof is valid for multidimensional version of the scheme, though all the numerical illustrations, on classical benchmark problems, are all one dimensional because we have an easy access to the exact solution for comparison. We conclude by explanning that the method is general and can be used in a different setting as the specific one used here, for example finite volume, of discontinuous Galerkin methods.

## Authors

• 19 publications
• 1 publication
03/28/2021

### Arbitrary Lagrangian-Eulerian hybridizable discontinuous Galerkin methods for fluid-structure interaction

We present a novel (high-order) hybridizable discontinuous Galerkin (HDG...
09/21/2019

### Single-Step Arbitrary Lagrangian-Eulerian Discontinuous Galerkin Method for 1-D Euler Equations

We propose an explicit, single step discontinuous Galerkin (DG) method o...
10/13/2019

### A Conservative Finite Element Method for the Incompressible Euler Equations with Variable Density

We construct a finite element discretization and time-stepping scheme fo...
09/15/2021

### An energy-based discontinuous Galerkin method for dynamic Euler-Bernoulli beam equations

In this paper, an energy-based discontinuous Galerkin method for dynamic...
02/07/2020

### An Eulerian-Lagrangian discontinuous Galerkin method for transport problems and its application to nonlinear dynamics

We propose a new Eulerian-Lagrangian (EL) discontinuous Galerkin (DG) me...
08/03/2021

### The Sod gasdynamics problem as a tool for benchmarking face flux construction in the finite volume method

The Finite Volume Method in Computational Fluid Dynamics to numerically ...
07/11/2019

### Capillary Rise -- A Computational Benchmark for Wetting Processes

Four different numerical approaches are compared for the rise of liquid ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The Euler equations of fluid dynamics are, in their conservative version,

 (1)

As usual, is the density, is the velocity, is the total energy, is the internal energy, and is the pressure. The system is closed by an equation of state, . The simplest is that of calorically perfect gas where

 p=eγ−1

where the ratio of specific heats is a constant.

When the solution is smooth, the system (1) can be equivalently written in nonconservative form as:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂ρ∂t+div(ρu)=0,∂u∂t+(u⋅∇)u+∇pρ=0,∂e∂t+u⋅∇e+(e+p)divu=0. (2)

When the solution is not smooth, the form (2) is meaningless because the differential operators are no longer defined. This is why the form (1) is preferred, in particular in its weak form, see [16]. This fact has a very strong implication for the design of numerical schemes applied to (1): the Lax-Wendroff theorem implies and guaranties that a suitable numerical approximation should be written in term of flux.

However, the form (2) is better suited for engineering purposes, since one has a direct access to the velocity and the internal energy. Hence a rather natural question is to ask how to discretise the Euler equations directly from (2), and still have convergence to the correct weak solutions, at least formally.

One obvious way to achieve this, is to start from a conservation approximation of (1), and by simple algebraic manipulations which amount to multiply the numerical scheme by approximations of

 ⎛⎜ ⎜ ⎜⎝100uρ0u22ρu1⎞⎟ ⎟ ⎟⎠. (3)

lead to a scheme directly working on the primitive variable. This ”new” scheme is equivalent to the original one.

This is not exactly the question we want to address here. We are interested in designing locally conservative approximations of (2) for which the thermodynamic variables are approximated in while the velocity is globally continuous. This can be seen as an Eulerian version of the Lagrangian schemes designed in [6, 11] or [14] and the related works by these authors. A similar question has been addressed by Herbin and co-authors, see for example [15, 19, 17] in the finite volume context. In this references, the authors describe a class of numerical schemes where the thermodynamic variables and the velocity are piecewise constant but logically described on a staggered mesh. They show the convergence towards the weak solution. The scheme can be partially implicit, so that in the low Mach number limit the scheme ”degenerates” to a Mac-type scheme, see [18]. The schemes are second order in time and space.

In this paper, we describe a different technique which allows to reach an arbitrary order of accuracy, both in time and space. Before describing this method, we shortly review the class of residual distribution schemes that will be the main tool we use, thanks to some reinterpretation. Then we describe the scheme, and explain why it is locally conservative. Last, we show a variant of the Lax-Wendroff theorem that is adapted to our setting.

## 2 A high-order nonconservative approach

We have in mind a numerical approximations where the variables are piece-wise polynomial in simplex. We also assume that the velocity is globally continuous, in contrast to Discontinuous Galerkin(DG)–like approximations. This constraint is motivated by the choice that we want to extend the technique of [6], where a Petrov Galerkin technique, inspired from [7] and the reference therein. If nothing special is done, we need to invert a mass matrix. This can be cumbersome, and even impossible if we want to extend the techniques of [5] because the equivalent of the mass matrix changes at every time step. This is why we use a particular time stepping that we describe now. It relies on series of Euler forward type of discretisation. This is indeed the essential point: if one prefers to forget the globally continuous methods, rely on a DG–like approach, and use a Strong Stability Preserving (SSP) Runge-Kutta approach, one can extend our correction technique and build schemes that converge to a weak solution of the problem, starting from (2).

### 2.1 Time stepping approach

Consider a hyperbolic system in the form

 ∂U∂t+L(U)=0. (4)

For time discretization, we want to get a high order accurate approximation. To do so, we will use the Deferred Correction (DeC) approach. The aim of DeC schemes is to avoid implicit methods, without losing the high order of accuracy of a scheme. The high order method that we want to approximate will be denoted by . To use the DeC procedure, we also need another method, which is easy and fast to be solved with low order of accuracy . The DeC algorithm is providing an iterative procedure that wants to approximate the solution of the scheme in the following way:

 L1(U(1))=0 (5)
 L1(U(k))=L1(U(k−1))−L2(U(k−1)),withk=2,..,K, (6)

where is the number of iterations that we compute. We need as many iterations as the order of accuracy that we want to reach. We know from [12]:

###### Proposition 2.1.

Let and be two operators defined on , which depend on the discretization scale such that

• is coercive with respect to a norm, i.e., independent of , such that for any we have that

 α1∣∣∣∣U−V∣∣∣∣≤∣∣∣∣L1(U)−L1(V)∣∣∣∣,
• is Lipschitz with constant uniformly with respect to , i.e., for any

 ∣∣∣∣(L1(U)−L2(U))−(L1(V)−L2(V))∣∣∣∣≤α2Δ∣∣∣∣U−V∣∣∣∣.

We also assume that there exists a unique such that Then, if the DeC is converging to and after iterations the error is smaller than .

Let us show how to use this for solving the problem (1) or (2) on domain , . For now we forget the question of local conservation and

 L(U)=⎛⎜⎝div(ρu)div(ρu⊗u+pI)div((E+p)u)⎞⎟⎠

for the conservative form and we define

 L(U)=⎛⎜ ⎜ ⎜ ⎜⎝div(ρu)(u⋅∇)u+∇pρu⋅∇e+(e+p)divu⎞⎟ ⎟ ⎟ ⎟⎠

for the non–conservative form. We will describe the procedure in two steps. First, we consider the case of a scalar problem, and then we look at (1) or (2). The reason is that in (2) not all the variables play the same role, contrarily to (1), and it is easier to start with a system with one variable. We will proceed by forgetting the question of local conservation.

### 2.2 Scalar case

#### 2.2.1 Trial space

We consider a triangulation of made of non overlapping simplex that are generically denoted by . We assume the triangulation to be conformal, and define

 Vh={v∈L2(Ω) such that for any K,v|K∈Pk(Rd)}⊂L2(Ω)

where as usual, is the set of polynomials in of degree less or equal to . We also define

 Wh=Vh∩C0(Ω).

In each element

, a polynomial is defined by a set of degrees of freedom, for example the Lagrange points. We denote by

a generic degree of freedom. Here, for reasons that will be more clear later, we expand the polynomials in term of Bézier polynomials.

• One dimensional elements: In the element , we consider the barycentric coordinates

 λ1(x)=xi+1−xxi+1−xi,λ2(x)=x−xixi+1−xi=1−λ1.

If , the restriction of is defined in the element as follows, if , then the Bézier form vanishes. We describe the two families of Bézier form we will need:

• Linear: The degrees of freedom are the vertices, so

 φ(1)i=λ1,φ(1)i+1=λ2, and% σ=xi or xi+1 here.
• Quadratic: the degrees of freedom are identified to the vertices , and the mid-points .

 φ(2)σ(x)=⎧⎪ ⎪⎨⎪ ⎪⎩λ21 if % σ=xi2λ1λ2 if σ=xi+1/2λ22 if σ=xi+1.
• Multidimensional elements: we only describe the 2D cases, with triangles, but similar things are obtained for quadrangles, or 3D simplex. A triangle is made of 3 vertices denoted by , and . The barycentric coordinates with respect to the vertices are denotes by , and .

• Linear: the degrees of freedom are the vertices and and .

• Quadratic: The degrees of freedom are the 3 vertices , and , and the midpoint of the edges:

 σ4=σ1+σ22,σ5=σ2+σ32,σ6=σ2+σ12.

The Bézier polynomials are:

 φσi=Λ2i for i=1,2,3,φσ4=2Λ1Λ2,φσ5=2Λ3Λ2,φσ6=2Λ1Λ3.

Then, considering or , for any , we expand as

 u|K=∑σ∈KuσφKσ

where is any of the linear, quadratic (or more) function defined above. If then we have the expansion

 u=∑K∑σ∈KuKσφKσ,

and if , we can expand as

 u=∑σuσφσ.

With some abuse of language, we will use the second notation throughout this paper, depending if we see per element or more globally.

#### 2.2.2 Test space

As we said earlier, we rely on a Petrov-Galerkin approach. This means that the test functions will belong to a finite dimensional subspace of that can also be described by the degrees of freedom : we can identify functions of that are indexed by the and span this space. We denote them by . For example, in the SUPG method, we define in each by: for ,

 Ξσ(x)=φσ(x)+hK(∇UL(U)τK)⋅∇φσ(x).

Here is the diameter of and is a positive matrix. In [5, 13, 3, 21, 2], examples are given where depends on the solution, in order to get stability. In all the examples we are considering, the support of is that of .

#### 2.2.3 Description of the time discretisation

The first step is to define . Let’s integrate the system (4):

 ∫Ω∫tn+1tnΞσ(∂U∂t+L(U))dtdx=0,∫tn+1tn∫ΩΞσ(∂U∂t+L(U))dtdx=0.

Using the locality of the basis function, we have

 ∫Ω∫tn+1tnΞσ(∂U∂t+L(U))dtdx=∑K,σ∈K∫K∫tn+1tnΞσ(∂U∂t+L(U))dtdx,

and

 ∫Ω∫tn+1tnΞσ(∂U∂t+L(U))dtdx=∑K,σ∈K∫K∫tn+1tnΞσ(∂U∂t+L(U))dtdx,

so we will define the restriction of on . Then we introduce sub-time steps between and : and for any degree of freedom the corresponding approximations for . Of course, is the approximation at time

. This being done, we will consider Lagrange interpolation in time and will integrate exactly in time. More details can be found in

[7]. We describe here two examples as illustrations.

In the first one, we take , so that there is no intermediate points, and we simply get the Crank Nicholson scheme:

 ∫K∫tn+1tnΞσ(∂U∂t+L(U))dtdx=∫KΞσ(Un+1−Un)dx+Δt2(∫KΞσL(Un)dx+∫KΞσL(Un+1)dx) (7)

If we look for third order in time, we consider the mid point time step, and get

 (8)

The definition of is somewhat formal, we replace it by some approximation to be defined later, the only constraint is that we have the relation

 ∑σ∈KΦKσ(U)=∫K˜L(U)dx. (9)

The precise definition of the term depends on wether or .

• If , then and the integration is done via quadrature formula of order at least equal to the polynomial degree,

• If , then may be discontinuous across the faces of , we only consider a consistent approximation. If is in conservation form, then we set

 ∫K˜L(U)dx:=∫∂K^fndγ

where is a numerical flux. If instead, is not in conservation form, for example

 L(U)=d∑j=1aj(U)∂U∂xj

and one of the is not a jacobian, we simply set

 ∫K˜L(U)dx:=∫K(d∑j=1aj(U)∂U∂xj)dx

which is evaluated by means of quadrature formula.

We set, for (7), and define

as a vector indexed by the degrees of freedom, and the components are:

 L2σ(U1,U0)=∑K,σ∈K(∫KΞσ(U1−U0)dx+Δt2(ΦKσ(U0)+ΦKσ(U1)L2σ,K(U1,U2)) (10a) and for (8), U(tn)≈U0,U(tn+1/2)≈U1,U(tn+1)≈U2 and define L2σ by L2(U1,U2,U0)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∑K,σ∈K(L2,Iσ,K∫KΞσ(U1−U0)dx+Δt(524ΦKσ(U0)+13ΦKσ(U1)−124ΦKσ(U2))∑K,σ∈K(∫KΞσ(U2−U0)dx+Δt(16ΦKσ(U0)+13ΦKσ(U1)+16ΦKσ(U2)L2,IIσ,K) (10b) So again we have a ”vector” of ”components” (L2,Iσ,K,L2,IIσ,K)T.

The operator is defined by replacing the temporal terms by

 ∫Kφ(Un+ξ−Un)dx

and then by ”lumping” the mass matrix: this is why we use Bézier approximation since we are sure that the lumped mass is non zero because it is

 Cσ=∫Ωφσdx>0.

We set (it does not depend on ) and set, for (7)

 L1σ,K(U1,U0)=CK(U1,σ−U0,σ)+Δt∑K,σ∈KΦKσ(U0)=0 (11a) and similarly, for (8), L1σ,K(U1,U2,U0)=⎛⎜ ⎜ ⎜⎝CK(U1,σ−U0,σ)+Δt2∑K,σ∈KΦKσ(U0)=0CKσ(U2,σ−U0,σ)+Δt∑K,σ∈KΦKσ(U0)=0⎞⎟ ⎟ ⎟⎠ (11b)

The generalisation to formally higher order temporal schemes is obvious. Last, one can show that the two operators satisfies the requirements of the lemma 2.1.

The DeC iteration is then defined by

 L1(U(k+1))=L1(U(k))−L2(U(k)).

After simplifications, this gives the following iteration: knowing , we get by, for all , we have for (7),

 Cσ(U(k+1)1,σ−U(k)1,σ)+∑K,σ∈K[∫KΞσ(U(k)1−U0)dx+Δt2(ΦKσ(U0)+ΦKσ(U(k)1))] (12)

and for (8),

 Cσ(U(k+1)1,σ−U(k)1,σ)+∑K,σ∈K[∫KΞσ(U(k)1−U0)dx+Δt(524ΦKσ(U0)+13ΦKσ(U(k)1)−124ΦKσ(U(k)2))]Cσ(U(k+1)2,σ−U(k)2,σ)+∑K,σ∈K[∫KΞσ(U(k)2−U0)dx+Δt(16ΦKσ(U0)+13ΦKσ(U(k)1)+16ΦKσ(U(k)2))] (13)

For now, this is all we need.

### 2.3 Case of system (2).

We describe the residuals, and develop the method first for the case of the second order in time, since it is simpler. The generalisation to higher order is straightforward and done in the paragraph that follows.

We assume that the computational domain is covered by non-overlapping simplices The velocity field belongs to a kinematic space of finite dimension; it has a basis denoted by , where is the set of kinematic degrees of freedom with the total degrees of freedom given by . The thermodynamic quantities such as the internal energy, the density and the pressure belong to a thermodynamic space ; this space is also finite dimensional and its basis is . The set is the set of thermodynamical degrees of freedom with the total degrees of freedom . The kinematic space is formed by the quadratic Bernstein elements, while the thermodynamic space has a piecewise-linear basis. The velocity field is approximated by

where the are the linear/quadratic Bézier, and the density, the pressure and the internal energy, are given by

 ρ(x,t)=∑σE∈DEρσE(t)φσE(x),p(x,t)=∑σE∈DEpσE(t)φσE(x),e(x,t)=∑σE∈DEeσE(t)φσE(x).

where the are the piecewise constant/linear per elements functions. Note that the degrees of freedom for the velocity are assumed to be globally continuous, so in , while the thermodynamical ones are discontinuous across the boundary of elements, so in .

We can rewrite the Euler equations (2) in the following form

 ∂U∂t+A⋅∇U=0.

The only thing to do is to describe how the method of the previous section adapt to this case, and this amounts to describing the general structure residuals, that is how (9) is written. Since the velocity is globaly continuous, we write

 ∫K˜L(U)dx=∫K(u⊗u+∇pρ)dx

where and . Since we are on , these are simply polynomials, and the integration is carried out by numerical quadrature.

For the density, is in conservation form and we use a numerical flux :

 ∫K˜L(U)dx=∫∂K^fndγ.

Last, for the internal energy, we write

 ∫K˜L(U)dx=∫K(u⋅∇e+(e+p) div u)dx

and again we use quadrature formula.

In the numerical section, we will describe the residuals that we use.

## 3 A discussion on conservation

### 3.1 A set of sufficient conditions to achieve convergence to a weak solution

Again, to simplify the notations, we focus on the second order case, but the extension to the more general case is straightforward. In the appendix 6, we show a Lax Wendroff theorem for this type of discretisation. What we do here is to show how to go from the system in non–conservative form to the one in conservation form.

There is nothing to do for the density because it is already in conservation form and the standard proof [16] applies. There is no need to repeat it here since the proof we give for the momentum and the total energy, modulo some complications, is essentially the same.

Let us first look at the momentum. Considering a test function , we define the value of at time at the centroid of , and consider the following approximation of that we still denote by :

 ψ(x,t)=∑KψnK1K for t∈[tn,tn+1[.

Then we consider

 (14)

Then, we can write, introducing and ,

 ∫Kρ(p+1)(u(p+1)−u(p))dx=∑σV∈KΔuσV∫Kρ(p+1)φσVdx

and

 ∫Ku(p)(ρ(p+1)−ρ(p))dx=∑σE∈KΔρσE∫Ku(p)φσEdx.

Hence, (14) is rewritten as:

 (15)

where we have set for simplicity

 ωρ,p+1,KσV:=∫Kρ(p+1)φσVdx|CσV|% and ωu,p,KσV:=∫Ku(p)φσEdx|CσE|.

For the velocity, we have:

 |CσV|(u(p+1)σV−u(p)σV)+Δtn∑K,σV∈KΦuσV,K=0

where, for the scheme (7),

 ΦuσV,K=12(ΦuσV,K(U(p))+ΦuσV,K(U(0)))

and for the density, we have

 |CσE|(ρ(p+1)σE−ρ(p)σE)+Δtn∑K,σE∈KΦρσE,K=0.

We note that the sum reduces to one term, hence

 |CσE|(ρ(p+1)σE−ρ(p)σE)+ΔtnΦρσE,K=0

where again

 ΦρσE=12(ΦρσV(U(p))+ΦρσV(U(0)))

and is the element such that

Using these relations in (15), we get

 (16)

We see that

 ∑KψnK[∑σV∈Kωρ,p+1,KσV{∑K′,σV∈K′ΦuσV,K′}+∑σE∈Kωρ,p+1,KσVΦρσE,K]=∑σV∈K[∑K′,σV∈K′∑K,σV∈K∩K′ωu,p,KσVψnKΦuσV,K′]I+∑K∑σE∈KψnKωu,p,KσVΦρσE

Then, can be written as:

 I=ψnσV∑K′,σV∈K′{∑K,σV∈K∩K′ωρ,p+1,KσV}ΦuσV,K′+[∑K′,σV∈K′∑K,σV∈K∩K′ωρ,p+1,KσV(ψnK−ψnσV)ΦuσV,K′]:=DσV

so that, setting

 ∑KψnK[∑σVωρ,p+1,KσV{∑K′,σV∈K′ΦuσV,K′}+∑σE∈Kωu,p,KσVΦρσE,K]=∑K∑σV∈KψnσVωρ,p+1,KσVΦuσV,K+∑σVDσV+∑KψnK∑σE∈Kωu,p,KσVΦρσE,K=∑KψnK[∑σV∈Kωρ,p+1,KσVΦuσV,K+∑σE∈Kωu,p,KσVΦρσE,K]+∑K(∑σV∈K(ψnσV−ψnK)ωu,p,KσVΦuσV,K:=FK+∑σV∈KDσV)

so that we get the master equation:

 (17a) with FK=∑σV∈K(ψnσV−ψnK)ωu,p,KσVΦuσV,KDK=∑K′,σV∈K′[∑K,σV∈K∩K′ωu,p,KσV(ψnK−ψnσV)ΦuσV,K′]ωρ,p+1σV=∑K,σV∈K∫Kρ(p+1)φσVdx|CσV|,ωu,p,KσV=∫Ku(p)φσEdx|CσE| (17b)

Then let us look at the total energy. The first thing to do is to remark that (with similar notations as before) that

 Δρu2=u(p+1)⋅Δρu+ρ(p)u(p)⋅Δu

which combined to

 Δρu=ρ(p+1)Δu+u(p)Δρ

gives:

 Δρu2=(ρ(p+1)u(p+1)+ρ(p)u(p))⋅Δu+u(p+1)⋅u(p)Δρ.

To simplify, we will set

 ˜m=ρ(p+1)u(p+1)+ρ(p)u(p)2,˜q2=u(p+1)⋅u(p).

Using these relations, we see that

 ∑KψK∫KΔEdx=∑KψK(∫KΔedx+∫K~m⋅Δudx+12