3.2.19. Porous-Fluid Coupling Algorithm

This section provides a brief overview of the current porous/fluid coupling algorithm, as it is intended to be used in simulations of composite fires using Fuego to model the fluid region and coupling to Aria to model the porous region.. This is a loosely-coupled algorithm, relying on framework interpolation transfers of nodal fields between the porous region and the low-Mach fluid region and region-region Picard loops to converge the overall problem within a timestep.

Note that the shorthand is adopted where the porous region is described as region $A$ and the low-Mach free fluid region is described as region $B$ , with the interface between them referred to as $\Gamma_{AB}$ and other boundaries not a part of this interface are referred to as $\Gamma \backslash \Gamma_{AB}$ .

3.2.19.1. Fluid Flow

3.2.19.1.1. Bulk Equations

paragraph{Porous Continuity Equation} The porous region contains a condensed phase (the solid skeleton of the porous system) and a gas phase occupying the pores of the condensed phase. The condensed phase is not discussed explicitly in this description, although it interacts with the gas phase through things like its permeability and porosity, and its decomposition which can produce gas-phase mass through chemical source terms.

The porous gas-phase continuity equation within a porous region, to be solved for the gas-phase pressure $p_g$ , is

(3.534) $\frac{ \partial (\bar{\psi} \rho_g) }{ \partial t } + \frac{ \partial (\rho_g u_{j,g}) }{ \partial x_j } = \dot{\omega}'''_{fg},$

where $\bar{\psi}$ is the mixture-averaged condensed-phase porosity, $\rho_g$ is the gas-phase density, and $u_{j,g}$ is the gas-phase velocity vector computed from Darcy’s approximation as

(3.535) $u_{j,g} = -\frac{\bar{K}}{\mu_g} \left( \frac{\partial p_g}{\partial x_j} + \rho_g g_j \right),$

where $\bar{K}$ is the mixture-averaged condensed-phase permeability, $\mu_g$ is the gas-phase viscosity, and $g_j$ is the gravity vector. The term $\dot{\omega}'''_{fg}$ represents the formation rate of gas-phase mass from the condensed phase.

Multiplying (3.534) by an arbitrary test function $w$ and integrating over the domain $\Omega$ while integrating the advection term by parts, yields the variational form of the continuity equation that is solved for $p_g$ using the Galerkin finite element method,

(3.536) $\int\limits_\Omega w \left( \frac{ \partial (\bar{\psi} \rho_g) }{\partial t} - \dot{\omega}'''_{fg} \right) d\Omega - \int\limits_\Omega \frac{\partial w}{\partial x_j} \rho_g u_{j,g} d\Omega + \int\limits_\Gamma w \rho_g u_{j,g} n_j d\Gamma = 0,$

where $n_j$ is the boundary surface normal. The boundary flux term is then split into contributions on the interface between regions $A$ and $B$ and off the interface so that they may be treated separately. The continuity equation then takes the form

(3.537) $\int\limits_\Omega w \left( \frac{ \partial (\bar{\psi} \rho_g) }{\partial t} - \dot{\omega}'''_{fg} \right) d\Omega - \int\limits_\Omega \frac{\partial w}{\partial x_j} \rho_g u_{j,g} d\Omega$

(3.538) $+ \int\limits_{\Gamma \backslash \Gamma_{AB}} w \rho_g u_{j,g} n_j d\Gamma + \int\limits_{\Gamma_{AB}} w F_A d\Gamma = 0,$

where $F_A$ is the imposed flux on the porous side ( $A$ ) of the $\Gamma_{AB}$ interface. A detailed description of the coupling boundary flux is given in Section Coupling Boundary Conditions

Low-Mach Continuity Equation

The continuity equation within the low-Mach fluid region, to be solved for the pressure $p$ , is

(3.539) $\frac{ \partial \rho }{ \partial t } + \frac{ \partial \rho u_j }{ \partial x_j } = S,$

where $\rho$ is the fluid density, $u_j$ is the fluid velocity, and $S$ is a generic mass volumetric source term. Integrating (3.539) over a CVFEM control volume and using the Gauss divergence theorem on the advection and diffusive flux terms, yields the integral form of the continuity equation that is solved,

(3.540) $\int\limits_\Omega \left( \frac{ \partial \rho }{\partial t} - S \right) d\Omega + \int\limits_\Gamma \rho u_j n_j d\Gamma = 0.$

Similar to the porous continuity equation, the boundary flux term is split into contributions both on and off the $\Gamma_{AB}$ interface, yielding

(3.541) $\int\limits_\Omega \left( \frac{ \partial \rho }{\partial t} - S \right) d\Omega + \int\limits_{\Gamma \backslash \Gamma_{AB}} \rho u_j n_j d\Gamma + \int\limits_{\Gamma_{AB}} F_B d\Gamma = 0.$

The interface coupling flux is described in Coupling Boundary Conditions.

Low-Mach Momentum Equation

The momentum equation within the low-Mach fluid region, to be solved for the velocity $u_i$ , is

(3.542) $\frac{ \partial \rho u_i }{ \partial t } + \frac{ \partial \rho u_j u_i }{ \partial x_j } = \frac{ \partial \sigma_{ij} }{ \partial x_j } + \rho g_i$

where the Cauchy stress tensor is given by

(3.543) $\sigma_{ij} = \tau_{ij} - p \delta_{ij}$

in terms of the viscous stress tensor

(3.544) $\tau_{ij} = \mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) - \frac{2}{3} \mu \frac{\partial u_k}{\partial x_k} \delta_{ij}.$

Integrating (3.542) over a CVFEM control volume and using the Gauss divergence theorem on the advection and diffusive flux terms, yields the integral form of the momentum equation that is solved,

(3.545) $\int\limits_\Omega \frac{ \partial \rho u_i }{\partial t} d\Omega + \int\limits_\Gamma \rho u_j u_i n_j d\Gamma - \int\limits_\Gamma \sigma_{ij} n_j d\Gamma - \int\limits_\Omega \rho g_i d\Omega = 0.$

Multiplying this equation by an arbitrary test function $w$ , integrating the advection and stress terms by parts, and splitting the boundary flux terms into on-interface and off-interface contributions yields

(3.546) $\int\limits_\Omega \frac{ \partial \rho u_i }{\partial t} d\Omega - \int\limits_\Omega \rho g_i d\Omega + \int\limits_{\Gamma \backslash \Gamma_{AB}} \rho u_j u_i n_j d\Gamma - \int\limits_{\Gamma \backslash \Gamma_{AB}} \sigma_{ij} n_j d\Gamma$

(3.547) $+ \int\limits_{\Gamma_{AB}} \rho u_j u_i n_j d\Gamma - \int\limits_{\Gamma_{AB}} \sigma_{ij} n_j d\Gamma = 0.$

3.2.19.1.2. Coupling Boundary Conditions

Coupling between the porous and fluid regions is achieved using an interface flux that is imposed as a Robin-style boundary condition. This approach has been used successfully in the past for coupling incompressible Darcy and Stokes flows [80]. Here we generalize the coupling for compressible fluids and Navier-Stokes flow.

The fluxes applied to the porous and fluid continuity equations at the interface $\Gamma_{AB}$ are

(3.548) $F_A = \dot{m}_B \cdot \hat{n} + \beta (p_A - p_B)$

(3.549) $F_B = \dot{m}_A \cdot \hat{n} + \beta (p_B - p_A),$

where $\dot{m}_A = \rho_g \vec{u_g}$ , $\dot{m}_B = \rho \vec{u}$ , and the free constant $\beta$ is computed as

(3.550) $\beta = c \frac{\bar{K} \rho_g}{\mu_g h}$

with $h$ being a measure of the mesh size adjacent to the interface, and $c$ a user-specified scaling coefficient. The same value of $\beta$ is used on both sides of the interface because that results in excellent mass conservation even on coarse meshes. If a different value of $\beta$ is used on each side the method is still convergent but worse mass conservation is observed when solving on under-resolved meshes. Some attempts have been made to use an averaged penalty coefficient of the form

(3.551) $\beta = \frac{\beta_A + \beta_B}{2}$

(3.552) $\beta_A = \frac{\bar{K} \rho_g}{\mu_g h}$

(3.553) $\beta_B = \frac{\tau^c}{h}, \text{ or } \beta_B = \frac{\mu}{h P_{ref}},$

however they resulted in an impractically large number of Picard iterations to converge for some test problems.

A distinguishing condition BC for velocity is applied to the low-mach momentum equation in the form

(3.554) $u_j - (u_{j,n}^D + u_{j,t}^D) = 0,$

where $u_{j,n}^D$ is the imposed normal component of velocity and $u_{j,t}^D$ is the imposed tangential component of velocity. The normal component is computed directly from the continuity flux at the interface,

(3.555) $u_{j,n}^D = \frac{F_B}{\rho} n_j.$

The tangential component is based on a variation of the classical Beavers-Joseph-Saffman condition [81, 82] for the slip velocity which has been extended to non-planar surfaces in multidimensional flow [83], which defines a provisional model velocity

(3.556) $u_j^{\rm BJS} = - \frac{\sqrt{\bar{K}}}{\alpha \mu} \left( n_i \tau_{ij} \right)$

where $\bar{K}$ is the permeability of the porous region at the interface, $\mu$ is the viscosity of the local fluid at the interface, $\tau_{ij}$ is the viscous stress tensor of the fluid at the interface, and $\alpha$ is a dimensionless model parameter that is a function of the microstructure of the porous material, which has been found to have typical values near $0.1$ [82]. The tangential component of this vector quantity is used as the tangential component of the distinguishing condition velocity, and is computed as

(3.557) $u_{j,t}^D = u_j^{\rm BJS} - \left( u_k^{\rm BJS} n_k \right) n_j.$

3.2.19.2. Enthalpy Transport

3.2.19.2.1. Bulk Equations

Porous Gas-Phase Enthalpy Equation

The gas-phase enthalpy equation within a porous region, to be solved for the gas-phase temperature $T_g$ , is

(3.558) $\frac{ \partial (\bar{\psi} \rho_g h_g) }{ \partial t } + \frac{ \partial (\rho_g u_{j,g} h_g) }{ \partial x_j } = - \frac{ \partial q_j^{h,g}}{ \partial x_j} + \left( \frac{\partial (\bar{\psi} p_g ) }{\partial t} + u_{j,g} \frac{\partial p_g}{\partial x_j} \right)$

(3.559) $+ h_{cv}\left( \bar{T} - T_g \right) + \sum_k \left( \dot{\omega}'''_{s,fk} - \dot{\omega}'''_{s,dk} \right) h_{k,g}$

where $h_g$ is the mixture-averaged gas-phase enthalpy, $h_{cv}$ is the volumetric heat transfer coefficient, $\bar{T}$ is the porous condensed-phase temperature, $\left( \dot{\omega}'''_{s,fk} - \dot{\omega}'''_{s,dk} \right)$ is the formation and destruction of gas-phase species due to heterogeneous reactions, and $h_{k,g}$ is the gas-phase enthalpy of chemical species $k$ . The gas-phase energy diffusive flux vector $q_j^{h,g}$ is modeled as

(3.560) $q_j^{h,g} = -\bar{\psi} \rho_g D_g \frac{\partial h_g}{\partial x_j},$

where $D_g$ is the mixture-averaged gas-phase mass diffusivity.

Note that, in (3.559), there is some concern that the pressure spatial derivative term, $u_{j,g} \frac{\partial p_g}{\partial x_j}$ , is incorrect. A crude re-derivation of this equation indicates that its form should instead be $\frac{u_{j,g}}{\bar{\psi}} \frac{\partial ( \bar{\psi} p_g ) }{\partial x_j}$ . A more formal re-derivation from first principles is required to decide conclusively on the correct form of this term, so it is left in its current form for now. Additionally, the diffusive flux vector is also of concern since the current form was derived under the assumption of constant specific heat, equal species mass diffusivities, and unity Lewis number. These assumptions may not be valid in future simulations, meaning that this term should possibly be returned to the standard Fick’s law version that includes a contribution due to enthalpy transport by differential diffusion of chemical species. Again, this term is left in its current form for the present work.

Multiplying (3.559) by an arbitrary test function $w$ and integrating over the domain $\Omega$ while integrating the advection and diffusion terms by parts, yields the variational form of the enthalpy equation that is solved for $h_g$ using the Galerkin finite element method,

(3.561) $\int\limits_\Omega w \left( \frac{ \partial (\bar{\psi} \rho_g h_g) }{ \partial t } - \left( \frac{\partial (\bar{\psi} p_g ) }{\partial t} + u_{j,g} \frac{\partial p_g}{\partial x_j} \right) - h_{cv}\left( \bar{T} - T_g \right) - \sum_k \left( \dot{\omega}'''_{s,fk} - \dot{\omega}'''_{s,dk} \right) h_{g,k} \right) d\Omega$

(3.562) $- \int\limits_\Omega \frac{\partial w}{\partial x_j} \left( \rho_g u_{j,g} h_g \right) d\Omega + \int\limits_\Gamma w \left( \rho_g u_{j,g} h_g \right) n_j d\Gamma$

(3.563) $- \int\limits_\Omega \frac{\partial w}{\partial x_j} q_j^{h,g} d\Omega + \int\limits_\Gamma w q_j^{h,g} n_j d\Gamma = 0.$

The boundary flux terms are then split into contributions on the interface between regions $A$ and $B$ and off the interface so that they may be treated separately. The enthalpy equation then takes the form

(3.564) $\int\limits_\Omega w \left( \frac{ \partial (\bar{\psi} \rho_g h_g) }{ \partial t } - \left( \frac{\partial (\bar{\psi} p_g ) }{\partial t} + u_{j,g} \frac{\partial p_g}{\partial x_j} \right) - h_{cv}\left( \bar{T} - T_g \right) - \sum_k \left( \dot{\omega}'''_{s,fk} - \dot{\omega}'''_{s,dk} \right) h_{g,k} \right) d\Omega$

(3.565) $- \int\limits_\Omega \frac{\partial w}{\partial x_j} \left( \rho_g u_{j,g} h_g \right) d\Omega - \int\limits_\Omega \frac{\partial w}{\partial x_j} q_j^{h,g} d\Omega$

(3.566) $+ \int\limits_{\Gamma \backslash \Gamma_{AB}} w \left( \rho_g u_{j,g} h_g \right) n_j d\Gamma + \int\limits_{\Gamma \backslash \Gamma_{AB}} w q_j^{h,g} n_j d\Gamma + \int\limits_{\Gamma_{AB}} w J_A^H d\Gamma = 0.$

where $J_A^H$ is the imposed flux on the porous side ( $A$ ) of the $\Gamma_{AB}$ interface. A detailed description of the coupling boundary flux is given in Coupling Boundary Conditions.

Low-Mach Enthalpy Equation

The enthalpy equation within the low-Mach fluid region, to be solved for the fluid temperature $T$ , is

(3.567) $\frac{ \partial (\rho h) }{ \partial t } + \frac{ \partial (\rho u_j h) }{ \partial x_j } = - \frac{ \partial q_j^h}{ \partial x_j} - \frac{ \partial q_j^r}{ \partial x_j} + \left( \frac{\partial p}{\partial t} + u_j \frac{\partial p}{\partial x_j} \right) + \tau_{ij} \frac{\partial u_i}{\partial x_j}$

where $h$ is the mixture-averaged fluid enthalpy, $q_j^r$ is a source term due to radiation absorption and emission, and $p$ is the fluid pressure. The diffusive flux vector is given by

(3.568) $q_j^h = -\lambda \frac{\partial T}{\partial x_j} + \sum_k \rho h_k Y_k \hat{u}_{jk},$

where $\lambda$ is the mixture thermal conductivity, $h_k$ is the enthalpy of species $k$ , $Y_k$ is the mass fraction of species $k$ , and $\hat{u}_{jk}$ is the diffusion velocity of species $k$ in the $j$ direction.

Integrating (3.567) over a CVFEM control volume and using the Gauss divergence theorem on the advective and diffusive flux terms, yields the integral form of the enthalpy equation to be solved,

(3.569) $\int\limits_\Omega \frac{ \partial (\rho h) }{ \partial t } d\Omega + \int\limits_\Omega \frac{ \partial q_j^r}{ \partial x_j} d\Omega - \int\limits_\Omega \left( \frac{\partial p}{\partial t} + u_j \frac{\partial p}{\partial x_j} \right) d\Omega - \int\limits_\Omega \tau_{ij} \frac{\partial u_i}{\partial x_j} d\Omega$

(3.570) $+ \int\limits_\Gamma \left( \rho u_j h \right) n_j d\Gamma + \int\limits_\Gamma q_j^h n_j d\Gamma = 0.$

The boundary flux terms are then split into contributions on the interface between regions $A$ and $B$ and off the interface so that they may be treated separately. The enthalpy equation then takes the form

(3.571) $\int\limits_\Omega \frac{ \partial (\rho h) }{ \partial t } d\Omega + \int\limits_\Omega \frac{ \partial q_j^r}{ \partial x_j} d\Omega - \int\limits_\Omega \left( \frac{\partial p}{\partial t} + u_j \frac{\partial p}{\partial x_j} \right) d\Omega - \int\limits_\Omega \tau_{ij} \frac{\partial u_i}{\partial x_j} d\Omega$

(3.572) $+ \int\limits_{\Gamma \backslash \Gamma_{AB}} \left( \rho u_j h \right) n_j d\Gamma + \int\limits_{\Gamma \backslash \Gamma_{AB}} q_j^h n_j d\Gamma + \int\limits_{\Gamma_{AB}} w J_B^h d\Gamma = 0,$

where $J_B^h$ is the imposed flux on the fluid side ( $B$ ) of the $\Gamma_{AB}$ interface. A detailed description of the coupling boundary condition is given in Coupling Boundary Conditions.

3.2.19.2.2. Coupling Boundary Conditions

Coupling enthalpy transport between the porous and fluid regions is complicated by the use of a two temperature model in the porous region.

To resolve this complication the energy flux applied to the fluid region has a diffusive/conductive component from the gas phase in the porous region, an advective component from the gas phase in the porous region, a convective component from the condensed phase in the porous region, and a penalty coefficient to enforce temperature continuity between the porous gas phase and the fluid. This takes the form

(3.573) $J_B^h = J_{A,g}^\mathrm{diff} + J_{A,g}^\mathrm{adv} + J_{A,c}^\mathrm{conv} + \overline{\left(\frac{\lambda}{h}\right)} (T_f - T_g),$

where $J_{A,g}^\mathrm{diff}$ is the diffusive energy transport from the porous gas phase, $J_{A,g}^\mathrm{adv}$ is the advective energy transport from the porous gas phase, $J_{A,c}^\mathrm{conv}$ is the convective energy transport from the porous condensed phase, and $\overline{\left(\frac{\lambda}{h}\right)}$ is the averaged thermal conductivity / mesh size between the porous and fluid regions. As with the flow coupling boundary conditions this same penalty coefficient is used in both regions to get the best energy conservation on coarse meshes.

The advective energy transport component takes the form

(3.574) $J_{A,g}^\mathrm{adv} = F_B h_{AB},$

where $h_{AB}$ is the upwinded interface enthalpy (i.e. it is either $h_A$ or $h_{B}$ depending on the direction of $F_B$ ). The convective component from the condensed phase has the form

(3.575) $J_{A,c}^\mathrm{conv} = (1-\bar{\psi}) \frac{h_{cv}}{a_s} (T_f - T_c),$

where $h_{cv}$ is the volumetric heat transfer coefficient of the porous region and $a_s$ is the specific surface area ( $m^2/m^3$ ) of the porous medium. This formulation of the convective component assumes that the convective heat transfer between the condensed phase and the free fluid is consistent with the convective heat transfer in the bulk of the porous medium that results in the volumetric heat transfer term of the bulk equations.

The coupling back to the porous region is derived based on the assumption that

(3.576) $J_{A,g}^h + J_{A,c}^h = J_B^\mathrm{diff} + J_B^\mathrm{adv},$

that is, the fluid region applies advective and diffusive energy transport components to the porous region as a whole. The flux applied to the condensed phase is assumed to be the same as the convective flux component it applies to the free fluid,

(3.577) $J_{A,c}^h = (1-\bar{\psi}) \frac{h_{cv}}{a_s} (T_c - T_f).$

The flux applied to the porous gas phase is then given by

(3.578) $J_{A,g}^h = J_B^\mathrm{diff} + J_B^\mathrm{adv} - J_{A,c}^h + \overline{\left(\frac{\lambda}{h}\right)} (T_g - T_f),$

where the advective component is computed in the same manner as is done for the advective flux applied to the fluid region,

(3.579) $J_{B}^\mathrm{adv} = F_A h_{AB}.$

3.2.19.3. Species Transport

3.2.19.3.1. Bulk Equations

Porous Gas-Phase Species Equation

The gas-phase species equation within a porous region, to be solved for the gas-phase mass fraction $Y_{k,g}$ of species $k$ , is

(3.580) $\frac{ \partial (\bar{\psi} \rho_g Y_{k,g}) }{ \partial t } + \frac{ \partial (\rho_g u_{j,g} Y_{k,g}) }{ \partial x_j } = - \frac{ \partial q_{kj}^{Y,g}}{ \partial x_j} + \left( \dot{\omega}'''_{s,fk} - \dot{\omega}'''_{s,dk} \right) + \left( \dot{\omega}'''_{g,fk} - \dot{\omega}'''_{g,dk} \right)$

where $\left( \dot{\omega}'''_{s,fk} - \dot{\omega}'''_{s,dk} \right)$ is the formation and destruction of gas-phase species due to heterogeneous reactions, and $\left( \dot{\omega}'''_{g,fk} - \dot{\omega}'''_{g,dk} \right)$ is the formation and destruction of gas-phase species due to homogeneous reactions. The gas-phase species diffusion flux vector $q_{kj}^{Y,g}$ is modeled as

(3.581) $q_{kj}^{Y,g} = -\bar{\psi} \rho_g D_{k,g} \frac{\partial Y_{k,g}}{\partial x_j},$

where $D_{k,g}$ is the gas-phase mass diffusivity for species $k$ . Note that if the mass diffusivities are not equal for all species, then an additional correction is required to maintain mass conservation.

Multiplying (3.580) by an arbitrary test function $w$ and integrating over the domain $\Omega$ while integrating the advection and diffusion terms by parts, yields the variational form of the species equation that is solved for $Y_{k,g}$ using the Galerkin finite element method,

(3.582) $\int\limits_\Omega w \left( \frac{ \partial (\bar{\psi} \rho_g Y_{k,g}) }{ \partial t } - \left( \dot{\omega}'''_{s,fk} - \dot{\omega}'''_{s,dk} \right) - \left( \dot{\omega}'''_{g,fk} - \dot{\omega}'''_{g,dk} \right) \right) d\Omega$

(3.583) $- \int\limits_\Omega \frac{\partial w}{\partial x_j} \left( \rho_g u_{j,g} Y_{k,g} \right) d\Omega + \int\limits_\Gamma w \left( \rho_g u_{j,g} Y_{k,g}\right) n_j d\Gamma$

(3.584) $- \int\limits_\Omega \frac{\partial w}{\partial x_j} q_{kj}^{Y,g} d\Omega + \int\limits_\Gamma w q_{kj}^{Y,g} n_j d\Gamma = 0.$

The boundary flux terms are then split into contributions on the interface between regions $A$ and $B$ and off the interface so that they may be treated separately. The species equation then takes the form

(3.585) $\int\limits_\Omega w \left( \frac{ \partial (\bar{\psi} \rho_g Y_{k,g}) }{ \partial t } - \left( \dot{\omega}'''_{s,fk} - \dot{\omega}'''_{s,dk} \right) - \left( \dot{\omega}'''_{g,fk} - \dot{\omega}'''_{g,dk} \right) \right) d\Omega$

(3.586) $- \int\limits_\Omega \frac{\partial w}{\partial x_j} \left( \rho_g u_{j,g} Y_{k,g} \right) d\Omega - \int\limits_\Omega \frac{\partial w}{\partial x_j} q_{kj}^{Y,g} d\Omega$

(3.587) $+ \int\limits_{\Gamma \backslash \Gamma_{AB}} w \left( \rho_g u_{j,g} Y_{k,g}\right) n_j d\Gamma + \int\limits_{\Gamma \backslash \Gamma_{AB}} w q_{kj}^{Y,g} n_j d\Gamma + \int\limits_{\Gamma_{AB}} w J_A^{Y_k} d\Gamma = 0.$

where $J_A^{Y_k}$ is the imposed flux on the porous side ( $A$ ) of the $\Gamma_{AB}$ interface. A detailed description of the coupling boundary flux is given in Coupling Boundary Conditions.

Low-Mach Species Equation

The species equation within the low-Mach fluid region, to be solved for the mass fraction $Y_k$ for species $k$ , is

(3.588) $\frac{ \partial (\rho Y_k) }{ \partial t } + \frac{ \partial (\rho u_j Y_k) }{ \partial x_j } = - \frac{ \partial q_{kj}^Y}{ \partial x_j} + \dot{\omega}'''_k$

where $\dot{\omega}'''_k$ is the volumetric mass formation rate if species $Y_k$ , and the diffusive flux vector is given by

(3.589) $q_{kj}^Y = -\rho \hat{u}_{j,k} Y_k,$

with $\hat{u}_{j,k}$ being the species diffusion velocity. Several forms for this velocity are possible, with the simplest being

(3.590) $\hat{u}_{j,k} = -D \frac{1}{Y_k} \frac{\partial Y_k}{\partial x_j}$

for equal mass diffusivities $D$ for all species. A more complex form is needed for unequal mass diffusivities, which is not presented here.

Integrating (3.588) over a CVFEM control volume and using the Gauss divergence theorem on the advective and diffusive flux terms yields the integral form of the species equation that is solved,

(3.591) $\int\limits_\Omega \frac{ \partial (\rho Y_k) }{ \partial t } d\Omega - \int\limits_\Omega \dot{\omega}'''_k d\Omega + \int\limits_\Gamma \left( \rho u_j Y_k \right) n_j d\Gamma + \int\limits_\Gamma q_{kj}^Y n_j d\Gamma = 0.$

The boundary flux terms are then split into contributions on the interface between regions $A$ and $B$ and off the interface so that they may be treated separately. The species equation then takes the form

(3.592) $\int\limits_\Omega \frac{ \partial (\rho Y_k) }{ \partial t } d\Omega - \int\limits_\Omega \dot{\omega}'''_k d\Omega + \int\limits_{\Gamma \backslash \Gamma_{AB}} \left( \rho u_j Y_k \right) n_j d\Gamma + \int\limits_{\Gamma \backslash \Gamma_{AB}} q_{kj}^Y n_j d\Gamma + \int\limits_{\Gamma_{AB}} w J_B^{Y_k} d\Gamma = 0.$

where $J_B^{Y_k}$ is the imposed flux on the fluid side ( $B$ ) of the $\Gamma_{AB}$ interface. A detailed description of the coupling boundary condition is given in Coupling Boundary Conditions.

3.2.19.3.2. Coupling Boundary Conditions

Coupling species transport across the porous-fluid interface is relatively simple compared to enthalpy transport. As with the flow coupling Robin style boundary conditions are applied on both the porous and fluid regions, but with both diffusive and advective flux components.

For the flux of a species $k$ this takes the form:

(3.593) $J_{A}^{Y_k} &= J_B^\mathrm{diff} + F_A \rho_g Y_{k,AB} + \overline{\left(\frac{D_k \rho}{h}\right)}(Y_{k,A} - Y_{k,B})$

(3.594) $J_{B}^{Y_k} &= J_A^\mathrm{diff} + F_B \rho Y_{k,AB} + \overline{\left(\frac{D_k \rho}{h}\right)}(Y_{k,B} - Y_{k,A})$

where $Y_{k,AB}$ is the upwinded interface mass fraction, equivalent to $h_{AB}$ from the enthalpy coupling. Once again the same penalty coefficient is used on each side in order to get good mass conservation even on coarse meshes.