A nonlinear Stokes-Biot model for the interaction of a non-Newtonian fluid with poroelastic media I: well-posedness of the model.

# A nonlinear Stokes-Biot model for the interaction of a non-Newtonian fluid with poroelastic media I: well-posedness of the model.

Ilona Ambartsumyan Department of Mathematics, University of Pittsburgh, Pittsburgh, PA 15260, USA; {ila6@pitt.edu, tqn4@pitt.edu, yotov@math.pitt.edu}; partially supported by NSF grant DMS 1418947 and DOE grant DE-FG02-04ER25618.    Vincent J. Ervin22footnotemark: 2  Department of Mathematical Sciences, Clemson University, Clemson, SC 29634-0975, USA; {vjervin@clemson.edu}.    Truong Nguyen11footnotemark: 1     Ivan Yotov11footnotemark: 1
July 14, 2019
###### Abstract

We propose and analyze a model for solving the coupled problem arising in the interaction of a free fluid with a poroelastic medium. The flow in the fluid region is described by the Stokes equations and in the poroelastic medium by the quasi-static Biot model. The focus of the model is on quasi-Newtonian fluids that exhibit a shear-thinning property. We establish existence and uniqueness of the solution for two alternative formulations of the proposed model.

## 1 Introduction

The interaction of a free fluid with a deformable porous medium is a challenging multiphysics problem that has a wide range of applications, including processes arising in gas and oil extraction from naturally or hydraulically fractured reservoirs, designing industrial filters, and blood-vessel interactions. The free fluid region can be modeled by the Stokes or the Navier-Stokes equations, while the flow through the deformable porous medium is modeled by the quasi-static Biot system of poroelasticity [4]. The two regions are coupled via dynamic and kinematic interface conditions, including balance of forces, continuity of normal velocity, and a no slip or slip with friction tangential velocity condition. These multiphysics models exhibit features of coupled Stokes-Darcy flows and fluid-structure interaction (FSI). There is extensive literature on modeling these separate couplings, see e.g. [15, 25, 33] for Stokes-Darcy flows and [20, 19, 22] for FSI. More recently there has been growing interest in modeling Stokes-Biot couplings, which can be referred to as fluid-poroelastic structure interatction (FPSI). The well-posedness of the mathematical model is studied in [36]. A variational multiscale stabilized finite element method for the Navier-Stokes-Biot problem is developed in [2]. In [9] a non-iterative operator-splitting method is developed for the Navier-Stokes-Biot model with pressure Darcy formulation. The well posedness of a related model is studied in [12]. The Stokes-Biot problem with a mixed Darcy formulation is studied in [8] and [1] using Nitsche’s method and a Lagrange multiplier, respectively, to impose the continuity of normal velocity on the interface. An optimization-based iterative algorithm with Neumann control is proposed in [13]. A reduced-dimension fracture model coupling Biot and an averaged Brinkman equation is developed in [10]. Alternative fracture models using the Reynolds lubrication equation coupled with Biot have also been studied, see e.g. [23, 28].

All of the above mentioned works are based on Newtonian fluids. In this paper we develop FPSI with non-Newtonian fluids, which, to the best of our knowledge, has not been studied in the literature. In many applications the fluid exhibits properties that cannot be captured by a Newtonian fluid assumption. For instance, during water flooding in oil extraction, polymeric solutions are often added to the aqueous phase to increase its viscosity, resulting in a more stable displacement of oil by the injected water [27]. In hydraulic fracturing, proppant particles are mixed with polymers to maintain high permeability of the fractured media [26]. In blood flow simulations of small vessels or for patients with a cardiovascular disease, where the arterial geometry has been altered to include regions of recirculation, one needs to consider models that can capture the sheer-thinning property of the blood [24].

In this work we use nonlinear Stokes equations to model the free fluid in the flow region and the Biot model for the fluid in the poroelastic region. The later is based on a linear stress-strain constitutive relationship and a nonlinear Darcy flow. The coupling conditions between the two subdomians include mass conservation, conservation of momentum and the Beavers-Joseph-Saffman slip with friction condition. We focus on fluids that possess the sheer thinning property, i.e., the viscosity decreases under shear strain, which is typical for polymer solutions and blood. Viscosity models for such non-Newtonian fluids include the Power law, the Cross model and the Carreau model [14, 30, 27, 31]. The Power law model is popular because it only contains two parameters, and it is possible to derive analytical solutions in various flow conditions [5]. On the other hand, it implies that in the flow region the viscosity goes to infinity if the deformation goes to zero, which may not be representative in certain applications and also significantly complicates the theoretical aspects of the problem. The Cross and Carreau models were deduced empirically as improvements of the Power law model. In both of these models the viscosity is strictly greater than zero and bounded, but knowledge of three parameters is required. We assume that the viscosity in each subdomain satisfies one such model, with dependence on the magnitude of the deformation tensor and the magnitude of Darcy velocity in the fluid and poroelastic regions, respectively. We further assume that along the interface the fluid viscosity is a function of the fluid and structure interface velocities.

Our model is built on the models presented in [17] and [1]. In [17], the coupling of generalized nonlinear Stokes and Darcy equations is studied, while a linear Stokes-Biot model is considered in [1]. We allow for both unbounded viscosity models, such as the Power law, as well as bounded models such as the Cross and Carreau models. In the former case, the analysis is done in an appropriate Sobolev space setting, using spaces such as , where is the viscosity shear thinning parameter. In the latter case, the analysis reduces to the Hilbert space setting. Nonlinear Stokes-Darcy models with bounded viscosity have been studied in [18, 16, 11], while the unbounded case is considered in [17].

Following the approach in [1], we enforce the continuity of normal velocity on the interface through the use of a Lagrange multiplier. The resulting weak formulation is a nonlinear time-dependent system, which is difficult to analyze, due to to the presence of the time derivative of the diplacement in some non-coercive terms. We consider an alternative mixed elasticity formulation with the structure velocity and elastic stress as primary variables, see also [36]. In this case we obtain a system with a degenerate evolution in time operator and a nonlinear saddle-point type spatial operator. The structure of the problem is similar to the one analyzed in [37], see also [6] in the linear case. However, the analysis in [37] is restricted to the Hilbert space setting and needs to be extended to the Sobolev space setting. Furthermore, the analysis in [37] is for monotone operators, see [38], and as a result requires certain right hand side terms to be zero, while in typical applications these terms may not be zero. Here we explore the coercivity of the operators to reformulate the problem as a parabolic-type system for the pressure and stress in the poroelastic region. We show well posedness for this system for general source terms and that the solution satisfies the original formulation. We also prove that the solution to the original formulation is unique and provide a stability bound. In part two of this sequence of papers, we consider a semidiscrete finite element approximation of the system, carry out stability and error analysis, and present numerical experiments.

The rest of the paper is organized as follows. In Section 2 we introduce the governing equations. Section 3 is devoted to the weak formulation, upon which we base the numerical method, and an alternative formulation, which is needed for the purpose of the analysis. In Section 4 we prove the well-posedness of the alternative and original formulations.

## 2 Problem set up

Let be a Lipschitz domain, which is subdivided into two non-overlapping and possibly non-connected regions: fluid region and poroelastic region . Let denote the (nonempty) interface between these regions and let and denote the external parts of the boundary . We denote by and he unit normal vectors which point outward from and , respectively, noting that on . Let be the velocity-pressure pairs in , , , and let be the displacement in .

We assume that the flow in is governed by the nonlinear generalized Stokes equations with homogeneous boundary conditions on :

 −∇⋅{\boldmathσ}f(uf,pf)=ff,∇⋅uf=qfin Ωf×(0,T],uf=0on Γf×(0,T], (2.1)

where and denote the deformation and the stress tensors, respectively:

 D(uf)=12(∇uf+∇uTf),{\boldmathσ}f(uf,pf)=−pfI+2ν(D(uf))D(uf),

where stands for the identity operator. We consider a generalized Newtonian fluid with the viscosity dependent on the magnitude of the deformation tensor, in particular shear-thinning fluids with a decreasing function of . We consider the following models [14, 30], where , , and are constants:

Carreau model:

 ν(D(uf))=ν∞+(ν0−ν∞)/(1+Kf|D(uf)|2)(2−r)/2, (2.2)

Cross model:

 ν(D(uf))=ν∞+(ν0−ν∞)/(1+Kf|D(uf)|2−r), (2.3)

Power law model:

 ν(D(uf))=Kf|D(uf)|r−2. (2.4)

In turn, in we consider the quasi-static Biot system [4]

 −∇⋅{\boldmathσ}p({\boldmathη% }p,pp)=fpin Ωp×(0,T], (2.5) νeff(up)κ−1up+∇pp=0,∂∂t(s0pp+αp∇⋅{\boldmathη}p)+∇⋅up=qpin Ωp×(0,T], (2.6) up⋅np=0 on ΓNp×(0,T],pp=0 on ΓDp×(0,T],{\boldmath% η}p=0 on Γp×(0,T], (2.7)

where and are the elasticity and poroelasticity stress tensors, respectively,

 {\boldmathσ}e({\boldmathη}p)=λp(∇⋅{\boldmathη}p)I+2μpD({% \boldmathη}p),{\boldmathσ}p({\boldmathη}p,pp)={\boldmathσ}e({\boldmathη}p)−αpppI, (2.8)

is the Biot-Willis constant, , are the Lamè coefficients, is a storage coefficient, is a scalar uniformly positive and bounded permeability function, and . To avoid the issue with restricting the mean value of the pressure, we assume that . We further assume that . We note that even though the analysis of our formulation is valid for a symmetric and positive definite permeability tensor, we restrict it to , due to assumptions made in the derivations of some of the viscosity functions suitable for modeling non-Newtonian flow in porous media. In particular, we consider the following two models for the effective viscosity in [27, 31], where , , and are constants:

Cross model:

 νeff(up)=ν∞+(ν0−ν∞)/(1+Kp|up|2−r), (2.9)

Power law model:

 νeff(up)=Kp(|up|/(√κmc))r−2, (2.10)

where is a constant that depends on the internal structure of the porous media.

Following [36, 2], the interface conditions on the fluid-poroelasticity interface , are mass conservation, balance of normal stress, the Beavers-Joseph-Saffman (BJS) slip with friction condition [3, 34], and conservation of momentum:

 uf⋅nf+(∂{\boldmathη}p∂t+up)⋅np=0on Γfp, (2.11) −({\boldmathσ}fnf)⋅nf=ppon Γfp, (2.12) −({\boldmathσ}fnf)⋅tf,j=νIαBJS√κ−1(uf−∂{% \boldmathη}p∂t)⋅tf,jon Γfp, (2.13) {\boldmathσ}fnf=−{\boldmathσ}pnpon Γfp, (2.14)

where , , is an orthogonal system of unit tangent vectors on and is an experimentally determined friction coefficient. We note that the continuity of flux takes into account the normal velocity of the solid skeleton, while the BJS condition accounts for its tangential velocity. We assume that along the interface the fluid viscosity is a function of the magnitude of the tangential component of the slip velocity given by the Cross model (2.9) or the Power law model (2.10), where . For the rest of the paper we will write , or keeping in mind that these are nonlinear functions as defined above.

The above system of equations is complemented by a set of initial conditions:

 pp(0,x)=pp,0(x),{\boldmathη}p(0,x)={\boldmathη}p,0(x) in Ωp.

In the following, we make use of the usual notation for Lebesgue spaces , Sobolev spaces and Hilbert spaces . For a set , the inner product is denoted by for scalar, vector and tensor valued functions. For a section of a subdomain boundary we write for the inner product (or duality pairing). We also denote by a generic positive constant independent of the discretization parameters.

Adopting the approach from [17, 18], we assume that the viscosity functions satisfy one of the two sets of assumptions (A1)–(A2) or (B1)–(B2) below. Let and let be given by . For , let satisfy, for constants and ,

 (G(x+h)−G(x))⋅h ≥C1|h|2, (A1) |G(x+h)−G(x)| ≤C2|h|, (A2)

or

 (G(x+h)−G(x))⋅h ≥C3|h|2c+|x|2−r+|x+h|2−r, (B1) |G(x+h)−G(x)| ≤C4|h|c+|x|2−r+|x+h|2−r, (B2)

with the convention that if , and if and . From (B1)–(B2) it follows that there exist constants such that for [35]

 (G(s)−G(t),s−t)G≥C5⎛⎝∫Ω|G(s)−G(t)||s−t|dx+∥s−t∥2Lr(G)c+∥s∥2−rLr(G)+∥t∥2−rLr(G)⎞⎠, (2.15) (G(s)−G(t),w)G≤C6∥∥∥|s−t|c+|s|+|t|∥∥∥2−rrL∞(G)(|G(s)−G(t)|,|s−t|)1/r′G∥w∥Lr(G). (2.16)
###### Remark 2.1.

It is shown in [16] that conditions (A1)–(A2) are satisfied for given in the Carreau model (2.2) with , in which case . A similar argument can be applied to show that (A1)–(A2) hold for the Cross model, with given in (2.3) for Stokes and given in (2.9) for Darcy, in the case of . Furthurmore, it is shown in [35] that conditions (B1)–(B2) with hold in the case of the Carreau model (2.2) with , and that conditions (B1)–(B2) with hold for the Power law model (2.4) and (2.10).

## 3 Variational formulation

We will consider two cases when defining the functional spaces, depending on which set of assumptions holds. In the case (B1)–(B2), we consider Sobolev spaces. For a given let be its conjugate, satisfying . Let

 Vf ={vf∈W1,r(Ωf)d:vf=0 on Γf}, Wf=Lr′(Ωf), (3.1)

with the corresponding norms

 ∥vf∥Vf=∥vf∥(W1,r(Ωf))d,∥wf∥Wf=∥wf∥Lr′(Ωf).

With , let

 Vp ={vp∈Lr(div;Ωp):vp⋅np=0 on ΓNp}, Wp=Lr′(Ωp), Xp ={{\boldmathξ}p∈H1(Ωp)d:{% \boldmathξ}p=0 on Γp}. (3.2)

with norms

 ∥vp∥rVp =∥vp∥r(Lr(Ωp))d+∥∇⋅vp∥r(Lr(Ωp)), ∥wp∥Wp=∥wp∥Lr′(Ωp), ∥{\boldmathη}p∥Xp =∥{\boldmathη}p∥(H1(Ωp))d.

In the case of (A1)–(A2), we consider Hilbert spaces, with the above definitions replaced by

 Vf ={vf∈H1(Ωf)d:vf=0% on Γf}, Wf=L2(Ωf), (3.3) Vp ={vp∈H(div;Ωp):vp⋅np=0 on ΓNp}, Wp=L2(Ωp). (3.4)

The global spaces are products of the subdomain spaces. For simplicity we assume that each region consists of a single subdomain.

###### Remark 3.1.

For simplicity of the presentation, for the rest of the paper we focus on the case (B1)–(B2), which is the technically more challenging case. The arguments apply directly to the case (A1)–(A2).

### 3.1 Lagrange multiplier formulation

To derive the weak formulation, we multiply (2.1)–(2.6) by appropriate test functions and integrate each equation over the corresponding region, utilizing the boundary and interface conditions (2.11)–(2.14). Integration by parts in the first equation in (2.1), (2.5), and the first equation in (2.6) leads to the Stokes, Darcy and the elasticity functionals

the bilinear forms

 b⋆(⋅,⋅):V⋆×W⋆⟶R,b⋆(v,w):=−(∇⋅v,w)Ω⋆,⋆=f,p,

and the interface term

 IΓfp=−⟨{\boldmathσ}fnf,vf⟩Γfp−⟨{\boldmathσ}pnp,{% \boldmathξ}p⟩Γfp+⟨pp,vp⋅np⟩Γfp.

This term is incorporated into the weak formulation by introducing a Lagrange multiplier which has a meaning of normal stress/Darcy pressure on the interface:

 λ=−({\boldmathσ}fnf)⋅nf=pp,on Γfp.

With introduced, we have, using (2.12), (2.13) and (2.14),

 IΓfp=aBJS(uf,∂t{\boldmath% η}p;vf,{\boldmathξ}p)+bΓ(vf,vp,{\boldmathξ}p;λ),

where

 aBJS(uf,{\boldmathη}p;vf,{\boldmathξ}p) =d−1∑j=1⟨νIαBJS√κ−1(uf−{\boldmathη}p)⋅tf,j,(vf−{\boldmathξ}p)⋅tf,j⟩Γfp, bΓ(vf,vp,{\boldmathξ}p;μ) =⟨vf⋅nf+({\boldmathξ}p+vp)⋅np,μ⟩Γfp.

For the term to be well-defined, we choose the Lagrange multiplier space as . It is shown in [17] that in the case , if , then can be identified with a functional in . Furthermore, for , , and for , . Therefore, with , the integrals in are well-defined.

The variational formulation reads: given , , , , and , , find, for , , such that for all , , , , , and ,

 af(uf,vf)+adp(up,vp)+aep({\boldmathη}p,{\boldmathξ}p)+aBJS(uf,∂t{\boldmathη}p;vf,{% \boldmathξ}p)+bf(vf,pf)+bp(vp,pp) +αpbp({\boldmathξ}p,pp)+bΓ(vf,vp,{\boldmathξ}p;λ)=(ff,vf)Ωf+(fp,{\boldmathξ}p)Ωp, (3.5) (s0∂tpp,wp)Ωp−αpbp(∂t{\boldmathη}p,wp)−bp(up,wp)−bf(uf,wf)=(qf,wf)Ωf+(qp,wp)Ωp, (3.6) bΓ(uf,up,∂t{% \boldmathη}p;μ)=0. (3.7)

Note that is well-defined, since for , we have that and .

Although related models have been analyzed previously, e.g. the non-Newtonian Stokes-Darcy model was investigated in [17] and the Newtonian dynamic Stokes-Biot model was studied in [36], the well posedness of (3.5)–(3.7) has not been established in the literature. Analyzing this formulation directly is difficult, due to the presence of in several non-coercive terms. Instead, we analyze an alternative formulation and show that the two formulations are equivalent.

### 3.2 Alternative formulation

Our goal is to obtain a system of evolutionary saddle point type, which fits the general framework studied in [37]. Following the approach from [36], we do this by considering a mixed elasticity formulation with the structure velocity and elastic stress as primary variables. Recall that the elasticity stress tensor is connected to the displacement through the relation [7]:

 A{\boldmathσ}e=D({\boldmathη}p), (3.8)

where is a symmetric and positive definite compliance tensor. In the isotropic case has the form

 A{\boldmathσ}e=12μp({% \boldmathσ}e−λp2μp+dλptr({\boldmathσ}e)I),   with  A−1{% \boldmathσ}e=2μp{\boldmathσ}e+λptr({\boldmathσ}e)I. (3.9)

The space for the elastic stress is with the norm .

The derivation of the alternative variational formulation differs from the original one in the way the equilibrium equation (2.5) is handled. As before, we multiply it by a test function and integrate by parts. However, instead of using the constitutive relation of the first equation in (2.8), we use only the second equation in (2.8), resulting in

 ∫Ωp({\boldmathσ}e:D(vs)−αppp∇⋅vs)dx−∫Γfp{\boldmathσ}pnp⋅vsds=∫Ωpfp⋅vsdx.

We eliminate the displacement from the system by differentiating (3.8) in time and introducing a new variable , which has a meaning of structure velocity. Multiplication by a test function gives

 ∫Ωp(A∂t{\boldmathσ}e:\boldmathτe−D(us):\boldmathτe)dx=0.

The rest of the equations are handled in the same way as in the original weak formulation, resulting in the same Stokes and Darcy functionals, and , respectively, and the same interface term . Defining the bilinear forms and ,

 bs(vs,{\boldmathτ}e):=(D(vs),{\boldmathτ}e)Ωp,asp({% \boldmathσ}e,{\boldmathτ}e):=(A{\boldmathσ}e,{\boldmathτ}e)Ωp,

we obtain the following weak formulation: given , , , , and , , for , find , such that for all , , , , , , ,

 af(uf,vf)+adp(up,vp)+aBJS(uf,us;vf,vs)+bf(vf,pf)+bp(vp,pp) +αpbp(vs,pp)+bs(vs,{\boldmathσ}e)+bΓ(vf,vp,vs;λ)=(ff,vf)Ωf+(fp,vs)Ωp, (3.10) (s0∂tpp,wp)Ωp+asp(∂t{\boldmathσ}e,{\boldmathτ}e)−αpbp(us,wp)−bp(up,wp)−bs(us,{\boldmathτ}e)−bf(uf,wf) =(qf,wf)Ωf+(qp,wp)Ωp, (3.11) bΓ(uf,up,us;μ)=0. (3.12)

We can write (3.10)–(3.12) in an operator notation as a degenerate evolution problem in a mixed form:

 =f(t) in Q′, (3.13) ∂∂tE2s(t)−Bq(t)+Cs(t) =g(t) in S′, (3.14)

where we define , the space of generalized displacement variables, as

 Q={q=(vp,vs,vf)∈Vp×Xp×Vf},

and, similarly, the space , consisting of generalized stress variables, as

 S={s=(wp,\boldmathτe,wf,μ)∈Wp×{\boldmathΣ}e×Wf×Λ}.

The spaces and are equipped with norms:

 ∥q∥Q =∥vp∥Vp+∥vs∥Xp+∥vf∥Vf, ∥s∥S =∥wp∥Wp+∥\boldmathτe∥{% \boldmathΣ}e+∥wf∥Wf+∥μ∥Λ.

The operators , and the functionals , are defined as follows:

 A=⎛⎜ ⎜⎝νeffκ−1000αBJSγ′tνI√κ−1γt−αBJSγ′tνI√κ−1γt0−αBJSγ′tνI√