Energy Conservation in Conservation Laws with Source Terms

Question

I'm wondering if anyone can help me understand energy conservation when using conservation law methods (i.e. Riemann solver, High-Resolution Wave-Propagation Methods) with the addition of source terms:

$ \frac{\partial \vec{Q}}{\partial t} + \nabla \cdot \vec{F} = \vec{\psi} $,

where, $\vec{Q}$, and $\vec{F}$ are the conserved quantities and their corresponding fluxes. $\vec{\psi}$ are the non-hyperbolic source terms.

A 1D example of my Euler equation setup with an external gravitational field is:

$ \frac{\partial}{\partial t}\left[ \rho \right] + \frac{\partial}{\partial z}\left[ \rho w \right] = 0$

$ \frac{\partial}{\partial t}\left[ \rho w \right] + \frac{\partial}{\partial z}\left[ \rho w^2 + p \right] = -\rho g$

$\frac{\partial}{\partial t}\left[ e \right] + \frac{\partial}{\partial z}\left[ (e+p)w\right] = -\rho g w$

where, $t$ is time, $z$ is the spatial coordinate, $\rho$ is density, $w$ is fluid velocity, $e$ is the energy (internal + kinetic), $g$ is the acceleration due to gravity in the $z$ direction, and $p$ is the pressure. The system is closed by defining an equation of state:

$p = (\gamma -1)e_{i} $,

where $e_{i} = e - \frac{1}{2}\rho w^2$ is the internal energy, and $\gamma$ is the ratio of specific heats.

The system of equations is solved using a wave-propagation, Riemann solver with an F-Wave approach using the PyClaw software (http://www.clawpack.org/doc/pyclaw/). This means the non-hyperbolic gravity source terms are included in the flux differencing (rather than solved using some fractional step source splitting technique). The advantage of including the gravity source terms in the flux differencing is that in an equilibrium case, the source terms will exactly cancel the pressure gradient and maintain numerical equilibrium. This has proven to work very well.

My question is regarding the conservation of energy, particularly the total energy that also includes gravitational energy, $ e_{total} = e_i + e_{KE} + e_{gravitational}$. Is there any expectation that the total energy $e_{total}$ should be conserved as precisely as the "conserved" quantity energy $e$ when the hyperbolic terms are modified with source terms? In a test problem where an equilibrium case is perturbed, I see excellent conservation of $e$, but not great conservation of the total energy $e_{total}$. I'm not sure what to make of it.

I can provide a more detailed explanation of the equilibrium setup, and the perturbed case if that is helpful.

Thanks, Wes

Answer 1

I have not studied this particular system before, and I'm sure that someone who has could say much more than I will.

I don't think there is any reason to expect that discretizing the equations in the form you have them will lead to highly accurate conservation of total energy. In this form, conservation of total energy depends on cancellation of certain flux terms and source terms, but since they are discretized in different ways they will only cancel up to truncation errors.

One way to get conservation of total energy, up to roundoff, is to rewrite the equations in conservative form. I worked this out quickly but I think it's correct. Since

$$\frac{\partial e}{\partial t} = \frac{\partial e_{total}}{\partial t} + \frac{\partial}{\partial z}(\rho g w z) - \rho g w,$$

we have

$$\frac{\partial e_{total}}{\partial t} + \frac{\partial}{\partial z}((e+p+\rho g z)w) = 0.$$

If you discretize this conservation law in place of your current energy equation, you will get numerical conservation of total energy, up to roundoff. In Clawpack this means that you will need to write a new Riemann solver; I suspect some Riemann solver for this form has been worked out in the literature. Also, you should use interface values (not cell-center values) of $z$ for the flux term that depends on it (so that the fluxes at interfaces cancel).

Answer 2

I am just looking at Käppeli and Mishra, "Well-balanced schemes for the Euler equations with gravitation". They make one key observation -- use the enthalpy to reconstruct values at cell boundaries -- and derive from this a formulation where the truncation errors of flux and source terms cancel.

The paper is surprisingly easy to read; kudos to the authors.