1 Introduction
1.1 Scope
We address the construction of appropriate monolithic solvers for multiphysics fluidporomechanical 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 multidomain 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, nonconforming, 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); RuizBaier et al. (2021) (and starting from the BiotStokes 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 RuizBaier (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 StokesDarcy advanced in Boon et al. (2021b); Chidyagwai et al. (2016); Holter et al. (2021); Karper et al. (2009). Another advantage of the threefield 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 RuizBaier 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 displacementbased formulations. Other flow regimes that are challenging include lowstorage 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 RuizBaier (2016); Boon et al. (2021a) and use parameterweighted norms to achieve robustness. However, as we will see, combining proper preconditioners for Stokes and Biot singlephysics 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 StokesDarcy 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 BiotStokes 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 wellposedness of the system using a global infsup 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 nonconforming 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 subboundary 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 vectorvalued 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
(1.1) 
supplied with the norm
where denotes the extensionbyzero operator
Furthermore, the restriction of to , defined as
belongs to the dual of . Its norm is
(1.2) 
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
(1.3) 
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,
(1.4) 
and therefore the restriction of a distribution to the interface, , is in the dual space and its norm is
(1.5) 
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 nonviscous 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 BiotStokes equations arising after a backward Euler semidiscretization in time, with time step , read
(2.1a)  
(2.1b)  
(2.1c)  
(2.1d)  
(2.1e) 
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 BiotWillis poroelastic coefficient. Here we have used the total pressure , as an additional unknown Lee et al. (2017); Oyarzúa and RuizBaier (2016).We furthermore supply boundary conditions as follows
(2.2a)  
(2.2b)  
(2.2c)  
(2.2d) 
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 socalled BeaversJosephSaffman condition for tangential fluid forces Bukač et al. (2015); Showalter (2005), which in the present setting reduce to
(2.3a)  
(2.3b)  
(2.3c)  
(2.3d) 
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 welldefined 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
(2.4) 
Consequently, we have the following weak form for the BiotStokes coupling: Find such that
(2.5a)  
(2.5b)  
(2.5c)  
(2.5d)  
(2.5e) 
where
System (2.5) differs from that analyzed in RuizBaier 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 lefthand side of (2.5)) of the form
(2.6) 
where the dependence on the model parameters is clearly identified. In particular, the interface coupling terms on the first offdiagonal blocks depend on the inverse of permeability.
We note that can be regarded as defining a perturbed saddlepoint problem, with
where the composing blocks are defined as
(2.7a)  
(2.7b) 
In turn, wellposedness of the BiotStokes system (2.5) in the product space (2.4) can be established using the abstract BrezziBraess 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 threefield 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 singlephysics or subphysics 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 BiotStokes formulation (2.5) defined on the subdomains , with boundary conditions such that the left edge of is a noslip boundary while the top and bottom edges will form . Similarly, the top and bottom edges on the Biot side are considered stressfree while the right edge is clamped.
Based on the wellposedness 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
(2.8) 
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 threefield 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 velocitydisplacement space we will also investigate the (blockdiagonal) preconditioner
(2.9) 
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 BeaversJosephSaffman 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 lowestorder TaylorHood 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 
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 Wellposedness of the BiotStokes system
Let us group the variables as and and introduce the weighted norm
(3.1) 
with
(3.2a)  
(3.2b)  
(3.2c) 
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.
Proof.
The operator norm is defined as
and condition (3.3a) states the continuity of . To show this, we first use the CauchySchwarz inequality to derive
(3.4) 
It therefore remains to show that is continuous. Another application of the CauchySchwarz 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 BanachNečasBabuška (BNB) theorem (see, e.g., Ern and Guermond (2004)). In particular, we aim to prove that
(3.5) 
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
(3.6) 
Next, the infsup condition proven in Lemma 3.1, below, allows us to construct such that
(3.7) 
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 CauchySchwarz and Young’s inequality. Setting gives us
(3.8) 
Lemma 3.1.
There exists a such that
Proof.
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 wellposedness 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 Stokesextension of into by setting up another auxiliary Stokes problem (and still in the Biot domain): Find a pair
Comments
There are no comments yet.