Two-dimensional simulation of fluid–structure interaction using lattice-Boltzmann methods

Computers and Structures 79 (2001) 2031±2037

www.elsevier.com/locate/compstruc

Two-dimensional simulation of ¯uid±structure interaction using lattice-Boltzmann methods M. Krafczyk *, J. T olke, E. Rank, M. Schulz Lehrstuhl f ur Bauinformatik, Technische Universitat Munchen, Arcisstrasse 21, D-80290 Munich, Germany

Abstract The development of ¯ow instabilities due to high Reynolds number ¯ow in arti®cial heart-value geometries inducing high strain rates and stresses often leads to hemolysis and related highly undesired e€ects. Geometric and functional optimization of arti®cial heart valves is therefore mandatory. In addition to experimental work in this ®eld it is meanwhile possible to obtain increasing insight into ¯ow dynamics by computer simulation of re®ned model problems. Here we present two-dimensional simulation results of the coupled ¯uid±structure problem de®ned by a model geometry of an arti®cial heart value with moving lea¯ets exposed to a channel ¯ow driven by transient boundary conditions representing a physiologically relevant regime. A modi®ed lattice-Boltzmann approach is used to solve the coupled problem. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: CFD; Fluid±structure interaction; Heart valve; Lattice-Boltzmann methods

1. Introduction Many important aspects of hemodynamic aortic blood ¯ow have been investigated in experimental studies (see the literature cited in Ref. [1]). Unfortunately these studies sometimes lack spatial resolution, have problems analyzing transient phenomena or partially do not allow the detailed observation of interesting variables crucial for extended theoretical modeling. These drawbacks and limitations have suggested that numerical methods could be of signi®cant help when analyzing the dynamics of ¯owing blood. Numerical simulation has indeed led to a better understanding of blood ¯ow behavior and has allowed medical hemodynamicists to gain an improved understanding of the behavior of the containing structures (vessels and heart)

*

Corresponding author. E-mail addresses: [email protected] (M. Krafczyk), [email protected] (J. T olke), [email protected] (E. Rank), [email protected] (M. Schulz).

as well as the contained ¯ow (blood). Both of these structures can now be investigated in more detail, especially considering the in¯uence of transient e€ects. An important malfunction of these processes can be identi®ed as the pathological state known as thrombosis resulting from high shear stresses or increased coagulation due to the blood stagnancy in contact with the artery walls. Basically both phenomena will produce coagulated masses of blood, which ¯ow through the cardiovascular system until they reach a small crosssection vessel, causing the dangerous obstruction. Many long time problems in valve hemodynamics are due to the intrinsic non-linear weakly turbulent ¯ow dynamics and the resulting shear stresses. Yet an improvement of valve design requires a detailed understanding of these dynamic system properties and demands substantial further experimental research. As in vivo experimental investigations are extremely dicult, major process has to be obtained by in vitro measurements and numerical simulations. In the ®eld of computational and numerical analysis, several works can be found in the technical literature investigating steady (see e.g. the literature cited in Ref. [1]) and turbulent blood ¯ows. Most of the simulations described in the literature used FEM, FV or

0045-7949/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 1 ) 0 0 0 5 0 - 5

2032

M. Krafczyk et al. / Computers and Structures 79 (2001) 2031±2037

FD schemes to discretize the equations under investigation, usually the Navier±Stokes equation coupled with a continuity equation. Ideally one would like to analyze a fully transient 3D model including moving lea¯ets for a peak Reynolds number of about 1000±5000. Unfortunately such studies are still lacking. In a recent work [1] we analyzed the transient 3D ¯ow for ®xed heart-valve lea¯ets and a weakly turbulent ¯ow regime. To our knowledge the present study for the ®rst time investigates the coupled ¯uid±structure system of moving leaflets driven by a physiological blood ¯ow during a complete heart beat cycle. 2. The method: lattice-Boltzmann models for complex ¯ow phenomena For decades the scienti®c community has struggled to develop and optimize suitable algorithms like ®nite volume, ®nite element, spectral or ®nite di€erence methods to simulate high Reynolds number Navier± Stokes problems. Apart from these mainstream approaches the last years have seen further attempts to model ¯ow problems on particle-based methods such as the so-called lattice±gas method. An LG system consists of a large ensemble of idealized particles living on the nodes of a grid of given symmetry and undergoes two kinds of processes: collision and propagation. The collision is usually implemented as a set of prede®ned rules mapping pre- to post-collision states while obeying particle and momentum conservation. Frisch et al. [2] showed that on macroscopic scales with length scales large compared to the node distance and time scales large compared to the collisional relaxation time, the dynamics of such a particle ensemble is adequately described by the incompressible Navier±Stokes equations. Thus it is possible to simulate incompressible Navier± Stokes problems using LG methods. Further investigations in the last years have shown that LG methods are a valuable extension of `standard' discretization methods for ¯ow problems, especially for porous media and multi-phase ¯ow. On the other hand, they are not well suited for high Reynolds number ¯ow because (in addition to other problems) the minimum viscosity (implicitly de®ned by the collision rules) is quite large leading to a prohibitively huge number of degrees of freedom compared to state-of-the-art methods. Having studied the properties of LG methods there have been numerous attempts to improve the original approach leading to a variety of so-called lattice-Boltzmann (LB) methods, one of which has been used in this work. In order to see the di€erence between LB methods and well-known Navier±Stokes discretizations we will give a short theoretical overview over the ®rst technique. The reader interested in detailed information is referred to the cited literature.

The Boltzmann equation of ‡~ vrf ˆ X ot

…1†

used in statistical physics describes the dynamics of a continuous normalized particle distribution function f …~ x;~ v; t†, which is the probability to ®nd a particle with microscopic velocity ~ v…~ x; t†. The collision operator X contains details of the physical interaction between particles and can be of arbitrary complexity. The ®rst step in order to construct an LB model is to discretize the microscopic velocity space with a discrete set of vectors ~ ei , resulting in a so-called LB equation ofi …~ x; t† ‡~ ei rfi …~ x; t† ˆ Xi …ffj …~ x; t†jj 2 f1; . . . ; ngg†; ot i 2 f1; . . . ; ng

…2†

which in fact is a system of n ®rst-order PDEs coupled via Xi as will be seen below. It is known from statistical physics that for a collision operator of minimum complexity to generate ¯uid behavior the so-called single time relaxation approximation (STRA) from can be introduced 1 …fi s

Xi ˆ

fieq †

…3†

which is readily interpreted as a source term resulting from the deviation of fi from an equilibrium function fieq which still has to be de®ned. s is a microscopic relaxation time. The introduction of the minimum set of macroscopic ¯uid variables (i.e. velocity ~ u and density q) is done by de®ning qˆ

n X iˆ1

fi ;

~ uˆ

n 1X ei fi~ q iˆ1

…4†

and furthermore making an ansatz for the equilibrium distribution fieq of the form u† ‡ Ci …~ u†2 ‡ Di~ u2 Š fieq ˆ q‰Ai ‡ Bi …~ ei~ ei~

…5†

with constants Ai , Bi , Ci and Di still to be de®ned. The LB equation (2) can now be discretized in space and time, e.g. via a ®rst-order FD scheme resulting in fi …~ xk ‡~ ei dt; t ‡ dt† ˆ fi …~ xk ; t† ‡ Xi …~ xk ; t†

…6†

where k ˆ 1; . . . ; m and m is the number of grid nodes. This set of m  n equations can be interpreted as a system of evolution equations for a set of particle distribution functions fi , evaluated at nodes ~ xk of a uniform lattice with lattice spacing dh mimicking particle collision at these nodes and particle propagation to neighboring nodes located at ~ xk ‡~ ei dt during dt. Usually dh and dt are set of unity.

M. Krafczyk et al. / Computers and Structures 79 (2001) 2031±2037

linear convergence in time with respect to the exact solution of the corresponding Navier±Stokes problem.

Demanding that n X

Xi ˆ

n X

iˆ1

ei ˆ 0 Xi~

…7†

iˆ1

implying a conservation of mass and momentum, it can be shown that a multi-scale Chapman±Enskog expansion of the kinetic moments of the set of fi allows to choose appropriate coecients in Eq. (5) in order to obtain a set of equations for the dynamics of the macroscopic ¯uid properties. These equations can (in the low Mach number limit j~ uj  cs ) with a speed of sound cs ) be proven to be the incompressible Navier±Stokes equation o…~ u† ‡~ ur~ uˆ ot

1 rp1 ‡ mr2 …q~ u† q0

…8†

where p1 is the pressure and q0 is the conserved initial mean density, the continuity equation op ‡ r…q~ u† ˆ 0 ot

…9†

and an equation of state p1 ˆ c2s q

…10†

For the simulations presented here we used a 2D square grid with nodes connected to their neighboring nodes by grid vectors ~ ei dt …dt ˆ 1† given by the n ˆ 9 columns of the matrix (in units of dh)   0 1 0 1 0 1 1 1 1 N^ ˆ 0 0 1 0 1 1 1 1 1 where f1 represents the probability distribution of rest particles. The speed of sound p for  the 2D square grid can be computed to be cs ˆ …1= 3†…dh=dt† and the kinematic viscosity m is related to the relaxation time s via mˆ

2s

1 dh2 ; 6 dt

…s > 0:5†

2033

…11†

when using the so-called D2Q9 model [3], where the equilibrium distributions are de®ned as    eia ua ua ub eia eib fieq ˆ tp q 1 ‡ 2 ‡ 2 d ab cs 2cs c2s 4 tp ˆ ; i ˆ 1 9 …12† 1 tp ˆ ; i 2 f2; . . . ; 5g 9 1 tp ˆ ; i 2 f6; . . . ; 9g 36 a and b denote Cartesian coordinates and Einstein's summation convention is used. Numerical experiments have shown (e.g. Ref. [4]) that under certain assumptions such an LB scheme can be tuned to give second-order convergence in space and

2.1. Fluid±structure coupling Several basic schemes have been using in CFD to simulate coupled ¯uid±structure problems. They di€er mainly in the way the discretizations for the ¯uid and the solid part are coupled. Using independant grids for each problem part requires mapping between the primary variables of these meshes whereas the use of arbitrary Lagrangean Eulerian (ALE) methods (e.g. Ref. [5]) requires the deformation of the ¯uid mesh to adapt to the contours of the solid(s) implying interpolation and problems when topological changes are occurring during the course of a simulations (e.g. in inelastic collisions of particles). The scheme being used in this work follows a di€erent concept and is based on the work of Ladd [6,7]. In the following the basic ideas of this approach are sketched. Using a ®xed mesh for the coupling problem implies that nodes at a certain time belonging to a solid region may become part of the ¯ow domain after a later timestep and vice versa. Instead of extrapolating the corresponding structural/¯uid variables for these nodes from neighboring nodes which did not experience a `domain change' in the last timestep(s), one can rede®ne the problem in the following way: Let there be two ¯uid domains being separated by a virtual rigid wall which inhibits exchange of mass between the two domains (e.g. a balloon with a rigid hull). If the computation of the force acting onto the solid object is just taking into account the shear and pressure induced forces of the outer ¯uid region, the rigid body motion is independant of the dynamics of the inner ¯uid. The virtual rigid wall implies a strong velocity coupling between neighboring nodes of the inner and outer ¯uids. Thus a node changing from the inner to the outer domain can be expected to represent a `reasonably' close ¯ow velocity to that of its (outer domain) nodal neighbors, serving as a good `initial value' for the explicit time stepping procedure (Eq. (6)) used in this work. This approach can only be expected to give reasonable results, when the solids velocity is small compared to the lattice distance times the timestep size used. Now we will see how this concept can be incorporated into an LB scheme: In standard LB computations with ®xed geometries no-slip conditions …~ u ˆ~ 0† on wall nodes are applied by swapping all anti-parallel distributions fi and fi0 for which ~ ei ‡~ ei0 ˆ ~ 0 in every timestep (the so-called bounce-back rule). For moving geometric objects a wall node is no longer subject to zero velocity but to no-slip BCs in the moving coordinate system of the solid. The velocity of a node representing a part of the edge of a geometric object of mass M is determined by the translational and angular velocity (with respect to

2034

M. Krafczyk et al. / Computers and Structures 79 (2001) 2031±2037

a ®xed or moving rotation center ~ R0 ) of that object ~; x ~† resulting from the interaction between the ¯uid …U and the object on a larger scale as will be explained below. Thus a ®rst necessary ingredient for a two-way ¯uid±structure coupling is the need to set the velocity at a boundary node (by modifying all of his distributions fi ) to ~‡x ~  …~ ~ ub ˆ U r ‡ 12! ei

~ R0 †

…13†

without modifying the local stress tensor given by  X 1 ^ 1 Sab ˆ eia eib …fi fieq † …14† 2s i This is achieved by an additional term modifying the post-collisional bounce-back procedure for inner and outer boundary nodes (see Fig. 1) fi …~ ei † r ‡~ ei ; t ‡ 1† ˆ fi0 …~ r ‡~ ei ; t‡ † ‡ 2Bi q…~ ub~ ei † fi0 …~ r; t ‡ 1† ˆ fi …~ r; t‡ † 2Bi q…~ ub~

…15†

Here Bi was de®ned in Eq. (5) and t‡ indicates the postcollision state of f. ~ acting on the outer solid surface can The total force K be obtained by summing up all the force contributions (from all relevant links) of the corresponding outer wall nodes j given by   X 1 1 ~ kj ~ ei ; t ‡ ˆ rj ‡ ~ 2…fi …~ r; t‡ † fi0 …~ r ‡~ ei ; t‡ † 2 2 all links i 2Bi q…~ ei †~ ub~ ei †

…16†

where the factor …1=2†~ ei Eqs. (13) and (16) comes from the assumption that the virtual wall is situated exactly between an outer and an inner boundary node. Finally a simple time averaging is done to minimize oscillations due to so-called staggered moments ! X  1 1 ~ ~ ˆ 1 K…t ~ K…t† kj ~ rj ‡ ~ 1† ‡ ei ; t ‡ …17† 2 2 2 j which in turn allows to compute the new velocity of the solid object (being updated every two timesteps)

~ …t ‡ 1† ˆ U ~ …t U ~…t ‡ 1† ˆ x ~…t x

~ 1† ‡ 2M 1 K…t†; T …t† 1† ‡ 2I^ 1~

…18†

where ~ T represents the torque and I^ the moment of inertia tensor of the solid object. For M ! 1 the whole procedure reduces to the trivial case of non-moving solid objects. The code has been veri®ed computing various test problems similar to those in Refs. [6,7] including nontrivial properties like drag coecients for various Reynolds numbers. Details of these computations will be subject to a future publication [9].

3. Description of the problem setup and simulation details The simulation geometry (sketching a 2D cut through a CarboMedics [10] type bi-lea¯et aortic heart valve being located into a section of an artery) is shown in Fig. 2. Two pivoted lea¯ets share one rotational degree of freedom. The extreme positions of the lea¯ets are in e€ect, when the artery is closed or when the angle between the lea¯ets and the symmetry axis goes below a value /0 . The real life diameter H of the tube is taken to be 20 mm, the length W is 60 mm. The blood ¯ow enters the tube along the x-axis with a fully developed parabolic velocity pro®le of a time dependant maximum amplitude shown in Fig. 3 (from Ref. [8]). The out¯ow condition is set to p ˆ p0 and uy ˆ 0 at the outlet, the initial ¯ow ®eld is resting, the same holds for the lea¯ets, which are closed at t ˆ 0. For a physiologically realistic scenario we use upeak ˆ 230 …mm=s†;

m ˆ 3:3 …mm2 =s†;

T0 ˆ

67 …s† 60 …19†

resulting in a peak Reynolds number of about 1400. As we use a uniform grid, the number of nodes for the spatial discretization is quite large (2000  700 nodes) allowing a reasonable resolution of the lea¯et geometry through all stages of the cardiac cycle being simulated.

Fig. 1. Inner and outer wall nodes with corresponding links.

M. Krafczyk et al. / Computers and Structures 79 (2001) 2031±2037

2035

Fig. 2. Sketch of the simulated geometry.

Fig. 3. Transient in¯ow BC imposed by the heart beat used for the simulation.

The temporal resolution of the explicit time integration scheme was 4:5  105 timesteps per cardiac cycle.

4. Results Although the structure of the emerging ¯ow patterns is far from trivial, the simulations here lack spontaneous symmetry breaking known from the 3D case using ®xed lea¯ets [1]. As the symmetry breaking in the ®rst study was induced by 3D ¯ow e€ects, the pertained symmetry of the ¯ow ®eld in the present 2D study appears reasonable. Fig. 4 gives an impression of the ¯ow ®eld structure via the visualization of the magnitude of vorticity for times t 2 f0:0889; 0:222; 0:356; 0:467; 0:622g ‰T0 Š. Additionally, we measured the magnitude of rxy ˆ mq……oux =oy† ‡ …ouy =ox†† and found shear stress peaks in

Fig. 4. Visualizations of the ¯ow ®eld via vorticity magnitude.

the vicinity of the front and back end of the lea¯ets of about 10±15 (N/m2 ) exempli®ed by a contour plot in Fig. 5. The resulting ¯ow dynamics could not be

2036

M. Krafczyk et al. / Computers and Structures 79 (2001) 2031±2037

Fig. 5. Contour plot of jjrxy jj at t ˆ 0:156T0 .

compared to actual experiments in detail, but the resulting shear stresses were of the same order of magnitude as observed in experiments of similar in vitro ¯ow problems (e.g. Refs. [11,12]) and are well below the critical shear stress range to cause necrosis of red blood cells or lethal erythrocyte and thrombocyte damages. An interesting feature of the simulation results is the ¯attening of the lea¯ets due to vortex±lea¯et interactions which can clearly be identi®ed when watching the dynamics evolve during a movie [13].

5. Discussion Our primary goal was to test the suitability of LB methods for a coupled ¯uid±structure ¯ow problem in bioengineering. Despite the neglection of e€ects like artery elasticity e€ects or non-Newtonian stress±strain relationships, the computation presented here using the LB method gave reasonable results for the prediction of velocity- and stress-®elds compared to experimental results. The extension of this work for 3D problems is pricipically straightforward (see e.g. Ref. [15]), but at the present stage poses tremendous CPU time requirements for Reynolds numbers beyond 1000. The main reason is the very high number of timesteps necessary for a stable explicit integration. A detailed understanding of the mechanisms leading to the onset of undesired instabilities is subject of ongoing research. It should be mentioned that the timestepping procedure for a standard LB implementation without ¯uid±structure coupling is much more stable implying substantially smaller computational requirements.

The direct simulation of thrombose generation due to hemodynamical e€ects based on a minimal model to our knowledge ®rst was undertaken in Ref. [14]. Apart from important biological and biochemical modeling problems still being open in this ®eld we feel that at least from a numerical point of view it now seems possible to consider tackling the direct simulation of the genesis of thrombic mass in a system of moving arti®cial heartvalve lea¯ets combining and extending this work and Ref. [14]. Yet we feel that this would require a lot of additional e€orts and the joint cooperation of physicians, engineers and mathematicians. The use of a uniform grid decreases the possibilities of the present algorithm to model the detailed in¯uence of lea¯ets with varying thicknesses and other geometric details which are nevertheless important for the accurate prediction of local peak stress ®elds e.g. in the moment of valve closure. It should be pointed out, that more advanced discretizations could be used to solve Eq. (2). Only fully adaptive schemes controlled by a posteriori error estimation [16] could provide discretizations with better global resolution than the uniform grid used in this study and are subject of ongoing research. References [1] Krafczyk M, Schulz M, Rank E, Cerrolaza M. Analysis of 3D transient blood ¯ow passing through an arti®cial aortic valve by lattice-Boltzmann methods. J Biomech 1998;31: 453±62. [2] Frisch U, D'Humieres D, Hasslacher B, Lallemand P, Pomeau Y, Rivet JP. Lattice gas hydrodynamics in two and three dimensions. Complex Syst 1987;1:649.

M. Krafczyk et al. / Computers and Structures 79 (2001) 2031±2037 [3] Qian YH, D'Humieres D, Lallemand P. Lattice BGK models for Navier±Stokes equation. Europhys Lett 1992;17(6):479±84. [4] Noble DR, Chen S, Georgiadis JG, Buckius RO. A consistent hydrodynamic boundary condition for the LB method. Phys Fluids 1995;7:203. [5] Wall WA, Ramm E. Fluid±structure interaction based upon a stabilized (ALE) ®nite element method. In: Idelson SR, Onate E, Dvorkin EN, editors. Computational Mechanics. Spain: CIMNE; 1998. [6] Ladd AJC. Numerical simulations of particulate suspensions via a discretized Boltzmann equation, Part I. Theoretical foundation. J Fluid Mech 1994;271: 285. [7] Ladd AJC. Numerical simulations of particulate suspensions via a discretized Boltzmann equation, Part II. Numerical results. J Fluid Mech 1994;271:311. [8] Stettler JC, Niederer P, Anliker M, Casty M. Theoretical analysis of arterial hemodynamics. . .. Ann Biomed Eng 1981;9.

2037

[9] Krafczyk M. Gitter-Boltzmann Methoden: von der Theorie zur Anwendung, Professorial Thesis, LS Bauinformatik, Tu Munich, 2001 [in German]. [10] http://www.carbomedics.com. [11] Lee SC, Chandran KB. Numerical simulation of instantaneous back¯ow through central clearance of. . .. Med Biol Eng Comp 1995;33:257. [12] Sakhaeimanesh AA, Morsi Y, Tansley G. In: der Sloten J, et al., editors. Proc 10th Conf Eur Soc Biomech, Leuven, Belgium. 1996. [13] http://www.inf.bauwesen.tu-muenchen.de/start.htm. [14] Krafczyk M, Berrios R. Simulating biological ¯ows using lattice-Boltzmann methods. In: Prado O, Rao M, Cerrolaza M, editors. Simulation con Methods Numericos. London: SVMNI; 1998. [15] Qi D. Non-spheric colloidal suspensions in three-dimensional space. Int J Mod Phys C 1997;8(4):985±97. [16] Becker R, Johnson C, Rannacher R. Adaptive error control for multigrid FE methods. Computing 1995;55: 271±88.