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

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.

## Authors

• 22 publications
• 1 publication
12/06/2021

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

Invariant finite-difference schemes are considered for one-dimensional m...
11/22/2021

### Memory erasure with finite-sized spin reservoir

Landauer's erasure principle puts a fundamental constraint on the amount...
09/05/2019

### A Harten's Multiresolution Framework for Subdivision Schemes

Harten's Multiresolution framework has been applied in different context...
11/09/2019

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

This work proposes and investigates a new model of the rotating rigid bo...
10/07/2021

### Lagrangian Neural Network with Differentiable Symmetries and Relational Inductive Bias

Realistic models of physical world rely on differentiable symmetries tha...
06/22/2022

### Locally conservative and flux consistent iterative methods

Conservation and consistency are fundamental properties of discretizatio...
04/21/2022

### 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

 ∂u∂t+div f(u)=0, (1)

defined on a space-time domain , where

 u=(ρ,ρv,E) and f=(ρv,ρv⊗v+pId ,(E+p)v)=(f1,…,fd),d=1,2,3.

The total energy writes

 E=ρe+12ρv2,

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:

 ∂S∂t+ div (ρS)=0(≤0 in % general.) (2)

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

 Tds=d(eρ)−pρ2dρ.

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

 ∂ρv∂t+ div (ρv⊗v)+∇p=0,

we write, in 2D,

 ∂ρv∂t+∂f1∂x+∂f2∂y=0,f1=(ρv2x+pρvxvy)=(f11f21),f2=(ρvxvyρv2y+p)=(f12f22),

We note that .

 ∂ρx∧v∂t=x∧∂ρv∂t=−x∧(∂f1∂x1+∂f2∂x2)=−x1(∂f21∂x1+∂f22∂x2)+x2(∂f11∂x1+∂f12∂x2)=−∂∂x1(f21x1−x2f11)−∂∂x2(x1f22−f12x2)+(f21−f12),

Hence we have:

 ∂J∂t+∂∂x1(f21x1−x2f11)+∂∂x2(x1f22−f12x2)=0. (3)

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

and since

 f32=f23=ρu2u3,f13=f31=ρu3u1,f21=f12=ρu1u2,

we obtain the additional conservation relation:

 ∂J1∂t+∂∂x1(x2f31−x3f21)+∂∂x2(x2f32−x3f22)+∂∂x3(x2f33−x3f23)=0,∂J2∂t+∂∂x1(x3f11−x1f31)+∂∂x2(x3f12−x1f32)+∂∂x3(x3f13−x1f33)=0,∂J3∂t+∂∂x1(x1f21−x2f11)+∂∂x2(x1f22−x2f12)+∂∂x3(x1f23−x2f13)=0. (4)

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

 |Cσ|=∫Ωφσdx>0. (5)

This is not necessarily true for Lagrange polynomials.

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

 u=∑σuσφσ,

Then the solution at time can be written as follows

 un=∑σunσφσ.

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)

 div f(u)=0. (6)

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))

 ΦK(u)=∫∂Kf(u)⋅ndx(=∫K div f(u)dx), (7)
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

 ΦK(u)=∑σ∈KΦKσ, ∀K∈Ω, (8)
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))

 ∑K,σ∈KΦKσ(u)=0, ∀σ∈Ω. (9)

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

 ∑K,σ∈KΦKσ,x(u)+∑Γ,σ∈γΦΓσ,x(u)=0, (10)

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

 ∑K,σ∈KΦKσ,x(u)=∫∂Kf(u).n, (11a) and for the boundary condition (here Dirichlet in weak form) ∑Γ,σ∈ΓΦΓσ,x(u)=∫Γ(Fn(u,g)−f(u).n). (11b)

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

 (12)

where

 ∑σ∈KβKσ(u)=Id(=1%inthescalarcase).

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

 ∑σ∈KξKσ=0.

Then we define with, for ,

 ξσ(x)=ξKσ.

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

 −∫Ω∇v⋅f(u)dx+∫∂ΩvFndγ+∑K∫Kξ(u,v)div f(u)dx=0,

where

 ξ(x)=∑σvσξσ.

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

 ∫stn∫Ω(v+ξ(v))∂u∂tdxdt−∫stn∫Ω∇v⋅f(u)dxdt+∫stn∫∂ΩvFndγdt+∑K∫stn∫Kξ(v) div f(u)dxdt=0.

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

 ∫Ω(v+ξ(v))(u(x,s)−u(x,tn))dx−∫stn∫Ω∇v⋅f(u)dxdt+∫stn∫∂ΩvFndγdt+∑K∫stn∫Kξ(v) div f(u)dxdt=0,

and then, for any ,

This suggest to introduce the space-time fluctuations

 ΦKσ(u,s)=∫K(φσ+ξσ)(u(x,s)−u(x,tn))dxΦKσ,t(u,s)+∫stnΦKσ,x(u,s)dt (13)

where

 ΦKσ,x(u,s)=−∫K∇φσ⋅f(u)dx+∫Kξσdiv f(u)dx+∫∂Kφσf(u)⋅ndγ.

### 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

 tn=t(0)

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

 I0(s)=u(0) if s∈[tn,tn+1].

Another choice could have been made

 I0(s)=u(l) if s∈[t(l),t(l+1)[.

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 ,

 [L1(U;un)]σ,(l):=|Cσ|uσ,(l)−|Cσ|uσ,(0)+∑K,σ∈K∫t(l)t(0)I0(ΦKσ,x(U),s)ds, (14)

where

 I0(ΦKσ,x(U),s)=I0(ΦKσ,x(u(0)),ΦKσ,x(u(1)),…,ΦKσ,x(u(M)),s),

i.e. here

 [L1(U;un)]σ,(l) =|Cσ|uσ,(l)−|Cσ|uσ,(0)+(t(l)−t(0))∑K,σ∈KI0(ΦKσ,x(u(0)))

This is a first order explicit approximation in time.

• A high order approximation: for any and ,

 [L2(U;un)]σ,(l):=∑K,σ∈K∫t(l)t(0)IM(ΦKσ(U),s)ds, (15)

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

 ∫t(l)t(0)IM(ΦKσ,x(U))ds=M∑k=0θlkΦKσ,x(U). (16)

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

 [L2(U;un)]σ,(l)=0, (17)

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

 (18)

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

 |Cσ|(U(p+1)σ−U(p)σ)=−∑K,σ∈KΦKσ(U(p)), (19)

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,

 α1∥U−V∥≤∥L1(U)−L1(V)∥,
• is uniformly Lipschitz continuous with Lipschitz constant , i.e., there exists independent of , such that for any U and V,

 ∥(L1(U)−L2(U))−(L1(V)−L2(V))∥≤α2Δ∥U−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 :

 L1−L2=O(Δt). (20)

## 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

 |Cσ|(xσ∧(m(p+1)σ−m(p)σ))+∑K,σ∈Kxσ∧ΦKm,σ(U(p))+∑Γ,σ∈Γxσ∧ΦΓm,σ(U(p))=0,

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:

 ∑σ∈Kxσ∧ΦKm,σ(U(p))=∫K(J(p)−J(0))dx+Δt∫∂KI2(G(U(p)))⋅ndγ:=ΦKJ (21)
 ∑σ∈Γxσ∧ΦΓm,σ(U(p))=Δt∫∂K(^G−I2(G(U(p)))⋅n)dγ:=ΦΓJ (22)

with the flux for the angular momentum and

 G(U(p))=12(G(U(p))+G(U(0)))

or

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

 (23)
 ∑σ∈ΓΦΓm,σ=Δt∫ΓFm(U(p))⋅ndγ (24)

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

 ∫Kf(x)dx≈|K|3∑σ∈Kfσ,

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:

 ∑σ∈KrKσ=0∑σ∈Kxσ∧rKσ=ΦKJ−∑σ∈Kxσ∧ΦKm,σ(U(p)):=ΨK (25)

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

 ∑σ∈ΓrΓσ=0∑σ∈Γxσ∧rΓσ=ΦΓJ−∑σ∈Γxσ∧ΦΓm,σ(U(p)):=ΨΓ (26)

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

 (x2−x1)∧r2+(x3−x1)∧r3=Ψ∈R

We define

 r2=r(x3−x1),r3=r(x1−x2),r∈R

and get

 r(det(x2−x1,x3−x1)+det(x3−x1,x1−x2))=Ψ,

i.e. since , we have:

 r=Ψ4|T|,r2=r(x3−x1),r3=r(x1−x2),r1=r(x2−x3). (27)

### 3.2 Solution for tetrahedrons

Again, , so that

 4∑j=2(xj−x1)∧rj=Ψ.

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

 r2=α2x13+β2x14r3=α3x12+β3x14r4=α4x12+β4x13

This gives:

 Ψ=α2x21∧x13+β2x21∧x14+α3x31∧x12+β3x31∧x14+α4x41∧x12+β4x41∧x13

Then, since ,

 Ψ⋅x21=β3det(x21,x31,x14)+β4det(x21,x41,x13)Ψ⋅x31=β2det(x31,x21,x14)+α4det(x31,x41,x12)Ψ⋅x41=α2det(x41,x21,x13)+α3det(x41,x31,x12)

Now,

 3|T|=det(x21,x31,x14)=−det(x21,x41,x13)=−det(x31,x21,x14)=det(x31,x41,x12)=det(x41,x21,x13)=−det(x41,x31,x12)

so that we have:

 Ψ⋅x21=3|T|(β3−β4)Ψ⋅x31=3|T|(−β2+α4)Ψ⋅x41=3|T|(α2−α3).

This suggests to assume , , , i.e

 β3=−β4=Ψ⋅x216|T|α4=−β2=Ψ⋅x316|T|α2=−α3=Ψ⋅x416|T|

Remembering that , we see that

 r2=16|T|(Ψ⋅x41x13−Ψ⋅x31x14)=16|T|Ψ∧(x14∧x31)r3=16|T|(Ψ⋅x21x14−Ψ⋅x41x12)=16|T|Ψ∧(x21∧x14)r4=16|T|(Ψ⋅x31x12−Ψ⋅x21x13)=16|T|Ψ∧(x31∧x12).

Last,

 x14∧x31+x21∧x14+x31∧x12=x14∧x31+x14∧x12+x31∧x12=x14∧x31+(x14+x31)∧x12=x14∧x31+x34∧x12=(x13+x34)∧x31+x34∧x12=x34∧x31+x34∧x12=x34∧x32

and then we have:

 (28)

## 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

 U(p+1)σ=U(p)σ+δU(p)σ.

Then

 U(p+1)=∑σU(p+1)σφσ.

The angular momentum, integrated, is

We introduce the notations

 yσ=1|Cσ|∫ΩxBσdx,
 zKσ=1|K|∫KxBσdx,

so that we can rewrite the total kinetic momentum as:

 ∫Ωx∧m(p+1)dx=∑σ|Cσ|yσ∧m(p+1)σ,

and then we can write the update:

 ∑σ|Cσ|yσ∧m(p+1)σ=∑σ|Cσ|yσ∧m(p)σ−∑σyσ∧(∑K,σ∈KΦKm,σ)=∑σ|Cσ|yσ∧m(p)σ−∑σyσ∧(∑K,σ∈KΦKm,σ+∑Γ,σ∈ΓΦΓm,σ)=∑σ|Cσ|yσ∧m(p)σ−∑