Observation

When using a forward-Euler method for the time integration of the momentum equation for an inviscid-flow, it appears that the kinetic energy of the flow grows unbounded in time, regardless of the timestep size.

Problem Statement

Estimate the change in total kinetic energy when using forward-Euler to integrate the Euler momentum equations in a periodic box.

Approach

To solve this problem, we need to do the following:

1. Formulate a FE semi-discrete formula for the Euler equations.
2. Construct the local kinetic energy (i.e. pointwise) by taking the dot product of the semi-discrete velocity field.
3. Using this information, construct an equation for the change in total kinetic energy with time, e.g. ${\displaystyle \frac{\Delta K}{\Delta t}} $.
4. Analyze the right-hand-side (RHS) of this equation in the context of a periodic box. If the RHS > 0, then the kinetic energy will increase.

Show Me the Math

It is best to use tensor notation for this problem. We start with the Euler equations

$$\begin{equation} \frac{\partial u_{i}}{\partial x_{i}}=0 \end{equation}$$

$$\begin{equation} \frac{\partial u_{i}}{\partial t}=-u_{j}\frac{\partial u_{i}}{\partial x_{j}}-\frac{\partial p}{\partial x_{i}}\equiv F_{i} \end{equation}$$

where $i\in\{1,2,3\}$ denotes the $i^{\text{th}}$ spatial direction.

Using forward-euler time integration, we construct the semi-discrete formula

$$\begin{equation} \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}+\mathcal{O}\left(\Delta t\right)=F_{i}^{n};\quad F_{i}^{n}\equiv u_{j}^{n}\frac{\partial u_{i}^{n}}{\partial x_{j}}+\frac{\partial p^{n}}{\partial x_{i}} \end{equation}$$

where $t^{n+1}=t^{n}+\Delta t$ and $\Delta t$ is the fixed timestep size. For simplity of the exposition, we drop the “order” notation and set

$$\begin{equation} u_{i}^{n+1}=u_{i}^{n}-\Delta t\, F_{i}^{n}\label{eq:ke-euler-discrete} \end{equation}$$

Now, in a staggered arrangement, the kinetic energy lives at cell centers while the velocity components are stored on staggered volumes. For simpliciy, we assume a uniform, staggered grid. Consequently, the staggered-volumes coincide with the faces of the cell centered ones. Then, one calculates the local kinetic energy, $k$, as

$$\begin{equation} k=\tfrac{1}{2}u_{i}u_{i}=\tfrac{1}{2}(u_{1}^{2}+u_{2}^{2}+u_{3}^{2}) \end{equation}$$

where Einstein summation is implied on repeated inidces. Note that the local kinetic energy is defined as the pointwise kinetic energy – the kinetic energy in every cell. Now, we can formulate the following sub-problem: given a discrete velocity field $u_{i}^{n}$ such that $\frac{\partial u_{i}^{n}}{\partial x_{i}}=0$ (discretely), estimate the kinetic energy, $ k^{n+1}$, at the next timestep. In other words, use $\eqref{eq:ke-euler-discrete}$ to calculate $k^{n+1}$.

We start by multiplying $\eqref{eq:ke-euler-discrete}$ through with $u_{i}^{n+1}$, we have

$$\begin{equation} u_{i}^{n+1}u_{i}^{n+1}=u_{i}^{n+1}\left(u_{i}^{n}-\Delta t\, F_{i}^{n}\right) \end{equation}$$

Or, using $k=\frac{1}{2}u_{i}u_{i}$, we have

$$\begin{equation} 2k^{n+1}=\left(u_{i}^{n}-\Delta t\, F_{i}^{n}\right)\left(u_{i}^{n}-\Delta t\, F_{i}^{n}\right) \end{equation}$$

Expanding the right-hand-side, we get

$$\begin{equation} 2k^{n+1}=u_{i}^{n}u_{i}^{n}-2\Delta t\, u_{i}^{n}F_{i}^{n}+\Delta t^{2}F_{i}^{n}F_{i}^{n} \end{equation}$$
or
$$\begin{equation} 2k^{n+1}=2k^{n}-2\Delta t\, u_{i}^{n}F_{i}^{n}+\Delta t^{2}F_{i}^{n}F_{i}^{n}\label{eq:kn1-0} \end{equation}$$
We now rearrange $\eqref{eq:kn1-0}$ as follows
$$\begin{equation} \frac{k^{n+1}-k^{n}}{\Delta t}=-u_{i}^{n}F_{i}^{n}+\frac{1}{2}\Delta tF_{i}^{n}F_{i}^{n} \end{equation}$$
or
$$\begin{equation} \frac{\Delta k}{\Delta t}=-u_{i}^{n}F_{i}^{n}+\frac{1}{2}\Delta tF_{i}^{n}F_{i}^{n}\label{eq:dkdt-0} \end{equation}$$

This equation tells us how much the implied kinetic energy changes with time. It is an implied kinetic energy equation given that it was constructed from the velocity field not from the transported kinetic energy.

We now focus our attention on the first term in the RHS of $\eqref{eq:dkdt-0}$, i.e. $u_{i}^{n}F_{i}^{n}$. Substituting for $F_{i}^{n}$, we have

$$\begin{equation} u_{i}^{n}F_{i}^{n}=u_{i}^{n}\left(u_{j}^{n}\frac{\partial u_{i}^{n}}{\partial x_{j}}+\frac{\partial p^{n}}{\partial x_{i}}\right)=u_{i}^{n}u_{j}^{n}\frac{\partial u_{i}^{n}}{\partial x_{j}}+u_{i}^{n}\frac{\partial p^{n}}{\partial x_{i}}\label{eq:ui-fn} \end{equation}$$

The purpose now is to try to convert this term into a divergence form. We will see why later. Starting with the first term $u_{i}^{n}u_{j}^{n}\frac{\partial u_{i}^{n}}{\partial x_{j}}$, we know that

$$\begin{equation} \frac{\partial u_{i}^{n}u_{i}^{n}u_{j}^{n}}{\partial x_{j}}=u_{i}^{n}u_{i}^{n}\frac{\partial u_{j}^{n}}{\partial x_{j}}+u_{j}^{n}\frac{\partial u_{i}^{n}u_{i}^{n}}{\partial x_{j}}=u_{i}^{n}u_{i}^{n}\frac{\partial u_{j}^{n}}{\partial x_{j}}+2u_{i}^{n}u_{j}^{n}\frac{\partial u_{i}^{n}}{\partial x_{j}} \end{equation}$$

But, assuming that continuity is satified discretely at every cell, we can set $\frac{\partial u_{j}^{n}}{\partial x_{j}}=0$. Then, we recover

$$\begin{equation} u_{i}^{n}u_{j}^{n}\frac{\partial u_{i}^{n}}{\partial x_{j}}=\frac{1}{2}\frac{\partial u_{i}^{n}u_{i}^{n}u_{j}^{n}}{\partial x_{j}}=\frac{\partial k^{n}u_{j}^{n}}{\partial x_{j}} \end{equation}$$
Lastly, the second term in $\eqref{eq:ui-fn}$ is easily replaced by
$$\begin{equation} u_{i}^{n}\frac{\partial p^{n}}{\partial x_{i}}=\frac{\partial u_{i}^{n}p^{n}}{\partial x_{i}}-p^{n}\frac{\partial u_{i}^{n}}{\partial x_{i}} \end{equation}$$

Again, by assuming that continuity is satisfied discretely, we can drop the last term in the previous equation and use

$$\begin{equation} u_{i}^{n}\frac{\partial p^{n}}{\partial x_{i}}=\frac{\partial u_{i}^{n}p^{n}}{\partial x_{i}} \end{equation}$$

Upon substitution of the modified terms back into the implied kinetic energy equation $\eqref{eq:dkdt-0}$, we have

$$\begin{equation} \frac{\Delta k}{\Delta t}=-\frac{\partial k^{n}u_{j}^{n}}{\partial x_{j}}-\frac{\partial u_{i}^{n}p^{n}}{\partial x_{i}}+\frac{1}{2}\Delta tF_{i}^{n}F_{i}^{n} \end{equation}$$

In principle, for an inviscid flow in a period box, there should be no production or dissipation of total kinetic energy. The total kinetic energy is defined as the volumetric integral of $k$ over the domain, i.e.

$$\begin{equation} \mathrm{K}=\int_{V}k\,\mathrm{d}V \end{equation}$$

Since there is no production or dissipation in our problem, we expect that $\text{K}=const$ or

$$\begin{equation} \frac{\partial}{\partial t}\int_{V}k\,\mathrm{d}V=0 \end{equation}$$

Looking at our semi-discrete equation, upon integration over the periodic box, we have

$$\begin{equation} \frac{\Delta}{\Delta t}\int_{V}k\mathrm{d}V=-\int_{V}\frac{\partial k^{n}u_{j}^{n}}{\partial x_{j}}\mathrm{d}V-\int_{V}\frac{\partial u_{i}^{n}p^{n}}{\partial x_{i}}\mathrm{d}V+\frac{1}{2}\Delta t\int_{V}F_{i}^{n}F_{i}^{n}\mathrm{d}V \end{equation}$$

Being periodic, we have, for any vector $\mathbf{v}$,

$$\begin{equation} \int_{V}\nabla\cdot\mathbf{v}\mathrm{d}V=\int_{V}\frac{\partial v_{i}}{\partial x_{i}}\mathrm{d}V=\int_{\mathcal{S}}\mathbf{v}\cdot\mathbf{n}\mathrm{d}\mathcal{S}=0 \end{equation}$$

Then, all terms that are written in divergence form vanish identically. Hence, one is left with,

$$\begin{equation} \frac{\Delta\mathrm{K}}{\Delta t}=\frac{1}{2}\Delta t\int_{V}F_{i}^{n}F_{i}^{n}\mathrm{d}V\geq0 \end{equation}$$

and therefore, the total kinetic energy grows unboundedly with time when using a forward-Euler scheme on the Euler equations, in a periodic box.