A Variational Finite Element Discretization of Compressible Flow

Abstract

We present a finite element variational integrator for compressible flows. The numerical scheme is derived by discretizing, in a structure-preserving way, the Lie group formulation of fluid dynamics on diffeomorphism groups and the associated variational principles. Given a triangulation on the fluid domain, the discrete group of diffeomorphisms is defined as a certain subgroup of the group of linear isomorphisms of a finite element space of functions. In this setting, discrete vector fields correspond to a certain subspace of the Lie algebra of this group. This subspace is shown to be isomorphic to a Raviart–Thomas finite element space. The resulting finite element discretization corresponds to a weak form of the compressible fluid equation that does not seem to have been used in the finite element literature. It extends previous work done on incompressible flows and at the lowest order on compressible flows. We illustrate the conservation properties of the scheme with some numerical simulations.

Appendices

Euler–Poincaré Variational Principle

In this Appendix we first recall the Euler–Poincaré principle for invariant Euler-Lagrange systems on Lie groups. This general setting underlies the Lie group description of incompressible flows recalled in Sect. 2.1 due to [2], in which case the Lie group is $G={\text {Diff}}_\mathrm{vol}(\Omega )$. It also underlies the semidiscrete setting, in which case the Lie group is $G=G_h$. In this situation, however, a nonholonomic constraint needs to be considered, see Appendix B. Then, we describe the extension of this setting that is needed to formulate the variational formulation of compressible flow and its discretization.

1.1 Euler–Poincaré Variational Principle for Incompressible Flows

Let G be a Lie group and let $L:TG \rightarrow {\mathbb {R}}$ be a Lagrangian defined on the tangent bundle TG of G. The associated equations of evolution, given by the Euler-Lagrange equations, arise as the critical curve condition for the Hamilton principle

$$\begin{aligned} \delta \int _0^T L(g(t), \dot{g}(t))\mathrm{d} t=0, \end{aligned}$$

(74)

for arbitrary variations $\delta g$ with $\delta g(0)= \delta g(T)=0$.

If we assume that L is G-invariant, i.e., $L(gh, \dot{g}h)=L(g, \dot{g})$, for all $h\in G$, then L induces a function $\ell :{\mathfrak {g}}\rightarrow {\mathbb {R}}$ on the Lie algebra ${\mathfrak {g}}$ of G, defined by $\ell (u)=L(g,\dot{g})$, with $u=\dot{g} g^{-1}\in {\mathfrak {g}}$. In this case, the equations of motion can be expressed exclusively in terms of u and $\ell$ and are obtained by rewriting the variational principle (74) in terms of $\ell$ and u(t). One gets

$$\begin{aligned} \delta \int _0^T \ell (u(t))\mathrm{d}t=0,\quad \text {for}~\delta u = \partial _t v + [v,u], \end{aligned}$$

(75)

where $v(t)\in {\mathfrak {g}}$ is an arbitrary curve with $v(0)=v(T)=0$. The form of the variation ${\delta u}$ in (75) is obtained by a direct computation using $u= \dot{g} g^{-1}$ and defining $v= {\delta g} g^{-1}$.

In order to formulate the equations associated with (75) one needs to select an appropriate space in nondegenerate duality with ${\mathfrak {g}}$ denoted ${\mathfrak {g}}^*$ (the usual dual space in finite dimensions). We shall denote by $\langle \!\langle \,, \rangle \!\rangle : {\mathfrak {g}}^* \times {\mathfrak {g}}\rightarrow {\mathbb {R}}$ the associated nondegenerate duality pairing. From (75) one directly obtains the equation

$$\begin{aligned} \Big \langle \Big \langle \partial _t \frac{\delta \ell }{\delta u}, v\Big \rangle \Big \rangle + \Big \langle \Big \langle \frac{\delta \ell }{\delta u}, [u,v]\Big \rangle \Big \rangle =0,\quad \forall \; v\in {\mathfrak {g}}. \end{aligned}$$

(76)

In (76), the functional derivative $\frac{\delta \ell }{\delta u}\in {\mathfrak {g}}^*$ of $\ell$ is defined in terms of the duality pairing as

$$\begin{aligned} \Big \langle \Big \langle \frac{\delta \ell }{\delta u}, \delta u \Big \rangle \Big \rangle = \left. \frac{d}{d\epsilon }\right| _{\epsilon =0}\ell (u+\epsilon \delta u). \end{aligned}$$

In finite dimensions, and under appropriate choices for the functional spaces in infinite dimensions, (76) is equivalent to the Euler–Poincaré equation

$$\begin{aligned} \partial _t \frac{\delta \ell }{\delta u}+ {\text {ad}}^*_u\frac{\delta \ell }{\delta u} =0, \end{aligned}$$

where the coadjoint operator ${\text {ad}}^*_u: {\mathfrak {g}}^* \rightarrow {\mathfrak {g}}^*$ is defined by $\langle \!\langle {\text {ad}}^*_um,v\rangle \!\rangle =\langle \!\langle m,[u,v]\rangle \!\rangle$.

For incompressible flows, without describing the functional analytic setting for simplicity, we have $G={\text {Diff}}_\mathrm{vol}(\Omega )$ and ${\mathfrak {g}}= {\mathfrak {X}}_\mathrm{vol}(\Omega )$ the Lie algebra of divergence free vector fields parallel to the boundary. We can choose ${\mathfrak {g}}^*= {\mathfrak {X}}_\mathrm{vol}(\Omega )$ with duality pairing $\langle \!\langle \,, \rangle \!\rangle$ given by the $L^2$ inner product. A direct computation gives the coadjoint operator ${\text {ad}}^*_u m= {\mathbf {P}}(u \cdot \nabla m + \nabla u ^{\mathsf {T}} m)$, where ${\mathbf {P}}$ is the Leray-Hodge projector onto ${\mathfrak {X}}_\mathrm{vol}(\Omega )$. One directly checks that in this case (75) yields the Euler equations (3) for incompressible flows.

1.2 Euler–Poincaré Variational Principle for Compressible Flows

The general setting underlying the variational formulation for compressible fluids starts exactly as before, namely, a system whose evolution is given by the Euler-Lagrange equations for a Lagrangian defined on the tangent bundle of a Lie group G. The main difference is that the Lagrangian depends parametrically on some element $a_0\in V$ of a vector space (the reference mass density $\varrho _0$ in the case of the barotropic compressible fluid, the reference mass and entropy densities $\varrho _0$ and $S_0$ for the general compressible fluid) on which G acts by representation, and, in addition, L is invariant only under the subgroup of G that keeps $a_0$ fixed. If we denote by $L(g, \dot{g}, a_0)$ this Lagrangian and by $a\in V\mapsto a\cdot g \in V$ the representation of G on V, the reduced Lagrangian is defined by $\ell (u, a)= L(g, \dot{g}, a_0)$, where $u= \dot{g}g^{-1}$, $a= a_0\cdot g^{-1}$.

The Hamilton principle now yields the variational formulation

$$\begin{aligned} \delta \int _0^T \ell (u(t),a(t))\mathrm{d}t=0,\quad \text {for}~\delta u = \partial _t v + [v,u] \text { and}~\delta a= - a\cdot v, \end{aligned}$$

(77)

where $v(t)\in {\mathfrak {g}}$ is an arbitrary curve with $v(0)=v(T)=0$. The form of the variation $\delta u$ in (77) is the same as before, while the expression for ${\delta a}$ is obtained from the relation $a= a_0\cdot g^{-1}$.

From (77) and with respect to the choice of a spaces ${\mathfrak {g}}^*$ and $V^*$ in nondegenerate duality with ${\mathfrak {g}}$ and V, with duality pairings $\langle \!\langle \,,\rangle \!\rangle$ and $\langle \,,\rangle _V$, one directly obtains the equations

$$\begin{aligned} \Big \langle \Big \langle \partial _t \frac{\delta \ell }{\delta u}, v\Big \rangle \Big \rangle + \Big \langle \Big \langle \frac{\delta \ell }{\delta u}, [u,v]\Big \rangle \Big \rangle + \Big \langle \frac{\delta \ell }{\delta a},a\cdot v \Big \rangle _V=0,\quad \forall \; v\in {\mathfrak {g}}. \end{aligned}$$

(78)

The continuity equation

$$\begin{aligned} \partial _t a + a \cdot u=0 \end{aligned}$$

arises from the definition $a(t)= a_0 \cdot g(t)^{-1}$. In a similar way with above, (78) now yields the Euler–Poincaré equations

$$\begin{aligned} \partial _t \frac{\delta \ell }{\delta u}+ {\text {ad}}^*_u\frac{\delta \ell }{\delta u} =\frac{\delta \ell }{\delta a}\diamond a, \end{aligned}$$

(79)

where $\frac{\delta \ell }{\delta a}\diamond a\in {\mathfrak {g}}^*$ is defined by $\left\langle \!\left\langle \frac{\delta \ell }{\delta a}\diamond a,v\right\rangle \!\right\rangle = -\Big \langle \frac{\delta \ell }{\delta a},a\cdot v \Big \rangle _V$, for all $v\in {\mathfrak {g}}$. We refer to [15] for a detailed exposition.

For the compressible fluid, in the continuous case we have $G= {\text {Diff}}(\Omega )$ and ${\mathfrak {g}}={\mathfrak {X}}(\Omega )$ the Lie algebra of vector fields on $\Omega$ with vanishing normal component to the boundary. We choose to identify ${\mathfrak {g}}^*$ with ${\mathfrak {g}}$ via the $L^2$ duality pairing. Consider the Lagrangian (65) of the general compressible fluid. Using the expressions ${\text {ad}}^*_u m = u \cdot \nabla m + \nabla u^{\mathsf {T}} m + m {\text {div}} u$, $\frac{\delta \ell }{\delta u}= \rho u$, $\frac{\delta \ell }{\delta \rho }= \frac{1}{2}|u|^2 - e(\rho ) - \rho \frac{\partial e}{\partial \rho }+\eta \frac{\partial e}{\partial \eta }-\phi$, $\frac{\delta \ell }{\delta s}= - \frac{\partial e}{\partial \eta }$, $\frac{\delta \ell }{\delta \rho }\diamond \rho = \rho \nabla \frac{\delta \ell }{\delta \rho }$, and $\frac{\delta \ell }{\delta s}\diamond s= s \nabla \frac{\delta \ell }{\delta s}$, one directly obtains

$$\begin{aligned} \rho (\partial _tu + u \cdot \nabla u)= - \nabla p - \rho \nabla \phi \end{aligned}$$

from (79), with $p= \rho ^2\frac{\partial e}{\partial \rho }$. For the semidiscrete case, one uses a nonholonomic version of the Euler–Poincaré equations (79), reviewed in the next paragraph (Fig. 3).

Fig. 3

Relative errors in the energy $E(t) = \int _\Omega \left( \frac{1}{2} \rho |u|^2 + \rho e(\rho ) + \rho \phi \right) \, \mathrm{d}x$ during the Rayleigh–Taylor instability simulation

Remarks on the Nonholonomic Euler–Poincaré Variational Formulation

Hamilton’s principle can be extended to the case in which the system under consideration is subject to a constraint, given by a distribution on the configuration manifold, i.e., a vector subbundle of the tangent bundle. This is known as the Lagrange-d’Alembert principle and, for a system on a Lie group G and constraint ${\Delta _G} \subset TG$, it is given by the same critical condition (74) but only with respect to variations satisfying the constraint, i.e., ${\delta g}\in {\Delta _G}$.

In the G-invariant setting recalled in Sect. A.1 it is assumed that the constraint ${\Delta _G}$ is also G-invariant and thus induces a subspace $\Delta \subset {\mathfrak {g}}$ of the Lie algebra. In the more general setting of Sect. A.2, one can allow ${\Delta _G}$ to be only $G_{a_0}$-invariant, although for the situation of interest in this paper, ${\Delta _G}$ is also G-invariant.

The Lagrange-d’Alembert principle yields now the Euler–Poincaré–d’Alembert principle (77) in which we have the additional constraint $u(t)\in \Delta$ on the solution and $v(t)\in \Delta$ on the variations, so that (78) becomes

$$\begin{aligned} \Big \langle \Big \langle \partial _t \frac{\delta \ell }{\delta u}, v\Big \rangle \Big \rangle + \Big \langle \Big \langle \frac{\delta \ell }{\delta u}, [u,v]\Big \rangle \Big \rangle + \Big \langle \frac{\delta \ell }{\delta a},a\cdot v \Big \rangle _V=0,\quad \text {for all}~v\in \Delta , \text {where}~u\in \Delta . \end{aligned}$$

(80)

In presence of the nonholonomic constraint, (79) becomes

$$\begin{aligned} \partial _t \frac{\delta \ell }{\delta u}+ {\text {ad}}^*_u\frac{\delta \ell }{\delta u} -\frac{\delta \ell }{\delta a}\diamond a\in \Delta ^\circ ,\qquad u \in \Delta , \end{aligned}$$

(81)

where $\Delta ^\circ = \{m \in {\mathfrak {g}}^* \mid \langle \!\langle m,u\rangle \!\rangle =0, \;\forall \; u\in \Delta \}$.

There are two important remarks concerning (80) and (81) that play an important role for the variational discretization carried out in this paper. First, we note that although the solution belongs to the constraint, i.e., $u\in \Delta$, the equations depend on the expression of the Lagrangian $\ell$ on a larger space, namely, on $\Delta + [\Delta , \Delta ]$. It is not enough to have its expression only on $\Delta$. This is a main characteristic of nonholonomic mechanics. Second, a sufficient condition to get a solvable differential equation is that the map $u \in \Delta \mapsto \frac{\delta \ell }{\delta u } \in {\mathfrak {g}}^*/\Delta ^\circ$ is a diffeomorphism for all a.

Polynomials

Below we prove two facts about polynomials that are used in the proof of Proposition 3.3. We denote by $H_r(K)$ the space of homogeneous polynomials of degree r on a simplex K. To distinguish powers from indices, we denote coordinates by $x_1,x_2,\dots ,x_n$ rather than $x^1,x^2,\dots ,x^n$ in this section.

Lemma C.1

Let K be a simplex of dimension $n \ge 1$. For every integer $r \ge 0$,

$$\begin{aligned} \left\{ \sum _{i=1}^N p_i q_i \mid N \in {\mathbb {N}}, \, p_i, q_i \in P_r(K), \, i=1,2,\dots ,N \right\} = P_{2r}(K). \end{aligned}$$

Proof

This follows from the fact that every monomial in $P_{2r}(K)$ can be written as a product of two monomials in $P_r(K)$. $\square$

Lemma C.2

Let K be a simplex of dimension $n \in \{2,3\}$. For every integer $r \ge 0$,

$$\begin{aligned} \left\{ \sum _{i=1}^N p_i \nabla q_i \mid N \in {\mathbb {N}}, \, p_i, q_i \in P_r(K), \, i=1,2,\dots ,N \right\} = P_{2r-1}(K)^n. \end{aligned}$$

Proof

We proceed by induction.

Denote $Q_r(K) = \left\{ \sum _{i=1}^N p_i \nabla q_i \mid N \in {\mathbb {N}}, \, p_i, q_i \in P_r(K), \, i=1,2,\dots ,N \right\}$. By inductive hypothesis, $Q_r(K)$ contains $Q_{r-1}(K)=P_{2r-3}(K)^n$. It also contains $H_{2r-2}(K)^n$. Indeed, if $f e_k \in H_{2r-2}(K)^n$ with $k \in \{1,2,\dots ,n\}$ and f a monomial, then $f=x_j g$ for some $g \in H_{2r-3}(K)$ and some $j \in \{1,2,\dots ,n\}$, so $fe_k = \sum _i (x_j p_i) \nabla q_i \in Q_r(K)$ for some $p_i,q_i \in P_{r-1}(K)$ by inductive hypothesis. Thus, $Q_r(K)$ contains $P_{2r-2}(K)^n$.

Next we show that $Q_r(K)$ contains every $u \in H_{2r-1}(K)^n$. Without loss of generality, we may assume $u = fe_1$ with $f \in H_{2r-1}(K)$ a monomial. When $n=2$, the only such vector fields are $u = x_1^a x_2^{2r-1-a} e_1$, $a=0,1,\dots ,2r-1$, which can be expressed as

$$\begin{aligned} x_1^a x_2^{2r-1-a} e_1 = {\left\{ \begin{array}{ll} \frac{1}{r} x_1^{a-r+1} x_2^{2r-1-a} \nabla (x_1^r), &{} \text{ if } a \ge r-1, \\ \frac{1}{a+1} x_2^r \nabla ( x_1^{a+1} x_2^{r-1-a} ) - \frac{r-1-a}{(a+1)r} x_1^{a+1} x_2^{r-1-a} \nabla (x_2^r), &{} \text{ if } a < r-1. \end{array}\right. } \end{aligned}$$

The case $n=3$ is handled similarly by considering the vector fields $f e_1$ with

$$\begin{aligned} f \in \{x_1^a x_2^b x_3^{2r-1-a-b} \mid a,b \ge 0, \, a+b \le 2r-1 \}. \end{aligned}$$

$\square$