Incompressible flows in elastic domains: an immersed boundary method approach

Incompressible flows in elastic domains: an immersed boundary method approach

Applied Mathematical Modelling 29 (2005) 35–54 www.elsevier.com/locate/apm Incompressible flows in elastic domains: an immersed boundary method approa...

814KB Sizes 0 Downloads 52 Views

Applied Mathematical Modelling 29 (2005) 35–54 www.elsevier.com/locate/apm

Incompressible flows 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 flows in two-dimensional domains bounded by elastic boundaries. It presents the basic intermediate steps involved in the derivation of a solution methodology, from a scientific 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, briefly 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 modified to introduce the inflow and outflow conditions into the problem through the use of sources and sinks. The methodology is applied to simulate two problems: the steady flow between two parallel plates, for which the exact solution is known and can be used to validate the approach, and the periodic flow 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 fluid 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 flux, 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-floating 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 flexible 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 first approach to this problem is the study of a laminar flow in one of the tube sections along its length, where sources and sinks can be employed to introduce inflow and outflow conditions. This simplified version of the problem has motivated the study of flows in two-dimensional flexible domains. Problems such as this, involving the interaction between a fluid flow 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 field, 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 flow 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 fluid 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 flows 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 inflow and outflow conditions. Section 3 presents two problems which serve to illustrate the approach: the steady flow between two parallel plates, used to validate the methodology, and the periodic flow in a winding channel, a more difficult 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 flow [10], the peristaltic pumping of solid particles [8], the fluid dynamics of the inner ear [3], the aquatic animal locomotion [4,7,9], the three-dimensional incompressible flows in collapsible tubes [26], the biofilm process [6], the two-dimensional flow past a cylinder [12], and the arteriolar blood flow and mass transport [2]. Next, a mathematical model and a computational scheme for a two-dimensional incompressible flow in an arbitrary, elastic channel-like domain will be presented. 2.1. Mathematical model Assuming that the walls delimiting the channel are infinitely thin, massless, and completely immersed in the same fluid which flows in their interior, the immersed boundary method models the dynamical interaction of the fluid 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 coefficient, are both constant, and u(x, t), p(x, t) are, respectively, the velocity of the fluid 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 inflow and outflow conditions, vanishing only outside of the regions where there is fluid inflow or fluid outflow. In the absence of any other force to drive the flow, 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 defined 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 defined on the Eulerian domain X, while the fluid-boundary interaction equation (4) is defined 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 inflow and outflow 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 fixed 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 configuration, 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 flow 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 flexibility allowed for the channel wall; for example, the bigger j is adopted the stiffer the channel wall becomes. Fðs; tÞ ¼ F B ðs; tÞ þ F T ðs; tÞ ¼

2.1.2. Modelling the inflow and outflow conditions Following the ideas of Peskin [16], and of Arthurs et al. [2], inflow and outflow 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 effect on the flow 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 flux given by the volume rate of flow (the source strength). If, for the sake of simplicity, periodic boundary conditions are adopted for the velocity field, 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 flux 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 ‘‘inflow’’ (from the source), and ‘‘O’’ for ‘‘outflow’’ (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 flux 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 flow and mass transport in arterioles. To model the flux 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 difference 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 fluid is neither created nor destroyed inside X, as it should be. 2.2. Computational scheme The basic idea is to solve fluid Eqs. (1) and (2) on a fixed 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 defined 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 finite 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 difference operator X sþ1 X s 1 : 2Ds The finite difference operators employed to approximate in the x-direction the first derivatives are the forward, backward and centered differences, 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 defined for the y-direction, their notation being given by Dþ y , Dy , and Dy . In terms of these finite difference operators, the discretizations of the Gradient and of the Divergence differential 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 coefficients of ‘‘correction terms’’ which incorporate into the flow effects of the mth source-sink pair. Hence, (25) and (26) give un+1 and pn+1 as the corrected incompressible flow by taking into account inflow and outflow 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 fluxes Qnþ1 m , and pm , 1 6 m 6 M, are independent from each other, (29) and (30) follow the coefficients  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 efficiently through the application of Fast Poisson Solvers, such as the Fast Fourier Transform (FFT). To finish the computation of (25) and (26), it is still missing to determine the fluxes Qnþ1 m . If one considers the inner product between two scalar functions defined 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 fluxes 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 define the solution fields un+1 and pn+1. An explicit form can be obtained by first 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 configuration Xn: F n ¼ D0s ½T n sn  jðX n X 0 Þ: 2. Spread the force Fn from the boundary to the surrounding fluid 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 fluxes, Qnþ1 m , 1 6 m 6 M, solving the linear system (38). 6. Use the flux values to find 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 fluid 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 flow between two parallel plates Consider the laminar, incompressible flow of a viscous fluid between two infinite 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 flow 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 flow 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 flux Q* as a function of the inflow pressure, p(x1), the outflow pressure, p(x2), and R, the resistance to the flow imposed by the channel formed by the plates, which in its turn, depends on the viscosity of the fluid l, on the channel half-width a, and on L, the distance between the two reference points x1 and x2, where inflow and outflow pressures are known. The mass flux 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 profile through the channel, as a function of the twodimensional steady state flux 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 flow takes place on an infinite domain. However, its numerical simulation can only take place on a finite 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 effect of the pressure gradient responsible for driving the flow is modelled through the introduction of one source-sink pair into X (M = 1). The functions wI1 ðxÞ and wO 1 ðxÞ, defining 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 fixed and is chosen independently from the mesh spacing h. I and O are the compact support centers of wI1 and wO 1 , respectively. The flux 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 fluid 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 fluid, 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 fluid is initially at rest, u(x,0) = (0, 0) cm s 1. 3.2. Periodic flow in a winding channel The methodology presented is extremely versatile. To illustrate it further, a periodic flux 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 modifications 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 flow between two parallel plates The numerical simulations for the steady state flow 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 flux 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 differences 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 flux 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 different situations and for several others different 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 first 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 fluid that flows 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 profiles 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 flow 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 profiles (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 refined the numerical velocity profile gradually approaches the theoretical parabolic profile. 4.2. Periodic flow in a winding channel Given a small  > 0, as it is shown in the Fig. 7, the flux function Qs(t) switches its sign in the intervals (0.25 , 0.25 + ) and (0.75 , 0.75 + ) and, hence, the flow 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, finally, 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 flow 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 flux 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 flows 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. Briefly, the scientific motivation was mentioned: the numerical simulation of the performance of ventricular assist devices, a bioengineering problem with a high social impact. A two-dimensional simplification 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 modified to incorporate inflow/outflow conditions, and a singular forcing term was included in the momentum equation to take into account the elastic stresses acting on the elastic boundaries. Inflow and outflow 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 flow 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 flow 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. Affeld 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 final 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. Affeld, 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 flow 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 Scientific Computation, vol. 3, 1998. [6] R. Dillon, L.J. Fauci, A. Fogelson, D. Gaver III, Modeling biofilm processes using the immersed boundary method, J. Comp. Physiol. 129 (1996) 57. [7] L.J. Fauci, Interaction of oscillating filaments––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 flow in the heart: (II) contractile fibers, 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 Microfilms #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 flow 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 flow in the heart: (I) immersed elastic fibers in a viscous incompressible fluid, J. Comp. Physiol. 81 (1989) 372–405. [19] C.S. Peskin, D.M. McQueen, Computational biofluid 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 Cliffs, NJ, 1996. [21] C.S. Peskin, D.M. McQueen, A general method for the computer simulation of biological systems interacting with fluids, 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 flows 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 Microfilms #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.