Conservative scheme compatible with some other conservation laws: conservation of the local angular momentum

by   Rémi Abgrall, et al.

We are interested in building schemes for the compressible Euler equations that are also locally conserving the angular momentum. We present a general framework, describe a few examples of schemes and show results. These schemes can be of arbitrary order.



page 14

page 15

page 22

page 23

page 24

page 25

page 26


Invariant finite-difference schemes with conservation laws preservation for one-dimensional MHD equations

Invariant finite-difference schemes are considered for one-dimensional m...

Memory erasure with finite-sized spin reservoir

Landauer's erasure principle puts a fundamental constraint on the amount...

A Harten's Multiresolution Framework for Subdivision Schemes

Harten's Multiresolution framework has been applied in different context...

The rotating rigid body model based on a non-twisting frame

This work proposes and investigates a new model of the rotating rigid bo...

Lagrangian Neural Network with Differentiable Symmetries and Relational Inductive Bias

Realistic models of physical world rely on differentiable symmetries tha...

Locally conservative and flux consistent iterative methods

Conservation and consistency are fundamental properties of discretizatio...

A novel formulation for the evolution of relativistic rotating stars

We present a new formulation to construct numerically equilibrium config...
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

We are interested in the solution of the Euler equations for compressible flows


defined on a space-time domain , where

The total energy writes

and the pressure is given by the equation of state . This function satisfies the usual convexity requirements. Equation (1) is complemented by initial and boundary conditions.

In addition to (1), the solution is known to satisfy other conservation relation. In the case of smooth flows, it is known that satisfies:


where and is the specific entropy defined by the Gibbs relation:

In order to fully define the entropy, we need a complete equation of state, from which one can deduce from standard thermodynamical relations.

There are indeed other relations, for example the preservation of the angular momentum. Defining , and from

we write, in 2D,

We note that .

Hence we have:


In 3 dimensions, we have a similar results: a direct computation of the time evolution of the angular momentum equation give:

and since

we obtain the additional conservation relation:


The weak form of the relations (3) and (4) are also satisfied by . To see this, in the 2D case, one has to take the test functions and where has compact support and is in space-time domain, and subtract to get the weak form of (3). The same method also works in 3D.

Remark 1.1.

The PDEs (3) and (4) are translational invariant. This means that if satisfies (3) or (4) in the frame of coordinates , then will satisfy (3) or (4) in the frame of coordinates .

Developing schemes that locally preserves the angular momentum has recently attracted some attention, see e.g [9, 10] and the reference therein. In [5, 7, 1, 4] has been developed a simple technique that enables to correct an initial scheme so that an additional conservation is satisfied. This has been applied to some non conservative form of the Euler equation, also to satisfy entropy conservation when this is relevant, and to various schemes: residual distribution schemes and discontinuous Galerkin ones. The goal of this paper is to show how one can extend the method to the problem of angular preservation, for second and higher order schemes, in a fully discrete manner, both in time and space.

The paper is organised as follows. In section 2 we briefly introduce a high order residual distribution scheme for the unsteady version of (1). In section 3 we construct a second order scheme for (1) that preserves locally the angular momentum; then we show this is possible for triangles and tetrahedrons. In section 4 we provide a possible solution to conservation of momentum for higher than second order schemes. In section 5 we study the validity and effectiveness of the proposed strategy, considering three different problems. In section 6 we describe how to extend this to discontinuous Galerkin schemes. Finally, section 7 provide some conclusions and future perspectives.

2 Numerical scheme

2.1 Some generalities

The goal is to construct a scheme for (1) that preserves locally the angular momentum. For this we use a high order residual distribution scheme for the unsteady version of (1) that we briefly describe now, see [6, 2] for more details.

We are given a tessellation of by non overlapping simplex denoted generically by , and more precisely . We assume the mesh to be conformal which means that is either empty, reduce to a full face or a full edge or a vertex. From this, we consider . We are given a basis of , where the

are the degrees of freedom (DOFs). The restriction of

to any is assumed here to be a Bézier polynomial, see appendix C for the notations we use. The reason is that for any , we have


This is not necessarily true for Lagrange polynomials.

We are looking for an approximation of the solution of the form

Then the solution at time can be written as follows

The coefficients , are chosen by a numerical method and we use the following technique to calculate the coefficients in the approximation.

2.2 Residual distribution scheme for steady problems

First, we consider a steady version of system (1)

(a) Compute the total residual.
(b) Compute the nodal residuals.
(c) Collect all the nodal residual contributions.
Figure 1: Illustration of the main steps of the residual distribution approach.

The main steps of the residual distribution approach can be summarized as follows

  1. For any element , compute a fluctuation term (total residual) (see figure 1(a))

  2. For every DOF within the element K, define the nodal residuals as the contribution to the fluctuation term (see figure 1(b)) such that

  3. The resulting scheme is obtained for each DOF by collecting all the nodal residual contributions from all elements K surrounding a node (see figure 1(c))


We can write a similar formulation for the boundary conditions. For any DOF we can split (6) into the internal and boundary contributions


where and are the residuals corresponding to the spatial discretization and is an edge on the boundary of . If we suppose the Dirichlet boundary condition on , for any K and , and satisfy the following conservation relations

and for the boundary condition (here Dirichlet in weak form)

In the appendix A, we provide several examples of such fluctuations. For now, and in order to simplify the discussion on the approximation of the unsteady problem (1), we introduce a ”variational” form of the fluctuation. In all the known cases, we can write as



The distribution coefficients can be scalar in the scalar case, and matrices in the system case.

Then, we can rewrite the scheme in a Petrov-Galerkin fashion:

where . To avoid confusion, we will write as , thus removing the functional dependency of this term with respect to the solution . We remove it, but we do not forget it!

We note that the functions satisfies

Then we define with, for ,

Defining , we can formally rewrite the scheme as: for any , find such that


Here the functional dependence of with respect to exists but is implicit.

2.3 Residual distribution scheme for unsteady problems

Using the results presented above, it becomes possible to describe the unsteady version of the scheme. In order to get a consistent approximation, we simply ”multiply” (1) with the Petrov-Galerkin test function , and integrate in time. Note that will not depend on time. Between and , we get

Since is independent of time, this can be equivalently rewritten as:

and then, for any ,

This suggest to introduce the space-time fluctuations



2.4 Iterative timestepping method

In order to integrate the system, we proceed in two steps. First we introduce sub-time steps in the interval : we subdivide the interval into sub-intervals obtained from a partition

In all our examples, this partition will be regular: , but other choices could have been made, especially for very high order accuracy in time. We construct an approximation of at times , denoted by . Then we introduce

the Lagrange interpolant (in time) of degree

defined from this subdivision. We also introduce the piece-wise constant interpolant defined as

Another choice could have been made

for . The notation

represents the vector

i.e the vector of all the approximations for the sub-steps. Note that and . We need residuals, computed for time and that satisfies, for any ,

We set .

Then we introduce two approximations of (1):

  • A first order approximation in time: for any and ,



    i.e. here

    This is a first order explicit approximation in time.

  • A high order approximation: for any and ,


    After having performed exact integration in time to obtain the approximation for every sub-steps, the time integration can be written in the form


Ideally, we would like to solve for each and each ,


but this is very difficult in general and the resulting scheme derived by operator is implicit. Instead we use a defect correction method and proceed within the time interval as follows:

  1. Set ,

  2. For any , define by


Since is explicit, can be obtained explicitly. In our case, this amounts to a multi-step method where each step writes as


The scheme (18) is completely explicit. One has the following result, see [6]:

Proposition 2.1.

If two operators and depending on the discretization scale , are such that:

  • There exists a unique such that

  • is coercive, i.e., there exists independent of , such that for any U and V,

  • is uniformly Lipschitz continuous with Lipschitz constant , i.e., there exists independent of , such that for any U and V,

    Then if the defect correction method is convergent, and after p iterations the error is smaller than

Therefore if there is a unique solution of (17), then after steps, , so that only steps are needed. Of course this is true thanks to the conditions (5), (14), (15) and the condition that is a first order approximation of :


3 Angular momentum preservation: second order case

It is clear that the residuals depends where the order in time. We can compute the variation of the angular momentum from (19). For any , we have

so that

where here is the momentum component of the residual at the DOF in the element and is the momentum component of the residual at the DOF in the boundary . An easy condition for having local conservation is:


with the flux for the angular momentum and


depending on how the system (1) has been discretized.

Of course, in general, the relation (21) and (22) cannot be satisfied. In order to be satisfied, as well as keeping


with a similar definition of the momentum flux, we will present a perturbation of the momentum residuals. At this level, it is important to provide the quadrature formula. For (21) and (23), we consider

so it is exact for (23) and only approximate for (21), but second order.

To achieve this, following [7], we introduce a perturbation of the momentum residual, , such that the new momentum residual is . We must have:


Since (22) is not true, as above, we introduce a vectorial correction of the momentum residual . We need:


can be derived in the same way.

In dimension and for an element with DOFs, we have equations for unknowns. Since is always larger than , the system has a priori solutions. However, finding a close form formula is element dependent. In the following, we show this is possible for triangles and tetrahedrons. Once these solutions are described, we have a scheme that locally preserves the angular momentum and globally conserves

up to boundary terms.

Remark 3.1 (Translation invariance).

The angular momentum is defined after a frame has been defined, and if one makes a translation of vector , the scheme is translational invariant as in the continuous case.

3.1 Solution for triangular elements

We first have , so

We define

and get

i.e. since , we have:


3.2 Solution for tetrahedrons

Again, , so that

In order to simplify the notations, we introduce . We are looking for

This gives:

Then, since ,


so that we have:

This suggests to assume , , , i.e

Remembering that , we see that


and then we have:


4 Angular momentum preservation: the high order case

The idea is similar, the key point is to characterize the quadrature formula that describes . The additional difficulty is that if we make the geometrical identification of the DOFs with the Greville points. Here to simplify the notations, denotes the residual at evaluated for the momentum, it does not contain the contribution for the density or the energy. In this section we always use this short hand notation, except at the end.

We start again from (19) that we rewrite as


The angular momentum, integrated, is

We introduce the notations

so that we can rewrite the total kinetic momentum as:

and then we can write the update: