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).
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.
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.
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
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, andis 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.
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)|
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 .
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:
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.
Next, the inf-sup condition proven in Lemma 3.1, below, allows us to construct such that
where we have also used Cauchy-Schwarz and Young’s inequality. Setting gives us
There exists a such that
The proof follows similarly to (Boon et al., 2021b, Section 3). Let be given. We proceed in five steps.
With a given total pressure in the Biot domain, we set up an auxiliary Stokes problem: Find that weakly satisfy
in 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
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