Discontinuous Galerkin method for a three-dimensional coupled fluid-poroelastic model with applications to brain fluid mechanics

TL;DR

This work extends a polytopal discontinuous Galerkin method to a 3D fluid–poroelastic coupling problem featuring multi-network poroelasticity (MPE) and (Navier-)Stokes CSF flow, with Beavers-Joseph-Saffman interface conditions. It provides a stable, optimally convergent semi-discrete scheme and a fully discrete time-stepping strategy (Newmark for elasticity, θ-method for the rest), along with a rigorous a priori error analysis. The study compares Stokes vs. Navier–Stokes for CSF and quantifies the impact of BJS friction on interface dynamics, displacement, and pressure in an idealized brain geometry. It also demonstrates notable computational advantages of polyhedral meshes over traditional hexahedral meshes for domains with fine geometric features, offering practical insights for brain fluid-poromechanics simulations and paving the way for more personalized, efficient simulations.

Abstract

The modeling of the interaction between a poroelastic medium and a fluid in a hollow cavity is crucial for understanding, e.g., the multiphysics flow of blood and Cerebrospinal Fluid (CSF) in the brain, the supply of blood by the coronary arteries in heart perfusion, or the interaction between groundwater and rivers or lakes. In particular, the cerebral tissue's elasticity and its perfusion by blood and interstitial CSF can be described by Multi-compartment Poroelasticity (MPE), while CSF flow in the brain ventricles can be modeled by the (Navier-)Stokes equations, the overall system resulting in a coupled MPE-(Navier-)Stokes system. The aim of this paper is three-fold. First, we aim to extend a recently presented discontinuous Galerkin method on polytopal grids (PolyDG) to incorporate three-dimensional geometries and physiological interface conditions. Regarding the latter, we consider the Beavers-Joseph-Saffman (BJS) conditions at the interface, modeling the friction between the fluid and the porous medium. Second, we analyze the computational efficiency of the proposed method on a domain with small geometrical features, thus demonstrating the advantages of employing polyhedral meshes. Finally, by a comparative numerical investigation, we assess the fluid-dynamics effects of the BJS conditions and of employing either Stokes or Navier-Stokes equations to model the CSF flow. The semidiscrete numerical scheme for the coupled problem is stable and optimally convergent. Temporal discretization is obtained using Newmark's $β$-method for the elastic wave equation and the $θ$-method for the remaining equations. The theoretical error estimates are verified by numerical simulations on a test case with a manufactured solution, and a numerical investigation is carried out on a three-dimensional geometry to assess the effects of interface conditions and fluid inertia on the system.

Figures (10)

• Figure 1: Computational domain: poroelastic region $\Omega_{\rm el}$ and fluid region $\Omega_{\rm f}$, interface $\Sigma$ between them (blue), and external boundaries $\Gamma_\text{out}$ (red), $\Gamma_\text{w}$ (grey), $\Gamma_\text{D}$ (light grey) and $\Gamma_\text{N}$ (green).
• Figure 2: Polygonal elements sharing an internal face (left) or a face on the interface $\Sigma$ (right).
• Figure 3: Verification test of \ref{['sec:conv']}: computed relative errors in the energy norm \ref{['eq:normcoerc']} versus $h$ for different polynomial degrees $m=1,2,3$ (log-log scale).
• Figure 4: Test case of \ref{['sec:polyvsstd']}. Left: computational domain with small inclusions (magenta holes). Center: hexahedral mesh $\mathscr T_{26816}$ colored according to the agglomerated polyhedral mesh $\mathscr T_{45}$. Right: clip and zoom of $\mathscr T_{45}$ close to some inclusions (in white).
• Figure 5: Test case of \ref{['sec:polyvsstd']}: computational costs using a hexahedral mesh of $N=26816$ elements and two different polyhedral meshes of $N=45$ or $N=91$ elements, for polynomial degrees $m=1,2,3$ (on the hexahedral mesh) and $m=1,2,3,4,5,6$ (on the polyhedral meshes). Left: convergence errors $E_{L^2}, E_{H^1}$. Right: computational time for assembling the linear system and its solution.

Theorems & Definitions (8)

• Remark 1: Application to brain fluid-poromechanics
• Remark 2: Skew-symmetry of the advection form $N_{\rm f}$
• Remark 3: Derivation of the interface forms $\mathfrak J$ and $\mathfrak G$
• Theorem 1: Stability estimate
• Remark 4
• Theorem 2: A priori error estimate
• Remark 5
• Lemma 1