Fluid–Poroelastic Structure Interaction (FPSI)

• Fluid–poroelastic structure interaction (FPSI) is a coupled multiphysics problem that integrates free fluid dynamics with deformation and Darcy filtration in porous media.
• FPSI formulations utilize Stokes/Navier–Stokes equations in fluid regions and Biot equations in poroelastic media, enforcing mass conservation and stress balance via tailored interface laws.
• Applications of FPSI span fractured rock and multilayer biological tissues, with numerical strategies including monolithic, partitioned, and Schur-complement methods ensuring stable simulations.

Fluid–poroelastic structure interaction (FPSI) denotes a class of coupled multiphysics problems in which a free fluid or a fracture-confined fluid interacts with a deformable, fluid-saturated porous solid. In the canonical interface setting, a free-fluid region $\Omega_f$ is coupled to a poroelastic region $\Omega_p$ across an interface $\Gamma_{fp}$; in fractured-media variants, a lower-dimensional fracture $\Sigma$ exchanges mass and momentum with a poroelastic matrix, with the fracture aperture determined by deformation (Ambartsumyan et al., 2017, Hanowski et al., 2016). FPSI therefore combines elements of fluid–structure interaction, Stokes–Darcy transmission, Darcy filtration, and Biot poromechanics. Representative models use Stokes or Navier–Stokes equations for the free fluid, Biot equations for the porous skeleton and pore pressure, and interface laws expressing mass conservation, stress balance, and Beavers–Joseph–Saffman slip, or, in fractured materials, aperture-dependent Darcy flow and pressure-induced crack opening (Bukac et al., 2014).

1. Model classes and physical scope

Representative FPSI formulations fall into several recurrent classes.

Variant	Governing subsystems	Representative papers
Stokes–Biot interface FPSI	Stokes in $\Omega_f$, quasi-static or dynamic Biot in $\Omega_p$, interface transmission on $\Gamma_{fp}$	(Ambartsumyan et al., 2017, Li et al., 2020, Dalal et al., 2024)
ALE/Navier–Stokes–Biot FPSI	ALE Navier–Stokes in a moving fluid domain and Biot in a deforming porous domain	(Guo et al., 2024)
Fractured poroelastic media	Bulk poroelasticity plus lower-dimensional fracture Darcy flow with aperture coupling	(Hanowski et al., 2016)
Thin or multilayer poroelastic structures	Stokes flow coupled to poroelastic plates, membranes, or multilayered walls	(Scharf et al., 14 Jul 2025, Brandt et al., 17 May 2026, Bukac et al., 2013)

In the standard Stokes–Biot formulation, the free fluid satisfies

$$-\nabla\cdot \sigma_f(u_f,p_f)=f_f,\qquad \nabla\cdot u_f=q_f \quad \text{in }\Omega_f,$$

with

$$\sigma_f(u_f,p_f)=-p_f I + 2\mu D(u_f),\qquad D(u_f)=\frac12(\nabla u_f+\nabla u_f^T),$$

while the poroelastic medium satisfies

$$-\nabla\cdot \sigma_p(\eta_p,p_p)=f_p,\qquad \mu K^{-1}u_p+\nabla p_p=0,$$

$\Omega_p$0

with

$\Omega_p$1

(Ambartsumyan et al., 2017, Houédanou, 2020).

Several extensions alter this baseline in essential ways. Fully inertial formulations retain solid and fluid inertia and analyze a hyperbolic–parabolic Biot–Stokes filtration system in 3D (Avalos et al., 2024). Moving-domain formulations place the fluid in ALE coordinates and couple it to a Biot solid on a reference configuration (Guo et al., 2024). Fracture-based FPSI treats the fracture as a $\Omega_p$2-dimensional manifold $\Omega_p$3, with tangential Darcy flow and an aperture $\Omega_p$4 that is itself an unknown determined by deformation (Hanowski et al., 2016). Large-deformation variational models also exist, in which a poroelastic body deforms inside an incompressible Navier–Stokes fluid and the coupling is mediated by drag and shared energetic structure rather than only by small-strain transmission laws (Benesova et al., 2021).

A common misconception is that FPSI is merely a minor variant of classical FSI. The Stokes–Biot literature instead treats it as a prototype problem that “shares common traits with fluid-structure interaction” and also “resembles the Stokes-Darcy coupling,” with the porous subdomain carrying its own flow, storage, and stress response (Bukac et al., 2014).

2. Interface laws and coupling mechanisms

The defining feature of FPSI is the interface law. In the Stokes–Biot setting, three conditions recur throughout the literature. The first is mass conservation across $\Omega_p$5: $\Omega_p$6 The second is normal and total stress balance: $\Omega_p$7 The third is the Beavers–Joseph–Saffman slip-with-friction condition: $\Omega_p$8 (Houédanou, 2020, Li et al., 2020, Ambartsumyan et al., 2017).

These relations encode different physical channels. Mass conservation couples free-fluid motion to solid velocity and Darcy seepage. The stress conditions transmit fluid traction to the poroelastic skeleton and identify pore pressure as the normal interfacial load. The BJS term models tangential slip with friction, interpolating between no-slip-type behavior and permeable-interface slip.

In fractured poroelastic media, the interface law becomes aperture-dependent and lower-dimensional. The fracture pressure $\Omega_p$9 and tangential flux $\Gamma_{fp}$0 satisfy

$\Gamma_{fp}$1

$\Gamma_{fp}$2

with coupling conditions

$\Gamma_{fp}$3

and mechanical loading

$\Gamma_{fp}$4

The aperture is not prescribed but determined by

$\Gamma_{fp}$5

(Hanowski et al., 2016).

This aperture law is the central nonlinear feedback in fracture FPSI: deformation determines transmissibility, transmissibility determines fracture pressure, and fracture pressure loads the matrix. Near a crack tip, the aperture behaves asymptotically like

$\Gamma_{fp}$6

so $\Gamma_{fp}$7 cannot be bounded away from zero (Hanowski et al., 2016).

In direct-contact nonlinear moving-domain FPSI, additional geometric terms appear. One representative interface set is

$\Gamma_{fp}$8

$\Gamma_{fp}$9

(Kuan et al., 25 Aug 2025). This suggests that, beyond the small-strain Stokes–Biot regime, FPSI interface laws may simultaneously encode kinematics, filtration, friction, and nonlinear geometric pressure balance.

3. Weak formulations and well-posedness theory

FPSI weak formulations are diverse because the coupled fields have mismatched regularity. One influential strategy is the mixed variational Stokes–Biot formulation with an interface pressure

$\Sigma$0

used as a Lagrange multiplier enforcing the mass-conservation constraint

$\Sigma$1

(Houédanou, 2020). Closely related formulations introduce two interface multipliers, one for the normal traction/pressure and one for the structure-velocity trace, because mixed elasticity variables do not admit a classical trace of the structure velocity (Li et al., 2020). More recent monolithic formulations split the interface transmission into three Lagrange multipliers,

$\Sigma$2

with different interface regularity spaces for normal traction, tangential traction, and Darcy flux (Castro et al., 9 Dec 2025, Castro et al., 20 Jan 2026).

Weak imposition by Nitsche’s method is a second major paradigm. Here the interface laws are recast as Nitsche interface operators, with a parameter $\Sigma$3 selecting symmetric, incomplete, or skew-symmetric variants (Bukac et al., 2014). This is attractive for partitioning because time-lagging of the interface operators produces loosely coupled subproblems.

The analytical theory reflects the model class. For the mixed Stokes–Biot formulation based on weakly symmetric elasticity and mixed Darcy flow, existence and uniqueness are proved for the continuous weak formulation, followed by stable and parameter-robust semi-discrete mixed finite element error estimates (Li et al., 2020). For inertial 3D Biot–Stokes filtration with BJS conditions, a semigroup approach yields $\Sigma$4-semigroup generation of the dynamics operator via Lumer–Phillips, with strong and generalized solutions, and weak solutions obtained by density; the analysis extends to the degenerate case $\Sigma$5 (Avalos et al., 2024).

Fracture FPSI requires a different functional setting because $\Sigma$6 near crack tips. The coupled fluid–fluid weak problem is posed in weighted spaces such as $\Sigma$7, $\Sigma$8, $\Sigma$9, and $\Omega_f$0, with $\Omega_f$1 and $\Omega_f$2 as weights; for fixed $\Omega_f$3, the fluid subproblem is linear, symmetric, and elliptic in this weighted setting, and existence and uniqueness follow from Lax–Milgram under a curvature-dependent coercivity condition (Hanowski et al., 2016).

The largest unresolved regularity obstacle appears in nonlinear moving-domain bulk-bulk FPSI. In the direct-contact setting where the interface is determined by the trace of a bulk Biot displacement, finite-energy regularity only gives $\Omega_f$4, which is not enough to define the geometry pointwise in 2D. The consequence is that no weak-solution framework was available for the original unregularized direct-contact problem; existence results are instead proved for regularized interface problems, either by smoothing problematic Biot displacement terms (Kuan et al., 2023) or by spatial convolution at scale $\Omega_f$5, which defines regularized moving domains and yields a weak solution for each fixed $\Omega_f$6 (Kuan et al., 25 Aug 2025).

4. Discretization strategies and algorithmic architectures

FPSI numerics span monolithic, partitioned, and domain-decomposition architectures. In monolithic mixed finite element methods, the Stokes, Darcy, and elasticity variables are discretized simultaneously. Typical choices include stable Stokes pairs such as MINI, Taylor–Hood, or Crouzeix–Raviart; stable mixed Darcy pairs such as Raviart–Thomas or Brezzi–Douglas–Marini; conforming displacement spaces; and nonmatching interface multiplier spaces induced by Darcy normal traces or elasticity stress traces (Houédanou, 2020, Ambartsumyan et al., 2017, Li et al., 2020). These formulations are particularly attractive in locking-prone regimes because stress and Darcy flux are primary unknowns and interface constraints can be imposed weakly on nonmatching meshes (Li et al., 2020).

Fracture FPSI uses a distinct discretization. The fracture pressure $\Omega_f$7 is approximated with standard first-order Lagrange finite elements on a $\Omega_f$8-dimensional fracture mesh, while the bulk displacement and pressure use XFEM enrichments that capture both discontinuity across the fracture and crack-tip singularities. The displacement space contains Heaviside enrichment $\Omega_f$9 and crack-tip functions $\Omega_p$0; the pressure space uses analogous enrichments with $\Omega_p$1 and $\Omega_p$2 (Hanowski et al., 2016).

Partitioned methods are organized around interface reformulations. Nitsche-based splitting lags interface operators in time, producing loosely coupled schemes whose fluid and Biot subproblems are solved independently; additional stabilization operators are introduced to control time-lagging residuals (Bukac et al., 2014). Robin–Robin methods recast the transmission conditions as Robin conditions with interface parameters and, in newer formulations, auxiliary interface variables representing Robin data. This yields non-iterative decoupled Stokes and Biot solves at each step (Dalal et al., 2024). Fully parallelizable Robin–Robin or Nitsche-type schemes go further by evaluating all interface data at the previous time level so that the fluid and structure systems can be assembled and solved independently and simultaneously, with no sub-iterations (Guo et al., 2024, Wang et al., 22 Mar 2026, He et al., 8 Apr 2026).

Higher-order explicit splitting has also been developed. One representative scheme uses BDF2 time stepping together with Adams–Bashforth extrapolation of interface data through a Robin reformulation, producing a second-order, fully explicit partitioned algorithm for fixed-domain Stokes–Biot FPSI (Wang et al., 18 Jun 2026). At the opposite end of the algorithmic spectrum are hybridized or Schur-complement formulations: non-overlapping domain decomposition reduces the global problem to an interface system for Lagrange multipliers solved by GMRES, with local Stokes or Biot subproblems handled in parallel (Ghumman et al., 15 Jun 2026), while Schur-complement partitioned methods based on interface stress and Darcy-flux multipliers decouple the fully discrete Stokes and Biot solves after solving an interface equation (Castro et al., 20 Jan 2026).

5. Stability, error control, locking, and linear solvers

The stability properties of FPSI time splittings are sharply scheme-dependent. Some recent partitioned methods are proved unconditionally stable. For the linearized fixed-domain Stokes–Biot model, a fully parallelizable loosely coupled scheme in a moving-domain framework satisfies

$\Omega_p$3

so no parameter-dependent CFL restriction is required (Guo et al., 2024). A locking-aware Robin–Robin scheme based on auxiliary variables $\Omega_p$4 and $\Omega_p$5 is likewise proved unconditionally stable and robust with respect to extreme poroelastic parameters (He et al., 8 Apr 2026). The non-iterative Robin–Robin method of a different line of work is also unconditionally stable, with a discrete energy satisfying

$\Omega_p$6

in the analyzed regime (Dalal et al., 2024).

Other splittings remain conditionally stable. The Nitsche time-lagging strategy requires a CFL-like restriction $\Omega_p$7 for general stability and, in some estimates, $\Omega_p$8 (Bukac et al., 2014). The multilayered Stokes–Biot–plate splitting is proved stable only under timestep conditions involving $\Omega_p$9, $\Gamma_{fp}$0, permeability, density, and coupling coefficients, with terms scaling like

$\Gamma_{fp}$1

(Scharf et al., 14 Jul 2025). The second-order BDF2–AB2 explicit splitting requires a parabolic CFL restriction

$\Gamma_{fp}$2

(Wang et al., 18 Jun 2026). Older arterial multilayer FPSI splitting likewise admits a conditional stability bound with a timestep restriction linear in mesh size (Bukac et al., 2013).

Error control also spans a priori and a posteriori regimes. For Stokes/Biot FPSI with a Lagrange multiplier interface formulation, residual-based a posteriori indicators are built from element, jump, and interface residuals; the global estimator

$\Gamma_{fp}$3

is proved reliable and efficient (Houédanou, 2020). For explicit fully discrete splitting, an energy-dissipation framework together with Ritz-type projections yields unconditional a priori estimates of the form

$\Gamma_{fp}$4

(Wang et al., 22 Mar 2026), while the second-order explicit scheme achieves

$\Gamma_{fp}$5

in the relevant norms under the CFL restriction (Wang et al., 18 Jun 2026).

Locking and parameter robustness are central concerns. Mixed elasticity formulations were introduced precisely to avoid storativity-related and Poisson locking, with stability and error constants independent of $\Gamma_{fp}$6 and the compliance lower bound $\Gamma_{fp}$7 (Li et al., 2020). The four-field reformulation with total-pressure-type variable $\Gamma_{fp}$8 is designed to be locking-free as $\Gamma_{fp}$9 and to suppress early-time pressure oscillations when $$-\nabla\cdot \sigma_f(u_f,p_f)=f_f,\qquad \nabla\cdot u_f=q_f \quad \text{in }\Omega_f,$$0, permeability is small, and the time step is small (He et al., 8 Apr 2026).

Solver technology mirrors the saddle-point and interface nature of the equations. Partitioned Nitsche discretizations have been used effectively as preconditioners for monolithic systems, with numerical evidence of almost mesh-independent GMRES behavior (Bukac et al., 2014). Domain decomposition for mixed Stokes–Biot yields a nonsymmetric but positive definite interface operator, fields-of-value estimates for GMRES, and subdomain solves that are naturally parallelizable (Ghumman et al., 15 Jun 2026). Schur-complement Lagrange multiplier methods reduce the coupled problem to interface unknowns approximating normal stress, tangential stress, and Darcy flux; matrix-free preconditioned Krylov solvers then recover decoupled fluid and poroelastic solves at each time step (Castro et al., 20 Jan 2026).

6. Specialized geometries, applications, and current directions

FPSI now covers a wide range of geometries beyond the classical fixed interface. In fractured rock, the coupling of poroelastic deformation, matrix Darcy flow, and lower-dimensional fracture Darcy flow provides a model for deformation–flow interaction with internal fracture endpoints and crack-tip singularities (Hanowski et al., 2016). In multilayer biological structures, the coupled system may include a thin poroelastic or elastic interface layer together with a thick Biot medium, as in artery-like wall models with a Koiter membrane plus a poroelastic wall or a thin Biot plate plus a thick Biot layer (Bukac et al., 2013, Scharf et al., 14 Jul 2025). In thin-interface reductions, a fully averaged poroelastic Kirchhoff plate replaces a volumetric Biot plate so that both elastodynamics and pressure live on a codimension-one interface, simplifying analysis and numerics while retaining BJS coupling and normal-flux continuity (Brandt et al., 17 May 2026).

Moving-domain FPSI has also diversified. ALE Navier–Stokes/Biot schemes have been demonstrated on 2D and 3D problems, including channel flow with poroelastic obstacles and a microfluidic chip with cylindrical poroelastic pillars, with fluid and structure solved in parallel on non-conformal meshes (Guo et al., 2024). Large-deformation variational models address an incompressible Navier–Stokes fluid surrounding a fully saturated poroelastic material that may deform largely, with weak-solution existence shown up to first contact time (Benesova et al., 2021). Nonlinear geometric direct-contact FPSI, where fluid and Biot domains occupy regions of the same spatial dimension separated by an interface defined by the trace of the Biot displacement, is substantially harder: existing existence theories rely on regularized moving interfaces and, in one line of work, a thin-plate insertion followed by a limit $$-\nabla\cdot \sigma_f(u_f,p_f)=f_f,\qquad \nabla\cdot u_f=q_f \quad \text{in }\Omega_f,$$1 to obtain regularized-interface weak solutions (Kuan et al., 25 Aug 2025).

Several comparative points are now clear. First, weak imposition of interface conditions by Nitsche or Lagrange multipliers is not merely a discretization convenience; it is often required by trace regularity or by nonmatching meshes (Ambartsumyan et al., 2017, Li et al., 2020). Second, partitioning is not synonymous with unconditional stability or preserved accuracy: some schemes are unconditionally stable, whereas others incur splitting defects of size $$-\nabla\cdot \sigma_f(u_f,p_f)=f_f,\qquad \nabla\cdot u_f=q_f \quad \text{in }\Omega_f,$$2, $$-\nabla\cdot \sigma_f(u_f,p_f)=f_f,\qquad \nabla\cdot u_f=q_f \quad \text{in }\Omega_f,$$3, or $$-\nabla\cdot \sigma_f(u_f,p_f)=f_f,\qquad \nabla\cdot u_f=q_f \quad \text{in }\Omega_f,$$4, and may degrade the approximation unless additional stabilization or timestep restrictions are imposed (Bukac et al., 2014, Wang et al., 18 Jun 2026). Third, direct-contact nonlinear moving-domain FPSI remains analytically incomplete at the unregularized weak-solution level; the current existence theory is regularized, and weak-classical consistency has been established only under stronger assumptions, including poroviscoelasticity and the existence of a classical solution of the original problem (Kuan et al., 2023).

Taken together, these developments show FPSI to be a mathematically heterogeneous field unified by a small set of structural principles: conservation of normal mass flux, balance of total traction, interfacial tangential slip or no-slip, and reciprocal coupling between pore pressure, Darcy filtration, and solid deformation. The resulting equations support several rigorously analyzed formulations—mixed, Nitsche, multiplier-based, semigroup, variational, XFEM, Robin–Robin, and Schur-complement—and these formulations are tailored to distinct regimes: fixed or moving domains, thin or bulk porous structures, fracture interfaces, multilayer tissues, and locking-prone poroelastic parameter ranges (Brandt et al., 17 May 2026, He et al., 8 Apr 2026, Hanowski et al., 2016).