Coupled Flow and Mechanics in a 3D Porous Media with Line Sources

In this paper, we consider the numerical approximation of the quasi-static, linear Biot model in a 3D domain Ω when the right-hand side of the mass balance equation is concentrated on a 1D line source δ_Λ. This model is of interest in the context of medicine, where it can be used to model flow and deformation through vascularized tissue. The model itself is challenging to approximate as the line source induces the pressure and flux solution to be singular. To overcome this, we here combine two methods: (i) a fixed-stress splitting scheme to decouple the flow and mechanics equations and (ii) a singularity removal method for the pressure and flux variables. The singularity removal is based on a splitting of the solution into a lower regularity term capturing the solution singularities, and a higher regularity term denoted the remainder. With this in hand, the flow equations can now be reformulated so that they are posed with respect to the remainder terms. The reformulated system is then approximated using the fixed-stress splitting scheme. Finally, we conclude by showing the results for a test case relevant for flow through vascularized tissue. This numerical method is found to converge optimally.



There are no comments yet.


page 7


The fixed-stress splitting scheme for Biot's equations as a modified Richardson iteration: Implications for optimal convergence

The fixed-stress splitting scheme is a popular method for iteratively so...

Efficient Solvers for Nonstandard Models for Flow and Transport in Unsaturated Porous Media

We study several iterative methods for fully coupled flow and reactive t...

Regularity and finite element approximation for two-dimensional elliptic equations with line Dirac sources

We study the elliptic equation with a line Dirac delta function as the s...

Physics-informed deep learning for flow and deformation in poroelastic media

A physics-informed neural network is presented for poroelastic problems ...

A meshless method for the solution of incompressible flow equations

This article reports on the efficiency of a co-located diffuse approxima...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The coupling of mechanics and flow in porous media is relevant for a wide range of applications, occurring for instance in geophysics bause ; erlend and medicine simula ; Nagashima1987BiomechanicsOH . We put forward a model relevant for simulating perfusion, i.e., blood flow, and deformation in vascularized tissue. This problem is of high interest in the context of medicine, as clinical measurements of perfusion provide important indicators for e.g. Alzheimer’s disease Chen1977 ; IturriaMedina2016 , stroke Markus353 and cancer Gillies1999 . Moreover, both the tissue and blood vessels are elastic, and these properties constitute another valuable clinical indicator. Vascular compliance, as one example, is reduced in cases of vascular dementia, but not in cases Alzheimer’s disease brain5 .

We consider the fully coupled quasi-static, linear Biot system biot1 , modeling a poroelastic media when the source term in the flow equation is concentrated on . Let denote a bounded open 3D domain with smooth boundary and a collection of straight line segments embedded in . The equations on the space-time domain read:

Find such that:


where denotes the displacement,

the (linear) strain tensor,

the pressure, the Darcy’s flux, the Biot coefficient, the permeability tensor divided by the fluid viscosity, the fluid density,

the gravity vector,

the Biot modulus, and are Lamé parameters, the contribution from body forces and a source term . We assume to be strictly positive and uniformly bounded and , . Additionally, is assumed scalar-valued, as required by the singularity removal method. For simplicity, we use homogeneous boundary conditions , on and initial conditions , in .

The system eq:h1-eq:h3 is made non-standard by a generalized Dirac line source of intensity in the right-hand side of eq:mass. The line source is defined mathematically as

Physically, it is introduced to model the mass exchange between the vascular network and the surrounding tissue. This exchange occurs through the capillary blood vessels, which have radii ranging from 5 to 10 micrometer. These blood vessels are too small to be captured as 3D objects in a mesh; instead, they are typically reduced to being one-dimensional line segments, see e.g. 1dmodel1 ; koppl2017 ; Laurino2019 ; ingeborg3 ; vidotto2018 . The system eq:h1-eq:h3 would then be on the same form as the one considered in simula , with the exception that the exchange term is here concentrated on the 1D vascular network.

We take and focus our attention on the challenges introduced by the line source . The line source induces and to be singular, i.e., they both diverge to infinity on . Consequently, one has but ; the solution is then not regular enough to fit the analytic framework of source22 ; florin . They both prove global convergence for the fixed-stress splitting scheme applied to Biot’s equations. Moreover, these singularities are expensive to resolve numerically, making the solution highly challenging to approximate.

In order to tackle this issue, we here combine two strategies: (i) a fixed-stress splitting scheme that decouples the mechanics equation eq:h1 from the flow equations eq:mass-eq:h3, and (ii) a singularity removal method for the flow equations. For an introduction to the fixed-stress splitting scheme, we refer to the works of Mikelić et al. source22 and Both et al. florin ; for an introduction to the singularity removal method, we refer to our earlier work ingeborg ; ingeborg2 .

2 Mathematical Model and Discretization

In this section, we begin by introducing a splitting method that decomposes and into higher and lower regularity terms. Here, the lower regularity terms are given explicitly. The model eq:h1-eq:h3 can be reformulated so that it is given with respect to the higher regularity terms; we refer to this as the singularity removal method. Next, we show how this model can readily be approximated by the fixed-stress splitting scheme.

2.1 Singularity Removal Method

For the sake of notational simplicity, we assume to be constant; a spatially varying could be handled as shown in (ingeborg, , Sect. 3.3). Let denote the endpoints of the line segment . From (ingeborg, , Sect. 3.2), we have a function defined as


with , , and as the normalized tangent vector of . Centrally, this function solves in the appropriate weak sense; i.e., we have:

Having this function in hand, we next formulate the following splitting ansatz:


where and . The terms and capture the singular part of the solution, and are explicitly given via the function . This allows and

to enjoy higher regularity and improved approximation properties. Assume for the moment that the solution

is given. Inserting the splitting eq:phelp into eq:mass-eq:h3 one finds the following reformulated flow equation:

Find such that:


where . Here, eq:h3-reform is straightforward to obtain. For eq:mass-reform, we used that

In the last line we used the product rule to obtain along with the relation .

The value of the reformulation lies in the fact that can now be expected to belong to . To see this, note that by assumption. can be shown by straightforward calculation; it follows that . Finally, one can show that ; for verification of this, we refer to the calculations in (ingeborg2, , Sect. 4.2) along with the embedding .

Let now eq:h1 and eq:mass-reform-eq:h3-reform denote the reformulated Biot equation. As , this system fits the analytic framework of florin .

2.2 Fixed-Stress Splitting Scheme

Next, we show how the reformulated Biot equations eq:h1 and eq:mass-reform-eq:h3-reform can be approximated via the fixed-stress splitting scheme from florin . Let be the triangularization of the domain with mesh size . We let be a partition of the time interval with and define a constant time step size for . To discretize the system, we employ backward Euler for time and a finite element method for space. The solutions are approximated with linear piecewise polynomials, constant piecewise polynomials and lowest-order Raviart-Thomas spaces for the displacement, pressure and flux, respectively. The discrete spaces are given by:

with and as the linear and constant piecewise polynomials.

Take now to be the -inner product and to be the initial values of the solution. We assume the solution of the displacement, pressure and flux is known for the previous time step. The time-discretization of eq:h1 and eq:mass-reform-eq:h3-reform then reads:

Given , find such that

for all .

The idea of the fixed-stress splitting scheme is to decouple the flow and mechanics equation while keeping an artificial volumetric stress constant. Here, is referred to as the drained bulk modulus. We consider the theoretically optimal tuning parameter with florin .
We define a sequence , . Let denote the current iteration step and denote the previous iteration step. Then initialize , and by , and , respectively. The algorithm iterates until a stopping criterion is reached. The full scheme reads:

Step 1: Given . Find such that :


Step 2: Update the full pressure and flux solutions: and .
Step 3: Given . Find , such that :


3 Numerical Results

In this section, we provide numerical convergence results for a test case using parameters relevant for flow through vascularized tissue. Let the medium of consideration be an isotropic, homogeneous porous medium and a positive scalar quantity. We let be the solutions at iteration step and the solutions at the previous iteration step . The procedure stops when reaching the following criterion:

where are given tolerances (see Table 1).

Symbol: Quantity: Value: Reference:
Permeability divided mm mPa s (Weiss2017, , Table 1)
by the fluid viscosity
Tuning parameter mPa (simula, , Table 6)
Biot modulus mPa (Guo2017, , Table 2)
Biot coefficient
Poisson’s ratio
Gravitational vector mm s biotbook
Final time s
Time step s
Absolute error tolerance
Relative error tolerance
Table 1: Material parameters used to solve the Biot’s equations with lower dimensional source terms eq:h1-eq:h3 in Sect. 3. There is a wide rage of parameters used in literature. The ones represented here are a sample of representative parameters.

Let be a cube discretized by tetrahedrons. The numerical results are obtained by the fixed-stress splitting scheme proposed in Sect. 2.2 and the programming platform FEniCS LoggMardalEtAl2012a . Convergence is tested against the following analytic solutions

where is a pulsative intensity function. We selected two points and to describe the line segment. Then computed the solutions using mixed-finite element formulations for the correction terms and , and the solution displacement is calculated with the conformal finite element formulation.

Table 2 shows the error and convergence rates obtained using the parameters listed in Table 1. The singularity removal based fixed-stress splitting scheme is seen to converge optimally for each variable , , and . The plots for this problem are illustrated in Fig. 1. The figure includes the plot of the full pressure, the magnitude of the full flux and the magnitude of the displacement, accordingly.


Figure 1: Left: Plot of the reconstructed pressure. Middle: Magnitude of the full flux. Right: Magnitude of the displacement. All plots are numerical solutions obtained by the fixed-stress splitting scheme eq:FSSSs2-eq:FSSS2 with one line source.
Table 2: Errors and convergence rates obtained solving eq:FSSSs2-eq:FSSS2 with analytical solutions found in this section. For reference, optimal convergence rates are listed in the bottom row.


  • (1) Bause, M., Radu, F.A., Köcher, U.: Space-time finite element approximation of the Biot poroelasticity system with iterative coupling. Computer Methods in Applied Mechanics and Engineering 320, 745 – 768 (2017)
  • (2) Biot, M.A.: General theory of three‐dimensional consolidation. Journal of Applied Physics 12(2), 155 – 164 (1941)
  • (3) Both, J.W., Borregales, M., Kumar, K., Nordbotten, J.M., Radu, F.A.: Robust fixed stress splitting for Biot’s equations in heterogeneous media. Applied Mathematics Letters 68, 101 – 108 (2017)
  • (4) Chen, Y., Wolk, D., Reddin, J., Korczykowski, M., Martinez, P., Musiek, E., Newberg, A., Julin, P., Arnold, S., Greenberg, J., Detre, J.: Voxel-level comparison of arterial spin-labeled perfusion mri and fdg-pet in alzheimer disease. Neurology 77(22), 1977–1985 (2011)
  • (5) D’Angelo, C., Quarteroni, A.: On the coupling of 1D and 3D diffusion-reaction equations: Application to tissue pefusion problems. Mathematical Models and Methods in Applied Sciences 18(8), 1481–1504 (2008)
  • (6) Dhoat, S., Ali, K., Bulpitt, C.J., Rajkumar, C.: Vascular compliance is reduced in vascular dementia and not in alzheimer’s disease. Age and Ageing 37(6), 653–659 (2008)
  • (7) Formaggia, L., Quarteroni, A., Veneziani, A.: Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System, vol. 1 (2009)
  • (8) Gillies, R.J., Schomack, P.A., Secomb, T.W., Raghunand, N.: Causes and effects of heterogeneous perfusion in tumors. Neoplasia 1(3), 197–207 (1999)
  • (9) Gjerde, I., Kumar, K., Nordbotten, J.M.: A singularity removal method for coupled 1d-3d flow models (2018). ArXiv:1812.03055 [math.AP]
  • (10) Gjerde, I., Kumar, K., Nordbotten, J.M.: A mixed approach to the poisson problem with line sources (2019). ArXiv:1910.11785 [math.AP]
  • (11) Gjerde, I., Kumar, K., Nordbotten, J.M., Wohlmuth, B.: Splitting method for elliptic equations with line sources. ESAIM: M2AN 53(5) (2019)
  • (12) Guo, L., Vardakis, J., Lassila, T., Mitolo, M., Ravikumar, N., Chou, D., Lange, M., Sarrami-Foroushani, A., Tully, B., Taylor, Z., Varma, S., Venneri, A., Frangi, A., Ventikos, Y.: Subject-specific multiporoelastic model for exploring the risk factors associated with the early stages of alzheimer’s disease. Interface Focus 8(1) (2017)
  • (13) Iturria-Medina, Y., , Sotero, R.C., Toussaint, P.J., Mateos-Pérez, J.M., Evans, A.C.: Early role of vascular dysregulation on late-onset alzheimer’s disease based on multifactorial data-driven analysis. Nature Communications 7(1) (2016)
  • (14) Köppl, T., Vidotto, E., Wohlmuth, B., Zunino, P.: Mathematical modeling, analysis and numerical approximation of second-order elliptic problems with inclusions. Mathematical Models and Methods in Applied Sciences 28(05), 953–978 (2018)
  • (15) Laurino, F., Zunino, P.: Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction. ESAIM: Mathematical Modelling and Numerical Analysis 53(6), 2047–2080 (2019)
  • (16) Lee, J., Piersanti, E., Mardal, K.A., Rognes, M.: A mixed finite element method for nearly incompressible multiple-network poroelasticity. SIAM Journal on Scientific Computing 41(2), A722–A747 (2019)
  • (17) Logg, A., Mardal, K.A., Wells, G.N., et al.: Automated Solution of Differential Equations by the Finite Element Method. Springer (2012)
  • (18) Markus, H.S.: Cerebral perfusion and stroke. Journal of Neurology, Neurosurgery & Psychiatry 75(3), 353–361 (2004)
  • (19) Mikelić, A., Wheeler, M.F.: Convergence of iterative coupling for coupled flow and geomechanics. Computational Geosciences 17(3), 455–461 (2013)
  • (20) Nagashima, T., Tamaki, N., Matsumoto, S., Horwitz, B., Seguchi, Y.: Biomechanics of hydrocephalus: a new theoretical model. Neurosurgery 21(6), 898–904 (1987)
  • (21) Storvik, E., Both, J.W., Kumar, K., Nordbotten, J.M., Radu, F.A.: On the optimization of the fixed-stress splitting for Biot’s equations. International Journal for Numerical Methods in Engineering 120(2), 179–194 (2019)
  • (22) Vidotto, E., Koch, T., Köppl, T., Helmig, R., Wohlmuth, B.: Hybrid models for simulating blood flow in microvascular networks. Multiscale Modeling & Simulation 17(3), 1076–1102 (2019)
  • (23) Weiss, C.: Finite element analysis for model parameters distributed on a hierarchy of geometric simplices. GEOPHYSICS 82(4), 1–52 (2017)