Derivation of a poroelastic flexural shell model A.M. was partially supported by the Programme Inter Carnot Fraunhofer from BMBF (Grant 01SF0804) and ANR.

# Derivation of a poroelastic flexural shell model††thanks: A.M. was partially supported by the Programme Inter Carnot Fraunhofer from BMBF (Grant 01SF0804) and ANR.

Andro Mikelić
Université de Lyon, CNRS UMR 5208,
Université Lyon 1, Institut Camille Jordan,
43, blvd. du 11 novembre 1918, 69622 Villeurbanne Cedex, France
E-mail: Andro.Mikelic@univ-lyon1.fr
Josip Tambača
Department of Mathematics
University of Zagreb
Bijenička 30, 10000 Zagreb, Croatia
August 2, 2019
###### Abstract

In this paper we investigate the limit behavior of the solution to quasi-static Biot’s equations in thin poroelastic flexural shells as the thickness of the shell tends to zero and extend the results obtained for the poroelastic plate by Marciniak-Czochra and Mikelić in [16]. We choose Terzaghi’s time corresponding to the shell thickness and obtain the strong convergence of the three-dimensional solid displacement, fluid pressure and total poroelastic stress to the solution of the new class of shell equations.

The derived bending equation is coupled with the pressure equation and it contains the bending moment due to the variation in pore pressure across the shell thickness. The effective pressure equation is parabolic only in the normal direction. As additional terms it contains the time derivative of the middle-surface flexural strain.

Derivation of the model presents an extension of the results on the derivation of classical linear elastic shells by Ciarlet and collaborators to the poroelastic shells case. The new technical points include determination of the strain matrix, independent of the vertical direction, in the limit of the rescaled strains and identification of the pressure equation. This term is not necessary to be determined in order to derive the classical flexural shell model.

Keywords: Thin poroelastic shell, Biot’s quasi-static equations, bending-flow coupling, higher order degenerate elliptic-parabolic systems, asymptotic methods.

AMS classcode MSC 35B25, MSC 74F10, MSC 74K25, MSC 74Q15, MSC 76S.

## 1 Introduction

A shell is a three dimensional body, defined by its middle surface and a neighborhood of a small dimension (the thickness) along the normals to it. The shell is said to be thin when the thickness is much smaller then the minimum of its two radii of the curvature and of the characteristic length of the middle surface .

The basic engineering theory for the bending of thin shells is known as Kirchoff-Love theory or Love’s first approximation. The equations were derived by the so called ”direct method” (see [18] and references therein) and not from the three dimensional equations. A derivation from the three dimensional, at the rigor of the continuum mechanics, is due to Novozhilov and we refer again to [18] for both linear and nonlinear models.

A different approach to deriving the shell equations is to suppose that the middle surface is given as , where be an open bounded and simply connected set with Lipschitz-continuous boundary and is a smooth injective immersion (that is and matrix is of rank two). The vectors , , are linearly independent for all and form the covariant basis of the tangent plane to the -surface .

Then the reference configuration of the shell is of the form , , where and

 r=r(y,x3)=X(y)+x3a3(y),a3(y)=a1(y)×a2(y)|a1(y)×a2(y)|. (1.1)

The associated equations of linearized three dimensional elasticity are then written in curvilinear coordinates with respect to . Their solutions represent the covariant components of the displacement field in the reference configuration .

Then Ciarlet and collaborators developed the asymptotic analysis approach where the normal direction variable was scaled by setting . This change of variables transforms the PDE to a singular perturbation problem in curvilinear coordinates on a fixed cylindrical domain. With such approach Ciarlet and collaborators have established the norm closeness between the solution of the original three dimensional elasticity equations and the Kirchhoff-Love two dimensional flexural and membrane shell equations, in the limit as For details, we refer to the articles [9] and [10] and to the books [6] and [7]. For the complete asymptotic expansion we refer to the review paper [13] by Dauge et al. Further generalizations to nonlinear shells exist and were obtained using convergence. We limit our discussion to the linear shells.

In the everyday life we frequently meet shells (and other low dimensional bodies) which are saturated by a fluid. Many living tissues are fluid-saturated thin bodies like bones, bladders, arteries and diaphragms and they are interpreted as poroelastic plates or shells. Furthermore, industrial filters are an example of poroelastic plates and shells.

Our goal is to extend the above mentioned theory to the poroelastic shells. In the case of the poroelastic plates, derivation of the mathematical model was undertaken in [16]. As in the case of plates, these are the shells consisting of an elastic skeleton (the solid phase) and pores saturated by a viscous fluid (the fluid phase). Interaction between the two phases leads to an overall or effective behavior described by the poroelasticity equations instead of the Navier elasticity equations, coupled with the mass conservation equation for the pressure field. The equations form Biot’s poroelasticity PDEs and can be found in [2], [3] and in the selection of Biot’s publications [23].

The effective linear Biot’s model corresponds to the homogenization of the complicated pore level fluid-structure interaction problem based on the continuum mechanics first principles, i.e. the Navier equations for the solid structure and by the Navier-Stokes equations for the flow. Small deformations are supposed and the interface between phases is linearized. The small parameter of the problem is the ratio between characteristic pore size and the domain size. If, in addition, we consider a periodic porous medium with connected fluid and solid phases, then the two-scale poroelasticity equations can be obtained using formal two-scale expansions in the small parameter, applied to the pore level fluid-structure equations. For details we refer to the book [20], the review [1] and the references therein.

Convergence of the homogenization process for a given frequency was obtained in [19], using the two-scale convergence technique. Convergence in space and time variables was proved by Mikelić and collaborators in [11] and [14]. The upscaling result was presented in detail in [16] and we avoid to repeat it here. We point out that the upscaled model depends on the particular time scale known as Terzaghi’s time . It is equal to the ratio between the viscosity times the characteristic length squared and the shear modulus of the solid structure multiplied by the permeability. If the characteristic time is much longer than than flow dominates vibrations and the acceleration and memory effects can be neglected. The model is then the quasi-static Biot model. For its derivation from the first principles using homogenization techniques see [17].

For the direct continuum mechanics approach to Biot’s equations, we refer to the monograph of Coussy [12]. Using a direct approach, a model for a spherical poroelastic shell is proposed in [22].

In this paper we follow the approach of Ciarlet, Lods and Miara, as presented in the textbooks [6] and [7], and rigorously develop equations for a poroelastic flexural shell.

Successful recent approaches to the derivation of linear and nonlinear shell models use the elastic energy functional. In our situation, presence of the flow makes the problem quasi-static, that is time-dependent, and non-symmetric. The equations for the effective solid skeleton displacement contain the pressure gradient and have the structure of a generalized Stokes system, with the velocity field replaced by displacement. Mass conservation equation is parabolic in the pressure and contains the time derivative of the volumetric strain.

We recall that the quasi-static Biot system is well-posed only if there is a relationship between Biot parameters multiplying the pressure gradient in the displacement system and the time derivative of the divergence of the displacement in the pressure equation. We were able to obtain in [16] the corresponding energy estimate. Similar estimates, for the equations in the curvilinear coordinates, will be obtained here.

In addition, there exists a major difference, with respect to the limit of the normalized term, compared to a classical derivation of the Kirchhoff-Love shell. In our poroelastic case, the limit also contains the pressure field. Furthermore, the pressure oscillations persist and we prove the regularity and uniqueness for the limit problem. As in the poroelastic plate case, it has a richer structure than the classical bending equation. We expect more complex time behavior in this model.

## 2 Geometry of Shells and Setting of the Problem

We study the deformation and the flow in a poroelastic shell , , where the injective mapping is given by (1.1), for and , diam . We recall the middle surface is the image by a smooth injective immersion of an open bounded and simply connected set , with Lipschitz-continuous boundary . We use the linearly independent vectors , , to form a covariant basis of the tangent plane to the -surface . The contravariant basis of the same plane is given by the vectors defined by We extend these bases to the basis of the whole space by the vector given in (1.1) (). Now we collect the local contravariant and covariant bases into the matrix functions

 (2.1)

The first fundamental form of the surface , or the metric tensor, in covariant or contravariant components are given respectively by

 aαβ=aα⋅aβ,aαβ=aα⋅aβ,α,β=1,2.

Note here that because of continuity of and compactness of , there are constants such that

 mcx⋅x≤Ac(y)x⋅x≤Mcx⋅x,x∈R3,y∈¯¯¯ωL. (2.2)

These estimates, with different constants, hold for as well, as it is the inverse of . The second fundamental form of the surface , also known as the curvature tensor, in covariant or mixed components are given respectively by

 bαβ=a3⋅∂βaα=−∂βa3⋅aα,bβα=2∑κ=1aβκbκα,α,β=1,2.

The Christoffel symbols are defined by

 Γκαβ=aκ⋅∂βaα=−∂βaκ⋅aα,α,β,κ=1,2.

We will sometime use for . The area element along is , where . By (2.2) it is uniformly positive, i.e., there is such that

 0

We also need the covariant derivatives which are defined by

 bκβ|α=∂αbκβ+2∑τ=1(Γκατbτβ−Γτβαbκτ),α,β,κ=1,2. (2.4)

In order to describe our results we also need the following differential operators:

 γαβ(\bf v)=12(∂αvβ+∂βvα)−2∑κ=1Γκαβvκ−bαβv3,α,β=1,2, (2.5) ραβ(\bf v)=∂αβv3−2∑κ=1Γκαβ∂κv3+2∑κ=1bκβ(∂αvκ−2∑τ=1Γτακvτ)+2∑κ=1bκα(∂βvκ−2∑τ=1Γτβκvτ) +2∑κ=1bκα|βvκ−2∑κ=1bκαbκβv3,α,β=1,2, (2.6) nαβ|β=∂βnαβ+2∑κ=1Γαβκnβκ+2∑κ=1Γββκnακ,α,β=1,2, nαβ|αβ=∂α(nαβ|β)+2∑κ=1Γκακnαβ|β,α,β=1,2,

defined for smooth vector fields and tensor fields n. The upper face (respectively lower face) of the shell is (respectively . is the lateral boundary, . We recall that the small parameter is the ratio between the shell thickness and the characteristic horizontal length is .

We note that Biot’s diphasic equations describe behavior of the system at so called Terzaghi’s time scale where is the characteristic domain size, is dynamic viscosity, is permeability and is the shear modulus. For the list of all parameters see Table 1.

Similarly as in [16], we chose as the characteristic length , which leads to the Taber-Terzaghi transversal time . Another possibility was to choose the longitudinal time scaling with . It would lead to different scaling in (2.9) and the dimensionless permeability coefficient in (3.3) would not be but . In the context of thermoelasticity, one has the same equations and Blanchard and Francfort rigorously derived in [5] the corresponding thermoelastic plate equations. We note that considering the longitudinal time scale yields the effective model where the pressure (i.e. the temperature in thermoelasticity) is decoupled from the flexion.

Then the quasi-static Biot equations for the poroelastic body take the following dimensional form:

 ~σ=2μe(~u)+(λ%div~u−α~p)I in ~ΩℓL, (2.7) − div ~σ=−μ△~u−(λ+μ)▽div ~u+α▽~p=0 in ~ΩℓL, (2.8) ∂∂t(βG~p+α % div ~u)−kη△~p=0 in % ~ΩℓL. (2.9)

Note that and is the stress tensor. All other quantities are defined in Table 1.

We impose a given contact force and a given normal flux at . At the lateral boundary we impose a zero displacement and a zero normal flux. Here is the outer unit normal at the boundary. At initial time we prescribe the initial pressure .

Our goal is to extend the Kirchhoff-Love shell justification by Ciarlet, Lods et al and by Dauge et al to the poroelastic case. We announce briefly the differential equations of the flexural poroelastic shell in dimensional form. Note that our mathematical result will be in the variational form and that differential form is only formal and written for reader’s comfort.

Effective dimensional equations:

The model is given in terms of which is the vector of components of the displacement of the middle surface of the shell in the contravariant basis and which is the pressure in the 3D shell. Let us denote the bending moment (contact couple) due to the variation in pore pressure across the plate thickness by

 \bf m=ℓ312~Cc(Ac\boldmathρ(\bf ueff))Ac+2μαλ+2μ∫ℓ/2−ℓ/2y3peffdy3Ac, (2.10)

where is given by (2.6) and is the elasticity tensor, usually appearing in the classical shell theories, given by

 ~CcE=2μλλ+2μtr(E)I+2μE,E∈M3×3sym.

Then the model in the differential formulation reads as follows:

 2∑β=1(−(nαβ+2∑κ=1bακmκβ)|β−2∑κ=1bακ(mκβ|β))=(P+ℓL)α+(P−ℓL)α in ωL,α=1,2, (2.11) 2∑α,β=1(mαβ|αβ−2∑κ=1bκαbκβmαβ−bαβnαβ)=(P+ℓL)3+(P−ℓL)3 in ωL, 12(∂αueffβ+∂βueffα)−2∑κ=1Γκαβueffκ−bαβueff3=0 in ωL,α,β=1,2, ueffi=0,i=1,2,3,∂ueff3∂\boldmathν=0 on ∂ωL,for % every t∈(0,T),
 (βG+α2λ+2μ)∂peff∂t−α2μλ+2μAc:\boldmathρ(∂\bf ueff∂t)y3−kη∂2peff∂(y3)2=0 (2.12) in (0,T)×ωL×(−ℓ/2,ℓ/2), kη∂peff∂y3=−VL, on (0,T)×ωL×({−ℓ/2}∪{ℓ/2}), peff=pℓL,ingiven at t=0.

Here are components of the contact force at in the covariant basis, , . Thus, the poroelastic flexural shell model in the differential formulation is given for unknowns and by equations (2.10), (2.11) and (2.12). The components of n are the contact forces, linked to the constraint , and being the Lagrange multipliers in the problem. The components of m are the contact couples. The first two equations in (2.11) can be found in the differential equation of the Koiter shell model (see [7, Theorem 7.1-3]). The third equation is the restriction of approximate inextensibility of the shell. The first equation in (2.12) is the evolution equation for the effective pressure with associated boundary and initial conditions in the remaining part of (2.12). Note also that the same model holds for the shell clamped only on the portion of the boundary, i.e., the boundary condition in the fourth equation in (2.11) holds for a subset of with positive measure.

In subsection 3.1 we present the dimensionless form of the problem. In subsection 3.2 we recall existence and uniqueness result of the smooth solution for the starting problem. Subsection 3.3 is consecrated to the introduction of the problem in curvilinear coordinates and the rescaled problem, posed on the domain . In subsection 3.4 we formulate the main convergence results. In Section 4 we study the a priori estimates for the family of solutions. Then in Section 5 we study convergence of the solutions to the rescaled problem, as . In Appendix we give properties of the metric and curvature tensors.

## 3 Problem setting in curvilinear coordinates and the main results

### 3.1 Dimensionless equations

We introduce the dimensionless unknowns and variable by setting

 β=βGμ;P=μUL;U~uε=~u;T=ηℓ2kμ;~λ=λμ; P~pε=~p;~yL=y;~x3L=x3;~rL=r;~XL=X;~tT=t;~σεμUL=~σ.

After dropping wiggles in the coordinates and in the time, the system (2.7)–(2.9) becomes

 −div~σε=−△~uε−~λ▽div ~uε+α▽~pε=0 in (0,T)×~Ωε, (3.1) ~σε=2e(~uε)+(~λ div ~uε−α~pε)I in (0,T)×~Ωε, (3.2) ∂∂t(β~pε+α div ~uε)−ε2△~pε=0 in (0,T)×~Ωε, (3.3)

where denotes the dimensionless displacement field and the dimensionless pressure. We study a shell with thickness and section . It is described by

 ~Ωε=r({(x1,x2,x3)∈ω×(−ε/2,ε/2)})=r(Ωε),

(respectively ) is the upper face (respectively the lower face) of the shell . is the lateral boundary, .

We suppose that a given dimensionless traction force is applied on and impose the shell is clamped on :

 ~σε\boldmathν=(2e(~uε)−α~pεI+~λ(div ~uε)I)% \boldmathν=ε3~P± on ~Σε+∪~Σε−, (3.4) ~uε=0, on ~Γε. (3.5)

For the pressure , at the lateral boundary we impose zero inflow/outflow flux:

 −▽~pε⋅\boldmathν=0. (3.6)

and at , we set

 −▽~pε⋅%\boldmath$ν$=±~V. (3.7)

Finally, we need an initial condition for at ,

 ~pε(x1,x2,x3,0)=ε~pin(x1,x2) in ~Ωε. (3.8)

Let . Then the weak formulation corresponding to (3.1)–(3.8) is given by

Find , such that it holds

 ∫~Ωε2 e(~uε):e(~v) dx+~λ∫~Ωε div ~uε div ~v dx−α∫~Ωε~pε div ~v dx (3.9) β∫~Ωε∂t~pε~q dx+∫~Ωεα%div∂t~uε~q dx+ε2∫~Ωε∇~pε⋅∇~q dx =ε2∫~Σε−~V~q ds−ε2∫~Σε+~V~q ds, for every ~q∈H1(~Ωε) and t∈(0,T), (3.10) ~pε|{t=0}=ε~pin, in ~Ωε. (3.11)

Note that for two matrices and the Frobenius scalar product is denoted by .

### 3.2 Existence and uniqueness for the ε-problem

In this subsection we recall the existence and uniqueness of a solution of the problem (3.9)-(3.11). We follow [16] and get

###### Proposition 1.

Let us suppose

 ~pin∈H20(~Ωε),P±∈H2(0,T;L2(ω;R3))and~V∈H1(0,T;L2(ω)),~V|{t=0}=0. (3.12)

Then problem (3.9)–(3.11) has a unique solution

### 3.3 Problem in Curvilinear Coordinates and the Scaled Problem

Our goal is to find the limits of the solutions of problem (3.9)–(3.11) when tends to zero. It is known from similar considerations made for classical shells that asymptotic behavior of the longitudinal and transverse displacements of the elastic body is different. The same effect is expected in the present setting. Therefore we need to consider asymptotic behavior of the local components of the displacement . It can be done in many ways, but in order to preserve some important properties of bilinear forms, such as positive definiteness and symmetry, we rewrite the equations in curvilinear coordinates defined by . Then we formulate equivalent problems posed on the domain independent of .

The covariant basis of the shell , which is the three-dimensional manifold parameterized by , is defined by

 gεi=∂ir:Ωε→R3,i=1,2,3.

Vectors are given by

 gε1=a1(y)+x3∂y1a3(y), gε2=a2(y)+x3∂y2a3(y), gε3=a3(y).

Vectors satisfying

 gj,ε⋅gεi=δij % on ¯¯¯¯Ωε,i,j=1,2,3,

where is the Kronecker symbol, form the contravariant basis on . The contravariant metric tensor , the covariant metric tensor and the Christoffel symbols of the shell are defined by

 gij,ε=gi,ε⋅gj,ε,gεij=gεi⋅gεj,Γi,εjk=gi,ε⋅∂jgεk on ¯¯¯¯Ωε,i,j,k=1,2,3.

We set

 Γi,ε=(Γi,εjk)j,k=1,…,3and~\boldmathγε(v)=12(∇v+∇vT)−3∑i=1viΓi,ε. (3.13)

Let . Until now we were using the canonical basis , for . Now the displacement is rewritten in the contravariant basis,

while for scalar fields we just change the coordinates

 ~pε∘r=pε,~q∘r=q,~V∘r=V,~pin∘r=pin,

on . The contact forces are rewritten in the covariant basis of the shell

 ~P±∘r=3∑i=1(P±)igεi on Σε±.

New vector functions are defined by

 uε=uεiei,v=viei,P±=(P±)iei.

Note that are not components of the physical displacement. They are just intermediate functions which will be used to reconstruct . The corresponding function space to is the space

 V(Ωε)={v∈H1(Ωε)3:v|Γε=0}.

Let and let

 CE=~λ(trE)I+2E,for allE∈M3×3sym. (3.14)

Then the problem (3.9)–(3.11) can be written by

 ∫ΩεC(Qε~\boldmathγε(uε)(Qε)T):(Qε~\boldmathγε(v)(Qε)T)√gεdy−α∫Ωεpεtr(Qε~\boldmathγε(v)(Qε)T)√gεdy (3.15) =ε3∫Σε+P+⋅v√gεds+ε3∫Σε−P−⋅v√gεds,\bf v∈V(Ωε), a.e. t∈[0,T], ∫Ωεβ∂pε∂tq√gεdy+∫Ωεα∂∂ttr(Qε~\boldmathγε(uε)(Qε)T)q√gεdy +ε2∫ΩεQε∇pε⋅Qε∇q√gεdy=ε2∫Σε−Vq√gεds−ε2∫Σε+Vq√gεds, q∈H1(Ωε), a.e. t∈[0,T], pε=εpin, for t=0.

This is the problem in curvilinear coordinates.

Problems for all and are posed on –dependent domains. In the sequel we follow the idea from Ciarlet, Destuynder [8] and rewrite (3.15) on the canonical domain independent of . As a consequence, the coefficients of the resulting weak formulation will depend on explicitly.

Let and let be defined by

 Rε(z)=(z1,z2,εz3),z∈Ω, ε∈(0,ε0).

By we denote the upper and lower face of . Let . To the functions , , , , , , , defined on we associate the functions , , , , , , , defined on by composition with . Let us also define

 V(Ω)={v=(v1,v2,v3)∈H1(Ω;R3):v|Γ=0}.

Then the problem (3.15) can be written as

 ε∫ΩC(Q(ε)\boldmathγε(u(ε))Q(ε)T):(Q(ε)% \boldmathγε(v)Q(ε)T)√g(ε)dz (3.16) −εα∫Ωp(ε)tr(Q(ε)\boldmathγ% ε(v)Q(ε)T)√g(ε)dz =ε3∫Σ+P+⋅v√g(ε)ds+ε3∫Σ−P−⋅v√g(ε)ds,v∈V(Ω),%a.e.t∈[0,T], ε∫Ωβ∂p(ε)∂tq√g(ε)dz+ε∫Ωα∂∂ttr(Q(ε)%\boldmath$γ$ε(u(ε))Q(ε)T)q√g(ε)dz +ε3∫ΩQ(ε)∇εp(ε)⋅Q(ε)∇εq√g(ε)dz =ε2∫Σ−Vq√g(ε)ds−ε2∫Σ+Vq√g(ε)ds,q∈H1(Ω),</