1 Introduction
We are interested in the solution of the Euler equations for compressible flows
(1) |
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:
(2) |
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:
(3) |
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:
(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.
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(5) |
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)
(6) |
![]() |
![]() |
![]() |
The main steps of the residual distribution approach can be summarized as follows
-
For any element , compute a fluctuation term (total residual) (see figure 1(a))
(7) -
For every DOF within the element K, define the nodal residuals as the contribution to the fluctuation term (see figure 1(b)) such that
(8) -
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))
(9)
We can write a similar formulation for the boundary conditions. For any DOF we can split (6) into the internal and boundary contributions
(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
(11a) | |||
and for the boundary condition (here Dirichlet in weak form) | |||
(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
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
where
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
(13) |
where
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 asAnother 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 ,
(14) where
i.e. here
This is a first order explicit approximation in time.
-
A high order approximation: for any and ,
(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
(16)
Ideally, we would like to solve for each and each ,
(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:
-
Set ,
-
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
(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,
-
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
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:
(21) |
(22) |
with the flux for the angular momentum and
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) |
(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
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:
(25) |
Since (22) is not true, as above, we introduce a vectorial correction of the momentum residual . We need:
(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
We define
and get
i.e. since , we have:
(27) |
3.2 Solution for tetrahedrons
Again, , so that
In order to simplify the notations, we introduce . We are looking for
This gives:
Then, since ,
Now,
so that we have:
This suggests to assume , , , i.e
Remembering that , we see that
Last,
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
Then
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: