- Email: [email protected]

Incompressible ﬂows in elastic domains: an immersed boundary method approach Santos Alberto Enriquez-Remigio, Alexandre Megiorin Roma

*

Instituto de Matema´tica e Estatı´stica, Universidade de Sa˜o Paulo, Caixa Postal 66281, Sa˜o Paulo, SP 05311-970, Brazil Received 1 September 2002; received in revised form 1 May 2004; accepted 19 July 2004 Available online 12 October 2004

Abstract This study illustrates how the immersed boundary method may be applied to perform the numerical simulation of incompressible ﬂows in two-dimensional domains bounded by elastic boundaries. It presents the basic intermediate steps involved in the derivation of a solution methodology, from a scientiﬁc motivation to the numerical results, which can be applied for both steady and transient problems, even when the boundaries have an arbitrary shape. Its motivation, brieﬂy presented, was borne in a bioengineering problem: the numerical simulation of the performance of ventricular assist devices. The mathematical model is composed by the Navier–Stokes equations, where the forcing term contains singular forces which arise from the elastic stresses acting on the boundaries. The incompressibility constraint is modiﬁed to introduce the inﬂow and outﬂow conditions into the problem through the use of sources and sinks. The methodology is applied to simulate two problems: the steady ﬂow between two parallel plates, for which the exact solution is known and can be used to validate the approach, and the periodic ﬂow in a winding channel, a transient problem in a non-trivial domain. 2004 Elsevier Inc. All rights reserved. Keywords: Incompressible viscous Navier–Stokes equations; Computational ﬂuid dynamics; Interface problems; Immersed boundary method

*

Corresponding author. Fax: +55 11 30916131. E-mail addresses: [email protected] (S.A. Enriquez-Remigio), [email protected] (A.M. Roma). URL: http://www.ime.usp.br/~roma.

0307-904X/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2004.07.007

36

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

1. Introduction Ventricular assist devices (VADs) are employed to maintain proper blood pressure and ﬂux, assisting blood circulation during cardiac cycles. Some of them, the paracorporeal VADs, can assist one or both ventricles for weeks or months, and others, the implantable ones, can assist the left ventricle for years! This explains why, nowadays, the use of these devices has increased as a therapeutic option for patients waiting for a cardiac transplantation [5]. For example, a pneumatically driven, free-ﬂoating membrane paracorporeal VAD developed at the Heart Institute (InCor), University of Sao Paulo Medical School (FM-USP) is shown in the Fig. 1. It is composed of a pumping unit (detail), and an impelling-control unit (not displayed). As can be observed from Fig. 1, the VAD pumping unit is connected to the patients heart by ﬂexible tubes, giving rise to relevant bioengineering questions [1], such as ‘‘Which is the most appropriate shape for a tube cross-section? Are tapers of any advantage, and if they are, where should they be placed?’’. The mathematical modelling and the numerical simulation of this problem in three-dimensions become quite involved. Nevertheless, if the tube is constrained to a plane, a ﬁrst approach to this problem is the study of a laminar ﬂow in one of the tube sections along its length, where sources and sinks can be employed to introduce inﬂow and outﬂow conditions. This simpliﬁed version of the problem has motivated the study of ﬂows in two-dimensional ﬂexible domains. Problems such as this, involving the interaction between a ﬂuid ﬂow and elastic interfaces, may appear in several branches of science such as engineering, physics, biology, and medicine. Regard-

Fig. 1. Schematic drawing of a paracorporeal VAD pumping unit (detail).

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

37

less the ﬁeld, as a rule, they share in common a high degree of complexity, often displaying intricate geometry or time-dependent elastic properties, turning the problem into a real challenge for applied scientists, from both the mathematical modelling and the numerical simulation points of view. In the early 70s, Peskin introduced a mathematical model and a computational method to study the ﬂow patterns of the blood around the mitral valve. Through years, Peskins method proved to be a powerful tool, applicable to general problems involving a transient incompressible viscous ﬂuid containing an immersed elastic interface, which may have time-dependent geometry or elastic properties, or both. This work aims to illustrate how Peskins method may be employed to perform the numerical simulation of incompressible ﬂows in two-dimensional domains bounded by elastic boundaries. Peskins immersed boundary method is described in Section 2, where the mathematical model and the computational scheme employed can be found, with a special attention dedicated to the inﬂow and outﬂow conditions. Section 3 presents two problems which serve to illustrate the approach: the steady ﬂow between two parallel plates, used to validate the methodology, and the periodic ﬂow in a winding channel, a more diﬃcult transient problem. The numerical results are included in Section 4, and the conclusion in Section 5.

2. Immersed boundary method Peskin [14,15] introduced a new method to perform numerical simulations of the hemodynamics around the heart valves, inside the cardiac chambers. Since he assumed, to simplify the problem, that all cardiac structures were equally immersed in blood, his method became known latter as the immersed boundary method. Some very intricate problems can be handled by Peskins method. As examples, one can mention besides the cardiac simulation problem [13,16–21,23], the platelet aggregation during the blood clotting [11], the suspension of particles in three-dimensional Stokes ﬂow [10], the peristaltic pumping of solid particles [8], the ﬂuid dynamics of the inner ear [3], the aquatic animal locomotion [4,7,9], the three-dimensional incompressible ﬂows in collapsible tubes [26], the bioﬁlm process [6], the two-dimensional ﬂow past a cylinder [12], and the arteriolar blood ﬂow and mass transport [2]. Next, a mathematical model and a computational scheme for a two-dimensional incompressible ﬂow in an arbitrary, elastic channel-like domain will be presented. 2.1. Mathematical model Assuming that the walls delimiting the channel are inﬁnitely thin, massless, and completely immersed in the same ﬂuid which ﬂows in their interior, the immersed boundary method models the dynamical interaction of the ﬂuid with the channel walls by the transient, viscous incompressible Navier–Stokes equations o uðx; tÞ þ ½uðx; tÞ ruðx; tÞ þ rpðx; tÞ ¼ lDuðx; tÞ þ f ðx; tÞ; q ot

ð1Þ

38

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

r uðx; tÞ ¼ wðx; tÞ;

ð2Þ

where q, the mass density, and l, the viscous coeﬃcient, are both constant, and u(x, t), p(x, t) are, respectively, the velocity of the ﬂuid and the hydrodynamical pressure at the point (x, t), for x 2 X, t P 0. The driving force, f (x, t), contains singular force terms arising from the elastic stresses acting on the channel walls. The function w(x, t), in the right hand side of (2), models the inﬂow and outﬂow conditions, vanishing only outside of the regions where there is ﬂuid inﬂow or ﬂuid outﬂow. In the absence of any other force to drive the ﬂow, f (x, t), in (1), is given by the singular elastic force distribution Z Fðs; tÞdðx Xðs; tÞÞ ds; ð3Þ f ðx; tÞ ¼ S

where d(Æ) is the Diracs delta, and F(s, t) is a force density deﬁned on the channel walls X(s, t), s 2 S, modelling the elastic behavior. From the no-slip condition, channel walls evolve in time accordingly to Z o : uðx; tÞdðx Xðs; tÞÞ dx; ð4Þ Xðs; tÞ ¼ uðXðs; tÞ; tÞ ¼ ot X where formally, the delta function properties were used. Note that the set of Eqs. (1)–(4) is a mixed Euler–Lagrangian formulation, where Navier– Stokes equations (1)–(3) are deﬁned on the Eulerian domain X, while the ﬂuid-boundary interaction equation (4) is deﬁned on the Lagrangian domain S. Next, a force density F, modelling the elastic behavior of the channel walls, and a scalar function w, modelling the inﬂow and outﬂow conditions, will be presented. 2.1.1. Modelling the elastic behavior of the channel walls If the channel wall points are assumed to be elastically connected to each other, applying Newtons second law, one obtains, recalling that the walls have negligible mass, the force F B ðs; tÞ ¼

o ðT sÞ; os

ð5Þ

where s is the unit tangent, s¼

oX=os koX=osk

ð6Þ

and T = T(koX/osk, s, t), the tension along the walls, is given by a generalized Hookes law. Further details on the derivation of (5) can be found in [18]. As a modelling option, yet another term will be added to FB to form the force density F. Tethers will be supposed to connect channel wall points to ﬁxed anchor points in X. A restoration force acting between these two kinds of points is given by F T ðs; tÞ ¼ j½Xðs; tÞ X 0 ðsÞ;

ð7Þ

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

39

where j is the elastic constant, and X0(s) is an anchor point position. Supposing that the channel walls have their equilibrium position at their initial conﬁguration, one has X0(s) = X(s, 0). In this case, the elastic force (7) drives a channel wall point back to its original position as the ﬂow pushes it away from it. Considering (5)–(7), the elastic behavior of the channel wall is modelled by o ðT sÞ j½Xðs; tÞ Xðs; 0Þ; ð8Þ os where the tension T and the constant j are channel dependent. By tuning T and j, one can control the ﬂexibility allowed for the channel wall; for example, the bigger j is adopted the stiﬀer the channel wall becomes. Fðs; tÞ ¼ F B ðs; tÞ þ F T ðs; tÞ ¼

2.1.2. Modelling the inﬂow and outﬂow conditions Following the ideas of Peskin [16], and of Arthurs et al. [2], inﬂow and outﬂow conditions can be introduced into the problem by employing sources and sinks placed conveniently in the domain. Due to the presence of sources and sinks, the scalar function w(x, t), in (2), must be non-zero on certain regions of the domain. To model the eﬀect on the ﬂow of only one source-sink pair, (2) can be written as : ð9Þ r u ¼ wðx; tÞ ¼ W1 ðxÞ Q1 ðtÞ; where W1(x) is a function specifying the location of the source and of the sink, and Q1(t) is the ﬂux given by the volume rate of ﬂow (the source strength). If, for the sake of simplicity, periodic boundary conditions are adopted for the velocity ﬁeld, the application of the Divergence Theorem to (9) leads to Z Z Z Z wðx; tÞ dx ¼ r u dx ¼ u n dS ¼ 0; ð10Þ Q1 ðtÞ W1 ðxÞ dx ¼ X

X

X

oX

where n is the unitary ‘‘exterior normal’’ to X. Neglecting the trivial case where the ﬂux Q1(t) is identical to zero, (10) implies that XW1(x) dx = 0, which imposes a constraint to the choice of the function W1(x). A particularly simple choice for W1(x) is : W1 ðxÞ ¼ wI1 ðxÞ wO ð11Þ 1 ðxÞ; where wI1 ðxÞ and wO 1 ðxÞ are two non-negative functions with compact support and integral equal to one. In the previous notation, the superscript ‘‘I’’ stands for ‘‘inﬂow’’ (from the source), and ‘‘O’’ for ‘‘outﬂow’’ (into the sink). This idea can readily be extended to M source-sink pairs, for any arbitrary positive integer M, by writing w(x, t) as wðx; tÞ ¼

M X

Wm ðxÞ Qm ðtÞ;

ð12Þ

m¼1 I O where Wm ðxÞ ¼ wIm ðxÞ wO m ðxÞ, with wm ðxÞ, wm ðxÞ non-negative functions with compact support and integral equal to one, and Qm(t) is the ﬂux related to the mth source-sink pair, for m = 1, 2, . . ., M. In (12), it is assumed that to every source there is a corresponding sink with same

40

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

strength in absolute value. A variation of (12) was employed in [2] to model blood ﬂow and mass transport in arterioles. To model the ﬂux through sources and sinks, a linear resistance model can be adopted IO : P m P m ðtÞ ; Qm ðtÞ ¼ Rm

ð13Þ

where the constants Pm and Rm, 1 6 m 6 M, are chosen to determine the properties of the mth source-sink pair, and P IO m ðtÞ is the average pressure diﬀerence between the pressures over the mth source and over its corresponding sink, Z IO pðx; tÞWm ðxÞ dx: ð14Þ P m ðtÞ ¼ X

Finally, note that integrating (12) on the domain X, one obtains Z Z Z M X I O wðx; tÞ dx ¼ wm ðxÞ dx wm ðxÞ dx Qm ðtÞ ¼ 0; X

m¼1

X

ð15Þ

X

which means that ﬂuid is neither created nor destroyed inside X, as it should be. 2.2. Computational scheme The basic idea is to solve ﬂuid Eqs. (1) and (2) on a ﬁxed uniform mesh and Eq. (4) on a moving mesh. Eulerian variables uij = (uij, vij) and pij, the numerical approximations of the velocity and of the pressure, are deﬁned on the computational mesh Xh ¼ fxij jxij ¼ ðx0 þ ih; y 0 þ jhÞ; 0 6 i; j 6 N g; where h is the mesh width (for convenience, kept the same both in x- and in y-directions). The discretization of the channel walls employs a moving computational mesh, represented as a ﬁnite collection of points Xs, 0 6 s 6 S, apart from each other by a distance Ds, usually taken as being h/2. Based on the work by Peskin and Printz [22], the discretization proposed for the equations of motion (1)–(4) is unþ1;0

unij ij q Dt

! ¼ f nij ;

ð16Þ

unþ1;1

unþ1;0 ij ij þ unij D0x unþ1;1 q ij Dt

unþ1;2

unþ1;1 ij ij þ vnij D0y unþ1;2 q ij Dt

!

nþ1;1 ¼ lDþ ; x Dx uij

ð17Þ

nþ1;2 ¼ lDþ ; y Dy uij

ð18Þ

!

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

41

the new velocity and pressure, unþ1 and pnþ1 ij ij , being obtained by solving the equations unþ1;2 ¼ unþ1 þ ij ij

Dt G pnþ1 ij ; q

ð19Þ

¼ wnþ1 D unþ1 ij ij :

ð20Þ

For Eqs. (3) and (4), one has the discretizations f nij

¼

S X

F ns dh ðxij X ns ÞDs;

ð21Þ

s¼0

¼ X ns þ Dt X nþ1 s

X

n 2 unþ1 ij dh ðxij X s Þ h :

ð22Þ

ij

The derivative with respect to s in (5) is approximated by the centered diﬀerence operator X sþ1 X s 1 : 2Ds The ﬁnite diﬀerence operators employed to approximate in the x-direction the ﬁrst derivatives are the forward, backward and centered diﬀerences, respectively D0s X s ¼

/iþ1;j /i;j /i;j /i 1;j /iþ1;j /i 1;j ; D

; D0x /i;j ¼ ; x /i;j ¼ h h 2h

0 these operators similarly deﬁned for the y-direction, their notation being given by Dþ y , Dy , and Dy . In terms of these ﬁnite diﬀerence operators, the discretizations of the Gradient and of the Divergence diﬀerential operators are, respectively, Dþ x /i;j ¼

G/i;j ¼ ðD0x /i;j ; D0y /i;j Þ

and

D ui;j ¼ D0x ui;j þ D0y vi;j :

The discrete Laplacian operator, L, is obtained by the product of the discrete Gradient operator by the discrete Divergence operator, D G/i;j ¼ L/i;j ¼

/iþ2;j þ /i 2;j þ /i;jþ2 þ /i;j 2 4/i;j

: 4h2 Observe that in (21) and (22), Diracs delta is approximated by the function dh. Motivated by a set of compatibility properties [16], the choice here is given by the product dh ðxi;j Þ ¼ d h ðxi Þ d h ðy j Þ; where d h ðzÞ ¼

(

ð23Þ

0:25 1 þ cos P2 z=h =h;

jzj < 2h;

0

jzj P 2h:

Other discretization alternatives can be found in [25,2].

ð24Þ

42

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

2.2.1. Numerical implementation of source-sink pairs As part of the numerical solution, one is required to solve on the computational mesh Xh the discrete equations (19) and (20), with wn+1 given by its full expression, (12). As in [2], employing the superposition principle, assume that the velocity and the pressure at time tn+1, un+1 and pn+1, can be written as the sums uþ unþ1 ¼

M X

um Qnþ1 m ;

ð25Þ

pm Qnþ1 m ;

ð26Þ

m¼1

pnþ1 ¼ pþ

M X m¼1

, where u p solve on Xh the incompressible problem with periodic boundary condition u unþ1;2 q ¼ G p; Dt D u¼0

ð27Þ

ð28Þ

and um , pm solve on the same grid the auxiliary problems um q ¼ G pm ; Dt D u m ¼ Wm ;

ð29Þ

ð30Þ

1 6 m 6 M, where periodic boundary conditions are also assumed. Observe that um and pm can be interpreted as coeﬃcients of ‘‘correction terms’’ which incorporate into the ﬂow eﬀects of the mth source-sink pair. Hence, (25) and (26) give un+1 and pn+1 as the corrected incompressible ﬂow by taking into account inﬂow and outﬂow conditions. To understand the origin of (29) and (30), substitute (25) and (26) into (19) to obtain M X q um m¼1

Dt

þ G pm

¼ 0; Qnþ1 m

ð31Þ

where (27) was employed, and (25) into (20) to obtain M X

ðD um Wm Þ Qnþ1 ¼ 0; m

ð32Þ

m¼1

where (28) and the full expression for wn+1 were both employed. Now, since the ﬂuxes Qnþ1 m , and pm , 1 6 m 6 M, are independent from each other, (29) and (30) follow the coeﬃcients um , immediately.

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

43

Note that the numerical solution of Eqs. (27), (28) and (29), (30), 1 6 m 6 M, is actually equivalent to solving Poisson equations with periodic boundary conditions. By substituting the divergence equation into the pressure equation, one obtains q q p m ¼ Wm ; ð33Þ D Gp ¼ D unþ1;2 and D G Dt Dt for 1 6 m 6 M, whose solutions can be obtained very eﬃciently through the application of Fast Poisson Solvers, such as the Fast Fourier Transform (FFT). To ﬁnish the computation of (25) and (26), it is still missing to determine the ﬂuxes Qnþ1 m . If one considers the inner product between two scalar functions deﬁned on the computational mesh Xh, f(x) and g(x), as : X hf ; gih ¼ fij gij h2 ; ð34Þ i;j

may be obtained from the discretization of the linear resistance model (13) then the ﬂuxes Qnþ1 m which can be written as ¼ Qnþ1 m

P m P IO;nþ1 P m hpnþ1 ; Wm ih m ¼ : Rm Rm

ð35Þ

In this case, Eqs. (25) and (26) implicitly deﬁne the solution ﬁelds un+1 and pn+1. An explicit form can be obtained by ﬁrst substituting (26) into (35), * + M X ¼ Pm pþ ; 1 6 m 6 M; pk Qnþ1 Rm Qnþ1 m k ; Wm k¼1

ð36Þ

h

which then, after some algebraic manipulation, can be rewritten as M X h pk ; Wm ih Qnþ1 þ Rm Qnþ1 ¼ P m hp; Wm ih ; k m

ð37Þ

1 6 m 6 M:

k¼1

Eq. (37) form a set of M linear equations in the unknowns Qnþ1 m , 1 6 m 6 M. In matrix notation, one has AQnþ1 ¼ b;

ð38Þ

where A is the matrix given by 2 6 6 6 6 6 6 6 6 6 6 4

h p1 ; W1 ih þ R1

h p2 ; W1 ih

h p 3 ; W1 i h

...

hpM ; W1 ih

h p1 ; W2 ih

h p2 ; W2 ih þ R2

h p 3 ; W2 i h

...

hpM ; W2 ih

h p1 ; W3 ih

h p2 ; W3 ih

h p3 ; W3 ih þ R3

...

hpM ; W3 ih

..

.. .

.. . h p 1 ; WM i h

h p2 ; WM ih

h p 3 ; WM i h

.

...

hpM ; WM ih þ RM

3 7 7 7 7 7 7; 7 7 7 7 5

44

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

2

P 1 h p ; W1 i h

6 6 P 2 h p ; W2 i h 6 6 6 P h p ; W3 i h 3 b¼6 6 6 .. 6 6 . 4

3 7 7 7 7 7 7 7 7 7 7 5

2

and

P M h p; WM ih

Qnþ1

Qnþ1 1

3

7 6 6 Qnþ1 7 6 2 7 7 6 6 nþ1 7 Q 6 ¼6 3 7 7: 7 6 6 .. 7 6 . 7 5 4 Qnþ1 M

2.3. Algorithm To obtain (un+1, Xn+1) from (un, Xn), one proceeds as follows: 1. Find the boundary force Fn explicitly from the boundary conﬁguration Xn: F n ¼ D0s ½T n sn jðX n X 0 Þ: 2. Spread the force Fn from the boundary to the surrounding ﬂuid mesh Xh: X F ns dh ðx X ns ÞDs; fn ¼ S

where x = (x, y) 2 Xh, dh(x) = dh(x)dh(y), with dh given by (24). 3. Solve (16)–(18) for the provisional velocity un+1,2. , 4. Find u p solving (27) and (28), and um , pm , for the mth source-sink pair, 1 6 m 6 M, solving (29) and (30). 5. Determine the source-sink ﬂuxes, Qnþ1 m , 1 6 m 6 M, solving the linear system (38). 6. Use the ﬂux values to ﬁnd the new velocity and pressure, un+1 and pn+1, performing the summations uþ unþ1 ¼

M X

um Qnþ1 m ;

m¼1

pþ pnþ1 ¼

M X

pm Qnþ1 m :

m¼1

7. Interpolate the new ﬂuid velocity to the wall points, updating their positions: X X nþ1 ¼ X ns þ Dt unþ1 dh ðx X ns Þh2 : s x

This completes the time step.

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

45

3. Application problems 3.1. Steady ﬂow between two parallel plates Consider the laminar, incompressible ﬂow of a viscous ﬂuid between two inﬁnite parallel plates located at y = a and y = +a, as shown in the Fig. 2. If one seeks a steady state solution in the form u(x, y) G (u(x, y), 0) and p G p(x), driven only by the pressure gradient, then the two-dimensional, viscous incompressible Navier–Stokes equations can be written as 2 ou op o u o2 u þ ; ð39Þ qu ¼ þ l ox ox ox2 oy 2 ou ¼ 0; ð40Þ ox with boundary conditions u(x, a) = u(x, +a) = 0, from the no-slip condition. In this context, (39) and (40) have a well known solution in closed form given by uðyÞ ¼

pðx2 Þ pðx1 Þ 2 ðy a2 Þ; 2lL

a 6 y 6 þa;

ð41Þ

with L = x2 x1 being the distance between two reference points x1 and x2, x1 < x2, where the pressures are known with the assumption that p(x1) > p(x2). The steady state volume rate of ﬂow through a cross-section S between the plates, S = {(x, y)jx = c, a 6 y 6 +a, c = constant}, is given by Z u ds; Q ¼ S

where ds = (1, 0) dy is the vector element of ‘‘area’’ normal to S. Hence,

Y Fixed wall +a p (x1 )

2

0

x1

p (x ) x2

X

–a Fixed wall

Fig. 2. Steady ﬂow between two parallel plates.

46

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

Q ¼

Z

þa

ðuðyÞ; 0Þ ð1; 0Þ dy ¼

a

Z

þa

a

pðx2 Þ pðx1 Þ 2 ðy a2 Þ dy; 2lL

from which one obtains, after integrating, Q ¼

pðx1 Þ pðx2 Þ R

ð42Þ

with 3lL : ð43Þ 2a3 The resulting expressions (42) and (43) are known as Poiseuille Law. It expresses the steady state ﬂux Q* as a function of the inﬂow pressure, p(x1), the outﬂow pressure, p(x2), and R, the resistance to the ﬂow imposed by the channel formed by the plates, which in its turn, depends on the viscosity of the ﬂuid l, on the channel half-width a, and on L, the distance between the two reference points x1 and x2, where inﬂow and outﬂow pressures are known. The mass ﬂux injected by the source, qQ1 , must be the same that goes through a cross-section at the center of the channel (conservation of the mass). Hence, employing (42) and (43), one obtains R¼

qQ1 ¼ q

pðx1 Þ pðx2 Þ 2a3 pðx1 Þ pðx2 Þ ; ¼q R1 L 3l

from which pðx1 Þ pðx2 Þ 3l ¼ Q1 3 : L 2a

ð44Þ

Substituting (44) into (41), one obtains uðyÞ ¼ Q1

3 2 ða y 2 Þ; 4a3

a 6 y 6 þa;

ð45Þ

which is the theoretical steady state velocity proﬁle through the channel, as a function of the twodimensional steady state ﬂux Q1 , and of the half-channel width. The exact solution (45) will be compared to the numerical solution. 3.1.1. Simulation parameters For the problem just described, the ﬂow takes place on an inﬁnite domain. However, its numerical simulation can only take place on a ﬁnite domain, which will be chosen as being the square X = [ 0.5, +0.5] · [ 0.5, +0.5] cm · cm. To prevent mass loss, semi-circular ‘‘caps’’ are employed at the ‘‘ends’’ of the channel which is formed by the parallel plates (Fig. 3). The eﬀect of the pressure gradient responsible for driving the ﬂow is modelled through the introduction of one source-sink pair into X (M = 1). The functions wI1 ðxÞ and wO 1 ðxÞ, deﬁning the location of the source and of the sink, were chosen as being wI1 ðxÞ ¼ dh ðx IÞ;

and

wO 1 ðxÞ ¼ dh ðx OÞ;

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

I

2h

2a

2h

2h

47

O 2h

L

Fig. 3. Idealized domain with semi-circular caps.

where dh (Æ) is the discrete Diracs delta (23),(24), whose compact support has size 2h ¼ 0:1455 cm. Note that h is kept ﬁxed and is chosen independently from the mesh spacing h. I and O are the compact support centers of wI1 and wO 1 , respectively. The ﬂux Q1(t) is given by the linear resistance model Q1 ðtÞ ¼

P 1 P IO 1 ðtÞ ; R1

ð46Þ

where P1 = 1.0 · 103 g s 2 and R1 = 1.0 · 105 g s 1 cm 3. The ﬂuid mass density by unit length is q = 1.0 g cm 2, the half-distance between the walls is a = 1.1 · 10 1 cm, the distance between the geometrical centers of the compact supports of

1 cm, and the viscosity is l = 1.0 · 10 2 g cm 1 s 1. wI1 ðxÞ and wO 1 ðxÞ is L = kO Ik = 5.82 · 10 The channel walls, immersed in the ﬂuid, have negligible mass and, by the mathematical modelling approach, exhibit elastic properties given by (8). As an option, it was assumed that FB G (0, 0) dyn cm 1, namely, the elastic forces arise only from tethered points. Since the exact solution derived is known only for a channel with rigid walls, the elastic constant j, in (7), was chosen to be large, j = 1.0 · 106 dyn cm 2. The mesh spacings employed were h = 1/32 cm, Ds = h/2 cm, and the ﬂuid is initially at rest, u(x,0) = (0, 0) cm s 1. 3.2. Periodic ﬂow in a winding channel The methodology presented is extremely versatile. To illustrate it further, a periodic ﬂux function will be employed, such that the linear resistance model (46) is now given by Q1 ðtÞ ¼ Qs ðtÞ aðtÞ P IO 1 ðtÞ; where the properties of the source-sink pair, P1 and R1, are allowed to change in time accordingly to P 1 =R1 ¼ Qs ðtÞ ¼ cosð2pt=T p Þ with period Tp = 1.0 s, and 1/R1 = a(t). The problem is transient. Moreover, the channel walls have a winding geometry as it will be seen in the next section. Except for the modiﬁcations explained above, all the simulation parameters to be employed are the same ones to be employed for the parallel plate problem.

48

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

4. Numerical results 4.1. Steady ﬂow between two parallel plates The numerical simulations for the steady state ﬂow between two parallel plates employed the parameters given in Section 3.1. The steady state was considered fully developed at time t = 1.0 s since changes in the ﬂux for latter time were less the 0.05% (Fig. 4). Table 1 presents the convergence ratios in the 2 and in the 1-norms of the diﬀerences between the numerical velocity uN, obtained at time t = 1.0 s on the cross-sections Sh = {(x, y)jx = 0, a 6 y 6 +a} ˙ Xh, h = 1/16, 1/32, 1/64, 1/128, and the exact velocity ue, obtained through (45), where the steady state ﬂux adopted was Q* = 0.01 (Fig. 4). The convergence ratios show that the immersed boundary method is of order one, as reported by other authors in several others diﬀerent situations and for several others diﬀerent applications [12,18,24].

1.0015 16x16 32x32 64x64 128x128

1.0010

1.0005

1.0000

0.9995

0.9990

0.9985

0.9980

0.9975

0.9970 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Fig. 4. Flux evolution in time for various mesh sizes.

Table 1 Convergence ratios, absolute errors in the middle of the channel at t = 1.0 s h N

1/16 16

Ratio

1/32 32

Ratio

1/64 64

Ratio

1/128 128

kuN uek2 kuN uek1

2.47 · 10 2 2.65 · 10 2

3.7 2.7

6.73 · 10 3 9.75 · 10 3

2.2 2.3

3.00 · 10 3 4.27 · 10 3

2.2 2.2

1.37 · 10 3 1.94 · 10 3

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

49

Table 2 Relative errors at time t = 1.0 s n kuN ue k2 kue k2 kuN ue k1 kue k1

16

32

64

128

0.40

0.11

0.04

0.02

0.38

0.14

0.06

0.028

The errors from Table 1, relative to the exact solution, can be found in the Table 2. Asymptoticaly, as expected, they decrease approximately by a factor of 2 as the computational meshes are doubled, meaning that the method behaves like a ﬁrst order method. The steady state reached is shown in Fig. 5. Note that, since the boundaries were supposed to be completely immersed in the same ﬂuid that ﬂows in its interior, although negligible, some spurious velocities are actually obtained outside of the channel as a result of the motion of the channel walls. Next, Fig. 6 shows the velocity proﬁles u(y) · y obtained at time t = 1.0 s on the cross-sections Sh = {(x, y)jx = 0, a 6 y 6 +a} ˙ Xh, h = 1/N, for N = 16, 32, 64, 128. The plots in solid lines

0.50

0.30

0.10

- 0.10

- 0.30

- 0.50

- 0.50

- 0.30

- 0.10

0.10

0.30

Fig. 5. Fluid ﬂow in the channel (N = 32, t = 1.0 s).

0.50

50

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

0.080

0.075

0.075 0.070 0.070 0.065 0.065 0.060 0.060 0.055

0.055

0.050

0.050

0.045 0.06

0.04

0.02

0

- 0.02

- 0.04

- 0.06

0.045 0.06

0.04

0.02

0

- 0.02

- 0.04

- 0.06

0.06

0.04

0.02

0

- 0.02

- 0.04

- 0.06

0.075

0.070 0.065

0.070

0.060 0.055

0.065 0.050 0.045

0.060

0.040 0.035

0.055

0.030 0.050

0.025 0.020

0.045

0.015 0.06

0.04

0.02

0

- 0.02

- 0.04

- 0.06

Fig. 6. Velocity proﬁles (u(y), y) obtained numerically on the cross-sections Sh = {(x, y)jx = 0, a 6 y 6 +a} ˙ Xh, h = 1/16, 1/32, 1/64, 1/128, are displayed in solid lines, while the ones obtained theoretically are displayed in dotted lines.

represent the numerical solution while the ones in dotted lines represent theoretical solution. As the mesh is reﬁned the numerical velocity proﬁle gradually approaches the theoretical parabolic proﬁle. 4.2. Periodic ﬂow in a winding channel Given a small > 0, as it is shown in the Fig. 7, the ﬂux function Qs(t) switches its sign in the intervals (0.25 , 0.25 + ) and (0.75 , 0.75 + ) and, hence, the ﬂow should switch its direction. It starts from left to right, t 2 (0.00, 0.25), then it comes to a halt at t = 0.25 s. Afterwards, it develops from right to left, t 2 (0.25, 0.75), and, ﬁnally, it switches back again, coming from left to right after halting once more at t = 0.75 s. In what follows, a was kept constant, a = 1/R1. Under these conditions, the obtained result is shown in Figs. 8 and 9, where the time evolution of the ﬂow can be seen at several instants of time.

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

51

1 0.8 0.6 0.4

Qs(t)

0.2 0

-0.2 -0.4 -0.6 -0.8 -1 0

0.25

0.5

0.75

1

t

Fig. 7. Periodic ﬂux function.

5. Conclusions This study illustrated how the immersed boundary method may be employed to obtain a methodology capable of performing the numerical simulation of viscous incompressible ﬂows in two-dimensional channel-like elastic domains. This methodology can handle both steady and non-steady problems, even when the elastic boundaries have an arbitrary shape. Brieﬂy, the scientiﬁc motivation was mentioned: the numerical simulation of the performance of ventricular assist devices, a bioengineering problem with a high social impact. A two-dimensional simpliﬁcation of this problem was the focus here. The basic mathematical model is given by the viscous incompressible Navier–Stokes equations, where the incompressibility constraint was modiﬁed to incorporate inﬂow/outﬂow conditions, and a singular forcing term was included in the momentum equation to take into account the elastic stresses acting on the elastic boundaries. Inﬂow and outﬂow conditions were brought into the problem through sources and sinks. A discretization for the equations of motion was presented, and the methodology applied to two problems, to the steady ﬂow between two parallel plates, for which the exact solution is known and could be used to validate the numerical approach for the steady state case, and to the periodic ﬂow in a winding channel, where the versatility of the approach was illustrated further by its application to a transient problem.

Acknowledgments For their many insights and suggestions during the elaboration of this study, the authors are in debt to Dr. K. Aﬀeld from Virchow Klinikum, Medizinische Fakulta¨t der Humboldt, Universita¨t

52

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

Fig. 8. Times t = 0.150, 0.200, 0.250, 0.300, 0.350, 0.450, 0.500, and 0.550 s (reading by column, from top to bottom, left to right).

zu Berlin, and to Dr. I. Cestari from InCor Bioengineering Division, FM-USP, and to Dr. L.C. de Castro Santos from IME-USP. We are also in debt to R.M. Yoshikawa for providing the original

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

53

Fig. 9. Times t = 0.650, 0.700, 0.750, 0.800, 0.850, and 1.00 s (reading by column, from top to bottom, left to right).

C code, from which the ﬁnal code could be prepared, and to M.M. Roma for the DAV schematic drawings. Financial support for this work was provided in part by CNPq, a Brazilian agency for research and technological development, grant #132499/1997-5. References [1] K. Aﬀeld, Virchow Klinikum, Medizinische Fakulta¨t der Humboldt––Universita¨t zu Berlin, Spandaver Damm 130, D-14050 Berlin, Germany, personal communication, 1996. [2] K.M. Arthurs, L.C. Moore, C.S. Peskin, E.B. Pitman, H.E. Layton, Modeling arteriolar ﬂow and mass transport using the immersed boundary method, J. Comp. Physiol. 147 (1998) 402–440. [3] R.P. Beyer, A computational model of the cochlea using the immersed boundary method, J. Comp. Physiol. 98 (1992) 145–162.

54

S.A. Enriquez-Remigio, A.M. Roma / Applied Mathematical Modelling 29 (2005) 35–54

[4] D.C. Bottino, An Immersed Boundary Model of Ameboid Deformation and Locomotion, Ph.D. Thesis, Tulane University, New Orleans, LA, 1998. [5] I.A. Cestari, S.A. Hayashida, L.F.P. Moreira, A. Bonı´sio, M. Maizatto, J.F. Iban˜ez, N.A.G. Stolf, A.A. Leirner, Avaliac¸a˜o do desempenho in vivo do dispositivo de assisteˆncia ventricular (dav) incor, in: Proceedings of the 11th IMACS World Congress on System Simulation and Scientiﬁc Computation, vol. 3, 1998. [6] R. Dillon, L.J. Fauci, A. Fogelson, D. Gaver III, Modeling bioﬁlm processes using the immersed boundary method, J. Comp. Physiol. 129 (1996) 57. [7] L.J. Fauci, Interaction of oscillating ﬁlaments––a computational study, J. Comp. Physiol. 86 (1990) 294–313. [8] L.J. Fauci, Peristaltic pumping of solid particles, Comput. Fluids 21 (1992) 583. [9] L.J. Fauci, C.S. Peskin, A computational model of aquatic animal locomotion, J. Comp. Physiol. 77 (1988) 85–108. [10] A. Fogelson, C.S. Peskin, A fast numerical method for solving the three-dimensional Stokes equations in the presence of suspended particles, J. Comp. Physiol. 79 (1988) 50. [11] A.L. Fogelson, A mathematical model and numerical method for studying platelet adhesion and aggregation during blood clotting, J. Comp. Physiol. 56 (1984) 111–134. [12] M.-C. Lai, Simulations of the Flow Past an Array of Circular Cylinders as a Test of the Immersed Boundary Method, Ph.D. Thesis, New York University, New York, 1998. [13] D.M. McQueen, C.S. Peskin, A three-dimensional method for blood ﬂow in the heart: (II) contractile ﬁbers, J. Comp. Physiol. 82 (2) (1989) 289–297. [14] C.S. Peskin, Flow Patterns Around Heart Valves: A Digital Computer Method for Solving the Equations of Motion, Ph.D. Thesis, Albert Einstein College of Medicine, Yeshiva University, University Microﬁlms #72–30, July 1972, New York, p. 378. [15] C.S. Peskin, Flow patterns around heart valves: a numerical method, J. Comp. Physiol. 10 (1972) 252–271. [16] C.S. Peskin, Numerical analysis of blood ﬂow in the heart, J. Comp. Physiol. 25 (1977) 220–252. [17] C.S. Peskin, Fiber architecture of the left ventricular wall: an asymptotic analysis, Commun. Pure Appl. Math. 42 (1989) 79–113. [18] C.S. Peskin, D.M. McQueen, A three-dimensional computational method for blood ﬂow in the heart: (I) immersed elastic ﬁbers in a viscous incompressible ﬂuid, J. Comp. Physiol. 81 (1989) 372–405. [19] C.S. Peskin, D.M. McQueen, Computational bioﬂuid dynamics, Contemp. Math. 141 (1993) 161–186. [20] C.S. Peskin, D.M. McQueen, in: H.G. Othmer, F.R. Adler, M.A. Lewis, J.C. Dallon (Eds.), Fluid Dynamics of the Heart and its Valves, Prentice-Hall, Inc., Englewood Cliﬀs, NJ, 1996. [21] C.S. Peskin, D.M. McQueen, A general method for the computer simulation of biological systems interacting with ﬂuids, in: SEB Symposium on Biological Fluid Dynamics, Leeds, England, July 5–8, 1994. [22] C.S. Peskin, B.F. Printz, Improved volume conservation in the computation of ﬂows with immersed elastic boundaries, J. Comp. Physiol. 105 (1993) 33–46. [23] B.F. Printz, Computer Modeling of Blood Flow through the Heart During the Complete Cardiac Cycle, Ph.D. Thesis, City University of New York, New York, 1992. [24] A.M. Roma, A Multilevel Self-adaptive Version of the Immersed Boundary Method, Ph.D. Thesis, Courant Institute of Mathematical Sciences, New York University, University Microﬁlms #9621828, New York, January 1996. [25] A.M. Roma, C.S. Peskin, M.J. Berger, An adaptive version of the immersed boundary method, J. Comp. Physiol. 153 (1999) 509–534. [26] M.E. Rosar, A Three-Dimensional Computer Model for Fluid Flow Through a Collapsible Tube, Ph.D. Thesis, Courant Institute of Mathematical Sciences, New York University, New York, June 1994.