18. Porous-Media Momentum Resistance

Porous resistance represents the momentum transferred from a fluid to an unresolved screen or solid matrix. In a pressure-based finite-volume method, this resistance enters at two related but distinct levels. It contributes a force to each cell momentum balance, and it also participates in the reconstruction of the provisional face flux from that momentum balance. A consistent treatment must address both levels.

The development follows the sequence by which the physical model enters the numerical method:

1. a continuum resistance law defines the momentum transfer and its dissipative character;

2. a finite-volume quadrature converts that force into cell and face quantities;

3. a nonlinear source split determines the matrix and right-hand side used by the momentum iteration;

4. a compatible face reconstruction carries the selected momentum response and source–pressure balance into continuity; and

5. the physical-time discretization supplies the required history and determines the damping of stiff transients.

The notation follows the finite-volume development in Chapter 4 and the time-integration development in Chapter 10. Thus \(P\) and \(N\) denote neighboring cells, \(f\) denotes their shared face, \(n\) labels a physical time level, \(k\) labels a nonlinear iteration at the new time level, and a star denotes a provisional predictor quantity. Porous-specific notation is introduced as it enters the derivation.

18.1. Continuum Model and Resistance Closures

With porous resistance included, the momentum equation may be written as

(18.1.1)\[\frac{\partial(\rho\vec{u})}{\partial t} +\nabla\cdot\left(\rho\vec{u}\otimes\vec{u}\right) =-\nabla p+\nabla\cdot\Matrix{\tau}+\vec{f}_b+\vec{f}_p,\]

where \(\rho\) is density, \(\vec{u}\) is velocity, \(p\) is pressure, \(\Matrix{\tau}\) is the viscous stress, \(\vec{f}_b\) contains other body forces, and \(\vec{f}_p\) is the porous force per unit volume. A local isotropic resistance law has the form

(18.1.2)\[\vec{f}_p=-C(\rho,\mu,q)\vec{u}, \qquad q=\lVert\vec{u}\rVert, \qquad C\geq 0.\]

Here \(\mu\) is dynamic viscosity and \(q\) is the speed. The coefficient \(C\) may contain both linear and nonlinear contributions. A representative Darcy–Forchheimer form is

(18.1.3)\[C(\rho,\mu,q)=\frac{\mu}{K}+\rho\beta_{\mathrm F} q,\]

where \(K\) is permeability and \(\beta_{\mathrm F}\) is a Forchheimer coefficient whose definition absorbs any conventional numerical factor. Other constitutive laws change the dependence of \(C\) on the state but not the source-linearization arguments below.

The sign of equation (18.1.2) makes the force dissipative:

(18.1.4)\[\vec{u}\cdot\vec{f}_p=-Cq^2\leq 0.\]

For steady one-dimensional flow through a layer of thickness \(L\), when pressure and porous resistance are the only streamwise forces,

(18.1.5)\[\frac{dp}{ds}=-C u_s, \qquad p_{\mathrm{in}}-p_{\mathrm{out}} =\int_0^L C u_s\,ds.\]

Here \(s\) is the streamwise coordinate and \(u_s\) is the velocity component in that direction. Larger resistance therefore produces a larger pressure loss at a given flow rate. For nonlinear resistance, increasing the flow rate can also increase \(C\) itself.

18.1.1. Examples of Isotropic Screen Closures

Two examples of scalar isotropic closures are the Armour–Cannon–Cady (ACC) correlation and a fitted quadratic correlation. ACC draws on woven-screen pressure-loss studies by Armour and Cannon and by Cady [AC68, Cad73] and has the coefficient

(18.1.1.1)\[C_{\mathrm{ACC}} = \frac{B}{B_{\mathrm{CFD}}}\frac{Q_{\mathrm{ACC}}}{\epsilon^2} \left( \alpha_{\mathrm{ACC}}\mu a_v^2 +\frac{\beta_{\mathrm{ACC}}\rho q}{D_p} \right).\]

Here \(\epsilon\) is void fraction, \(a_v\) is screen surface area per unit volume, \(D_p\) is a hydraulic pore diameter, \(Q_{\mathrm{ACC}}\) is a tortuosity factor, and \(\alpha_{\mathrm{ACC}}\), \(\beta_{\mathrm{ACC}}\), and \(B/B_{\mathrm{CFD}}\) are empirical correlation parameters. A purely quadratic fitted closure can be written as

(18.1.1.2)\[C_{\mathrm{fit}} =\frac{\rho K_{q,\mathrm{fit}}q}{2\epsilon^2}.\]

This fitted form is referred to here as numerically determined resistance (NDR), with \(K_{q,\mathrm{fit}}\) obtained from pressure-loss data. Both closures enter the momentum equation through \(\vec{f}_p=-C\vec{u}\) and inherit the dissipativity result only over parameter ranges for which \(C\geq0\).

The finite-volume development below assumes the positive scalar, isotropic resistance form in equation (18.1.2).

18.2. Spatial Finite-Volume Representation

Let \(P\) denote a control volume with volume \(\Omega_P\), centroid \(\vec{x}_P\), and integrated porous force

(18.2.1)\[\vec{S}_{p,P} \equiv \int_{\Omega_P}\vec{f}_p\,d\Omega.\]

Two useful discretizations differ in how this integral is constructed.

18.2.1. Standard Cell-Centered Quadrature

The standard approximation evaluates the resistance law at the cell center:

(18.2.1.1)\[\vec{S}_{p,P}^{\,\mathrm{std}} =\Omega_P\vec{f}_{p,P} =-\Lambda_P\vec{u}_P, \qquad \Lambda_P=C_P\Omega_P\geq 0.\]

Here \(\Lambda_P\) is the integrated scalar resistance coefficient. This is a local volume quadrature: the force direction is parallel and opposite to the cell-centered velocity.

18.2.2. Conservative Face-Distributed Quadrature

A face-based construction can tie the cell force more directly to the pressure loss represented between neighboring cells. Consider an internal face \(f\) between cells \(P\) and \(N\), with face centroid \(\vec{x}_f\) and area vector \(\vec{A}_f=A_f\hat{\vec{n}}_f\) directed from \(P\) to \(N\). Define the source-equivalent pressure increment

(18.2.2.1)\[\Delta p_{p,f} \approx \int_{\vec{x}_P}^{\vec{x}_N}\vec{f}_p\cdot d\vec{\ell} \approx \vec{f}_{p,P}\cdot(\vec{x}_f-\vec{x}_P) +\vec{f}_{p,N}\cdot(\vec{x}_N-\vec{x}_f).\]

The quantity \(\Delta p_{p,f}\) has units of pressure and approximates the pressure change associated with porous resistance along the segment joining the cell centers.

The force associated with this face-centered dual volume may be distributed between the adjacent cells according to a nonnegative geometric weight \(\eta_{P f}\):

(18.2.2.2)\[\vec{S}_{p,P\leftarrow f}^{\,c} =\eta_{P f}\Delta p_{p,f}\vec{A}_f, \qquad \vec{S}_{p,N\leftarrow f}^{\,c} =\left(1-\eta_{P f}\right)\Delta p_{p,f}\vec{A}_f, \qquad 0\leq\eta_{P f}\leq1.\]

For the line-segment quadrature in equation (18.2.2.1), a distance-based distribution is

(18.2.2.3)\[\eta_{P f} =\frac{\lVert\vec{x}_f-\vec{x}_P\rVert} {\lVert\vec{x}_f-\vec{x}_P\rVert+\lVert\vec{x}_N-\vec{x}_f\rVert}.\]

The cell source \(\vec{S}_{p,P}^{\,c}\) is the sum of its internal-face shares and the corresponding one-sided boundary-face contributions. The distribution is conservative because the two cell shares sum to the force associated with the shared dual volume. That force transfers momentum from the carrier fluid to the unresolved solid.

Because the conservative source is assembled from face-normal vectors, \(\vec{S}_{p,P}^{\,c}\) need not be exactly parallel to \(\vec{u}_P\) on a nonorthogonal mesh or where resistance changes between neighboring cells. This distinction matters when the source is linearized with a scalar momentum coefficient.

The standard and conservative constructions are alternative quadratures of the same integrated continuum force. The standard form is locally aligned with the cell velocity; the conservative form also supplies a face-path pressure increment and distributes its associated dual-volume force. Either spatial quadrature may be paired with explicit or implicit source placement.

18.3. Momentum-Source Linearization

Chapter 4 denotes the integrated pressure-gradient contribution by

(18.3.1)\[\vec{\mathcal{P}}_P \equiv \int_{\Omega_P}\nabla p\,d\Omega \approx \Omega_P(\nabla p)_P.\]

Let \({}^{(u)}\!a_{P,0}^{k}\) be the momentum diagonal before any porous contribution is made implicit, and let \(\vec{H}_P^k\) collect the nonpressure, nonporous terms on the right-hand side. One componentwise momentum row at nonlinear iteration \(k\) then has the vector shorthand

(18.3.2)\[{}^{(u)}\!a_{P,0}^{k}\vec{u}_P^{k+1} =\vec{H}_P^k+\vec{S}_{p,P} -\vec{\mathcal{P}}_P^k.\]

The coefficient \({}^{(u)}\!a_{P,0}^{k}\) contains the temporal, convective, and diffusive diagonal contributions. The vector \(\vec{H}_P^k\) contains neighbor terms, time-history terms, and other nonpressure contributions. Its superscript \(k\) indicates that its coefficients and explicit terms are assembled from the current nonlinear iterate; the implicit neighbor velocities inside \(\vec{H}_P^k\) belong to the same new predictor solve as \(\vec{u}_P^{k+1}\).

Here explicit and implicit describe source placement within the nonlinear momentum solve. Both forms evaluate the new-time porous force through the nonlinear iteration.

18.3.1. Standard Source: Explicit and Implicit Picard Forms

With the standard source of equation (18.2.1.1), an explicit source iteration is

(18.3.1.1)\[{}^{(u)}\!a_{P,0}^{k}\vec{u}_P^{k+1} =\vec{H}_P^k -\Lambda_P^k\vec{u}_P^k -\vec{\mathcal{P}}_P^k.\]

An implicit Picard treatment instead places the positive resistance coefficient on the left-hand side:

(18.3.1.2)\[\left({}^{(u)}\!a_{P,0}^{k}+\Lambda_P^k\right)\vec{u}_P^{k+1} =\vec{H}_P^k -\vec{\mathcal{P}}_P^k.\]

The coefficient \(C\), and hence \(\Lambda_P\), is evaluated from iteration \(k\) and held fixed while solving for \(\vec{u}^{k+1}\). This is a Picard linearization of the velocity-dependent resistance. After the nonlinear iteration converges, the iteration labels may be dropped and both rows reduce to

(18.3.1.3)\[{}^{(u)}\!a_{P,0}\vec{u}_P =\vec{H}_P -\Lambda_P\vec{u}_P -\vec{\mathcal{P}}_P.\]

The force and the implicit coefficient must come from the same constitutive law. The two placements then recover the same converged cell momentum residual through different iteration matrices.

18.3.2. Conservative Source and Deferred Remainder

For an explicitly treated conservative source,

(18.3.2.1)\[{}^{(u)}\!a_{P,0}^{k}\vec{u}_P^{k+1} =\vec{H}_P^k +\vec{S}_{p,P}^{\,c,k} -\vec{\mathcal{P}}_P^k.\]

Because \(\vec{S}_{p,P}^{\,c,k}\) need not be parallel to \(\vec{u}_P^k\), a scalar diagonal coefficient cannot represent the entire vector source. If \(\Lambda_P^k\) is placed on the diagonal, define the right-hand-side remainder

(18.3.2.2)\[\vec{R}_{p,P}^k =\vec{S}_{p,P}^{\,c,k} +\Lambda_P^k\vec{u}_P^k.\]

Using \(\vec{S}_{p,P}^{\,\mathrm{std},k}=-\Lambda_P^k\vec{u}_P^k\), the split can be written as

(18.3.2.3)\[\vec{R}_{p,P}^k =\vec{S}_{p,P}^{\,c,k}-\vec{S}_{p,P}^{\,\mathrm{std},k}, \qquad \vec{S}_{p,P}^{\,c,k} =\vec{S}_{p,P}^{\,\mathrm{std},k}+\vec{R}_{p,P}^k.\]

The corresponding implicit row is

(18.3.2.4)\[\left({}^{(u)}\!a_{P,0}^{k}+\Lambda_P^k\right)\vec{u}_P^{k+1} =\vec{H}_P^k +\vec{R}_{p,P}^k -\vec{\mathcal{P}}_P^k.\]

At the fixed point, the diagonal term and the remainder recover \(\vec{S}_{p,P}^{\,c}\) exactly. For the aligned case \(\vec{S}_{p,P}^{\,c}=-\Lambda_P\vec{u}_P\), the deferred remainder is zero.

Fig. 18.3.2.1 The scalar implicit coefficient represents the standard source \(\vec{S}_{p,P}^{\,\mathrm{std}}=-\Lambda_P\vec{u}_P\). The deferred remainder \(\vec{R}_{p,P}=\vec{S}_{p,P}^{\,c}-\vec{S}_{p,P}^{\,\mathrm{std}}\) closes the vector sum to the conservative source. Vector lengths are illustrative.

18.4. Face-Flux Reconstruction

In a collocated pressure-based method, the momentum diagonal determines the velocity response used to reconstruct the provisional face flux. Following the established face-flux notation, the volumetric flux and cell-to-cell pressure difference at an internal face are

(18.4.1)\[\dot V_f=\vec{u}_f\cdot\vec{A}_f, \qquad \Delta p_f=p_N-p_P.\]

The momentum diagonal and porous term on the right-hand side for each source placement are

(18.4.2)\[\begin{split}\begin{array}{c|c|c} \text{placement} &{}^{(u)}\!a_P^{k} &\vec{S}_{p,P}^{\,\mathrm{rhs},k}\\ \hline \text{explicit} &{}^{(u)}\!a_{P,0}^{k} &\vec{S}_{p,P}^k\\ \text{implicit} &{}^{(u)}\!a_{P,0}^{k}+\Lambda_P^k &\vec{S}_{p,P}^k+\Lambda_P^k\vec{u}_P^k . \end{array}\end{split}\]

This definition applies to either spatial quadrature. For the standard source, the implicit right-hand-side term is zero. For the conservative source, it is the deferred remainder \(\vec{R}_{p,P}\) from equation (18.3.2.2). The face reconstruction uses \({}^{(u)}\!a_{P,0}\) for explicit placement and \({}^{(u)}\!a_{P,0}+\Lambda_P\) for implicit placement. Its source–pressure balance uses the physical porous source; the deferred remainder belongs to the implicit cell row.

A local pressure–porous balance satisfies

(18.4.3)\[\vec{S}_{p,j}=\vec{\mathcal{P}}_j, \qquad j\in\{P,N\}, \qquad \Delta p_{p,f}=\Delta p_f.\]

At this local equilibrium, the paired pressure and porous contributions cancel within each source placement.

18.5. Time Integration and Resistance Stiffness

Source placement defines the nonlinear momentum row at a new physical time, while the time integrator defines the temporal residual and its history terms. Following Chapter 10, define the cell-integrated momentum

(18.5.1)\[\vec{Q}_P=(\rho\vec{u})_P\Omega_P=M_P\vec{u}_P, \qquad M_P=\rho_P\Omega_P,\]

and write the semi-discrete momentum balance as

(18.5.2)\[\frac{d\vec{Q}_P}{dt}+\vec{\mathcal{R}}_P=0, \qquad \vec{\mathcal{R}}_P =\vec{\mathcal{R}}_{P,0}-\vec{S}_{p,P}.\]

Here \(\vec{\mathcal{R}}_{P,0}\) is the nonporous spatial residual, and \(\vec{S}_{p,P}\) is the porous force on the right-hand side of the momentum balance.

18.5.1. BDF2 History and Face Response

Using the \(D_2\) notation of Chapter 10, BDF2 advances the conservative momentum according to

(18.5.1.1)\[D_2\vec{Q}_P^{\,n+1} +\vec{\mathcal{R}}_{P,0}^{\,n+1} -\vec{S}_{p,P}^{\,n+1} =\vec{0}.\]

For a constant timestep this is

(18.5.1.2)\[\frac{3\vec{Q}_P^{\,n+1}-4\vec{Q}_P^{\,n}+\vec{Q}_P^{\,n-1}} {2\Delta t} +\vec{\mathcal{R}}_{P,0}^{\,n+1} -\vec{S}_{p,P}^{\,n+1} =\vec{0}.\]

The porous force is evaluated at the new physical time; \(n\) and \(n-1\) enter through the momentum history. If density or control-volume size changes, each historical momentum \(\vec{Q}_P^{\,j}=M_P^{\,j}\vec{u}_P^{\,j}\) carries the mass from its own time level. For variable timesteps, \(D_2\) uses the coefficients defined in equation (10.3.2.6).

18.5.1.1. Time-Consistent Face Response

Chapter 4 separates the momentum diagonal into temporal and spatial pieces. With porous resistance, the complete predictor diagonal has the form

(18.5.1.3)\[\begin{split}{}^{(u)}\!a_P ={}^{(u)}\!a_{P,\mathrm{time}} +{}^{(u)}\!a_{P,\mathrm{spatial}} +\begin{cases} 0, & \text{explicit porous source},\\ \Lambda_P, & \text{implicit porous source}. \end{cases}\end{split}\]

For fixed \(M_P\) and constant \(\Delta t\), a directly assembled BDF2 derivative contributes \({}^{(u)}\!a_{P,\mathrm{time}}=3M_P/(2\Delta t)\). A factored BDF2 form distributes this temporal response between the predictor diagonal and a lagged correction. The predictor diagonal, current-iterate feedback, face response, and \(n\) and \(n-1\) history terms use the same selected factorization.

18.5.2. Crank–Nicolson Residual Blending

Let \(0\leq\beta\leq1\) blend backward Euler at \(\beta=0\) with Crank–Nicolson at \(\beta=1\). Using the residual convention of equation (18.5.2), the update is

(18.5.2.1)\[\frac{\vec{Q}_P^{\,n+1}-\vec{Q}_P^{\,n}}{\Delta t} +\left(1-\frac{\beta}{2}\right)\vec{\mathcal{R}}_P^{\,n+1} +\frac{\beta}{2}\vec{\mathcal{R}}_P^{\,n} =\vec{0}.\]

The residual at each time level contains the complete spatial balance. The old-time porous force therefore enters equation (18.5.2.1) with weight \(\beta/2\).

18.5.3. Resistance Stiffness Measures

The physical-time drag stiffness for one cell is

(18.5.3.1)\[z_P=\frac{\Lambda_P\Delta t}{M_P} =\frac{C_P\Delta t}{\rho_P}.\]

The ratio \(z_P\) compares the timestep with the local drag-decay time. A distinct iteration-level comparison is \(\Lambda_P/{}^{(u)}\!a_{P,0}\), the resistance relative to the nonporous momentum diagonal. If the backward-Euler temporal term dominates, then \({}^{(u)}\!a_{P,0}\approx M_P/\Delta t\) and the two ratios are approximately equal. Other time integrators change the current-time coefficient in the diagonal.

The ratio \(C\Delta t/\rho\) measures local drag stiffness. The coupled response also depends on convection, diffusion, pressure coupling, mesh geometry, nonlinear coefficient variation, and relaxation.

18.5.4. Nonlinear-Iteration Convergence

Hold \(\Lambda_P\) fixed and consider a local solve with positive nonporous diagonal \({}^{(u)}\!a_{P,0}>0\) and fixed right-hand side \(\vec{b}_P\). After applying velocity under-relaxation, an explicitly lagged source gives the recurrence

(18.5.4.1)\[\vec{u}_P^{k+1} =(1-\alpha_v)\vec{u}_P^k +\frac{\alpha_v}{{}^{(u)}\!a_{P,0}} \left(\vec{b}_P-\Lambda_P\vec{u}_P^k\right).\]

A perturbation \(\vec{e}_P^k\) from the fixed point therefore satisfies

(18.5.4.2)\[\vec{e}_P^{k+1} = \left[ 1-\alpha_v \left( 1+\frac{\Lambda_P}{{}^{(u)}\!a_{P,0}} \right) \right] \vec{e}_P^k.\]

Without under-relaxation, the factor is \(-\Lambda_P/{}^{(u)}\!a_{P,0}\), so the isolated iteration becomes increasingly oscillatory with resistance and ceases to be contractive when \(\Lambda_P/{}^{(u)}\!a_{P,0}\geq1\). More generally, it is contractive when

(18.5.4.3)\[0<\alpha_v< \frac{2}{1+\Lambda_P/{}^{(u)}\!a_{P,0}}.\]

When the temporal diagonal dominates, maintaining contraction with an explicit source therefore requires progressively smaller \(\Delta t\) or stronger under-relaxation.

Treating the constant resistance implicitly gives the recurrence

(18.5.4.4)\[\vec{u}_P^{k+1} =(1-\alpha_v)\vec{u}_P^k +\frac{\alpha_v}{{}^{(u)}\!a_{P,0}+\Lambda_P}\vec{b}_P.\]

The lagged linear feedback is removed, and the corresponding source-only error factor is \(1-\alpha_v\). Nonlinear variation of \(C\) can still require outer iterations because a Picard treatment lags the coefficient.

18.5.5. Physical-Time Damping

Consider the source-only linear model

(18.5.5.1)\[M_P\frac{d\vec{u}_P}{dt} =-\Lambda_P\vec{u}_P.\]

For one velocity component in cell \(P\), use the physical stiffness \(z_P=\Lambda_P\Delta t/M_P\). With resistance evaluated implicitly at the new physical time, equation (18.5.2.1) gives the amplification factor

(18.5.5.2)\[G(z_P,\beta) =\frac{1-\frac{\beta}{2}z_P} {1+\left(1-\frac{\beta}{2}\right)z_P}.\]

Pure Crank–Nicolson gives

(18.5.5.3)\[G(z_P,1)=\frac{1-z_P/2}{1+z_P/2}.\]

This factor remains bounded for every \(z_P\geq0\), but becomes negative for \(z_P>2\) and approaches \(-1\) as \(z_P\rightarrow\infty\). Pure Crank–Nicolson is therefore stable for this linear drag model but does not strongly damp a very stiff mode; the numerical solution can alternate in sign from one timestep to the next. Reducing \(\beta\) adds backward-Euler damping, and \(\beta=0\) gives \(G=1/(1+z_P)\).

This physical-time behavior is separate from the nonlinear contraction analyzed above. Implicit resistance can improve the nonlinear iteration, while the selected time integrator still determines damping, monotonicity, and temporal accuracy.

18.6. References

[AC68]

James C. Armour and Joseph N. Cannon. Fluid flow through woven screens. AIChE Journal, 14(3):415–420, 1968. doi:10.1002/aic.690140315.

[Cad73]

E. C. Cady. Study of thermodynamic vent and screen baffle integration for orbital storage and transfer of liquid hydrogen. Technical Report NASA-CR-134482, McDonnell Douglas Astronautics Company, 1973. URL: https://ntrs.nasa.gov/citations/19730025004.