3.2.4.12. Porous Flow Equations

In what follows we use a subscript $\alpha$ to denote a component index and $N_\alpha$ is the total number of components. We use a subscript $\beta$ to denote a phase and the total number of phases is $N_\beta$ . The subscript $s$ is used to refer to the solid phase of the porous skeleton; this phase is not included in $N_\beta$ . Summations over $\alpha$ and $\beta$ are assumed to range from one to $N_\alpha$ and $N_\beta$ , respectively, unless otherwise specified.

3.2.4.12.1. Generalized Equations

This write up follows, directly, the reports of [8], [9] and [10].

Porosity and Void Ratio

In modeling of porous materials, one often needs to know the relative volume of pore space ( $V_p$ ) and the solid porous skeleton ( $V_s$ ). The most common measure of this is the porosity, $\phi$ , which is defined as the ratio of the volume of pore space ( $V_p$ ) per unit total volume ( $V_t = V_p + V_s$ ), viz.,

$\phi = \frac{V_p}{V_t} = \frac{V_p}{V_p + V_s}.$

The void ratio ( $r$ ) is sometimes used as well and is defined as the ratio of the pore volume to the solid volume, viz.,

$r = \frac{V_p}{V_s}$

Clearly the porosity and void ratio are related,

$r &= \frac{\phi}{1 - \phi} \\ \phi &= \frac{r}{1 + r}$

Component Mass Balances

The mass balance for each component across all phases can be written

(3.42) $\pt{\massa} + \div\Fa = Q_\alpha \qquad \alpha = 1, \ldots, N_\alpha$

where $\massa$ is mass concentration of component $\alpha$ across all phases, $t$ is time, $\Fa$ is the total mass flux of component $\alpha$ relative to a fixed frame of reference and $Q_\alpha$ is the volumetric source of component $\alpha$ .

The mass concentration of component $\alpha$ can be written as

(3.43) $\massa = \phi\sum_\beta\Yab\denb S_\beta$

where $\phi$ is the interstitial volume fraction, a.k.a. the porosity, that ranges from zero to one, $\Y$ denotes mass fraction so that $\Yab$ is the mass fraction of component $\alpha$ in phase $\beta$ , $\denb$ is mixture density of phase $\beta$ and $S_\beta$ is saturation which is defined as the fraction of the interstitial pore space occupied by phase $\beta$ .

For completeness it is noted here that the mass fractions sum to one in each phase and that the saturations of all phases sum to unity,

(3.44) $\sum_\alpha \Yab = 1$

(3.45) $\sum_\beta S_\beta = 1$

The macroscopic mass flux in phase $\beta$ , $\denb\v_\beta$ is obtained from the extended Darcy’s Law,

(3.46) $\denb\v_\beta = \fb = -\frac{\denb \krb}{\viscb}\K\bcdot\left(\grad P_\beta + \denb\g\right)$

where $\krb$ is the relative permeability of phase $\beta$ and $\K$ is the permeability tensor which are both intrinsic properties of the porous material and so may vary spatially. The mixture viscosity of phase $\beta$ is $\viscb$ , $P_\beta$ is the pressure in phase $\beta$ and $\g$ is the gravitational acceleration.

The component mass flux, $\Fa$ can be written as the sum of the component mass fluxes across each phase,

(3.47) $\Fa = \sum_\beta \Fab.$

Within each phase, the component mass flux can be decomposed in its convective and diffusive components,

(3.48) $\Fab = \Yab\fb + \Jab$

where $\Jab$ is the diffusive flux of component $\alpha$ in phase $\beta$ . For example, for a gaseous, dilute binary mixture the diffusive flux may be modeled using Fick’s Law,

(3.49) $\Jab = -\denb \Dab \grad \Yab \qquad \mtext{(dilute binary gas mixture)}$

where $\Dab$ is the diffusion coefficient of component $\alpha$ in phase $\beta$ .

Porous Species

The porous species equation is identical to the mass balance equation, with all quantities being molar concentrations. The porous species equation is written as

(3.50) $\pt{c_\alpha} + \div \vector \Pi_\alpha = R_\alpha \qquad \alpha = 1, \ldots, N_\alpha,$

where $c_\alpha$ is the molar concentration, defined analogously to (3.43), $\Pi_\alpha$ is the molar flux, and $R_\alpha$ is the volumetric molar source of component $\alpha$ . The porous species equation exists to facilitate using expressions derived for use with the non porous-flow species equation. An analogous set of terms to those found in the mass balance equation can be derived on a molar concentration basis.

Energy Balance

In order to accommodate multicomponent transport problems we employ a enthalpy based energy conservation law,

(3.51) $\pt{e} + \div\q = Q_e$

where $e$ is the total internal energy, $\q$ is the total heat flux relative to a stationary frame of reference and $Q_e$ is the volumetric source of energy. The total internal energy is written as a sum of the energies of all components across all phases as well as the internal energy of the solid porous structure,

(3.52) $e = \left(1-\phi\right)\dens e_s + \phi \sum_\beta S_\beta \denb e_\beta$

where $\denb$ and $e_\beta$ are the mixture density and internal energy of phase $\beta$ . Generally, thermodynamic relations are required to relate these quantities to the component specific quantities.

The total heat flux can be decomposed into different modes of transport. Typically,

(3.53) $\q = -\lambda_T\grad T + \sum_\beta \sum_\alpha \hab \Fab$

where $\lambda_T$ is the effective thermal conductivity, $T$ is the temperature and $\hab$ is the partial enthalpy of component $\alpha$ in phase $\beta$ which is interpreted as the enthalpy of pure component $\alpha$ at the temperature and pressure of phase $\beta$ . Using equation (3.48), we can expand equation (3.53) such that

$\sum_\alpha \hab \Fab &= \sum_\alpha \hab \left(\Yab\fb + \Jab\right) \nonumber\\ &= \hb\fb + \sum_\alpha \hab\Jab$

where $\hb = \sum_\alpha\hab\Yab$ is the mixture enthalpy of phase $\beta$ . Thus, the total heat flux can be expressed as the sum of conductive, convective and diffusive contributions,

(3.54) $\q = -\lambda_T\grad T + \sum_\beta \denb \v_\beta \hb + \sum_\beta\sum_\alpha \hab\Jab$

3.2.4.12.2. Porous Flow Systems

To facilitate various flavors of single phase and multiphase flow in Aria, the porous flow system model was created to manage expression creation. In addition to the single phase model (Single Phase, Single Component, Isothermal Flow and Single Phase, Single Component, Non-isothermal Flow), there are several others available including: mixed single phase, CO2-salt-brine, two-phase-immiscible (Two-Phase Immiscible and Pressure-pressure Formulation for Two-phase Immiscible), Two-phase Air-Water, and Organic Material Decomposition (OMD). The basic purpose of the porous flow system model is to create the necessary expressions for each system. There are examples of each flow system in the regression tests, which are the best place to start when developing an input for a similar system. The input syntax for a porous flow system is given in Material Properties in addition to the required parameters for each porous flow system.

Single Phase, Single Component, Isothermal Flow

The simplest possible case is for a single phase flow ( $N_\beta=1$ , $S=1$ ) with a single component ( $N_\alpha=1$ , $Y=1$ ) at a fixed temperature. This is the single phase porous flow system. For a single component flow, there is no diffusive mass flux so $\J=\zerovector$ and $\F = \f$ where

(3.55) $\F = \f = \den\v = -\frac{\den \kr}{\mu}\K\bcdot\left(\grad P + \den\g\right)$

In this case, the mass balance equation is all that is needed to solve for the unknown pressure field $P$ .

(3.56) $\pt{d} + \div\f = Q$

Since the density is most likely a function of temperature, the constant uniform temperature will have to be specified. Finally, in this special case, equation (3.43) reduces to $d = \phi\den$ .

Single Phase, Single Component, Non-isothermal Flow

Building on the simple case developed in the section above (Single Phase, Single Component, Isothermal Flow), the next step is to accommodate a non-isothermal flow. Here, we still have a single phase flow ( $N_\beta=1$ , $S=1$ ) with a single component ( $N_\alpha=1$ , $Y=1$ ). For a single component flow, there is no diffusive mass flux so $\J=\zerovector$ and $\F = \f$ where

(3.57) $\F = \f = \den\v = -\frac{\den \kr}{\mu}\K\bcdot\left(\grad P + \den\g\right)$

Again, the mass balance equation is all that is needed to solve for the unknown pressure field $P$ .

(3.58) $\pt{d} + \div\f = Q$

Again, in this special case, equation (3.43) reduces to $d = \phi\den$ .

In order to solve for the unknown temperature, $T$ , we solve an energy transport equation in the enthalpy form as described in Energy Balance. Here, the energy flux in equation (3.53) reduces to

$\q = -\lambda_T\grad T + \den\v\h$

where $h$ is the pure component enthalpy which may be a function of $P$ and $T$ .

Two-Phase Immiscible

In this section, we define a formulation for two-phase flow of immiscible fluids, that is, without mass transfer between the phases. Common examples include oil/water, under pressure and temperature conditions where phase transitions are negligible, and systems involving non- aqueous phase liquids (NAPL) occurring in subsurface contamination by organic liquids, e.g., solvents, fuels, oils. If we disallow partitioning of components across phases, then we get the following system of mass balance equations for a wetting and non-wetting phase,

$\frac{\partial}{\partial t} (\rho_w \phi (1-S_n)) = \grad \cdot \left( \rho_w \frac{k_{rw}}{\mu_{w}} \K \cdot(\grad p - \rho_w \g)\right) + Q_w$

$\frac{\partial}{\partial t} (\rho_n \phi S_n) = \grad \cdot \left( \rho_n \frac{k_{rn}}{\mu_{n}} \K \cdot(\grad p + \grad p_c - \rho_n \g)\right) + Q_n$

In these equations, the subscript w represents wetting phase properties and the subscript n nonwetting phase properties. We have also incorporated the Darcy flux terms, including pressure and gravitational forces, as well as capillary pressure:

p_c = p_n - p_w

The symbol $p$ in the equations is used in place of wetting phase pressure, pw. Also, the pore space is assumed saturated with fluids,

S_w + S_n = 1

This constraint has been incorporated in the wetting phase mass balance. In addition to the density models described in the previous section, models for density and viscosity can be specified as constant, polynomial with respect to other primary or secondary variables, or by tabular functions, or a user-defined plug-in.

Both the capillary pressure and relative permeability are functions of phase saturations. Either of these can be specified by the user in a tabular form, as a function of either phase saturation. A plug-in can also be written by the user for either of these functions. Several commonly used functions are built-in to the code.

Pressure-pressure Formulation for Two-phase Immiscible

The $P_w-S_n$ formulation is generally known to have superior numerical behavior for $S_w \rightarrow 0$ , however, the capillary pressure in heterogeneous media creates a problem for a vertex-centered (e.g. FEM) numerical scheme. In FEM, because a node straddles a material interface, which separates two materials with different texture and therefore different wetting properties, the capillary pressure is continuous across the interface, but the saturation is not. So, a ( $P_w$ ,:math:S_n) formulation in FEM would require dealing with discontinuous DOF for $S_n$ . One alternative is to instead solve for two pressures. Two pairs of pressure-based DOFs available in Aria include ( $P_w,P_c$ ) and ( $P_w,P_n$ ). (Note that in single phase regions, where $S_w = 0$ or 1, then $P_w$ or $P_n$ are undefined. The formulation should work for the thermal battery problem so long as $S_w>0$ .

The following system of mass balance equations for a wetting and non-wetting phase describes an immiscible flow,

$\frac{\partial}{\partial t} (\rho_w \phi S_w) = \grad \cdot \left( \rho_w \frac{k_{rw}}{\mu_w}\K \cdot (\grad P_w - \rho_w \g)\right) + Q_w$

$\frac{\partial}{\partial t} (\rho_n \phi S_n) = \grad \cdot \left( \rho_n \frac{k_{rn}}{\mu_n}\K \cdot (\grad P_w + \grad P_c - \rho_n \g)\right) + Q_n$

In these equations, the subscript w represents wetting phase properties and the subscript n non-wetting phase properties. We have also incorporated the Darcy flux terms, including pressure and gravitational forces, as well as capillary pressure:

p_c = p_n - p_w

Also, the pore space is assumed saturated with fluids,

S_w + S_n = 1

The current implementation in Aria uses a van Genuchten function to model the relationship between the capillary pressure and the saturation. This relationship is inverted for saturation given the capillary pressure in the pressure-pressure formulation.

Two-phase Air-Water

This model describes the nonisothermal two-component, two-phase transport of water and air in a porous medium. The air is treated as a non-condensible gas, but it can dissolve in the liquid phase. The water can exist as water vapor or liquid. The model allows the fluids to partition between the two phases, depending on pressure and temperature conditions. Two single phase states (all liquid or all gas) and one two-phase state are allowed.

The component mass balance equations describing two-phase transport of water, air and energy are:

$\frac{\partial}{\partial t} ( \phi (S_l Y_{wl} \rho_l + S_g Y_{wg} \rho_g)) + \grad \cdot (Y_{wl} \rho_l \v_l + Y_{wg} \rho_g \v_g + \J_{wg}) + Q_w = 0$

$\frac{\partial}{\partial t} ( \phi (S_l Y_{al} \rho_l + S_g Y_{ag} \rho_g)) + \grad \cdot (Y_{al} \rho_l \v_l + Y_{ag} \rho_g \v_g + \J_{ag}) + Q_a = 0$

$\frac{\partial}{\partial t}\left[ (1-\phi)\rho_s e_s + \phi(\rho_l S_l e_l + \rho_g S_g e_g)\right] \\ + \grad \cdot (-\lambda_T \grad T + \rho_l \v_l h_l + \rho_g \v_g h_g + h_{wg} \J_{wg} + h_{ag} \J_{ag} ) \\ + Q_e = 0$

where

$\J_{ag} = - \rho_g D_{wa} \grad Y_{ag}$

is the gas-phase binary diffusion flux for water vapor through air. In these equations, subscript $l$ refers to the liquid phase and $g$ to the gas phase. In addition to previously defined variables: $Y_{\alpha\beta}$ is the mass fraction of component $\alpha$ ( $w$ for water, $a$ for air) in phase $\beta$ , $S_{\beta}$ is the saturation of phase $\beta$ , $\rho_{\beta}$ is the phase density, $e_{\beta}$ is the internal energy of phase $\beta$ , $\lambda_{T}$ is the effective thermal heat conduction, $h_\beta$ is the enthalpy of phase $\beta$ and $h_{\alpha g}$ is the component enthalpy in the gas phase.

The liquid and gas phase Darcy velocities are

$\v_l = -\frac{k_{rl}}{\mu_{l}} \K \cdot (\grad P_l + \rho_l \g)$

$\v_g = -\frac{k_{rg}}{\mu_g} \K \cdot (\grad P_g + \rho_g \g)$

In addition, the following constraints hold:

$S_l + S_g = 1 \; ; \quad Y_{w\beta} + Y_{a\beta} = 1 \; , \text{for } \beta = l \text{ or } g$

CO2-salt-brine

The CO2-salt-brine system is intended to model the interaction of CO2 gas, solid salt, and liquid brine. In this system there are three mass balance equations that must be solved for the saturations of each. This formulation describes a 3-component 2-phase flow for the water-CO2-NaCl system. The two phases are brine, denoted by subscript w, a liquid phase of water with salt and CO2 dissolved in it, and “gas” denoted by subscript g, a CO2-rich phase, with some water dissolved in it, but no salt. The formulation also deals with precipitation of salt from the liquid phase.

To deal with the precipitation of NaCl from the liquid to form solid salt, the fluid system is treated as a three-phase system, with two flowing phases composed of the brine and CO2 phases, and one non-flowing solid phase, which is the solid salt precipitate. Obviously, the formulation holds for other chloride salts, e.g. CaCl, KCl.

$\frac{\partial}{\partial t} \phi (S_w \rho_w x_w + S_g \rho_g y_w) + \grad \cdot (\rho_w x_w \v_w + \rho_g y_w \v_g + \J_w^{H2O} + \J_g^{H2O} ) = Q^{H2O}$

$\frac{\partial}{\partial t} \phi (S_w \rho_w x_{CO2} + S_g \rho_g y_{CO2}) + \grad \cdot (\rho_w x_{CO2} \v_w + \rho_g y_{CO2} \v_g + \J_w^{CO2} + \J_g^{CO2} = Q^{CO2}$

$\frac{\partial}{\partial t} \phi (S_w \rho_w x_{NaCl} + S_g \rho_s) + \grad \cdot (\rho_w x_{NaCl} \v_w + \J_w^{NaCl}) = Q^{NaCl}$

In the equations above, we have used the following nomenclature, $S_s$ is the saturation of solid salt precipitate, $S_w$ is the brine phase saturation, $S_g$ is the CO2-rich phase saturation, $x_w$ is the mass fraction of water in the brine phase, $x_{CO2}$ is the mass fraction of CO2 in the brine phase, $x_{NaCl}$ is the mass fraction of salt dissolved in the brine phase, $y_{CO2}$ is the mass fraction of CO2 in the CO2-rich (“gas”) phase, $y_w$ is the mass fraction of water in the gas phase, $\rho_s(T,p)$ is the density of solid precipitate, and $\J_{\beta}^{\alpha}$ is the diffusion flux of component $\alpha$ in phase $\beta$ which can also include dispersion. The subscript ‘w’ refers to the brine phase, ‘g’ refers to the CO2 phase and ‘s’ refers to the solid salt phase.

This model treats all the separate-phase CO2 as a single CO2-rich phase. The phase diagram for the CO2-H20 system displays a saturation line where liquid water, liquid CO2 and a vapor phase can co-exist. The saturated water vapor pressure on this line is small, so the vapor phase is mostly CO2.

$S_w + S_g + S_s = 1 \; ; \quad x_w + x_{CO2} + X_{NaCl} = 1 \; ; \quad y_w + y_{CO2} + y_{NaCl} = 1$

We also have the constraint on the over-all composition, which can replace either of the last two constraints above,

$z_w + z_{CO2} + z_{NaCl} = 1$

To solve for the pressure, in all compositional situations, we can solve a pressure equation, obtained by summing all component balances. Since the resulting equation is not independent, it replaces one of the three mass balance equations presented earlier.

$\frac{\partial}{\partial t} \phi (S_w \rho_w + S_g\rho_g + S_s\rho_s) + \grad \cdot(\rho_w\v_w + \rho_g \v_g) + Q^{H2O} + Q^{CO2} + Q^{NaCl} = 0$

In the numerical implementation, this equation could be obtained after all BCs have been applied to the discrete equations, by summing the three residual equations together.

Recall, the primary variables for this system are defined as a set of variables that allow one to compute all the remaining (secondary) variables appearing in the mass and energy balance equations. For this problem, the primary solution vector (assuming we will also solve an energy balance for $T$ ) is $[P,z_{CO2},\rho_{NaCl},T]$ , where

$z_{CO2} = \frac{S_w\rho_w x_{CO2} + S_g \rho_g y_{CO2}}{S_w \rho_w + S_g \rho_g} = (1-\nu) x_{CO2} + \nu y_{CO2}$

is the over-all mass fraction of CO2 in the flowing phases. We have introduced the fraction of gas, defined as,

$\nu = \frac{\rho_g S_g}{\rho_w S_w + \rho_g S_g}$

Similarly, the overall density of salt is defined by,

$\rho_{NaCl} = \rho_w S_w x_{NaCl} + \rho_s S_s = (1-\nu) x_{NaCl}F + \rho_s S_s$

where we have introduced, for later convenience, the following mass density,

$F = \rho_w S_w + \rho_g S_g$

Note that the accumulation term in the CO2 balance can be written as $\phi F z_{CO2}$ .

In this formulation, $P$ is solved for from the pressure equation, $z_{CO2}$ from the CO2 mass balance equation, $\rho_{NaCl}$ from the NaCl mass balance and $T$ from the energy balance equation. We assume the phase equilibrium for this system is done ahead of time.

The thermodynamics algorithm for this system is still in development and has not yet been completed.

Organic Material Decomposition (OMD)

Organic material decomposition (OMD) models currently define two-phase (solid and gas) problems in which the solid (or porous) species react to produce additional gas species that pressurize a system and volumetric heat transfer between phases. The notation used above does not include the porous solid as a phase, but since for OMD we are solving equations for the solid phase including reactions, we set $N_{\beta}=2$ , with it being understood that one of these phases is the solid phase. OMD problems will extend to include a liquid phase in the future. The OMD equations described here are based on the combustion model [11]. Further details about coupling boundary conditions between fluid-only and porous regions are provided in Coupling BCs for Organic Material Decomposition. There are separate mass balance equations for the solid phase and gas phase. The solid phase mass balance equation is used to solve for solid phase density:

(3.59) $\underbrace{\frac{\partial \bar{\rho}}{\partial t}}_{MASS} = \underbrace{-\dot{\omega}'''_{f,g} }_{SRC}.$

where $\bar{\rho}$ is the solid phase density and $\dot{\omega}'''_{f,g}$ is the rate of formation of gases due to heterogenous reaction kinetics. The mass balance equation for the gas phase solves for the pressure (from the gas phase density via the ideal gas law),

(3.60) $\underbrace{\frac{ \partial (\psi \rho_g) }{ \partial t }}_{MASS} + \underbrace{\frac{ \partial (\rho_g u_{j,g}) }{ \partial x_j }}_{ADV} = \underbrace{\dot{\omega}'''_{f,g}}_{SRC},$

where $\rho_g$ is the gas-phase density and $\psi$ is the mixture-averaged solid-phase porosity. The porosity is the void space not taken up by the volume fractions of each solid phase species, typically defined as a function of the volume fractions of each solid species,

$\psi=1 - \sum_{\alpha} X_{\alpha}.$

The volume fraction of each solid species is calculated as

$X_{\alpha}=\frac{\bar{\rho}Y_{\alpha}}{\rho_{s0,\alpha}},$

where $\rho_{s0,\alpha}$ is the solid (non-porous) density for each species. In equation (3.60), $u_{j,g}$ is the gas-phase velocity and gets computed from Darcy’s approximation as

(3.61) $\psi v_{j,g}=u_{j,g} = -\frac{\bar{K_{i j}}}{\mu_g} \left( \frac{\partial p_g}{\partial x_i} + \rho_g g_i \right),$

where $\bar{K_{i j}}$ is the mixture-averaged solid-phase permeability tensor defined as the normalized volume average of the solid phase species’ permeabilities,

$\bar{K_{i j}}=\frac{\sum_{\alpha} X_{\alpha} K_{\alpha j}}{\sum_{\alpha} X_{\alpha}}$

where

$K_{\alpha j}=K_{ref}\left( \frac{T}{T_{ref}} \right)^n\delta_{\alpha j}$

The intrinsic permeability is the product of the intrinsic permeability scaling and the intrinsic permeability. The reason it is separated in this way is that scalar models have access to many generic functions, whereas tensor entries can only be constants. A more complex function for the intrinsic permeability tensor would have to be written as a user subroutine.

In equation (3.61) the parameter $\mu_g$ is the gas-phase viscosity and $g_j$ is the gravity vector. The ideal gas law gets used for the gas density in the pressure form of equation (3.60)

$\rho_g=\frac{P\bar{M}}{RT_g}.$

Therefore the actual equation for pressure being solved here is

(3.62) $\frac{\partial}{\partial t} \left( \frac{ p_g \bar{M} \psi}{RT_g} \right) + \frac{\partial}{ \partial x_j } \left[ \frac{ p_g \bar{M}}{RT_g} \frac{\bar{K_{\alpha j}}}{\mu_g} \left( \frac{\partial p_g}{\partial x_j} + \frac{p_g\bar{M}}{RT_g} g_j \right) \right] = 0.$

where $\bar{M}$ is the mass averaged molecular weight,

$\bar{M}=\sum_{\alpha} Y_{\alpha} M_{\alpha},$

where $Y_{\alpha}$ are the gas phase mass fractions ans $M_{\alpha}$ are the molecular weights of the gaseous species.

The solid phase species equations solve for the solid phase mass fractions, $Y_{\alpha,s}$ , of species $\alpha$

(3.63) $\underbrace{\frac{ \partial (\bar{\rho} Y_{\alpha}) }{ \partial t }}_{MASS} = \underbrace{\left( \dot{\omega}'''_{f,\alpha} - \dot{\omega}'''_{d,\alpha} \right)}_{SRC}.$

where $\dot{\omega}'''_{f,\alpha}$ is the formation rate of species $\alpha$ and $\dot{\omega}'''_{d,\alpha}$ is its destruction rate. The gas-phase species equations solve for the gas-phase mass fraction $Y_{\alpha,g}$ of species $\alpha$

(3.64) $\underbrace{\frac{ \partial (\psi \rho_g Y_{\alpha ,g}) }{ \partial t }}_{MASS} + \underbrace{\frac{ \partial (\rho_g u_{j,g} Y_{\alpha ,g}) }{ \partial x_j }}_{ADV} + \underbrace{\frac{ \partial q_{\alpha j}^{Y,g}}{ \partial x_j}}_{DIFF} = \underbrace{\dot{\omega_{\alpha,g}}}_{SRC}$

No source terms are included in this model. The gas-phase species diffusion flux vector $q_{\alpha j}^{Y,g}$ is modeled as

$q_{\alpha j}^{Y,g} = -\psi \rho_g D_{\alpha,g} \frac{\partial Y_{\alpha,g}}{\partial x_j},$

where $D_{\alpha,g}$ is the gas-phase mass diffusivity for species $\alpha$ .

Enthalpy equations are used to solve for temperature, and OMD type problems allow for different temperature fields in the gas and solid phases, so these are described separately here.

The solid phase enthalpy equation for the solid phase temperature $\bar{T}$ is given by

(3.65) $\underbrace{\frac{ \partial ( \bar{\rho} c_p\bar{T}) }{ \partial t }}_{MASS} = \underbrace{-\frac{ \partial q_j^{h}}{ \partial x_j}}_{DIFF} + \underbrace{h_{cv}\left( T_g-\bar{T} \right)}_{SRC}$

There is no advection term here since the solid phase is stationary. The solid phase energy diffusive flux vector $q_j^{h,g}$ is given by Fourier’s law,

$q_j^{h,g} = -\kappa\frac{\partial \bar{T}}{\partial x_j},$

where $\kappa$ is the thermal conductivity.

The gas-phase enthalpy equation within a porous region, to be solved for the gas-phase enthalpy (and corresponding temperature) $h_g=c_p T_g$ , is

(3.66) $\underbrace{\frac{ \partial (\psi \rho_g h_g) }{ \partial t }}_{MASS} + \underbrace{\frac{ \partial (\rho_g u_{j,g} h_g) }{ \partial x_j }}_{ADV} = \underbrace{-\frac{ \partial q_j^{h,g}}{ \partial x_j}}_{DIFF} + \underbrace{\left( \frac{\partial \psi p_g }{\partial t} \right)}_{SRC} + \underbrace{h_{cv}\left( \bar{T} - T_g \right)}_{SRC}$

where $h_g$ is the mixture-averaged gas-phase enthalpy, $h_{cv}$ is the volumetric heat transfer coefficient which gets set outside of the solid/gas phases of the material definition, $\bar{T}$ is the porous solid-phase temperature and $T_g$ is the gas phase temperature. The gas-phase energy diffusive flux vector $q_j^{h,g}$ is modeled as

$q_j^{h,g} = -\psi \rho_g D_g \frac{\partial h_g}{\partial x_j},$

where $D_g$ is the mixture-averaged gas-phase mass diffusivity. The compressible pressure source term for the enthalpy equation has been derived in the context of Darcy flow. The volumetric transfer source accounts for the heat transfer between the solid and gas phases.

3.2.4.12.3. Porosity Models

There are several models for porosity in the porous flow capabilities in Aria. The following list provides a brief description for each.

Mesh Deforming

In the mesh (or solid) deforming porosity model, the expression for the porosity is given as

$\phi = 1.0 - \frac{b (1.0 - \phi_0)}{\text{det(F)}}$

where $\phi_0$ is the initial porosity, $b$ is the Biot coefficient, and $\text{F}$ is the mesh deformation gradient.

Rock Compressible

For the Rock Compressible porosity model, the porosity is given as

$\phi = \phi_0(1.0 + Cr(p - p_{ref}))$

where again, $\phi_0$ is the initial porosity, $Cr$ is a rock compressibility coefficient, $p$ is the pore pressure, and $p_{ref}$ is a reference pore pressure defined in the input parameters.

Coussy

For the Coussy porosity model, the expression for porosity is

$\phi = \phi_0 + \frac{1}{M}(p - p_{ref}) + b \epsilon_u$

where $M$ represents the formation compressibility, $\epsilon_u = \text{div}(\vector{u})$ is the divergence of the solid deformation vector, $b$ is the Biot coefficient, $p$ is the pore pressure, and $p_{ref}$ is a reference pore pressure.