# 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 and the low-Mach free fluid region is described as region , with the interface between them referred to as and other boundaries not a part of this interface are referred to as .

## 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 , is

(3.534)

where is the mixture-averaged condensed-phase porosity, is the gas-phase density, and is the gas-phase velocity vector computed from Darcy’s approximation as

(3.535)

where is the mixture-averaged condensed-phase permeability, is the gas-phase viscosity, and is the gravity vector. The term represents the formation rate of gas-phase mass from the condensed phase.

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

(3.536)

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

(3.537)

(3.538)

where is the imposed flux on the porous side () of the 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 , is

(3.539)

where is the fluid density, is the fluid velocity, and 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)

Similar to the porous continuity equation, the boundary flux term is split into contributions both on and off the interface, yielding

(3.541)

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 , is

(3.542)

where the Cauchy stress tensor is given by

(3.543)

in terms of the viscous stress tensor

(3.544)

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)

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

(3.546)

(3.547)

### 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 [81]. 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 are

(3.548)

(3.549)

where , , and the free constant is computed as

(3.550)

with being a measure of the mesh size adjacent to the interface, and a user-specified scaling coefficient. The same value of is used on both sides of the interface because that results in excellent mass conservation even on coarse meshes. If a different value of 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)

(3.552)

(3.553)

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)

where is the imposed normal component of velocity and is the imposed tangential component of velocity. The normal component is computed directly from the continuity flux at the interface,

(3.555)

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

(3.556)

where is the permeability of the porous region at the interface, is the viscosity of the local fluid at the interface, is the viscous stress tensor of the fluid at the interface, and 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 [83]. The tangential component of this vector quantity is used as the tangential component of the distinguishing condition velocity, and is computed as

(3.557)

## 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 , is

(3.558)

(3.559)

where is the mixture-averaged gas-phase enthalpy, is the volumetric heat transfer coefficient, is the porous condensed-phase temperature, is the formation and destruction of gas-phase species due to heterogeneous reactions, and is the gas-phase enthalpy of chemical species . The gas-phase energy diffusive flux vector is modeled as

(3.560)

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

Note that, in (3.559), there is some concern that the pressure spatial derivative term, , is incorrect. A crude re-derivation of this equation indicates that its form should instead be . 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 and integrating over the domain while integrating the advection and diffusion terms by parts, yields the variational form of the enthalpy equation that is solved for using the Galerkin finite element method,

(3.561)

(3.562)

(3.563)

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

(3.564)

(3.565)

(3.566)

where is the imposed flux on the porous side () of the 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 , is

(3.567)

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

(3.568)

where is the mixture thermal conductivity, is the enthalpy of species , is the mass fraction of species , and is the diffusion velocity of species in the 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)

(3.570)

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

(3.571)

(3.572)

where is the imposed flux on the fluid side () of the 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)

where is the diffusive energy transport from the porous gas phase, is the advective energy transport from the porous gas phase, is the convective energy transport from the porous condensed phase, and 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)

where is the upwinded interface enthalpy (i.e. it is either or depending on the direction of ). The convective component from the condensed phase has the form

(3.575)

where is the volumetric heat transfer coefficient of the porous region and is the specific surface area () 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)

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)

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

(3.578)

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

(3.579)

## 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 of species , is

(3.580)

where is the formation and destruction of gas-phase species due to heterogeneous reactions, and is the formation and destruction of gas-phase species due to homogeneous reactions. The gas-phase species diffusion flux vector is modeled as

(3.581)

where is the gas-phase mass diffusivity for species . 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 and integrating over the domain while integrating the advection and diffusion terms by parts, yields the variational form of the species equation that is solved for using the Galerkin finite element method,

(3.582)

(3.583)

(3.584)

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

(3.585)

(3.586)

(3.587)

where is the imposed flux on the porous side () of the 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 for species , is

(3.588)

where is the volumetric mass formation rate if species , and the diffusive flux vector is given by

(3.589)

with being the species diffusion velocity. Several forms for this velocity are possible, with the simplest being

(3.590)

for equal mass diffusivities 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)

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

(3.592)

where is the imposed flux on the fluid side () of the 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 this takes the form:

(3.593)

(3.594)

where is the upwinded interface mass fraction, equivalent to 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.