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{split}\begin{aligned} \frac{\partial(\rho X_k)}{\partial t} &=& -\nabla\cdot(\rho X_k\Ub) + \rho\omegadot_k,\label{enth:eq:species}\\ \frac{\partial\Ub}{\partial t} &=& -\Ub\cdot\nabla\Ub - \frac{1}{\rho}\nabla\pi - \frac{\rho-\rho_0}{\rho} g\eb_r,\label{eq:momentum}\end{aligned}\end{split}

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{split}\begin{aligned} \frac{\partial(\rho h)}{\partial t} + \nabla\cdot(\rho h\Ub) &=& \frac{Dp}{Dt} + \rho\Hnuc \nonumber \\ % \frac{\partial(\rho e)}{\partial t} + \nabla\cdot(\rho e \Ub) &=& - p\nabla\cdot \Ub + \rho\Hnuc \nonumber \end{aligned}\end{split}

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 \nabla \cdot \Ub.$$

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{D p}{Dt} = p_\rho \frac{D \rho}{Dt} + p_T \frac{D T}{Dt} + \sum_k p_{X_k} \frac{D X_k}{Dt}$

Then, by rearranging the terms, we get

(395)$\frac{D \rho}{Dt} = -\rho \nabla \cdot \Ub = \frac{1}{p_\rho} \left( \frac{D p}{Dt} - p_T \frac{D T}{Dt} - \sum_k p_{X_k} {\omegadot}_k \right) ,$

with $$p_\rho = \left.\partial p/\partial \rho\right|_{X_k,T}$$, $$p_{X_k} = \left.\partial p/\partial X_k \right|_{T,\rho,(X_j,j\ne k)}$$, and $$p_T = \left.\partial p/\partial T\right|_{\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{D h}{Dt} = \rho \left( h_p \frac{D p}{Dt} + c_p \frac{D T}{Dt} + \sum_k h_{X_k} \frac{D X_k}{Dt} \right) = \frac{D p}{D t} + \rho \Hnuc$

we can express $$DT/Dt$$ in terms of

(397)$\frac{DT}{Dt} = \frac{1}{\rho c_p} \left( (1 - \rho h_p) \frac{D p}{D t} - \sum_k \rho \xi_k \omegadot_k + \rho \Hnuc \right) , \label{eq:dTdt}$

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

We could then write,

(398)\begin{split}\begin{aligned} \nabla \cdot \Ub &=& \frac{1}{\rho p_\rho} \left( - \frac{D p}{D t} + \frac{p_T}{\rho c_p} \left( (1 - \rho h_p) \frac{D p}{D t} - \rho \sum_k \xi_k \omegadot_k + \rho \Hnuc \right) + \sum_k p_{X_k} \omegadot_k \right) \\ &=& \frac{1}{\rho p_\rho} \left( \frac{p_T}{\rho c_p}(1 - \rho h_p) - 1 \right) \frac{D p}{D t} + \frac{1}{\rho p_\rho} \left( \frac{p_T}{\rho c_p} (\rho \Hnuc - \rho \sum_k \xi_k \omegadot_k) + \sum_k p_{X_k} \omegadot_k \right) .\end{aligned}\end{split}

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 $$\Ub \cdot \nabla p_0$$, and the divergence constraint can be written

(399)$\nabla \cdot \Ub + \alpha \Ub \cdot \nabla p_0 = \frac{1}{\rho p_\rho} \left( \frac{p_T}{\rho c_p} \left( - \sum_k\rho \xi_k \omegadot_k + \rho \Hnuc \right) + \sum_k p_{X_k} \omegadot_k \right) \equiv \tilde{S} \label{eq:full_divu_constraint} ,$

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) . \label{eq:alphadef}$

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{D e}{Dt} = \rho \left( e_p \frac{D p}{Dt} + e_T \frac{D T}{Dt} + \sum_k e_{X_k} \frac{D X_k}{Dt} \right) = -p \nabla \cdot \Ub + \rho \Hnuc$

we can express $$DT/Dt$$ in terms of

(402)\begin{aligned} \frac{DT}{Dt} &=& \frac{1}{\rho e_T} \left( -p \nabla \cdot \Ub + \rho \Hnuc - \rho e_p \frac{D p}{Dt} - \rho \sum_k e_{X_k} \omegadot_k \right)\end{aligned}

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

(403)$-\rho \nabla \cdot \Ub = \frac{1}{p_\rho} \left( \frac{D p}{Dt} - \frac{p_T}{\rho e_T} \left( -p \nabla \cdot \Ub + \rho \Hnuc - \rho e_p \frac{D p}{Dt} - \rho \sum_k e_{X_k} \omegadot_k \right) - \sum_k p_{X_k} {\omegadot}_k \right) ,$

(404)$\left(-\rho - \frac{p p_T}{\rho e_T p_\rho} \right) \nabla \cdot \Ub = \frac{1}{p_\rho} \left( (1 + \frac{e_p p_T}{e_T}) \frac{D p}{Dt} - \frac{p_T}{\rho e_T} \left( \rho \Hnuc - \rho \sum_k e_{X_k} \omegadot_k \right) - \sum_k p_{X_k} {\omegadot}_k \right) ,$

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 $$\omegadot_k = \Hnuc = 0$$ for simplicity, then the constraint as derived using $$h$$ can be written

(405)$\nabla \cdot \Ub + \left( \frac{(1 - \rho h_p )p_T - \rho c_p}{\rho^2 c_p p_\rho} \right) \frac{D p_0}{D t} = 0$

and the constraint derived using $$e$$ can be written

(406)$\nabla \cdot \Ub + \left( \frac{\rho e_T + \rho e_p p_T}{\rho^2 e_T p_\rho + p p_T} \right) \frac{D p_0}{D t} = 0$

We note that if we evaluate both constraints for $$p = \rho R T,$$ 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)$\nabla \cdot \Ub + \frac{1}{\gamma p} \frac{D p_0}{D t} = 0$

Enthalpy vs Energy Equation

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

(408)$\frac{\partial(\rho h)}{\partial t} = -\nabla\cdot(\rho h\Ub) + \frac{Dp}{Dt} + \rho\Hnuc \label{eq:enthalpy}$

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{split}\begin{aligned} \frac{\partial(\rho h)}{\partial t} &=& -\nabla\cdot(\rho h\Ub) + \frac{Dp_0}{Dt} + \rho\Hnuc \nonumber \\ % \frac{\partial(\rho e)}{\partial t} + \frac{\partial p}{\partial t} &=& -\nabla\cdot(\rho e\Ub) -\nabla\cdot(p\Ub) + \frac{Dp_0}{Dt} + \rho\Hnuc \nonumber \\ % \frac{\partial(\rho e)}{\partial t} &=& -\nabla\cdot(\rho e\Ub) - p\nabla\cdot\Ub + \rho\Hnuc + \left \{ \frac{Dp_0}{Dt} - \frac{Dp}{Dt} \right \} \nonumber \end{aligned}\end{split}

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)$\frac{\partial(\rho e)}{\partial t} + \nabla\cdot(\rho e\Ub) = -p \; \nabla\cdot\Ub + \rho\Hnuc + \left \{ \frac{Dp_0}{Dt} - \frac{Dp}{Dt} \right \} \nonumber$

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

(411)$\frac{\partial(\rho e)}{\partial t} + \nabla\cdot(\rho e\Ub) = p_0 \nabla\cdot\Ub + \rho\Hnuc \nonumber$

The second equation subtracted from the first gives:

(412)$\frac{D (p_0 - p)}{Dt} - (p_0 - p) \nabla \cdot \Ub = 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{split}\begin{aligned} \frac{\partial(\rho h)}{\partial t} &=& -\nabla\cdot(\rho h\Ub) + \frac{Dp_0}{Dt} + \rho\Hnuc \nonumber \\ % \frac{\partial(\rho e)}{\partial t} + \frac{\partial p_0}{\partial t} &=& -\nabla\cdot(\rho e\Ub) -\nabla\cdot(p_0 \Ub) + \frac{Dp_0}{Dt} + \rho\Hnuc \nonumber \\ % \frac{\partial(\rho e)}{\partial t} &=& -\nabla\cdot(\rho e\Ub) - p_0 \nabla\cdot\Ub + \rho\Hnuc \nonumber\end{aligned}\end{split}

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{split}\begin{aligned} \frac{\partial\rho h}{\partial t} + \nabla\cdot\left(\rho h\Ub\right) &=& \frac{Dp_0}{Dt} + \rho H {} \nonumber\\ \frac{\gamma}{\gamma-1}\frac{\partial p}{\partial t} + \frac{\gamma}{\gamma-1}\nabla\cdot\left(p\Ub\right) &=& \frac{Dp_0}{Dt} + \rho H {} \nonumber\\ \frac{Dp}{Dt} &=& -p\nabla\cdot\Ub + \frac{\gamma}{\gamma-1}\left(\frac{Dp_0}{Dt}+\rho H\right) \end{aligned}\end{split}

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

(415)$\frac{Dp}{Dt} = -p\nabla\cdot\Ub - \left(\gamma-1\right)p_0\nabla\cdot\Ub + \rho H$

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

(416)$\nabla\cdot\Ub = \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}\left(1-\frac{p}{p_0}\right)\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}\left(1-\frac{p}{p_0}\right)\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

We believe that the original desire stems from our experience with smallscale combustion. There, stratification is not important and $$Dp_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)$\nabla \cdot(\beta_0\Ub) = \beta_0\left(S-\frac{1}{\overline{\Gamma_1} p_0} \frac{\partial p_0}{\partial t} - \frac{f}{\overline{\Gamma_1} p_0} \frac{p_0-p_\text{EOS}}{\dt}\right)$
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)$$?