Toward permeability of deformable porous media - Fraunhofer

O. Iliev, A. Mikelic, P. Popov

Fluid structure interaction problems in deformable porous media: Toward permeability of deformable porous media

Berichte des Fraunhofer ITWM, Nr. 65 (2004)

© Fraunhofer-Institut für Techno- und Wirtschaftsmathematik ITWM 2004

ISSN 1434-9973

Bericht 65 (2004)

Alle Rechte vorbehalten. Ohne ausdrückliche, schriftliche Gene h mi gung des Herausgebers ist es nicht gestattet, das Buch oder Teile daraus in irgendeiner Form durch Fotokopie, Mikrofi lm oder andere Verfahren zu reproduzieren oder in eine für Maschinen, insbesondere Daten ver ar be i- tungsanlagen, verwendbare Sprache zu übertragen. Dasselbe gilt für das Recht der öffentlichen Wiedergabe.

Warennamen werden ohne Gewährleistung der freien Verwendbarkeit benutzt.

Die Veröffentlichungen in der Berichtsreihe des Fraunhofer ITWM können bezogen werden über:

Fraunhofer-Institut für Techno- undWirtschaftsmathematik ITWMGottlieb-Daimler-Straße, Geb. 49

67663 KaiserslauternGermany

Telefon: +49 (0) 6 31/2 05-32 42Telefax: +49 (0) 6 31/2 05-41 39E-Mail: [email protected]: www.itwm.fraunhofer.de

Vorwort

Das Tätigkeitsfeld des Fraunhofer Instituts für Techno- und Wirt schafts ma the ma tik ITWM um fasst an wen dungs na he Grund la gen for schung, angewandte For schung so wie Be ra tung und kun den spe zi fi sche Lö sun gen auf allen Gebieten, die für Tech- no- und Wirt schafts ma the ma tik be deut sam sind.

In der Reihe »Berichte des Fraunhofer ITWM« soll die Arbeit des Instituts kon ti -nu ier lich ei ner interessierten Öf fent lich keit in Industrie, Wirtschaft und Wis sen -schaft vor ge stellt werden. Durch die enge Verzahnung mit dem Fachbereich Ma- the ma tik der Uni ver si tät Kaiserslautern sowie durch zahlreiche Kooperationen mit in ter na ti o na len Institutionen und Hochschulen in den Bereichen Ausbildung und For schung ist ein gro ßes Potenzial für Forschungsberichte vorhanden. In die Be- richt rei he sollen so wohl hervorragende Di plom- und Projektarbeiten und Dis ser -ta ti o nen als auch For schungs be rich te der Institutsmitarbeiter und In s ti tuts gäs te zu ak tu el len Fragen der Techno- und Wirtschaftsmathematik auf ge nom men werden.

Darüberhinaus bietet die Reihe ein Forum für die Berichterstattung über die zahl- rei chen Ko o pe ra ti ons pro jek te des Instituts mit Partnern aus Industrie und Wirt- schaft.

Berichterstattung heißt hier Dokumentation darüber, wie aktuelle Er geb nis se aus mathematischer For schungs- und Entwicklungsarbeit in industrielle An wen dun gen und Softwareprodukte transferiert wer den, und wie umgekehrt Probleme der Pra- xis neue interessante mathematische Fragestellungen ge ne rie ren.

Prof. Dr. Dieter Prätzel-WoltersInstitutsleiter

Kaiserslautern, im Juni 2001

FLUID STRUCTURE INTERACTION PROBLEMS INDEFORMABLE POROUS MEDIA: TOWARD PERMEABILITY

OF DEFORMABLE POROUS MEDIA

OLEG ILIEV, ANDRO MIKELIC, AND PETER POPOV

Abstract. In this work the problem of fluid flow in deformable porous mediais studied. First, the stationary fluid-structure interaction (FSI) problem isformulated in terms of incompressible Newtonian fluid and a linearized elasticsolid. The flow is assumed to be characterized by very low Reynolds numberand is described by the Stokes equations. The strains in the solid are smallallowing for the solid to be described by the Lame equations, but no restric-tions are applied on the magnitude of the displacements leading to stronglycoupled, nonlinear fluid-structure problem. The FSI problem is then solvednumerically by an iterative procedure which solves sequentially fluid and solidsubproblems. Each of the two subproblems is discretized by finite elementsand the fluid-structure coupling is reduced to an interface boundary condition.Several numerical examples are presented and the results from the numericalcomputations are used to perform permeability computations for different ge-ometries.

1. Introduction

Many applications require finding the effective properties of fluid flow in de-formable porous media. There problem has first been studied under the assump-tion that the fluid-structure interface displacements are small compared to themicroscale dimensions (pore size) [1, 12]. Applying homogenization leads to a lin-ear macroscopic Biot’s law [3]. Under less severe restrictions [9] have derived anonlinear macroscopic governing equations by allowing the interface displacementsto be of the same order as the pore size. The microscopic cell problems which ariseare formulated in terms of fluid flow around a rigid skeleton. This becomes possi-ble after a local rigid-body coordinate transformation of each unit cell after whichthe interface displacement become infinitely small in the new reference frame. Thefocus of this report is the case when it is not possible to perform a rigid body move-ment of the unit cell so that the local problems are reduced to the familiar Darcysetting. We use numerical techniques to solve fluid structure problems in channel-like geometries as well as in genuine 3-dimensional metallic foam structures. Basedon the numerical solution upscaling is then performed to obtain effective propertiesof the deformable porous medium.

Date: October 19, 2004.Key words and phrases. fluid-structure interaction, deformable porous media, upscaling, linear

elasticity, stokes, finite elements.

1

2 OLEG ILIEV, ANDRO MIKELIC, AND PETER POPOV

2. The Fluid-Structure Interaction (FSI) Problem

Before we present the fluid structure problem we begin with a brief summary ofthe notation used, the formulation of the fluid, and solid problems alone. Considera continuum body. The material points in the body are associated with points p inR

3. A body is defines as an open subset Ω ⊂ R3. We denote the reference configu-

ration by Ω0 which represents the body before the deformation has begun and thedeformed configuration by Ω. Further, we are concerned with stationary processesso the both the Lagrangian fields associated with the solid are time independent asare the Eulerian ones for the fluid. We will use p to denote material coordinatesand x to denote spatial ones.

2.1. Solid. Consider the solid first. It undergoes a continuous, invertible deforma-tion

x = x(p),

so the reference and deformed configurations are connected by Ω = x(Ω, t). Thedeformation gradient is denoted by F:

(1) F(p) = ∇x(p),

and the displacements u(p) are given by

(2) u(p) = x(p) − p.

The usual infinite small strain tensor is introduced:

(3) E(p) =12(∇u(p) + ∇u(p)T

)

and we denote by T(x) the Cauchy stress in the deformed configuration. Note,that it is spatial field. The field equations for elastic bodies are best formulated inmaterial description, so we need the first Piola-Kirchhoff stress tensor S(p) whichis related to T(x) by:

(4) S(p) = det(F(p))T(x(p))F−T (p).

Hyperelastic bodies in general are defined as materials for which S is the gradientof a potential1, that is there exists a function W (F) such that:

(5) S(F) =∂W (F)

∂F.

Then, given a body force b0 in the reference configuration, the boundary valueproblem for a solid is stated (in the reference configuration) as follows: Find u(p)such that

(6) ∇ · S + b0 = 0 in Ω0,

with Dirichlet

(7) u = u on ΓD0

and/or Neumann

(8) Sn0 = s on ΓN0

1Equivalently, hyperelastic materials can be defined as materials which produce non-negativework on a closed cycle, (Cf., e.g. [6, pg. 184-191])

FSI PROBLEMS IN DEFORMABLE POROUS MEDIA 3

boundary data, with the usual conditions ΓD0

⋂ΓN

0 = ∅ and ΓD0

⋃ΓN

0 = Γ0. Underthe assumption that S(I) = T(I) = 0 and the assumption that ∇u is small one canintroduce the fourth-order, linearized elasticity tensor

(9) C = DS(I)

It is a simple calculation to show that

(10) C = DT(I)

and

(11) S = C : E + o(∇u).

As a result, the balance of linear momentum (6) is linearized as

(12) −∇ · (C : E) = b0.

The relation (11) is known as Hooke’s law. Note that, for an isotropic material, theelasticity tensor C is necessarily expressed in terms of the two Lamme constants,λs and µs (cf. eg. [11, 6]), so the stress tensor reduces:

(13) S = λstr (E) I + 2µsE.

2.2. Newtonian fluid at low Reynolds number. Newtonian fluids are best de-scribed using spatial fields. For stationary problems (the spatial description of allinvolved quantities is time independent), one has a velocity v(x) and correspond-ingly, the symmetric part of the velocity gradient, namely, the stretching tensorD(x) given by:

(14) D(x) =12(∇v(x) + ∇v(x)T

)

By definition, a Newtonian fluid is one for which:

(15) T = −pI + 2µD

where µ is the absolute viscosity of the fluid2. The fluid must satisfies the conser-vation of mass, that is:

(16) ∇ · v = 0,

and conservation of linear momentum

(17) ρ(v · ∇)v = −∇p + µ∆v + b,

where b is a distributed body force (per unit volume) acting on the fluid. Withthe additional assumption that the term (v · ∇)v is small with respect to the restof the terms in 17 and can be disregarded, the balance of linear momentum takesthe form

(18) −µ∆v + ∇p = b.

The system (16,18) is known as the Stokes equations.

2To be precise, this is a definition of a Newtonian fluid which is also independent under changein observer. (Cf., e.g. [6, pg. 147-151])


(a) Reference configuration (b) Deformed configuration

Figure 1. Schematic of the fluid and structure domains.

2.3. Statement of the coupled FSI problem. Consider now the stationaryfluid-structure problem (Figure 1) in the deformed configuration Ω = Ωf ∪ Ωs,where the fluid occupies Ωf , the solid occupies Ωs and Ωf ∩ Ωs = ∅. The part ofthe boundary shared between the fluid and the solid is denoted by ΓI = ∂Ωf ∩∂Ωs.We also assume that the deformation x(p) is such that we don’t have contactproblems, break-up of the boundary and so on. We have two conditions on theinterface. First we need kinematic compatibility, that is, the velocity of the fluidon the interface should be equal to the velocity of the interface itself. This, for astationary problem, implies that

(19) v = 0 on ΓI .

Further, we also need continuity of tractions, namely (C.f., e.g. [9]):

(20) Tfn = Tsn on ΓI

where n = ns and ns is the outward normal to the solid domain. The stress Tf

in the fluid is given by equation (15). Further, we can express the Cauchy stressTs in terms of the Piola-Kirchhoff stress using equation (4), which together withHooke’s law (11) implies:

(21) −pn + 2µDn = det(F)−1(C : E)FTn on ΓI .

The kinematic conditionIn order to simplify notation, it is convenient to introduce introduce the sym-

metric part of gradient operator:

e(w) =12(∇w + (∇w)T

),

for some field w. Observe, that E = e(u) and D = e(v). The FSI problemtherefore consists of finding the interface between the two domains, a velocity,pressure and displacements which solve the Stokes (16), (18), and Lamme equations(12) respectively, and also satisfy the interface conditions (19) and (21). Moreformally, the FSI boundary-value problem is summarized bellow in terms of theunknowns ΓI , v, p and u: Find ΓI , v, p and u such that:

(22) ΓI =p + u(p)|∀p ∈ ΓI

0

,


(23)−µ∆v + ∇p=b in Ωf ,

∇ · v=0 in Ωf ,−∇ · (C : e(u))=b0 in Ωs

0,

(24) det(∇u + I)(−pI + 2µe(v)) (∇u + I)−T n0 = (C : e(u))n0 on ΓI0,

v satisfies the kinematic interface condition (19) and in addition v, p and u shouldalso satisfy any boundary conditions that might be specified on ∂Ω\ΓI. Equation(24) is the continuity of tractions (21) expressed on the reference position ΓI

0 of theinterface. Observe that the position of the interface is part of the boundary valueproblem, and the solid-fluid coupling term (24) makes it a nonlinear one.

2.4. Weak form of the Elasticity, Stokes and FSI problems. Consider aLipschitz domain Ω and let (·, ·)Ω be the usual inner product on L2(Ω) and as thereis no chance of confusion, also the inner product on

[L2(Ω)

]d, where d = 2, 3 isthe size of the spatial dimension. Let Hs

0(Ω), −1 ≤ s ≤ 1 be the Sobolev spacesand L2

0(Ω) be Hilbert space of functions in L2 having zero mean. For completedevelopment discussion on these subjects, including fractional Sobolev spaces, see[10]. Suppose that both Ωs

0 and Ωf are Lipschitz domains. To formulate theelasticity problem, one introduces the bilinear form

aΩ0(u,w) =∫

Ω0

(C : e(u)) : e(w)dx.

Let u0 ∈[H1/2(ΓD

0 )]d

be the Dirichlet data given on ΓD0 ⊂ ∂Ωs

0, s be the Neumann

data given on ΓN0 ⊂ ∂Ωs

0, and let b0 ∈[H−1(Ωf )

]d be the distributed body force.

The weak form of the linear elasticity problem is: Find u ∈[H1(Ωs

0)]d such that

(25) aΩs0(u,w)=(b0,w)Ωs

0+ (s,w)ΓN

0, ∀w ∈

[H1

D(Ωs0)]d

,

u=u0, on ΓD0 .

The stability of the weak problem is a standard subject that can be found forexample in [13]. In the case of the Stokes problem, assume, again for simplicity,that we are given a homogeneous Dirichlet problem. The weak form of the Stokesequation is: Find v ∈

[H1

0 (Ωf )]d, p ∈ L2

0(Ωf ) such that

(26) DΩf (v,w) − (p,∇ · w)Ωf =(b,w)Ωf , ∀w∈[H1

0 (Ωf )]d

,−(∇ · v, q)Ωf =0, ∀q ∈L2

0.

Here DΩf (v,w) is the vector Dirichlet form

DΩf (v,w) =∫

Ωf

µ∇v : ∇wdx.

For a complete discussion of the weak problem (26) the reader is referred to [5].Finally, let is introduce the form

(27) gΓI0(v,u, p,w) =

∫

ΓI0

det(∇u + I)(−pI + 2µe(v)) (∇u + I)−T n0

·wds.

One obvious way to restate the boundary value problem (23)-(24) in a weak form:Find the interface ΓI , the deformed configuration of the fluid domain Ωf , the dis-placements u ∈

[H1(Ωs

0)]d, velocity v ∈

[H1

0 (Ωf )]d and pressure p ∈ L2

0(Ωf ) such


that

(28)

DΩf (v,w) − (p,∇ ·w)Ωf =(b,w)Ωf , ∀w∈[H1

0 (Ωf )]d

,−(∇ · v, q)Ωf =0, ∀q ∈L2

0,

aΩs0(u,w)=(b0,w)Ωs

0+ gΓI

0(v,u, p,w), ∀w∈

[H1

D(Ωs0)]d

,

Γ=p + u(p)|∀p ∈ Γ0 .

3. Solution methods for the coupled FSI system

Upon introducing the finite-dimensional subspaces Uu, Uv and Up for the dis-placements, velocity and pressure respectively, the first three equations in (28) leadto the following nonlinear system of algebraic equations:

(29)

⎛

⎝A(u) CT (u) 0C(u) 0 0

0 0 K

⎞

⎠

⎛

⎝vpu

⎞

⎠ =

⎛

⎝f10

f2 + g(u,v,p)

⎞

⎠ ,

where the blocks A(u) and K correspond to the bilinear forms DΩf (·, ·) and aΩs0(·, ·),

while the blocks C(u) and C(u)T couple the velocity and pressure unknowns. Sincethe position of ΓI and hence Ωf depends on u, both A and C are functions of thedisplacement. The coupling between the fluid and the structure (21) appears onthe right hand side as the nonlinear function g(u,v,p) which corresponds to theform gΓI

0(·, ·, ·, ·). Note that we are evaluating the fluid tractions in the reference

configuration, i.e. on ΓI0.

3.1. Dirichlet-Neumann iterative scheme. One way to solve the BVP (23),(22) is to use an iterative scheme which successively solves separate problem on thetwo domains. Considering the following approach:

• Solve the Stokes equation in the fluid domain treating the solid as a rigidbody;

• Transfer the forces to the solid;• Calculate the displacement field in the solid and then update the fluid

domain.The second step needs some further clarification. In order to solve an elasticityproblem we need boundary conditions on the reference configuration, while thetractions computed from a fluid solution are given in the deformed configuration.For given velocities and pressure, the interface condition (21) specifies the tractionon the solid domain in the deformed configuration, so we have to use the definition(4) of the Piola-Kirchhoff stress tensor and convert the tractions to the referenceconfiguration. It is not difficult to check that the above algorithm corresponds toa fixed point iteration using the following linearization of of (29), starting withu0 = 0, v0 = 0, p0 = 0:

(30)

⎛

⎝A(uk) CT (uk) 0C(uk) 0 0

0 0 K

⎞

⎠

⎛

⎝vk+1

pk+1

uk+1

⎞

⎠ =

⎛

⎝f10

f2 + g(uk,vk+1,pk+1)

⎞

⎠ .

This iterative scheme can be expresses more explicitly in the following

Algorithm 3.1.1. (Dirichlet-Neumann domain decomposition method for the FSIproblem) Set u0 = 0. For k = 0, 1, ... until convergence do:


(1) Find vk, pk which satisfy the Stokes equations (16),(18) in Ωfk with the no-

slip boundary condition on the interface ΓIk and the appropriate boundary

conditions on ∂Ωfk \ ΓI

k.(2) Compute the traction tk = Tnk on the interface ΓI

k using equation (21).(3) Based on tk compute the tractions sk in the reference configuration of the

interface, i.e. ΓI0 using equation (4) and the current iterate for the displace-

ments uk.(4) Find uk+1 which satisfies the balance of linear momentum (12) in Ωs

k withSn0 = sk and the appropriate boundary data on ∂Ωs

0\ΓI0.

(5) Compute ΓIk+1 =

p + uk+1(p)|∀p ∈ ΓI

0

and Ωf

k+1:(6) Check convergence: ||uk+1−uk||ΓI

k+1< TOLERANCE ∗||uk+1||ΓI

k+1. The

norm is the discrete euclidian norm of the interface nodal values.

It is clear that if the interface converges to a fixed position we will have a veloc-ity and pressure field which satisfy the Stokes equation (16),(18), a displacementfield which satisfies the elasticity equations (12), and as a results of the convergedinterface, the interface condition (21) will also be satisfied.

3.2. Solution methods for the solid and fluid subproblems. At each itera-tions of algorithm 3.1.1 it is required to solve one solid and one fluid problem. Bothproblems are discretized using the FEM method. The solid problems is solvedby the Conjugate Gradient (CG) method (the two classical references are [7, 8])applied directly to the linear system of algebraic equations that results from thediscretization. The fluid problem is solved again by the CG method, but appliedto the Schur complement for the pressure variables. The next two sections describethe methodology in detail.

3.2.1. Elasticity problem. The Elasticity problem is solved by standard linear tri-angular elements. Given a triangulation T s

h of Ωs0, the approximation space for the

displacements is

(31) Uu =[

u ∈ C0(Ωs0)|u is linear on ∀τ ∈ T s

h

]d ⊂[H1(Ωs

0)]d

.

Denote by φji , i = 1...Nv, j = 1...d the nodal basis functions for the displacement

space Uu. The weak form (25) gives rise to the linear system

(32) Ku = b.

The stiffness matrix K is symmetric, positive definite and sparse, making it ideal forthe application of the CG method. In this work, the nodal unknowns are orderedby displacement, leading to the following partitioning into blocks of the stiffnessmatrix K:

(33) K =

⎛

⎜⎝

K11 · · · K1d

.... . .

...Kd1 · · · Kdd

⎞

⎟⎠

This ordering allows for a robust preconditioning of the linear system. Considerthe matrix KSDC = diag(K11, . . . ,Kdd) which contains only the diagonal blocks ofK. It can be shown that KSDC is an optimal preconditioner for K (see [4]), thatis, the condition number of the preconditioned matrix satisfies

(34) cond((KSDC)−1K) ≤ d − 1γ

1 − ν

1 − 2ν,


where γ is a constant independent of the mesh size3. In general the application of(KSDC)−1 at each CG iteration for (32) can be done by multigrid in linear time,leading to an optimal method. A different (and simpler) approach, selected here,is to use a MIC(0) factorization of KSDC which results in an O(h3/2) algorithm[4]. For the problems under consideration this proved to be adequate enough.

3.2.2. Stokes problem. The Stokes problem is solved using the LBB stable P2P1

(Taylor-Hood) element pair. Given a triangulation T fh of Ωf the approximation

spaces for the velocity and pressure are defined by

(35) Uv =[

v ∈ C0(Ωf )|v is quadratic polynomial on ∀τ ∈ T fh

]d

⊂[H1(Ωf )

]d

and

(36) Up =p ∈ C0(Ωf )|p is linear on ∀τ ∈ Th

⊂ H1(Ωf ) ⊂ L2(Ωf )

respectively. Denote again by φji , i = 1...Nv, j = 1...d the nodal basis functions

for the velocity space Uv and by ψi, i = 1...Np the nodal basis pressure space Up.Upon ordering the velocity unknowns first followed by the pressure, and applyingthe essential boundary conditions, the weak problem (26) results in a linear systemsof the block form:

(37)(

A CT

C 0

)(vp

)=(

fd

).

The block A = diag(A(1), ...,A(d)) has block-diagonal structure, each block givenby:

(38) A(l)ij = µ

∫

Ωf

d∑

k=1

φli,kφl

j,kdx, l = 1, ..., d,

and C = (C(1), ..., C(d)), where:

(39) C(l)ij = −

∫

Ωf

φli,lψjdx, l = 1, ..., d.

We only considered iterative method for the solution of (37). Since the linear systemis indefinite with zeros on the main diagonal, one can attempt to solve it by theMINRES method. However, preconditioning is not an obvious task, thereforewe chose the Schur complement method for the pressure [15, 2]. Since A is aninvertible matrix one can eliminate the first row of (37) to obtain

(40) CA−1CTp = CA−1f − d.

The Schur complement S = CA−1CT is a symmetric, positive definite system (seee.g. [14]) therefore we can solve the system (40) using the CG method. Now, ingeneral, S is a dense matrix and it is expensive to evaluate it explicitly. However,the CG algorithm only requires the computuation of the action of S on a vector, andthis can be done, provided that one can solve linear systems with A efficiently. Thevelocity block A is a block diagonal matrix, each block corresponding to a Laplacianstiffness matrix, and these can be inverted efficiently (in O(N) operations with

3It actually appears in Korn’s inequality, for details, see [4, 13]

Toward permeability of deformable porous media - Fraunhofer

Documents

Transcript of Toward permeability of deformable porous media - Fraunhofer