Parameter-robust methods for the Biot-Stokes interfacial coupling without Lagrange multipliers

In this paper we advance the analysis of discretizations for a fluid-structure interaction model of the monolithic coupling between the free flow of a viscous Newtonian fluid and a deformable porous medium separated by an interface. A five-field mixed-primal finite element scheme is proposed solving for Stokes velocity-pressure and Biot displacement-total pressure-fluid pressure. Adequate inf-sup conditions are derived, and one of the distinctive features of the formulation is that its stability is established robustly in all material parameters. We propose robust preconditioners for this perturbed saddle-point problem using appropriately weighted operators in fractional Sobolev and metric spaces at the interface. The performance is corroborated by several test cases, including the application to interfacial flow in the brain.



There are no comments yet.


page 3

page 19

page 20


The Biot-Stokes coupling using total pressure: formulation, analysis and application to interfacial flow in the eye

We consider a multiphysics model for the flow of Newtonian fluid coupled...

Discrete Adjoint Momentum-Weighted Interpolation Strategies

This Technical Note outlines an adjoint complement to a critical buildin...

Mixed methods for large-strain poroelasticity/chemotaxis models simulating the formation of myocardial oedema

In this paper we propose a novel coupled poroelasticity-diffusion model ...

The Darcy problem with porosity depending exponentially on the pressure

We consider the flow of a viscous incompressible fluid through a porous ...

Parameter-robust Multiphysics Algorithms for Biot Model with Application in Brain Edema Simulation

In this paper, we develop two parameter-robust numerical algorithms for ...

A novel pressure-free two-fluid model for one-dimensional incompressible multiphase flow

A novel pressure-free two-fluid model formulation is proposed for the si...
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

1.1 Scope

We address the construction of appropriate monolithic solvers for multiphysics fluid-poromechanical couplings interacting through an interface. Particular attention is payed to tracking parameter dependence of the continuous and discrete formulations so that the resulting numerical methods are robust with respect to typical scales in material constants spanning over many orders of magnitude. We adopt a multi-domain approach, where appropriate conditions for the coupling through the shared interface need to be imposed. We use the conditions proposed in Murad et al. (2001) (although, other forms and dedicated phenomena could be incorporated without much effort, such as fluid entry resistance Showalter (2005); Bergkamp et al. (2020)). The recent literature contains various numerical methods for (Navier-)Stokes/Biot interface formulations including mixed, double mixed, monolithic, segregated, conforming, non-conforming, and DG discretizations Ager et al. (2019); Ambartsumyan et al. (2018); Badia et al. (2009); Bukač et al. (2015); Caucao et al. (2020); Cesmelioglu (2017); Cesmelioglu and Chidyagwai (2020); Li and Yotov (2020); Taffetani et al. (2021); Wen and He (2020); Wen et al. (2021); Wilfrid (2020).

In Taffetani et al. (2021); Ruiz-Baier et al. (2021) (and starting from the Biot-Stokes equations advanced in Ambartsumyan et al. (2018, 2019)) the authors rewrite the poroelasticity equations using displacement, fluid pressure and total pressure (also as in the poromechanics formulations from Bürger et al. (2021); Kumar et al. (2020); Oyarzúa and Ruiz-Baier (2016)). Since fluid pressure in the poroelastic domain has sufficient regularity, no Lagrange multipliers are needed to enforce the coupling conditions, which resembles the different formulations for Stokes-Darcy advanced in Boon et al. (2021b); Chidyagwai et al. (2016); Holter et al. (2021); Karper et al. (2009). Another advantage of the three-field Biot formulation is its robustness with respect to the Lamé constants of the poroelastic structure. This robustness is of particular importance when we test the flow response to changes in the material properties of the skeleton and when the solid is nearly incompressible. The work Taffetani et al. (2021) focuses on the stability analysis and its precise implications on the asymptotics of the interface conditions when the permeability depends on porosity heterogeneity, whereas Ruiz-Baier et al. (2021) addresses the stability of the semi- and fully discrete problems, and the application to interfacial flow in the eye. Here, we extend these works by concentrating on deriving robust stability, on designing efficient block preconditioners (robust with respect to all material parameters) following the general theoretical formalism from Mardal and Winther (2011), and on the simulation of free flow interacting with interstitial flow in the brain. In such a context (and in the wider class of problems we consider in this paper), tissue permeability is of the order of 10m, and the incorporation of tangential interface transmission conditions usually involves terms that scale inversely proportional to the square root of permeability. Moreover, the solid is nearly incompressible, making the first Lamé parameter significantly larger than the other mechanical parameters and exhibiting volumetric locking for some types of displacement-based formulations. Other flow regimes that are challenging include low-storage cases Mardal et al. (2021). It is then important that the stability and convergence of the numerical approximations are preserved within the parameter ranges of interest.

Here we follow Lee et al. (2017, 2019); Oyarzúa and Ruiz-Baier (2016); Boon et al. (2021a) and use parameter-weighted norms to achieve robustness. However, as we will see, combining proper preconditioners for Stokes and Biot single-physics problems is not sufficient for the interface coupled problem. In fact, the condition number of the preconditioned system, although robust in mesh size, grows like the square root of the ratio between fluid viscosity and permeability. This phenomenon is demonstrated in Example 2.1, below. That is, the efficiency of seemingly natural preconditioners varies with the material parameters. In order to regain stability with respect to all parameters, we include both an additional fractional term involving the pressure and a metric term coupling the tangential fluid velocity and solid displacement at the interface, hereby increasing the regularity at the interface in a proper parameter dependent manner. This strategy draws inspiration from similar approaches employed in the design of robust solvers for Darcy and Stokes-Darcy couplings Bærland et al. (2020); Boon et al. (2021b); Holter et al. (2021).

1.2 Outline

We have organised the contents of this paper in the following manner. The remainder of this section contains preliminaries on notation and functional spaces to be used throughout the manuscript. Section 2 outlines the main details of the balance equations, stating typical interfacial and boundary conditions, and restricting the discussion to the steady Biot-Stokes coupled problem. There we also include the weak formulation and demonstrate that simple diagonal preconditioners based on standard norms do not lead to robustness over the whole parameter range. This issue is addressed in Section 3 where we show well-posedness of the system using a global inf-sup argument with parameter weighted operators in fractional spaces, which in turn assist in the design of robust solvers by operator preconditioning. Section 4 discusses finite element discretization of the coupled problem using both conforming and non-conforming elements; and it also contains numerical experiments demonstrating robustness of the fractional preconditioner and its feasibility for simplified simulations of interfacial flow in the brain.

1.3 Preliminaries

Let us consider a spatial domain , where , disjointly split into and . These subdomains respectively represent the region filled with an incompressible fluid and the elastic porous medium (a deformable solid matrix or an array of solid particles). We will denote by

the unit normal vector on the boundary

, and by the interface between the two subdomains, which is assumed sufficiently regular. We adopt the convention that on the normal vector points from to . We also define the boundaries and . The sub-boundary is further decomposed between and where we impose, respectively, no slip velocities and zero normal total stresses. Similarly, we split into and where we prescribe zero traction and clamped boundaries, respectively. For the analysis, the setup of trace spaces needs that and that , which can be satisfied if the interface meets the boundary at the Biot displacement and Stokes velocity boundaries (see Figure 1.1, left), where and . Our numerical tests will also include cases where the interface intersects boundaries and on the Stokes and Biot sides, respectively.

Figure 1.1: A configuration of subdomains and boundary partition for the Biot-Stokes coupled problem. (Left) Setup assumed for the analysis in Theorem 3.1. (Right) Example of another configuration investigated by some of the numerical experiments in Section 4.4.

For generic Sobolev spaces and a scalar , the weighted space refers to endowed with the norm . The intersection provided with the norm , is a Hilbert space Bærland et al. (2020); Bergh and Löfström (2012).

Vector fields and vector-valued spaces will be written in boldface. In addition, by we will denote the usual Lebesgue space of square integrable functions and denotes the usual Sobolev space with weak derivatives of order up to in , and denotes its vector counterpart. In addition, will denote the space of functions , , such that .

An (as well as ) inner product over a generic bounded domain is denoted as . The symbol will denote the pairing between the trace functional space and its dual , and we will also write to denote other, more general, duality pairings. Moreover, for , its normal and tangential traces defined by bounded surjective maps, will be denoted by and , respectively Holter et al. (2021). We also remark that, while is a subspace of as a set, the tangential trace of is in and that in our setting it is important to track the parameter dependence for the sake of robustness.

Pertaining to the poroelastic domain, and assuming momentarily that , we recall from (Gatica, 2014, Section 2.4.2) the definition of the space


supplied with the norm

where denotes the extension-by-zero operator

Furthermore, the restriction of to , defined as

belongs to the dual of . Its norm is


The boundary setup in Figure 1.1 is such that , and therefore a modification of (1.1)-(1.2) is required. With that purpose, we note that any can be continuously extended to in the space


where (see Figure 1.1, left), and

This treatment can be interpreted as replacing by the space of functions such that , where is a sufficiently regular trace function, positive on , and vanishing only on (see, for example, le Roux and Reddy (1993)).

The norm of the resulting extension is defined analogously as before,


and therefore the restriction of a distribution to the interface, , is in the dual space and its norm is


2 Governing equations and weak formulation

The momentum and mass balance equations for the flow in the fluid cavity are given by Stokes equations written in terms of fluid velocity and fluid pressure , whereas the non-viscous filtration flow through the porous skeleton can be described by Darcy’s law in terms of pressure head , and the porous matrix elastostatics are stated in terms of the solid displacement . The coupled Biot-Stokes equations arising after a backward Euler semi-discretization in time, with time step , read


where is a vector field of body loads, is the gravity acceleration, is the fluid viscosity,

is the strain rate tensor, and

is the infinitesimal strain tensor, are the density of the fluid and solid, respectively, are the first and second Lamé constants of the solid, is the heterogeneous tensor of permeabilities (satisfying a.e. in and for all ); is a source/sink term for the fluid pressure (which also includes pressures in the previous backward Euler time step); and are the total storage capacity and Biot-Willis poroelastic coefficient. Here we have used the total pressure , as an additional unknown Lee et al. (2017); Oyarzúa and Ruiz-Baier (2016).

We furthermore supply boundary conditions as follows


In order to close the system, we consider the classical transmission conditions on accounting for the continuity of normal fluxes, momentum balance, equilibrium of fluid normal stresses, and the so-called Beavers-Joseph-Saffman condition for tangential fluid forces Bukač et al. (2015); Showalter (2005), which in the present setting reduce to


where is the slip rate coefficient depending on the geometry of the domain, and we recall that the normal on the interface is understood as pointing from the fluid domain towards the porous structure , while , are orthonormal tangent vectors on , normal to .

We proceed to test (2.1a)-(2.1e) against suitable smooth functions and to integrate over the corresponding subdomain. The challenging model parameters are , , , , , and the magnitude of . Therefore, we will concentrate on the specific case where , and the model parameters (in particular ) are spatially constant. Following Karper et al. (2009), after applying integration by parts wherever adequate and using the transmission conditions (2.3a)-(2.3d), we arrive at the following remainder on the interface

which is well-defined and therefore no additional Lagrange multipliers are required to realize the coupling conditions. Also, in view of the character of the resulting variational forms in combination with the specification of boundary conditions (2.2), we define the Hilbert spaces

and the product space


Consequently, we have the following weak form for the Biot-Stokes coupling: Find such that



System (2.5) differs from that analyzed in Ruiz-Baier et al. (2021) in the ordering of the unknowns, and in that we obtain a symmetric multilinear formulation defined by a global operator (the coefficient matrix of the left-hand side of (2.5)) of the form


where the dependence on the model parameters is clearly identified. In particular, the interface coupling terms on the first off-diagonal blocks depend on the inverse of permeability.

We note that can be regarded as defining a perturbed saddle-point problem, with

where the composing blocks are defined as


In turn, well-posedness of the Biot-Stokes system (2.5) in the product space (2.4) can be established using the abstract Brezzi-Braess theory Braess (1996), after invoking separately the solvability and stability results for the Stokes subproblem Girault and Raviart (1986) and the Biot subproblem in the three-field total pressure formulation Lee et al. (2017). However, such a decoupled approach does not lead to stability independent of the material parameters and consequently, preconditioners based on the standard norms (also referred to as single-physics or sub-physics preconditioners) are not necessarily parameter robust. This issue is demonstrated next in Example 2.1.

Example 2.1 (Simple preconditioners using standard norms).

We consider the Biot-Stokes formulation (2.5) defined on the subdomains , with boundary conditions such that the left edge of is a no-slip boundary while the top and bottom edges will form . Similarly, the top and bottom edges on the Biot side are considered stress-free while the right edge is clamped.

Based on the well-posedness of (2.5) in the space (cf. (2.4)), we can readily consider a solution space with weighted inner product leading to the Riesz map (diagonal) preconditioner


We remark that the first and third blocks of together define a parameter robust preconditioner for the (standalone) Stokes problem and the remaining blocks form the robust three-field Biot preconditioner Lee et al. (2017). However, in the subproblem preconditioners are decoupled.

Alternatively, after observing that the operator in (2.7a) defines a norm over the velocity-displacement space we will also investigate the (block-diagonal) preconditioner


We note that in (2.9) the tangential components of the Stokes velocity and of the Biot displacement are coupled. In this sense, the preconditioner captures the interaction between the subsystems and, in particular, the coupling through the Beavers-Joseph-Saffman condition (2.3d).

To investigate the robustness of these preconditioners, we set typical physical parameters in (2.5) except for and , which are to be varied. Using a discretization in terms of the lowest-order Taylor-Hood elements (see more details in Section 4.1), we next consider the boundedness of the number of iterations of the preconditioned MinRes solver under mesh refinement and parameter variations. More precisely, using , (inverted by LU) we compare the number of iterations required for convergence determined by reducing the preconditioned residual norm by a factor . The initial vector is taken as random.

We report the results in Table 2.1. It can be seen that the number of MinRes iterations produced with the diagonal preconditioner is rather sensitive to variations in both and . In comparison, when fixing , the iterations appear to be more stable in when (2.9) is used. However, if is set, there is a clear deterioration of the performance for small values of also with the preconditioner .

Preconditioner (2.8) Preconditioner (2.9)
211 240 258 264 71 82 89 89
76 76 74 73 50 49 48 48
41 41 41 41 37 37 36 36
394 184 693 471 639
59 59 59 58 58 57 55 55
41 41 41 41 37 37 36 36
Table 2.1: Performance of preconditioners (2.8) and (2.9) for the Biot-Stokes system (2.5) in Example 2.1. Only the parameters and are varied away from 1. The lack of convergence after 750 MinRes iterations is indicated as –.

The improved performance of the preconditioner , which preserves the tangential coupling of the Stokes and Biot problems in (2.3d), over , where the components are decoupled, suggests to strengthening the coupling (involving the tangential traces) in order to obtain parameter robustness. However, from the point of view of the interface conditions (2.3a)-(2.3d), it is clear that the coupling in the normal direction is missing in .

With the above idea in mind, we proceed to establish well-posedness of (2.5) in the product space equipped with non-standard norms that include additional control at the interface reflecting/arising from the mass conservation condition (2.3a).

3 Well-posedness of the Biot-Stokes system

Let us group the variables as and and introduce the weighted norm




where the fractional norm is defined in (1.5).

In turn, we define the weighted product space as the space that contains all that are bounded in this norm. The subscript encodes the collection of weighting parameters . Moreover, the space allows for the natural decomposition:

Theorem 3.1.

Problem (2.5) is well-posed in the space equipped with the norm (3.1). In other words, the operator in (2.6) is a symmetric isomorphism satisfying


where are positive constants independent of .


The operator norm is defined as

and condition (3.3a) states the continuity of . To show this, we first use the Cauchy-Schwarz inequality to derive


It therefore remains to show that is continuous. Another application of the Cauchy-Schwarz inequality on the different terms together with a trace inequality provides this result.

In order to prove the second relation (3.3b), we aim to verify the assumptions of the Banach-Nečas-Babuška (BNB) theorem (see, e.g., Ern and Guermond (2004)). In particular, we aim to prove that


We do this by assuming that is given and by constructing an appropriate test function . Following, e.g., Anaya et al. (2021); Lee et al. (2017), we choose and , giving


Next, the inf-sup condition proven in Lemma 3.1, below, allows us to construct such that


We now use this test function, scaled by a constant to be chosen later, and using (3.2a)-(3.2c) along with (3.4) and (3.7), we derive

where we have also used Cauchy-Schwarz and Young’s inequality. Setting gives us


Finally, we take and put together (3) and (3.8) to arrive at

The combination of these two bounds shows that (3.5) holds. The BNB theorem now provides (3.3b). ∎

Lemma 3.1.

There exists a such that


The proof follows similarly to (Boon et al., 2021b, Section 3). Let be given. We proceed in five steps.

  1. With a given total pressure in the Biot domain, we set up an auxiliary Stokes problem: Find that weakly satisfy

    subject to the mixed boundary conditions

    By the well-posedness of this auxiliary problem (for a proof see, e.g., (Girault and Raviart, 1986, Chapter I)), the first component of the solution satisfies

  2. We next consider that the trace of is a distribution in , and focus on its restriction to the interface, , belonging to (cf., the end of Section 1).

    Let be the Riesz representative of , and consider a Stokes-extension of into by setting up another auxiliary Stokes problem (and still in the Biot domain): Find a pair