Porous Media Module

The porous module provides momentum resistance and local thermal equilibrium (LTE) models. porousScreens adds resistance without changing fluid storage. porousLTE includes fluid volume fraction, solid heat storage, and effective heat conduction. Neither model includes an interface slip law.

Momentum Resistance Setup

Load the module, define a material, and assign it to a region. This example applies isotropic resistance to cells with the volume tag porous_region in the VOG file:

loadModule: porous
{
  // ... standard flow and boundary settings
  porousScreens: <
    material=DarcyForchheimer(darcy=400,forchheimer=40),
    screens=[volumeTag(name=porous_region,material=material)]
  >
}

Material names are arbitrary. Each screens entry selects a material by name and a region by volume tag or geometry. With porousScreens, velocity is bulk (superficial) velocity, measured over the total flow area.

Inviscid cases (transport_model: none) require a nonzero initial velocity.

Source Treatment and Assembly

porositySourceOptions selects the source treatment and spatial assembly. The defaults are implicit and conservative:

porositySourceOptions: <implicit,conservative>

Source options

Option	Effect
implicit	Add resistance to the momentum diagonal.
explicit	Evaluate resistance using the current iterate and place it on the momentum right-hand side.
conservative	Assemble face contributions into adjacent cells.
standard	Evaluate the force at each cell center.

Choose implicit or explicit and conservative or standard; omitted choices retain their defaults.

Darcy–Forchheimer Model

DarcyForchheimer combines linear viscous drag and quadratic inertial drag. For scalar coefficients, the force per unit volume is

\[\vec S_p=-\left(\mu d+\tfrac12\rho f\lVert\vec u\rVert\right)\vec u\]

Here \(\rho\) is density, \(\mu\) is dynamic viscosity, and \(\vec u\) is velocity. The inputs darcy and forchheimer specify \(d\) and \(f\). Supply at least one resistance option; all supplied coefficients must be finite.

Resistance coefficients

Option	SI unit	Default	Meaning
darcy	1/m²	0	Inverse permeability \(d\geq0\).
permeability	m²	Unset	Permeability \(K>0\); sets \(d=1/K\). Mutually exclusive with darcy.
forchheimer	1/m	0	Inertial resistance \(f\geq0\).

Bare values use SI units; permeability=25 cm*cm equals darcy=400. For the convention \(\mu/K+\rho\beta_F\lVert\vec u\rVert\), use darcy=1/K and forchheimer=2*beta_F. The NDR model is recovered with darcy=0 and forchheimer=Kq/epsilon^2.

Directional Resistance

For directional resistance, supply three values for each coefficient, one per material axis. A scalar applies equally along all three axes. The force is

\[\vec S_p=-\left(\mu\mathbf D+\tfrac12\rho\lVert\vec u\rVert\mathbf F\right)\vec u\]

The entries of darcy and forchheimer are the principal values of \(\mathbf D\) and \(\mathbf F\). Both tensors share the material axes. Supply axis1 and axis2 together as finite, nonzero, perpendicular vectors. They are normalized; their cross product defines the third axis. The defaults are Cartesian x and y. Material axes are independent of the screen geometry’s axes.

Unequal principal coefficients require porousSGS and implicit treatment for stationary Cartesian incompressible flow. This example sets the first material axis at 30 degrees to x and includes the required solver settings:

porousScreens: <
  material=DarcyForchheimer(darcy=[500,50000,500],forchheimer=[2,10,2],
    axis1=[0.8660254,0.5,0],axis2=[-0.5,0.8660254,0]),
  screens=[volumeTag(name=porous_cylinder,material=material)]
>
flowCompressibility: incompressible
gridCoordinates: cartesian
porositySourceOptions: <implicit,conservative>
momentumEquationOptions: <linearSolver=porousSGS>
momentumInterpolationOptions: <form=old>
pressureBasedMethod: SIMPLE
timeIntegrator: BDF
inviscidFlux: SOU

FOU may replace SOU, and standard may replace conservative. The conservativePressure interpolation modifier is optional. Adjust relaxationFactor and maxIterations in momentumEquationOptions as needed for convergence.

Translational periodic boundaries are supported. Rotating frames, rotational periodicity, moving meshes, and interface boundaries are unsupported. With twoDimensions, resistance must not couple in-plane and out-of-plane velocity components. These restrictions also apply when porousSGS is selected for an isotropic material.

Fitting Pressure-Loss Measurements

tools/fit_porous_resistance.py fits nonnegative darcy and forchheimer coefficients. It requires Python, NumPy, and a CSV file with columns velocity_m_s and pressure_drop_pa. Include at least two distinct nonzero speed magnitudes. From the Stream source directory, run:

python3 tools/fit_porous_resistance.py measurements.csv \
  --length 0.2 --density 1.0 --viscosity 0.01 --output fit.json

--length is the porous-layer thickness \(L\); --density and --viscosity are fluid density and dynamic viscosity. Use SI units, superficial velocity, and \(\Delta p=p(0)-p(L)\), with positive velocity from 0 to \(L\). Exclude entrance, exit, and pipe losses. Fit each principal direction separately and specify material axes independently.

The command prints a material definition and writes fitted coefficients, measured and predicted pressure losses, and pressure-error estimates to fit.json.

Armour Cannon Cady (ACC) Mesh Model

The Armour–Cannon–Cady model [ArCa1968] [Cady1973] specifies screen resistance from empirical parameters. Its force per unit volume is

\[\vec S_p=-\frac{B}{B_{\mathrm{CFD}}}\frac{Q}{\epsilon^2} \left(\alpha\mu a^2+\frac{\beta\rho\lVert\vec u\rVert}{D_p}\right)\vec u\]

All parameters are required:

ACC parameters

Option	SI unit	Meaning
a	1/m	Screen surface area per unit volume.
B	m	Physical screen thickness.
BCFD	m	Thickness of the modeled porous region.
Dp	m	Hydraulic pore diameter.
Q	None	Mesh tortuosity factor.
alpha	None	Viscous resistance constant.
beta	None	Inertial resistance constant.
epsilon	None	Screen void fraction.

Use ACC(...) as a material definition in porousScreens:

mesh1=ACC(a=510235, B=8.9e-5, BCFD=1000.0,
  Dp=5e-6, Q=1.3, alpha=3.2, beta=0.19, epsilon=0.245)

Numerically Determined Resistance (NDR) Mesh Model

NDR specifies quadratic resistance:

\[\vec S_p=-\frac{\rho K_q\lVert\vec u\rVert}{2\epsilon^2}\vec u\]

Both parameters are required:

NDR parameters

Option	SI unit	Meaning
Kq	1/m	Quadratic resistance coefficient \(K_q\).
epsilon	None	Screen void fraction.

Enter Kq as a bare value in inverse meters:

mesh2=NDR(Kq=1742.782,epsilon=0.4)

Screen Geometries

Use circle, cylinder, or box in the screens list. Each entry selects a named material with material=....

Circular Screen

circle defines a planar circular screen.

Parameter	Unit	Description
radius	m	Radius of the circular screen.
center	m	Center point of the circular screen.
normal		Normal vector to the plane of the circular screen.

circle(radius=1.0, center=[0,0,0], normal=[0,1,0], material=material)

Cylindrical Screen

cylinder defines a cylindrical porous region.

Parameter	Unit	Description
radius	m	Radius of the cylindrical screen.
length	m	Length of the cylindrical screen.
center	m	Center point of the cylindrical screen.
axis		Axis vector along the length of the cylinder.

cylinder(radius=1.0, length=1.0, center=[0,0,0],
  axis=[0,1,0], material=material)

Note

For a cylinder aligned with a 2D mesh’s extrusion direction, extend its ends slightly beyond both extrusion planes. Coincident ends can prevent cells from being tagged. For a mesh spanning z = -0.5 m to 0.5 m, use center=[0,0,0], axis=[0,0,1], and length=1.1.

Box Screen

box defines a porous region by its center, orientation, and widths.

Parameter	Unit	Description
w1	m	Width of the box along axis1.
w2	m	Width of the box along axis2.
w3	m	Width of the box along axis3 (automatically determined by the cross product of axis1 and axis2).
center	m	Center point of the box screen.
axis1		First axis vector defining one edge of the box.
axis2		Second axis vector defining another edge of the box perpendicular to axis1.

box(center=[0,0,0], axis1=[1,0,0], axis2=[0,1,0],
  w1=1.0, w2=1.0, w3=1.0, material=material)

Local Thermal Equilibrium

LTE uses one temperature for the fluid and stationary solid. Select one isotropic material for the whole mesh or a named volume tag:

loadModule: porous
{
  porousLTE: <volumeTag=porous_region,porosity=.5,
    solidDensity=2500 kg/m/m/m,solidSpecificHeat=800 J/kg/K,
    effectiveConductivity=5 W/m/K,darcy=1e8,forchheimer=1000,
    heatSource=1e5 W/m/m/m>
  flowRegime: laminar
  flowCompressibility: compressible
  pressureBasedMethod: SIMPLE
  timeIntegrator: BDF
  inviscidFlux: SOU
  momentumInterpolationOptions: <form=old,conservativePressure>
  energyEquationOptions: <form=totalEnthalpy,linearSolver=SGS,
    relaxationFactor=.7,maxIterations=20>
  // ... EOS, grid, initial conditions, boundaries, and time settings
}

Cells outside volumeTag retain ordinary fluid properties. Omit volumeTag to fill the domain. A tag that selects no cells is an error.

LTE material options

Option	SI unit	Default	Meaning
porosity	—	Required	Fluid volume fraction, greater than zero and at most one.
solidDensity	kg/m³	Required	Positive density of the solid material.
solidSpecificHeat	J/(kg K)	Required	Positive, constant solid heat capacity.
effectiveConductivity	W/(m K)	Required	Nonnegative conductivity of the combined material, per total area.
darcy	1/m²	0	Nonnegative inverse permeability.
forchheimer	1/m	0	Nonnegative inertial resistance, using the one-half convention above.
heatSource	W/m³	0	Heat supplied per total volume; negative values remove heat.
volumeTag	—	Whole mesh	VOG volume tag occupied by the material.

Density remains the fluid EOS density. Velocity inputs and outputs are intrinsic fluid velocity; superficial velocity is porosity times intrinsic velocity. Resistance coefficients, specified mass or volume fluxes, and wall heat fluxes use total area. Do not multiply effectiveConductivity or heatSource by porosity.

LTE requires a single-species chemEOS fluid, total enthalpy, SIMPLE, BDF, FOU or SOU, old momentum interpolation, and a stationary Cartesian mesh. Use the ordinary compressible mass form of pressure correction and drhodt_predictor: zero. Species transport, turbulence, gravity, rotating frames, periodic or overset interfaces, total-pressure inlets, fixed-mass outlets, double flux, and PIMPLE are unsupported. Do not combine porousLTE with porousScreens. Keep material properties unchanged across a restart.

For conduction-dominated cases, use full energy relaxation and enough scalar-solver iterations to converge each time step.

Output Variables

For porousScreens, add screenID to plot_output to inspect regions. LTE provides porousLTEPorosity, porousLTECapacity (solid heat capacity per total volume), and porousLTEConductivity.

output/porous_lte_budget.dat is written at print_freq. Its columns are cycle, time, fluid mass, combined fluid/solid total energy, imposed bulk heat input, outward mass flow, and outward total-energy flow. Energy flow includes enthalpy transport, heat conduction, and viscous work. Standard volume integrals also include porosity and solid energy when LTE is active.

References

[ArCa1968]

10. 3. Armour, J. N. Cannon, Fluid flow through woven screens 1968.

[Cady1973]

5. 3. Cady, Study of Thermodynamic Vent and Screen Baffle Integration for Orbital Storage and Transfer of Liquid Hydrogen, 1973.