Notes on Enthalpy

Evolution Equations

The compressible and low Mach number formulations of the governing equations both share the unapproximated continuity and momentum equations shown here (where $p=p_0+\pi$).

(392)$$\begin{array}{r}\begin{array}{rlr}\frac{\partial (\rho X_k)}{\partial t}& =& -\mathrm{\nabla }\cdot (\rho X_k\mathbf{U})+\rho {\dot{\omega }}_k,\\ \frac{\partial \mathbf{U}}{\partial t}& =& -\mathbf{U}\cdot \mathrm{\nabla }\mathbf{U}-\frac{1}{\rho }\mathrm{\nabla }\pi -\frac{\rho -{\rho }_0}{\rho }g{\mathbf{e}}_r,\end{array}\end{array}$$

In the compressible formulation we complete the system with an energy equation as well as an equation of state; in the low Mach number formulation we can derive a constraint on the velocity by setting $p_{EOS}=p(\rho ,h,X_k)=p_0$ and differentiating the equation of state along particle paths. In this case adding the energy equation would over-constrain the system, but its solution can help us in providing a diagnostic capability for the solution. We can write the energy equation in terms of internal energy, $e$, or enthalpy, $h$; these two equations are analytically equivalent.

(393)$$\begin{array}{r}\begin{array}{rlr}\frac{\partial (\rho h)}{\partial t}+\mathrm{\nabla }\cdot (\rho h\mathbf{U})& =& \frac{Dp}{Dt}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}\\ \frac{\partial (\rho e)}{\partial t}+\mathrm{\nabla }\cdot (\rho e\mathbf{U})& =& -p\mathrm{\nabla }\cdot \mathbf{U}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}\end{array}\end{array}$$

In low Mach number combustion, $Dp/Dt=0,$ which makes the enthalpy equation preferable to the internal energy equation because we don’t need to evaluate $p\mathrm{\nabla }\cdot \mathbf{U}.$

Derivation of Velocity Constraint

Differentiating the equation of state, written in the form, $p=p(\rho ,T,X_k),$ along particle paths, we can write

(394)$$\frac{Dp}{Dt}=p_{\rho }\frac{D\rho }{Dt}+p_T\frac{DT}{Dt}+\sum _k{p}_{X_k}\frac{D{X}_k}{Dt}$$

Then, by rearranging the terms, we get

(395)$$\frac{D\rho }{Dt}=-\rho \mathrm{\nabla }\cdot \mathbf{U}=\frac{1}{p_{\rho }}(\frac{Dp}{Dt}-p_T\frac{DT}{Dt}-\sum _k{p}_{X_k}{\dot{\omega }}_k),$$

with $p_{\rho }={\partial p/\partial \rho |}_{X_k,T}$, $p_{X_k}={\partial p/\partial X_k|}_{T,\rho ,(X_j,j\ne k)}$, and $p_T={\partial p/\partial T|}_{\rho ,X_k}$.

Using the Enthalpy Equation

Now writing $h=h(p,T,X_k)$ and expanding $Dh/Dt$ and using the enthalpy evolution equation:

(396)$$\rho \frac{Dh}{Dt}=\rho (h_p\frac{Dp}{Dt}+c_p\frac{DT}{Dt}+\sum _k{h}_{X_k}\frac{D{X}_k}{Dt})=\frac{Dp}{Dt}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}$$

we can express $DT/Dt$ in terms of

(397)$$\frac{DT}{Dt}=\frac{1}{\rho c_p}((1-\rho h_p)\frac{Dp}{Dt}-\sum _k\rho {\xi }_k{\dot{\omega }}_k+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}),$$

where $c_p={\partial h/\partial T|}_{p,X_k}$ is the specific heat at constant pressure, ${\xi }_k={\partial h/\partial X_k|}_{p,T}$, and $h_p={\partial h/\partial p|}_{T,X_k}.$

We could then write,

(398)$$\begin{array}{r}\begin{array}{rlr}\mathrm{\nabla }\cdot \mathbf{U}& =& \frac{1}{\rho p_{\rho }}(-\frac{Dp}{Dt}+\frac{p_T}{\rho c_p}((1-\rho h_p)\frac{Dp}{Dt}-\rho \sum _k{\xi }_k{\dot{\omega }}_k+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}})+\sum _k{p}_{X_k}{\dot{\omega }}_k)\\ & =& \frac{1}{\rho p_{\rho }}(\frac{p_T}{\rho c_p}(1-\rho h_p)-1)\frac{Dp}{Dt}+\frac{1}{\rho p_{\rho }}(\frac{p_T}{\rho c_p}(\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}-\rho \sum _k{\xi }_k{\dot{\omega }}_k)+\sum _k{p}_{X_k}{\dot{\omega }}_k).\end{array}\end{array}$$

When we derived this expression we explicitly retained the dependence of $h$ on $p$, as shown by the presence of the $h_p$ term.

Then, replacing $p$ by $p_0(r)$, $Dp/Dt$ becomes $\mathbf{U}\cdot \mathrm{\nabla }{p}_0$, and the divergence constraint can be written

(399)$$\mathrm{\nabla }\cdot \mathbf{U}+\alpha \mathbf{U}\cdot \mathrm{\nabla }{p}_0=\frac{1}{\rho p_{\rho }}(\frac{p_T}{\rho c_p}(-\sum _k\rho {\xi }_k{\dot{\omega }}_k+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}})+\sum _k{p}_{X_k}{\dot{\omega }}_k)\equiv \tilde{S},$$

where we define

(400)$$\alpha (\rho ,T)\equiv -\left(\frac{(1-\rho h_p)p_T-\rho c_p}{{\rho }^2{c}_p{p}_{\rho }}\right).$$

Using the Energy Equation

Now writing $e=e(p,T,X_k)$ and expanding $De/Dt$ and using the energy evolution equation:

(401)$$\rho \frac{De}{Dt}=\rho (e_p\frac{Dp}{Dt}+e_T\frac{DT}{Dt}+\sum _k{e}_{X_k}\frac{D{X}_k}{Dt})=-p\mathrm{\nabla }\cdot \mathbf{U}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}$$

we can express $DT/Dt$ in terms of

(402)$$\begin{array}{rlr}\frac{DT}{Dt}& =& \frac{1}{\rho e_T}(-p\mathrm{\nabla }\cdot \mathbf{U}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}-\rho e_p\frac{Dp}{Dt}-\rho \sum _k{e}_{X_k}{\dot{\omega }}_k)\end{array}$$

where $e_T={\partial e/\partial T|}_{p,X_k}$ and $e_p={\partial e/\partial p|}_{T,X_k}.$ Then

(403)$$-\rho \mathrm{\nabla }\cdot \mathbf{U}=\frac{1}{p_{\rho }}(\frac{Dp}{Dt}-\frac{p_T}{\rho e_T}(-p\mathrm{\nabla }\cdot \mathbf{U}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}-\rho e_p\frac{Dp}{Dt}-\rho \sum _k{e}_{X_k}{\dot{\omega }}_k)-\sum _k{p}_{X_k}{\dot{\omega }}_k),$$

which leads to

(404)$$(-\rho -\frac{p{p}_T}{\rho e_T{p}_{\rho }})\mathrm{\nabla }\cdot \mathbf{U}=\frac{1}{p_{\rho }}((1+\frac{e_p{p}_T}{e_T})\frac{Dp}{Dt}-\frac{p_T}{\rho e_T}(\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}-\rho \sum _k{e}_{X_k}{\dot{\omega }}_k)-\sum _k{p}_{X_k}{\dot{\omega }}_k),$$

Note that we can replace $p$ by $p_0$ in the coefficient on the l.h.s. as well as on the r.h.s.

Comparison of Constraints

If we set ${\dot{\omega }}_k=H_{\mathrm{n}\mathrm{u}\mathrm{c}}=0$ for simplicity, then the constraint as derived using $h$ can be written

(405)$$\mathrm{\nabla }\cdot \mathbf{U}+\left(\frac{(1-\rho h_p)p_T-\rho c_p}{{\rho }^2{c}_p{p}_{\rho }}\right)\frac{D{p}_0}{Dt}=0$$

and the constraint derived using $e$ can be written

(406)$$\mathrm{\nabla }\cdot \mathbf{U}+\left(\frac{\rho e_T+\rho e_p{p}_T}{{\rho }^2{e}_T{p}_{\rho }+p{p}_T}\right)\frac{D{p}_0}{Dt}=0$$

We note that if we evaluate both constraints for $p=\rho RT,$ with $h=c_{p}T,$ $e=c_{v}T,$ $c_p=c_v+R$ and $\gamma =c_p/c_v,$ then both constraints reduce to

(407)$$\mathrm{\nabla }\cdot \mathbf{U}+\frac{1}{\gamma p}\frac{D{p}_0}{Dt}=0$$

Enthalpy vs Energy Equation

The full enthalpy equation, with no approximations, appears as:

(408)$$\frac{\partial (\rho h)}{\partial t}=-\mathrm{\nabla }\cdot (\rho h\mathbf{U})+\frac{Dp}{Dt}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}$$

Here, $h=e+p/\rho$ is the specific enthalpy, with $e$ the specific internal energy. In the low Mach number formulation, we replace $p$ with $p_0$ in the $Dp/Dt$ term, however, the definition of enthalpy implicitly contains a pressure. When calling the equation of state, we take $h$ and $\rho$ as inputs. The equation of state is expressed in terms of $T$ and $\rho$, so it iterates until it finds the $h$ that we desire. This $h$ will be of the form $h=e+p_{\mathrm{EOS}}/\rho$, where $p_{\mathrm{EOS}}$ is the pressure returned from the EOS. Note that $p_{\mathrm{EOS}}$ may not be equal to $p_0$—this may be what causes us to drive off of the constraint.

The mismatch between the pressure implicit in the definition of $h$ and $p_0$ can be seen by substituting $h=e+p/\rho$ into the enthalpy equation, where we replace $p$ with $p_0$ in the $Dp/Dt$ term:

(409)$$\begin{array}{r}\begin{array}{rlr}\frac{\partial (\rho h)}{\partial t}& =& -\mathrm{\nabla }\cdot (\rho h\mathbf{U})+\frac{D{p}_0}{Dt}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}\\ \frac{\partial (\rho e)}{\partial t}+\frac{\partial p}{\partial t}& =& -\mathrm{\nabla }\cdot (\rho e\mathbf{U})-\mathrm{\nabla }\cdot (p\mathbf{U})+\frac{D{p}_0}{Dt}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}\\ \frac{\partial (\rho e)}{\partial t}& =& -\mathrm{\nabla }\cdot (\rho e\mathbf{U})-p\mathrm{\nabla }\cdot \mathbf{U}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}+\{\frac{D{p}_0}{Dt}-\frac{Dp}{Dt}\}\end{array}\end{array}$$

However, if we solve the evolution equation for $e$ we would substitute $p_0$ for $p$ in this equation as well. Thus, we can pose the situation as the following. If we solve the evolution equation for $h$ then we effectively are solving

(410)$$\begin{array}{}& \frac{\partial (\rho e)}{\partial t}+\mathrm{\nabla }\cdot (\rho e\mathbf{U})=-p\mathrm{\nabla }\cdot \mathbf{U}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}+\{\frac{D{p}_0}{Dt}-\frac{Dp}{Dt}\}\end{array}$$

but if we solve the evolution equation for $e$ we are effectively solving

(411)$$\begin{array}{}& \frac{\partial (\rho e)}{\partial t}+\mathrm{\nabla }\cdot (\rho e\mathbf{U})=p_0\mathrm{\nabla }\cdot \mathbf{U}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}\end{array}$$

The second equation subtracted from the first gives:

(412)$$\frac{D(p_0-p)}{Dt}-(p_0-p)\mathrm{\nabla }\cdot \mathbf{U}=0,$$

but this equation is only true, in general, if $p=p_0$.

Suppose we solve the current enthalpy equation, but when we call the EOS, we subtract $p_0$ from $\rho h$ and then call the EOS with $e$ instead of $h$. This is equivalent to:

(413)$$\begin{array}{r}\begin{array}{rlr}\frac{\partial (\rho h)}{\partial t}& =& -\mathrm{\nabla }\cdot (\rho h\mathbf{U})+\frac{D{p}_0}{Dt}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}\\ \frac{\partial (\rho e)}{\partial t}+\frac{\partial p_0}{\partial t}& =& -\mathrm{\nabla }\cdot (\rho e\mathbf{U})-\mathrm{\nabla }\cdot (p_0\mathbf{U})+\frac{D{p}_0}{Dt}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}\\ \frac{\partial (\rho e)}{\partial t}& =& -\mathrm{\nabla }\cdot (\rho e\mathbf{U})-p_0\mathrm{\nabla }\cdot \mathbf{U}+\rho H_{\mathrm{n}\mathrm{u}\mathrm{c}}\end{array}\end{array}$$

which is identical to solving the energy equation with $p\to p_0$. This option is enabled in MAESTROeX via use_eos_e_instead_of_h = T.

Constant $\gamma$ Gas

Going back to the constant $\gamma$, ideal gas EOS, we can rewrite the enthalpy equation as a pressure evolution equation

(414)$$\begin{array}{r}\begin{array}{rlr}\frac{\partial \rho h}{\partial t}+\mathrm{\nabla }\cdot \left(\rho h\mathbf{U}\right)& =& \frac{D{p}_0}{Dt}+\rho H\\ \frac{\gamma }{\gamma -1}\frac{\partial p}{\partial t}+\frac{\gamma }{\gamma -1}\mathrm{\nabla }\cdot \left(p\mathbf{U}\right)& =& \frac{D{p}_0}{Dt}+\rho H\\ \frac{Dp}{Dt}& =& -p\mathrm{\nabla }\cdot \mathbf{U}+\frac{\gamma }{\gamma -1}(\frac{D{p}_0}{Dt}+\rho H)\end{array}\end{array}$$

Similarly, we can derive a pressure evolution equation from the energy equation

(415)$$\frac{Dp}{Dt}=-p\mathrm{\nabla }\cdot \mathbf{U}-(\gamma -1)p_0\mathrm{\nabla }\cdot \mathbf{U}+\rho H$$

Now, if we further make the assumption that $p_0$ is constant, $D{p}_0/Dt=0$, then the divergence constraint for such a gas reduces to

(416)$$\mathrm{\nabla }\cdot \mathbf{U}=\frac{\gamma -1}{\gamma p_0}\rho H.$$

Plugging this back into either of Eq.414 or Eq.415 gives

(417)$$\frac{Dp}{Dt}=\frac{\gamma -1}{\gamma }(1-\frac{p}{p_0})\rho H.$$

If $p_0$ is assumed constant and using div constraint for constant gamma, the difference between the enthalpy equation and the energy equation, Eq.412, can be rewritten as

(418)$$-\frac{Dp}{Dt}-\frac{\gamma }{\gamma -1}(1-\frac{p}{p_0})\rho H=0,$$

where the equality holds from Eq.417. In other words, for the constant $\gamma$ gas we have $p=p_0$ as expected.

Outstanding Questions

1. Why do we want to start with enthalpy instead of internal energy?

   We believe that the original desire stems from our experience with smallscale combustion. There, stratification is not important and $D{p}_0/Dt=0$, so the enthalpy equation becomes a conservation equation for $(\rho h)$.

2. Should we call the EOS with $h$ as is, or call the EOS with $e=h-p_0/\rho$?

3. When we stay on the constraint, i.e. $p_{EOS}=p_0$, then the equations for $e$ and for $h$ are identical. However, once we are off the constraint, do the terms in the current evolution equation for $h$ serve to drive us back to the constraint? Recall our current “volume discrepancy factor” acts as a source term which modifies the divergence constraint, which effectively modifies both $\rho$ and $T$ (or $e$ or $h$). The term in the enthalpy equation only modifies $\rho .$ Is this relevant and/or useful?? Recall that the current “volume discrepancy factor” takes the form of adding to the r.h.s. of the constraint:

   (419)$$\mathrm{\nabla }\cdot ({\beta }_0\mathbf{U})={\beta }_0(S-\frac{1}{\overline{{\mathrm{\Gamma }}_1}{p}_0}\frac{\partial p_0}{\partial t}-\frac{f}{\overline{{\mathrm{\Gamma }}_1}{p}_0}\frac{p_0-p_{\text{EOS}}}{\mathrm{\Delta }t})$$

4. Suppose we corrected the $h$ (or $e$ equation) by using the full $p$ instead of $p_0$? Would this be more or less consistent (one could imagine doing this as a correction after solving for $\pi$ earlier in the timestep).

5. In computing the thermodynamic coefficients in $S$ for the projection, don’t we need these to be in terms of $p_0$ instead of $p(h,\rho )$?