Simulation of a pulsatile total artificial heart: Development of a partitioned Fluid Structure Interaction model

Journal of Fluids and Structures 38 (2013) 187–204

Contents lists available at SciVerse ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Simulation of a pulsatile total artificial heart: Development of a partitioned Fluid Structure Interaction model ¨ Simon J. Sonntag n, Tim A.S. Kaufmann, Martin R. Busen, Marco Laumen, Torsten Linde, Thomas Schmitz-Rode, Ulrich Steinseifer Department of Cardiovascular Engineering, Institute of Applied Medical Engineering, RWTH Aachen University and University Hospital Aachen, Pauwelsstr. 20, 52074 Aachen, Germany

a r t i c l e i n f o

abstract

Article history: Received 25 April 2012 Accepted 24 November 2012 Available online 21 January 2013

Heart disease is one of the leading causes of death in the world. Due to a shortage in donor organs artificial hearts can be a bridge to transplantation or even serve as a destination therapy for patients with terminal heart insufficiency. A pusher plate driven pulsatile membrane pump, the Total Artificial Heart (TAH) ReinHeart, is currently under development at the Institute of Applied Medical Engineering of RWTH Aachen University. This paper presents the methodology of a fully coupled three-dimensional timedependent Fluid Structure Interaction (FSI) simulation of the TAH using a commercial partitioned block-Gauss–Seidel coupling package. Partitioned coupling of the incompressible fluid with the slender flexible membrane as well as a high fluid/structure density ratio of about unity led inherently to a deterioration of the stability (‘artificial added mass instability’). The objective was to conduct a stable simulation with high accuracy of the pumping process. In order to achieve stability, a combined resistance and pressure outlet boundary condition as well as the interface artificial compressibility method was applied. An analysis of the contact algorithm and turbulence condition is presented. Independence tests are performed for the structural and the fluid mesh, the time step size and the number of pulse cycles. Because of the large deformation of the fluid domain, a variable mesh stiffness depending on certain mesh properties was specified for the fluid elements. Adaptive remeshing was avoided. Different approaches for the mesh stiffness function are compared with respect to convergence, preservation of mesh topology and mesh quality. The resulting mesh aspect ratios, mesh expansion factors and mesh orthogonalities are evaluated in detail. The membrane motion and flow distribution of the coupled simulations are compared with a top-view recording and stereo Particle Image Velocimetry (PIV) measurements, respectively, of the actual pump. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Pulsatile total artificial heart Partitioned Fluid Structure Interaction simulation Artificial compressibility method Mesh deformation analysis Experimental validation

1. Introduction Cardiovascular Diseases (CVDs) are still the leading causes of death in the world (Mendis et al., 2011). Heart transplantation is often the last option for people suffering from end-stage heart failure when all other surgical or medical treatments have been exhausted. The main problem facing the field of transplantation is the severe shortage in donor organs. The need for transplantable hearts outweighs the number of organs available. According to Organ Procurement and

n

Corresponding author. Tel.: þ49 241 80 80540; fax: þ 49 241 80 82442. E-mail addresses: [email protected] (S.J. Sonntag), [email protected] (T.A.S. Kaufmann), ¨ [email protected] (M.R. Busen), [email protected] (M. Laumen), [email protected] (T. Linde), [email protected] (T. Schmitz-Rode), [email protected] (U. Steinseifer). 0889-9746/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jfluidstructs.2012.11.011

188

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

Transplantation Network (OPTN) and Scientific Registry of Transplant Recipients (SRTR) in the United States, about 10% of patients died between 2000 and 2009 whilst waiting for a donor heart (Organ Procurement and Transplantation Network, 2011). In 2009 the average waiting time was about eight months. Total Artificial Hearts (TAHs) are fully implantable prosthetic devices replacing the excised ventricles of the diseased native heart both anatomically and functionally. These devices are intended to act as a long-term bridge to transplantation (Copeland et al., 2004) or serve as a destination therapy in case transplantation is not possible. The double chamber pulsatile TAH ReinHeart is currently under development at the Institute of Applied Medical Engineering at the RWTH Aachen University (Finocchiaro et al., 2008). Each of the two pump chambers is compressed alternating by two linear driven plates with concurrent filling of one chamber and ejection from the other. In each chamber the plate pushes a flexible polyurethane membrane, which separates the blood from the motor compartment, into the chamber. As a consequence, the blood is ejected through the bileaflet mechanical outlet valves (St. Jude Medicals RegentTM) which are connected to the aorta and the pulmonary artery, respectively (Fig. 1a). This ejection phase is known as systole. Subsequently the pusher plate moves back to its initial position allowing blood to fill the chambers again through the bileaflet inlet valves (diastolic phase, Fig. 1b).

Fig. 1. (a) During systole a linear driven plate pushes a flexible membrane into the pump chamber. As a consequence, the blood is ejected through the outlet valves. (b) During diastole the pusher plate moves back to its initial position. Now the chamber is gradually filled with blood again through the inlet valves. The upper images show the left heart from the top view. The lower images show both heart chambers from the side view (cross section).

Fig. 2. (a) The structural components of the left heart chamber consist of the pusher plate, the housing and the membrane. The membrane has a uniform thickness of 0.6 mm apart from the marginal area. The pump chamber represents the fluid part of the FSI simulation and (b) side view of the CAD model with dimensions.

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

189

In this study the methodology of the Fluid Structure Interaction (FSI) simulation was developed for the membrane pumping process of the ReinHeart. Therefore, the left heart chamber was separated and the heart valves were neglected for simplicity reasons. The components of the CAD model are illustrated in Fig. 2. The objective was to conduct a stable and accurate simulation of the membrane motion. For this purpose, an extensive examination of the contact algorithm, mesh deformation, turbulence and coupling stability was indispensable. The resulting membrane motion and flow distribution of the FSI simulation are compared with a top-view recording and stereo Particle Image Velocimetry (PIV) measurements, respectively, of the actual ReinHeart. 2. Fluid structure interaction approach Most earlier simulations of pulsatile artificial hearts or pulsatile Ventricular Assist Devices (VADs) have disregarded FSI effects and used separate flow and structural simulations instead (Chesler and Kamm, 1998; Firouzi et al., 2005; Kiris et al., ¨ 1997; Konig et al., 1999; Okamoto et al., 2003; Sato et al., 2009). This was because of the high complexity and computational effort of coupled simulations. A fully coupled method was used by Shim et al. (2003) to investigate the flow in the flexible sac of the Korean artificial heart. The motion of the actuator was simplified by imposing a prescribed pressure on the outer surface of the sac. Yang et al. (2007) studied the flow and diaphragm motion in the pneumatically driven phoenix artificial heart. A commercial Computational Fluid Dynamics (CFD) solver was coupled with an in-house developed Finite Element Method (FEM) code. Avrahami et al. (2006) performed a numerical FSI analysis of the sac-type Berlin VAD. In these latter two studies, external pressure was also imposed directly on the flexible wall. An FSI simulation of a pulsatile ventricular assist device consisting of two fluid chambers separated by a flexible diaphragm was performed by Doyle et al. (2008) for an idealized Brunel University VAD, Zhang and Hisada (2001) for a cylindrical VAD and Moosavi and Fatouraee (2007) for the HeartSaver VAD. All of them considered both the driving fluid and blood dynamics. The study of Moosavi and Fatouraee (2007) was the most complete model among them, however, they only used a two dimensional model. To numerically solve Fluid Structure Interaction problems, mainly two approaches are used. The first one, the monolithic approach, solves the fluid and structure domain simultaneously by discretizing the problem into one system of equations using a single numerical scheme (e.g. FEM). This leads to a stable solution process for most of the cases since the mutual influence of the two fields are incorporated directly. However, for a large three dimensional engineering problem, such as the one presented in this report, a prohibitively amount of memory storage is needed. Moreover, since only one numerical scheme is applied for both physical domains, the applicability of the monolithic approach using commercial packages is limited (Hou et al., 2012). This is because structural problems are generally solved using the FEM and most commercial packages use the Finite Volume Method (FVM) to numerically treat CFD problems. In this study the partitioned approach was used. In this case the fluid and solid domain are separately treated with two distinct solvers. Information between the two systems is communicated across a domain interface. For an ‘implicit’ coupled scheme several sub-iterations (coupling iterations) are performed for each domain per time step until the data converges to the solution of the monolithic system. In the presented study, the fluid has a significant impact on the behavior of the membrane due to the flexibility of the structure. Therefore, a bidirectional data transfer (two-way FSI simulation) was required. As an individual solver for each domain is used, independent numerical schemes as well as more efficient and developed algorithms can be applied to solve the flow and structural equations. Furthermore, less memory storage is required compared to the monolithic approach. However, FSI simulations involving strong couplings such as in the case of the flexible membrane do generally not converge due to stability problems (Felippa et al., 2001; Causin et al., 2005). According to Hou et al. (2012) this hinders the use of the partitioned approach for broad FSI applications. This paper will show that with proper treatment the partitioned approach is well applicable for large engineering problems using commercial packages. The simulations were performed using the multiphysics capabilities of the commercial ANSYSTM 14.0 (Ansys Europe Ltd.) package. The finite volume CFD solver ANSYSTM CFX (Ansys Europe Ltd.) was coupled with the FEM based structural solver ANSYSTM Mechanical (Ansys Europe Ltd.). 3. Numerical set-up 3.1. Material and fluid modeling Structural domain: All the material properties of the solid bodies were taken from the materials data sheets. The density of the thermoplastic polyurethane (TPU) membrane was r ¼1120 kg/m3. An isotropic elastic material model was considered with a Young’s modulus of E¼ 5 MPa and a Poisson’s ratio of n ¼0.45. The stiffness behavior of the membrane was defined as being flexible. The housing and the pusher plate is made of AlCuMgPb Aluminum (linear isotropic) with a Young’s modulus of E¼72000 MPa, a Poisson’s ratio of n ¼0.31 and a density of r ¼2850 kg/m3. An uncoupled simulation with these two bodies designated as being flexible showed a negligible deformation ( o0.001 mm in maximum) of the pusher plate and housing. For this reason, a rigid stiffness behavior was specified for the pusher plate and housing. This reduces each rigid body to a single mass element at its centroid and only the faces in contact have to be meshed. The pump chamber is made of a MABS polymer (r ¼1080 kg/m3, E¼ 2000 MPa, n ¼0.45). In a separate structural simulation of the pump chamber with a defined internal pressure of 150 mmHg a maximal deformation of less than 0.2 mm was calculated.

190

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

Because of this insignificant deformation, the pump chamber was assumed to be stationary in the simulations of this study. Fluid domain: The numerical fluid model was chosen according to the Particle Image Velocimetry (PIV) measurements of the TAH in order to achieve an adequate comparison between the two systems. Thus, for the simulation a homogeneous, Newtonian water–glycerin mixture with a density of r ¼1135.9 kg/m3 and a dynamic viscosity of Z ¼ 3.4566 mPa s was utilized. For further hemodynamic evaluations of the blood pump the physical properties of blood and a non-Newtonian blood model based on Ballyk et al. (1994) will be used.

3.2. Boundary conditions For the sake of consistency, the boundary conditions were specified similar to the experiments. Structural domain: In contrast to other studies, the displacement was not directly applied to the membrane surface, but the pusher plate motion was incorporated. The simulations started with the pusher plate in its lowest position. The membrane was not pre-stressed at this moment. To adapt the displacement of the real pusher plate motion, a smooth periodic function (Fig. 3) was applied to the pusher plate. According to the experiments, a heart rate of 100 beats per minute was specified with a systole to diastole ratio of 1:1. The housing was constrained from all movement with a displacement constraint and the membrane was clamped at its end by a fixed support boundary condition. A contact was defined between the membrane and the housing and between the membrane and the pusher plate (Fig. 4, dotted line). The applied contact algorithm will be discussed in detail in the following section (Section 3.3). For the uncoupled structural simulations of this study a dummy pressure load of 120 mmHg and 15 mmHg were added during systole and diastole, respectively, to the upper surface of the membrane. Fluid domain: The negligence of the heart valves was compensated by switching the conditions at the inlet and outlet boundaries. In the experiments, a systolic afterload of 120 mmHg and an aortic diastolic pressure of 80 mmHg were set (normal physiological conditions) resulting in a peak cardiac output of about Qpe ¼ 20 l/min. Afterload is defined as the systemic vascular resistance to ventricular ejection. In order to model the systolic afterload, a fixed relative pressure

Fig. 3. Displacement profile of the pusher plate for two heart cycles with a maximum displacement of 18 mm. One cardiac cycle lasts for 0.6 s with an equal amount for diastole and systole. This corresponds to 100 beats per minute.

Fig. 4. Boundary conditions of the structural solver: the displacement (double arrow) was applied to the pusher plate. The end of the membrane and the housing were clamped. Contact was defined between the membrane and the pusher plate and between the membrane and the housing (dotted line). The fluid–solid interface was located at the upper surface of the membrane (arrows).

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

191

opening of Prelative ¼80 mmHg and a fluid resistance were combined. The resistance was calculated by using a loss coefficient based on the pressure difference of DP¼40 mmHg (pulse pressure) between systole and diastole. The loss coefficient has previously been applied to model the pressure drop due to cerebrovascular resistance in order to achieve physiological flow conditions in the aorta and its major branches (Benim et al., 2011; Kaufmann et al., 2012). The pressure drop is related to the loss coefficient according to the Darcy–Weisbach equation: 1 f rU 2n , 2

DP ¼

ð1Þ

where f is the non-dimensional loss coefficient, r is the density of the fluid and Un is the normal component of velocity at the outlet. The boundary static pressure Pstatic is specified using (Eq. (2a)) for outflow and (Eq. (2b)) for inflow conditions: Pstatic ¼ P relative þ DP,

ð2aÞ

Pstatic ¼ P relative DP:

ð2bÞ

The experimentally measured peak cardiac output of about Qpe ¼20 l/min was used to determine the loss coefficient empirically. Qpe is defined as Q pe ¼ U pe A ¼ U pe

pd2o 4

,

ð3Þ

where Upe is the flow normal velocity at peak systole, A is the cross-sectional area of the outlet tube and do the outlet diameter. Using a peak cardiac output of Qpe ¼20 l/min ( E3.33  10  4 m3/s) and an outlet diameter of do ¼2.46 cm, Eq. (3) yields an average peak systolic normal velocity of Upe E0.701 m/s at the outlet tube. Using a density of r ¼1135.9 kg/m3, a pressure drop of DP¼ 40 mmHg ( E5332.88 kg/(ms2)) and a normal velocity of Un ¼Upe ¼0.701 m/s the loss coefficient is calculated according to Eq. (1) as f¼

2DP U 2n

r

¼

2U5332:88 kg=m3  2  19:09: 1135:9 kg=m3 0:701 ms

ð4Þ

The inlet was closed during systole. During ventricular diastole, the preload of 15 mmHg was applied to the inlet by using a fixed pressure opening, while the outlet was closed. At the beginning of the simulation the relative pressure was ramped over several time steps to prevent a sudden pressure impact on the initially unstressed membrane. The initial values for velocity in the chamber were set to zero. No-slip boundary conditions were specified for the walls of the pump chamber. Coupling: A fluid–solid interface boundary condition was defined at the upper surface of the membrane in the structural domain (Fig. 4, upper arrows) and the matching surface at the bottom of the pump chamber in the fluid domain. 3.3. Contact algorithm Structural domain: The membrane was not bonded to the pusher plate. This ensured a passive filling of the pump chamber in order to prevent high negative pressures and suction of the atrium during diastole (Kwant et al., 2007). However, this resulted in highly nonlinear contact interaction between the structural components what made an extensive analysis of the contact algorithm necessary. The augmented Lagrange formulation (Laursen, 2002) was used to enforce the frictional contact conditions between the contact pairs. The Coulomb friction law (Laursen and Simo, 1993) was used. The coefficient of friction of polyurethane against aluminum in dry condition was taken as 0.5 according to the technical data sheet. By combining the advantageous features of the penalty method and the Lagrange multiplier method, the augmented Lagrange method has the benefit to be well conditioned, robust and to satisfy the constraints with finite penalties without introducing additional equations (Simo and Laursen, 1992). Furthermore, the augmented Lagrange method was best suited for the presented bending, buckling and sliding dominated problem. For the rigid-to-flexible contacts an asymmetric surface-to-surface behavior was defined with the rigid surfaces designated as the target surfaces and the deformable surfaces designated as the contact surfaces. The contact detection was located at the Gauss integration points, which provide a higher efficiency and accuracy for large deformation problems compared to other approaches (Bathe, 1996). At the beginning of each simulation any initial gaps and penetrations were adjusted between the initially touching elements. The contact searching region was set to a radius of 0.8 mm in order to allocate the right contact and target elements. To achieve a stable and robust convergence additional control of the contact problem behavior was necessary. This is possible in ANSYS with the normal contact stiffness factor FKN (ANSYS Mechanical APDL Contact Technology Guide, 2011). For the augmented Lagrange method the amount of penetration between the contact and target bodies depends on the normal contact stiffness (ANSYS Mechanical APDL Contact Technology Guide, 2011; Laursen, 2002). Its default value is based on the material properties and the size of the underlying elements. FKN acts as a scaling factor of the normal stiffness and therefore regularizes the amount of penetration once contact occurs. A high value of FKN leads to less penetration but can also lead to ill-conditioning of the stiffness matrix what prevents the convergence of the model. For bending dominated problems the normal stiffness is usually reduced to overcome such convergence difficulties while allowing some amount of penetration (ANSYS Mechanical APDL Contact Technology Guide, 2011). In order to find the right

192

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

Table 1 Mesh independence test of the structural domain: weighted mean relative difference (WMRD) between Mesh 1 and Mesh 3 and between Mesh 2 and Mesh 3 for the maximum equivalent stress and the maximum equivalent elastic strain.

Nr. of finite elements Nr. of nodes Maximum equivalent stress WMRD (%) Maximum equivalent elastic strain WMRD (%)

Mesh 1

Mesh 2

Mesh 3

42,117 141,310 5.79 12.52

83,504 273,988 1.58 2.34

185,885 615,785 / /

value with the lowest level of penetration while providing convergence a series of uncoupled structural simulations were conducted with a decreasing stiffness factor (FKNo1.0). With a penalty stiffness factor of 0.7 the penetration was acceptably small (o0.05 mm maximal penetration) and robust convergence under different operating environments was achieved. 3.4. Meshing and mesh independence tests Structural domain: The meshing of the structural parts was performed with ANSYSTM Meshing (Ansys Europe Ltd.). A hexahedral mesh with higher-order 20-node brick elements with uniform reduced integration was used for the membrane. The uniform reduced integration was utilized to prevent volumetric mesh locking effects. At least two layers of elements were used for the membrane to avoid hourglass modes. The contact surfaces of the membrane was meshed using 8-node quadratic quadrilateral contact elements and the target surfaces of the pusher plate and the housing were meshed using 8-node quadratic quadrilateral target elements. A mesh independence study was performed for the solid part in an uncoupled structural simulation. The weighted mean relative difference (WMRD) was applied in this study for all independence tests as a measure of the error. WMRD is defined as the mean difference divided by the mean values (Duffield et al., 2003): P 9g g 0 9 WMRD ¼ P  t t 0 t , ð5Þ t g t þ g t =2 where gt and g0t are scalar quantities of the simulation results for the time step t obtained with two different simulation set-ups. In this case, uncoupled simulations performed with different mesh densities (42,117, 83,504 and 185,885 elements) were compared. The maximum equivalent (von-Mises) stress and the maximum equivalent elastic strain were taken as the scalar quantities. The WMRD was calculated between the coarser meshes and the finest mesh (Table 1). The error between the coarsest and the finest mesh was significantly high for both quantities (5.79% and 12.52%). On the other hand, the second mesh was deemed adequate as there was no significant difference when the mesh density was increased to 185,885 elements (1.58% for maximum stress and 2.34% for maximum strain). Fluid domain: For the fluid domain a triangulated surface mesh and an unstructured tetrahedral volume mesh were generated and optimized using ICEM CFD (Ansys Europe Ltd.). A mesh independence test was carried out by performing FSI simulations with three different types of fluid mesh regarding their density (782,005, 1,806,108 and 4,696,852 elements, see Table 2). Mesh 2 and mesh 3 had a finer mesh at all surfaces of the pump chamber in comparison to the internal mesh, whereas mesh 1 was just refined in the area of the interface. Furthermore, for mesh 2 and mesh 3 a prismatic boundary layer with three and eight layers, respectively, was used to better resolve the flow near the walls. Table 2f illustrates the comparison of the instantaneous velocity profile at the outlet surface during systole (t¼ 0.8 s). Qualitatively, all the meshes share similar flow distributions. Table 2h and i show the WMRD of the maximum values of velocity and the arithmetic average of velocity, respectively, between mesh 1 and mesh 3 as well as between mesh 2 and mesh 3. The error for mesh 1 is significant (8.71% and 21.36%). Therefore, this mesh was not sufficient to accurately resolve the flow field. On the other hand, an error of 3.86% and 6.08% of the maximum values of velocity and of the arithmetic average of velocity, respectively, when using mesh 2 was deemed acceptable. Table 2j indicates the WMRD of the mass flow rate at the outlet. For both coarser meshes the error does not differ substantially and is below 0.4%. This exemplifies that the mass flow rate does not depend significantly on the mesh density. Moreover, the Wall Shear Stress (WSS) on the membrane was compared for the three meshes. Table 2g shows instantaneous WSS profiles for the three meshes at t ¼0.8 s. The computation without boundary layers, mesh 1, resulted in a low resolution and a significantly differing stress distribution. The WMRD between mesh 1 and mesh 3 was over 50% (Table 2k). Mesh 2 shows the same stress distribution compared to mesh 3, however, the WSS is significantly underestimated (Table 2g). The WMRD between mesh 2 and mesh 3 is still over 30% (Table 2k). Thus, an adequate mesh refinement in the boundary region is urgently necessary in order to accurately resolve the near-wall flow region. In conclusion, for an analysis of the flow distribution distal the chamber walls, mesh 2 is sufficient. For an analysis of the near-wall flow field mesh 3 should be used to make sure that the boundary layer is effectively treated. For comparison, Table 3 shows the used meshes of the fluid domain for similar cases. All of these studies used a coarser mesh without treatment of the wall proximity flow field. Zhang and Hisada (2001) admitted that the mesh density was insufficient to

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

193

Table 2 Mesh independence test of the fluid domain: weighted mean relative difference (WMRD) between Mesh 1 and Mesh 3 and between Mesh 2 and Mesh 3 for several flow quantities.

(a) (b) (c) (d) (e)

Nr. of finite volumes Nr. of nodes Nr. of prism. bd. layers Internal mesh size (mm) Surface mesh size (mm)

Mesh 1

Mesh 2

Mesh 3

782,005 155,188 0 2 0.5 at interface 2 at other walls

1,806,108 577,565 3 2

4,696,852 1,507,046 8 1

0.5

0.45

8.71 21.36 0.37 52.89

3.86 6.08 0.29 32.15

/ / / /

(f)

at outlet surface, t ¼0.8 s

(g)

at interface, t¼ 0.8 s

(h) (i) (j) (k)

Maximum Values of velocity WMRD (%) Arithmetic average of velocity WMRD (%) Mass flow rate at outlet WMRD (%) Maximum WSS at interface WMRD (%)

Table 3 Comparison of the fluid mesh density of previously performed studies. FE ¼finite element mesh, FV ¼ finite volume mesh, tet¼ tetrahedral, hex¼ hexahedral, n/s ¼not specified. Reference

Simulation type

Number of fluid elements/nodes

Mesh type

¨ Konig et al. (1999) Zhang and Hisada (2001) Shim et al. (2003) Avrahami et al. (2006) Okamoto et al. (2003) Yang et al. (2007) Moosavi and Fatouraee (2007) Doyle et al. (2008) Sato et al. (2009)

CFD FSI FSI FSI CFD FSI 2D FSI FSI CFD

46,720 elements 4296 elements 41,091 elements 16,250 nodes 100,000 elements 178,630 elements 37,299 elements n/s 117,750 elements

block FV FE tet FE hex FV tet FE tet FV 3-node triangular FE 4-node tet mini elements FE

accurately resolve all effects. Avrahami et al. (2006) used a coarser mesh for the FSI simulation compared to an uncoupled ¨ analysis to reduce the computation time. Konig et al. (1999), Moosavi and Fatouraee (2007) and Yang et al. (2007) claimed mesh independence. Yang et al. (2007), however, performed the independence test using the mass flow rate as the criterion. Our study has shown that other flow quantities should be used instead when performing a mesh independence test.

3.5. Mesh deformation Fluid domain: As a consequence of the structural interaction the mesh of the fluid domain is deformed. In ANSYSTM CFX (Ansys Europe Ltd.) this is obtained with a displacement diffusion model. On receiving the structural displacements of the membrane the volume nodes inside the fluid domain are diffused in location by solving the mesh displacement equations

194

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

(ANSYS CFX-Solver, 2011):   rU Gdisp rd ¼ 0,

ð6Þ 2

where d denotes the displacement relative to the previous mesh distribution and Gdisp (unit: [m /s]) is the mesh diffusivity. The latter parameter is also called mesh stiffness because it determines the stiffness of volume elements to change in shape and dimension. Adaptive remeshing of the fluid domain was avoided in this study because the remeshing procedure is computationally heavy. However, to maintain mesh quality and to prevent inversion of fluid elements in regions of large deformation, a function depending on the volume of the finite volumes (volcvol), the mesh aspect ratio and the distance from the nearest wall boundary (wall distance) was implemented for the mesh stiffness Gdisp:

Gdisp ¼

Aspect Ratio2 m3 1 m3 þ : 1=3 maxð0:000001 m, wall distanceÞ s s volcvol

ð7Þ

By increasing the stiffness of certain elements in relation to other elements a deformation of these is prevented. Eq. (7) causes small elements and elements which have already underwent significant changes in shape and dimension to become stiffer than the general mesh. Furthermore, by considering the wall distance in Eq. (7), the prismatic boundary layer mesh can be maintained. Because of the large deformations, a variable mesh stiffness was crucial in this study in order to preserve the mesh quality and the mesh distribution. Different approaches of the mesh stiffness function were compared and evaluated with respect to convergence, preservation of mesh topology and mesh quality (see Section 4.2). Eq. (7) was determined to be best suited for the presented problem. 3.6. Turbulence Fluid domain: Following Bachmann et al. (2000), the Reynolds number, which can be considered as a measure for the turbulence intensity of a flow, within the pumping chamber is defined as Re ¼

8Q r , pZdi

ð8Þ

where Q is the mean cardiac output, r the density and Z the dynamic viscosity of the fluid and di the diameter of the inlet duct. The TAH has a minimal cardiac output of Q¼4.5 l/min. Using Eq. (8) with Q¼ 4.5 l/min ( ¼7.5  10  5 m3/s), di ¼2.3 cm, r ¼1135.9 kg/m3 and Z ¼ 3.4566 mPa s ( ¼3.4566  10  3 kg/(m s)), the Reynolds number is about Re ¼2729. ¨ According to Konig et al. (1999), the sudden inflow into the pump chamber during diastole causes transition to turbulence even at Reynolds numbers below 800. Therefore, a laminar flow condition could not be assumed. The Reynolds Averaged Navier-Stokes (RANS) Shear Stress Transport (SST) k–o turbulence model by Menter (1994) with an enhanced wall function was chosen for this study. This model is well suited for such low Reynolds number flows without the need to use additional damping functions. 3.7. Structural and fluid solver Structural domain: ANSYSTM Mechanical uses the Finite Element Method (FEM) to discretize the underlying transient dynamic equations (Wilson, 2004) and the associated boundary conditions. Large deformation effects were included in the analysis taking stiffness changes into account due to geometric changes during deformation. To solve the resulting nonlinear discretized system of equations the full unsymmetric Newton–Raphson solution procedure with a sparse direct solver was applied. The unsymmetric approach is recommended for problems involving a coefficient of friction higher than 0.2 to achieve a better convergence behavior (ANSYS Mechanical APDL Contact Technology Guide, 2011). The direct solver ran in the in-core memory mode to achieve the best CPU performance. Iterative solvers turned out to be inefficient for this case. The FEM simulation was running on two cores of an 8-core Intels CoreTM i7-950 3.07 GHz CPU workstation with 24 GB of RAM. To prevent artificial oscillations of the membrane during the initial loading steps a Rayleigh damping for stiffness was applied for this period (Wilson, 2004). Oscillations occur as a consequence of the material’s elasticity. At the very least a small damping was required for numerical stability. However, the damping should be small enough to not affect the accuracy of the results. A series of uncoupled structural simulations were conducted with increasing damping. The optimal damping was achieved using a stiffness proportional Rayleigh damping coefficient of b ¼0.001 for the viscous damping matrix D¼ bK, where K is the system structural stiffness matrix. Fluid domain: ANSYSTM CFX uses the Finite Volume Method (FVM) to numerically solve the unsteady three-dimensional Navier–Stokes equations. To accommodate moving boundaries, the Arbitrary Lagrangian–Eulerian (ALE) formulation was employed (e.g., Donea et al. (2004)). Advection fluxes were evaluated by using the bounded high-resolution scheme based on Barth and Jesperson (1989). This approach automatically blends the 1st order upwind scheme and the 2nd order central differencing scheme depending on the local flow field. The second-order accurate and implicit backward Euler scheme

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

195

(ANSYS CFX-Solver, 2011) was utilized for the temporal discretization of the transient terms. The root mean square normalized values of the governing equation residuals were set to a target of 10  4. The CFD part was using six cores. 3.8. Fluid–solid coupling Coupling: Total forces and total mesh displacements are transferred across the fluid–solid interfaces. Because the driving force in this simulation is the pusher plate, it is more efficient to choose the mechanical solver to run first within each coupling iteration (ANSYS Mechanical APDL Coupled-Field Analysis Guide, 2011). In particular, in one coupling iteration the total mesh displacement data is calculated first in the ANSYSTM Mechanical solver and subsequently sent to the fluid solver. Before solving the transport equations, the mesh displacement equations (Eq. (6)) are solved in ANSYSTM CFX and the mesh coordinates are updated based on the resulting displacements. When the fluid solver has converged, the calculated forces on the interface are sent back to the structural solver. The next coupling iteration then commences. Because of the dissimilar meshes of the fluid and solid interface an interpolation of the transferred data is necessary. The globally conservative CPP interpolation method was used for the force transfer and the profile preserving NONC interpolation method was used for the displacement transfer (ANSYS Mechanical Coupled-Field, 2011) due to the rougher structural interface mesh compared to the fluid interface mesh. This way the force balance over the surface is preserved (ANSYS Mechanical APDL Coupled-Field Analysis Guide, 2011). To map a node from one grid to an element face on the other grid the bucket search algorithm (Jansen et al., 1992) was applied. Simulations were performed using double precision of calculations because of the large range in scales of mesh deformation and simulation quantities. The calculation proceeds to the next time step when the changes in transferred forces and displacements is smaller than the specified convergence target of 5  10  3 (block-Gauss–Seidel scheme). 3.9. Artificial compressibility method ¨ Coupling: It has previously been shown (Causin et al., 2005; Forster et al., 2007; Mok and Wall, 2001; Wall et al., 1999) that partitioned coupling of an incompressible flow with a slender flexible structure, like the membrane, can lead to a deterioration of the stability properties. This is known as the ‘artificial added mass instability’. The instability increases ¨ when the interaction between the fluid and the structure is strong due to a high fluid/structure mass density ratio (Forster et al., 2007), as was the case in this study (fluid/structure ratio of 1.01). As previously mentioned, in every time step data is transferred between the two solvers until equilibrium is reached. Despite the block-Gauss–Seidel scheme is implicit, the fluid iterations and the structural iterations are solved explicitly. This means that changes of displacements are not considered during the fluid iterations and forces acting on the interface are not adapted during the structural iterations. Because of the incompressibility of the fluid, disturbances of predicted interface displacements cannot be damped during the fluid iterations what results in an acceleration of the fluid. As a consequence large pressure gradients arise in the flow domain (Degroote et al., 2010). This load acts instantaneously as ¨ a tremendous artificial mass on the membrane what yields instabilities. Studies have shown (Forster et al., 2007) that the instability even occurs earlier by increasing the temporal accuracy, what is consistent with our observations. However, a small time step was necessary to accurately resolve the fluid flow and structural behavior and to reduce other numerical problems. One way to handle this instability is to introduce some kind of implicitness to the solution process of the two systems. ¨ One option is to heavily reduce the under-relaxation factors of the transferred data (Causin et al., 2005; Wuchner et al., 2007). These factors ensure that in each coupling iteration only a percentage of the calculated forces and displacements are transferred to the other solver. By damping down the solution process, this may stabilize the system, however, the computational efficiency is seriously limited due to the additional coupling iterations needed. Furthermore, in several cases it is not possible to achieve convergence of the interface loads even with very low under-relaxation factors. In this study, an artificial compressibility was applied to the fluid problem adjacent to the membrane in order to enhance stability (Degroote et al., 2010, 2008; Raback et al., 2001). The implicitness is introduced by adding a source term to the continuity equation of the interface boundary. By mitigating the incompressibility constraint at the interface, the pressure response is slowed down. This can be understood as introducing an indirect under-relaxation. However, the under-relaxation factors for all the transferred data were set to 1, so that no explicit under-relaxation was performed resulting in an efficient solution process. The effect of the artificial compressibility vanishes when the coupling iterations have converged (Degroote et al., 2010). Therefore, the accuracy of the solution is not affected, only the way in which convergence is achieved. 3.10. Time step independence test Structural domain: An uncoupled structural simulation was first performed to check for convergence of the structural part. In this analysis an initial time step size of 10  3 s was chosen. Subsequently, the time step size was decreased until convergence could be achieved. This was the case using a time step size below 5  10  4 s. Therefore, to keep the structural part of the coupled simulation stable a maximum coupling step size as low as 5  10  4 s was necessary.

196

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

Table 4 Time step independence test: weighted mean relative difference (WMRD) of several scalar quantities for FSI simulations performed with a time step size of 5  10  4 s and 1.25  10  4 s. Simulation quantity

WMRD (%)

Arithmetic average of velocity (fluid domain) Maximum values of velocity (fluid domain) Mass flow at outlet (fluid domain) Velocity magnitude at pump chamber roof (fluid domain) Maximum equivalent stress (structural domain) Maximum elastic equivalent strain (structural domain)

1.78 1.63 1.54 1.87 0.92 0.22

Table 5 Pulse cycle independence test: weighted mean relative difference (WMRD) of several scalar flow quantities for FSI simulations performed with different numbers of cardiac cycles (e.g. 2nd/3rd is the difference between the second and the third heart cycle). Simulation quantity

WMRD 2nd/3rd (%)

WMRD 2nd/4th (%)

WMRD 2nd/5th (%)

WMRD 3rd/4th (%)

WMRD 3rd/5th (%)

WMRD 4th/5th (%)

Arithmetic velocity average Maximum values of velocity Mass flow at outlet Velocity magnitude at roof

2.45 2.62 3.84 21.21

2.92 4.35 7.01 29.06

3.97 4.82 10.32 32.78

1.91 4.01 3.58 24.44

2.47 4.05 6.60 27.05

1.39 1.68 3.55 7.31

Coupling: To check for time step independence, FSI simulations with a time step size of 5  10  4 s (1200 time steps per cycle) and 1.25  10  4 s (4800 time steps per cycle) were performed. Fluid mesh 2 (Table 2) was used for this study. Using a finer mesh with such a low time step size of 1.25  10  4 s resulted in prohibitively high computation time. The computation time for two pump cycles with a time step size of 5  10  4 s by using fluid mesh 2 and fluid mesh 3 was about 4 days and 16 days, respectively. The weighted mean relative difference (Eq. (5)) was used as a measure of the error between the simulation performed with a time step size of 5  10  4 s and 1.25  10  4 s (Table 4). For the fluid domain the WMRD of the arithmetic average of velocity, the maximum values of velocity, the mass flow during systole at the outlet and the velocity magnitudes at the central point of the pump chamber roof were evaluated. The WMRD of the maximum values of equivalent (von-Mises) stress and equivalent elastic strain were taken as a measure for the structural domain. The error was less than 2% for all quantities. Therefore, using a time step size of 5  10  4 s was sufficient to provide accurate solutions. For comparison, Avrahami et al. (2006) reported time step independence using 1000 time steps per cycle. 3.11. Pulse cycle independence test Coupling: Previous studies of pulsatile blood pumps have shown that conventionally three cycles are sufficient to achieve time-periodic independence (Avrahami et al., 2006; Firouzi et al., 2005; Moosavi and Fatouraee, 2007). To verify a steady repeated pulsating flow domain, five cycles were simulated. Fluid mesh 2 (Table 2) was used again for this study. The weighted mean relative difference (WMRD) was taken as a measure to examine cycle-to-cycle variations. The same flow quantities as for the time independence test (Section 3.10) were compared (arithmetic average of velocity, maximum values of velocity, mass flow at outlet, velocity magnitude at pump chamber roof). The WMRD between the second and the third, the fourth and the fifth cycle was calculated, as well as the difference between the third and the fourth cycle, between the third and the fifth cycle and between the fourth and the fifth cycle (Table 5). The error is increasing when comparing the second cycle to the third, the fourth and the fifth cycle. Nevertheless, the maximum error of the arithmetic velocity average (3.97%) and the maximum values of velocity (4.82%) are still in an acceptable range. On the other hand, the maximum error of the mass flow at the outlet (10.32%) and the velocity magnitude at the chamber roof (32.78%) are significant. The third cycle does not show a significant improvement to the second cycle when looking at the difference to the fifth cycle. The comparison of the fourth and the fifth cycle delivers the best results (WMRDo7.4% for all flow quantities). Therefore, at least the fourth cycle should be used for analyses. However, when using the finest fluid mesh in this study, we only simulated two heart cycles to make the simulation computationally feasible. 4. Results 4.1. Coupling stability Partitioned block-Gauss–Seidel coupling of the incompressible fluid with the slender flexible membrane and a high fluid/structure mass density ratio of about unity led inherently to a deterioration of the stability. The combination of fixed pressure boundary condition and resistance at the outlet boundary during systole showed better results concerning

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

197

stability compared to a pressure boundary condition. Nevertheless, using the resistance outlet without additional stabilizing treatment convergence was only achieved using a time step size above 5  10  3 s. As the time step size was decreased, the coupling started to become instable again. However, due to the highly dynamic behavior of the membrane during diastole a step size as low as 5  10  4 s was at least necessary to achieve convergence of the structural part. Therefore, additional stabilizing treatment was needed for the presented simulation to obtain convergence. Using the artificial compressibility approach of Section 3.9, a fully stable simulation could be achieved, even at very low time steps. 4.2. Mesh deformation Mesh deformation is a very important aspect for solving FSI problems. Large deflections of the membrane were involved in this simulation which led to large mesh deformations inside the pump chamber. An end diastolic volume of EDV¼ 98.4 ml and an end systolic volume of ESV¼45.0 ml were calculated giving a stroke volume of SV¼EDV– ESV ¼53.4 ml, a mean cardiac output of Q¼SV  Heart rate¼5.34 l/min and an ejection fraction of EF ¼SV/EDV¼54.3%. The height from the upper side of the membrane at its centroid to the chamber ceiling is reduced by 95.2% (Fig. 5). As a consequence, an adequate deformation of the mesh was very important. As mentioned in Section 3.5, time-consuming adaptive remeshing was avoided. Instead, a variable mesh stiffness depending on certain mesh properties was specified for the fluid domain. Different approaches for the mesh stiffness function were compared and evaluated with respect to convergence, preservation of mesh topology and mesh quality. The finer the mesh of the fluid domain, the more difficult it is to preserve these conditions. Therefore, the finest fluid mesh (mesh 3, Table 2) was used for the calculations of this analysis. Furthermore, the first cardiac cycle was selected for the evaluations. Table 6 shows the resulting mesh at t¼0.4 s for simulations performed with two different mesh stiffness functions (the region of the shown mesh is illustrated in Fig. 5 with a rectangle). Without incorporating the wall distance in the mesh stiffness function (Table 6a), the prismatic boundary layer could not be preserved during diastole. The reason was that the deformed elements beyond the boundary layer had relatively a higher aspect ratio and lower volume compared to the prismatic elements. Therefore, the inner mesh became stiffer than the boundary region resulting in a severe deformation of the prismatic elements. However, not incorporating the aspect ratio and volume of the finite volumes (volcvol) led to divergence. On this account, the conventionally available mesh stiffness options in ANSYS CFX (Gdisp ¼(1/volcvol)n[m(2 þ 3n)/s], Gdisp ¼(1/wall distance)n[m(2 þ n)/s], with a specified exponent n) resulted in divergence during systole caused by excessive distortion of certain elements. On the  other hand, the mesh stiffness function Gdisp ¼ ðAspect Ratio2 =maxð0:000001 m, wall distanceÞ Þ m3 =s was appropriate for the diastole. However, not considering the volume of the finite volumes resulted in negative volume elements at end systole. Incorporating additional 1/volcvol [m5/s] resulted in very stiff elements in the inner mesh during diastole (due to the low value of volcvol) what caused deformation of the boundary layer again. An addition of 1/volcvol1/3 [m3/s] was determined adequate. For these reasons, mesh stiffness function (7) (Table 6b), which is considering the volume of the finite volumes, the wall distance

Fig. 5. Total mesh displacement in one cross section (dashed line in Fig. 9a) when the pusher plate is in its lowest (a, t ¼0 s) and highest (b, t ¼ 0.3 s) position. The rectangle illustrates the region of the mesh shown in Table 6.

198

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

Table 6 Comparison of the resulting fluid mesh at t ¼0.4 s for simulations performed with two different mesh stiffness functions (a and b). The region of the shown mesh is illustrated in Fig. 5 with a rectangle. Mess stiffness function (a)

Gdisp ¼

Aspect Ratio volcvol

(b)

Gdisp ¼

Aspect Ratio2 m3 maxð0:000001 m, wall distanceÞ s

2

Fluid mesh at t¼ 0.4 s

m5 s

þ

1 m3 1=3 s volcvol

and the mesh aspect ratio, was best suited for the presented problem. All the mesh properties stayed in an acceptable range and no excessive mesh distortion was observed. Furthermore, the prismatic boundary layer mesh was maintained. At the end of diastole (t¼0.6 s and t¼ 1.2 s) the mesh distribution almost returned to the initial configuration. In the following, the resulting mesh aspect ratios, mesh expansion factors and mesh orthogonalities when using mesh stiffness function (7) are presented. These measures are most relevant for mesh quality to the CFX-Solver (ANSYS CFX-Solver, 2011). Mesh aspect ratio: The aspect ratio is defined as the ratio of the maximum to minimum cell faces for all elements adjacent to a node. The maximum mesh aspect ratio remained smaller than 90 for the whole calculation. According to the ANSYSTM Users Guide (ANSYS CFX-Solver, 2011), the desirable limit is 100 and the acceptable limit is 1000 when running in double precision. Values above this limit lead to round-off errors and convergence problems. Fig. 6 shows a diagram of the ratio (%) of the fluid domain volume and the volume of elements with an aspect ratio420, 440 and 460 over time (t ¼0.2 s to t ¼0.4 s). For t ¼0.3 s, when the mesh deformation is maximal, the ratio of the fluid domain volume and the volume of elements with an aspect ratio 420,440 and 460 was 1.79%, 0.048% and 0.00036%, respectively. Therefore, only a very small number of highly stretched elements were present. Elements with an aspect ratio 420 occurred solely in the central area above the membrane, where the mesh is pressed together the most. In this region, the velocity gradients were low. The overall maximum velocity for elements with an aspect ratio 460 was 0.41 m/s. Furthermore, elements with an aspect ratio 440 were present only from about t ¼0.27 s to t ¼0.33 s. Mesh expansion: The mesh expansion factor is defined as the ratio of the largest to the smallest element volumes surrounding a node. The desirable range is o20 (ANSYS CFX-Solver, 2011). Mesh expansion factors exceeding this limit introduce discretization errors of transient and body force terms. A maximal mesh expansion factor of 53.4 could be observed. However, a mesh expansion factor 440 was only present for one element pair (1.42  10  7% of the whole domain). The maximal ratio of the fluid domain volume and the volume of elements with a mesh expansion factor 410, 415 and 420 was 0.046%, 0.0034% and 0.00045%, respectively. Consequently, no significant mesh expansion was present. The minimal volume of tetrahedral elements dropped from t ¼0 s to t ¼0.3 s by 35.8%. Mesh orthogonality: The mesh orthogonality is a measure of the conformity of the angles between adjacent element faces and the optimal angle (901 for quadrilateral faced elements and 601 for triangular faces elements). The orthogonality angle is defined as the area weighted average of 901-arccos(ns) (quadrilateral faced elements) and 601-arccos(ns) (triangular faced elements) for all integration point surfaces of a control volume (ANSYS CFX-Solver, 2011), where s is the

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

199

Fig. 6. Ratio (%) of the fluid domain volume and the volume of elements with an aspect ratio 420, 440 and 460 over time (t ¼0.2 s to t ¼ 0.4 s).

Fig. 7. Ratio (%) of the fluid domain volume and the volume of elements with an orthogonality angle o201, o151 and o101 over time (t¼0.2 s to t¼ 0.4 s).

vector that joins the centers of the adjacent control volume nodes and n is the associated face normal vector. The orthogonality angle should be 4 201 (ANSYS CFX-Solver, 2011). Discretization errors are introduced and increased when a high amount of elements below this range are present. This can result in poor convergence. An overall minimum orthogonality angle of 4.11 was measured. Fig. 7 shows a diagram of the ratio (%) of the fluid domain volume and the volume of elements with an orthogonality angle o201, o151 and o101 over time (t ¼0.2 s to t ¼0.4 s). Equal to the aspect ratio, elements with an orthogonality angle o201 occurred solely in the central area above the membrane. For t ¼0.3 s the ratio of the fluid domain volume and the volume of elements with an orthogonality angle o201, o151, o101and o71 was 1.89%, 0.91%, 0.0069% and 0.00011%, respectively. Therefore, some elements with angular deviation are present during end systole and early diastole. But the amount of elements with significant non-orthogonality is very low. 4.3. Membrane motion The accurate deformation of the membrane was crucial for this study. For this reason, the resulting membrane motion of the FSI simulation was compared to recordings of the actual pump. However, a detailed, quantitative validation of the

200

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

Fig. 8. Top view of the FSI simulation (left) and experimental set-up (right) with the pusher plate in the upper position (t ¼0.9 s). The circles highlight the dents of the membrane.

Fig. 9. Total mesh displacement of the membrane during early diastole (t¼ 0.98 s): (a) isometric view and (b) front view of one cross section (represented by a dashed line in a). A separation and buckling of the membrane distal the tubes can be observed.

membrane motion for the whole pump cycle was limited. This is because of the restricted field of view of the membrane in the experimental set-up. When viewed from the side, the shape of the membrane cannot be estimated due to the housing. Therefore, only a limited, qualitative comparison of the membrane motion was possible. An experimental study of a sacdriven blood pump with full optical access of the membrane has been performed by Hochareon et al. (2003). The finest fluid mesh (mesh 3, Table 2) and the second cardiac cycle were used for comparison. Fig. 8 shows a top view of both settings with the pusher plate in the upper most position (t¼0.9 s). A sinusoidal buckling of the membrane can be observed in both images. The amount and location of the dents are almost the same (highlighted with circles in Fig. 8). Furthermore, in the experiment we could see a sudden separation of the membrane from the pusher plate during early diastole due to the passive filling process. The membrane clung again back to the plate at about mid-to-late diastole. This separation (mainly occurring distal of the tubes) was also observed with the FSI simulations (Fig. 9).

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

201

Fig. 10. (a) Plane position of the FSI simulation and Particle Image Velocimetry (PIV) measurements. (b), (c), (d) Comparison of the FSI simulation (left) and PIV (right) velocities in one plane (a) for different time steps (b: t¼ 2.5 s, c: t ¼2.84 s and d: t ¼2.92 s).

202

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

4.4. Flow distribution Flow distribution is a very important aspect for a medical device when looking at blood flow related problems like presence of recirculation areas, blood damage and washout behavior. For this reason, the flow distribution was compared to stereo Particle Image Velocimetry (Stereo-PIV, e.g. Visscher et al., 2011) measurements in one plane (Fig. 10) in order to validate the results. Because the near-wall flow field cannot be compared using the PIV measurements, fluid mesh 2 (Table 2) was used for the FSI simulation. The fifth cardiac cycle was utilized for the evaluations. Areas containing no data can be found in the PIV measurements at the rear part of the outlet tube and close to the inlet. This is the result of shadowing effects. In the mid-to-late systolic phase the displacement of the membrane induces high velocities in the outlet tube (Fig. 10b). Distal to the outlet (back corner of the chamber) the velocities are relatively low. In this area the velocity is slightly overestimated in the simulation compared to the PIV measurements. A peak cardiac output of 16.97 l/min (t ¼2.5 s) was calculated using the FSI simulation. Thus, the experimentally measured peak cardiac output of about 20 l/min was slightly underestimated. Furthermore, a pressure of 116.2 mmHg was reached in the chamber at the peak systole of the simulation. At the end of systole and beginning of diastole, when the pusher plate has reached its maximum displacement, the velocity in the chamber decreases. This is due to the stagnation of the pusher plate. Because of the relaxation of the membrane deformation during mid-diastole, the flow velocity in the inlet tube becomes strong and the chamber is gradually filled with the fluid again. The comparison with the PIV measurements shows a significant difference at the inflow (Fig. 10c). Due to the bileaflet heart valve, three inlet jets are observed in the PIV measurement. On the other hand, a uniform inflow is present in the simulation. The inflow is also too strong because of the missing resistance of the valves. Therefore, the heart valves have a significant impact on the flow field. However, the resulting main vortexes in the pump chamber at mid-to-late diastole (Fig. 10c and d) as well as the velocity magnitudes are in good correlation. 5. Discussion and conclusion In this study a three-dimensional two-way Fluid Structure Interaction (FSI) simulation of the pumping process of the Aachen Total Artificial Heart (TAH) was successfully performed using a commercial partitioned block-Gauss–Seidel coupling package. Since structural problems are generally solved using the Finite Element Method (FEM) and most commercial packages use the Finite Volume Method (FVM) to numerically treat CFD problems, the applicability of the monolithic approach using commercial packages is very difficult. In contrast, in a partitioned simulation existing flow solvers and structural solvers can be used. However, the so-called ‘artificial added mass instability’ limits the general use of the partitioned approach for strongly coupled FSI problems. This paper presented that with proper treatment a stable partitioned FSI simulation can be achieved even for large engineering problems. The combination of fixed pressure boundary condition and resistance at the outlet boundary during systole improved the stability compared to a pressure boundary condition. This boundary condition is also physiologically more realistic than a simple outlet pressure. However, additional stabilizing treatment was needed for the presented simulation to obtain convergence. Time consuming underrelaxation was avoided. A much more efficient approach to achieve stability was to apply an artificial compressibility to the fluid domain adjacent to the membrane. The resulting membrane motion of the FSI simulation corresponded well with recordings of the actual pump. However, only a qualitative comparison of the membrane behavior was possible due to the restricted field of view in the experimental set-up. This restricted optical access is a limiting disadvantage of the experiment. On the other hand, by using the simulation results, the behavior and stress distribution of the membrane can be evaluated in detail. Mesh independence tests were performed both for the structural and fluid domain. The independence test for the fluid domain exemplified that the mass flow rate is not a suitable criteria for mesh-convergence of the simulation results, because this value does not depend significantly on the mesh density. Adaptive remeshing of the fluid mesh was not considered in this study because of its additional computational effort. However, an adequate mesh deformation is very important when large deflections of the solid part are involved. Such FSI cases often result in a termination of the calculation caused by negative volume elements. Adapting the mesh stiffness in critical areas turned out to be a good method to preserve the mesh quality and the mesh distribution without remeshing. On the other hand, the appropriate mesh stiffness function has to be determined for each FSI case individually depending on the specific mesh distribution and mesh displacement. Using a mesh stiffness function which is considering the volume of the finite volumes, the wall distance and the mesh aspect ratio, no significant mesh expansion as well as mesh stretching (bad aspect ratios) was present. Furthermore, the prismatic boundary layer was maintained. However, some elements with angular deviation were present during end systole and early diastole. But the amount of elements with significant non-orthogonality was very low. Moreover, in this region the velocity gradients were low. The mesh orthogonality was very critical concerning convergence. The conventional approaches for the mesh stiffness function available in ANSYS CFX resulted in very low orthogonality angles what led immediately to the termination of the simulation process. A time step independence test and pump cycle independence test was only performed using a coarse fluid mesh. A finer mesh should be used instead. However, using a finer fluid mesh and a low time step size of 1.25  10  4 s resulted in a prohibitively high computation time. This made the simulation within reasonable CPU time infeasible. Moreover, despite admitting errors of the mass flow rate and local velocity fields, when using the finest fluid mesh the second cardiac cycle

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

203

was utilized for the analyses of this study. The pulse cycle independence test (performed with a coarser mesh) showed that at least four cardiac cycles should be simulated in order to obtain a time periodic solution. Due to a Reynolds number above 2500 and a sudden expansion of the flow at the inlet tube during diastole a laminar ¨ condition cannot be assumed. However, despite admitting implications for accuracy, Konig et al. (1999), Yang et al. (2007), Avrahami et al. (2006) and Moosavi and Fatouraee (2007) assumed the flow in the pump chamber to be laminar in order to save computation time. In this study the Reynolds Averaged Navier-Stokes (RANS) Shear Stress Transport (SST) k–o turbulence model with an enhanced wall function was used. More advanced and scale resolving turbulence models, like the Large Eddy Simulation (LES) or the Delayed Detached Eddy Simulation (DDES), might be better suited for such an intermittent turbulent flow (Guilmineau et al., 2011; Mittal et al., 2001; Wang et al., 2012). However, these models require a considerably finer mesh and consequently prohibitively more CPU time and memory compared to RANS models. This is not yet applicable for the presented FSI case. In this study the bileaflet mechanical heart valves at the inlet and outlet were neglected. Instead a simple switching on/off boundary condition was applied. First comparisons with Particle Image Velocimetry (PIV) measurements showed a significant difference at the inflow. The main vortexes and the general flow characteristics in the pump chamber, however, are already in good correlation. Therefore, we can assume that the resulting membrane motion is not affected significantly by the negligence of the valves. Nevertheless, in order to completely validate the numerical model, the heart valves have to be considered. Separate FSI simulations of the heart valves using the Immersed Boundary Method (IBM) have already been performed successfully. However, the standalone FSI simulations (using a simple pipe geometry) are already computationally very time consuming. Therefore, incorporating a fully coupled FSI of both valves (inlet and outlet) in the TAH simulation will probably be beyond our computational capabilities. Besides the challenges of simulating the valve dynamics this is one reason why no simulations have yet been reported that incorporate two fully coupled passively moving valves in a 3D model of a pulsatile flow pump chamber simulation. Most of the studies used a switching open/ ¨ closed boundary condition (Konig et al., 1999; Shim et al., 2003; Yang et al., 2007), as was done in this study. Avrahami et al. (2006) performed an FSI simulation of a sac-type Ventricular Assist Device where the valves were modeled as either fully open or fully closed. This simplification at least provides a more realistic inflow profile than neglecting the contribution of the valves entirely. Kiris et al. (1997) considered the valve motion. However, the motion of the valves was predefined based on experimental data. Both options will also be considered. Moosavi and Fatouraee (2007) is the only reported study of a membrane-type assist device that included a fully coupled motion of the valves. Though, they only used a two dimensional model. The future objective is to perform design optimization of the TAH using the FSI simulations. The advanced engineering concepts and hemodynamic requirements of medical devices demand new tools to support the development process. Computational analysis promises lower cost and time compared to experimental approaches. A design optimization would address the localization of high shear stress and stagnation regions which may have adverse hemodynamic effects. The detection of areas prone to hemolysis and thrombus formation as well as a proper washout behavior is very important when designing these devices. However, due to the high computation time the feasibility of the FSI simulations as a design tool for the TAH is questionable. Despite trying to reduce the computation time where numerically possible and already making several simplifications, the computation time is still too high to effectively assist the design process. For this reason, further simplifications or assumptions are necessary to make, for example, Design of Experiments (DOE) practicable. This could address the assumption of a laminar flow condition, a coarser mesh density, the negligence of the pusher plate and housing, a predefined membrane motion, etc. Therefore, before performing a design optimization, a study is necessary to examine what simplifications have significant impact on the accuracy and what simplifications can be assumed without considerable loss of accuracy.

Acknowledgments This work was funded by the European Union (EU) and the state of North Rhine-Westphalia. References ANSYS CFX-Solver Modeling Guide, ANSYS Release 14.0, ANSYS Inc. ANSYS Mechanical APDL Contact Technology Guide, 2011. ANSYS Release 14.0, ANSYS Inc. ANSYS Mechanical APDL Coupled-Field Analysis Guide, 2011. ANSYS Release 14.0, ANSYS Inc. Avrahami, I., Rosenfeld, M., Raz, S., Einav, S., 2006. Numerical model of flow in a sac-type ventricular assist device. Artificial Organs 30, 529–538. Bachmann, C., Hugo, G., Rosenberg, G., Deutsch, S., Frontaine, A., Tarbell, J.M., 2000. Fluid dynamics of a pediatric ventricular assist device. Artificial Organs 24, 362–372. Ballyk, P.D., Steinman, D.A., Ethier, C.R., 1994. Simulation of non-Newtonian blood flow in an end-to-side anastomosis. Biorheology 31 (5), 565–586. Barth, T.J., Jesperson, D.C., 1989. The design and application of upwind schemes on unstructured meshes. AIAA Paper 89-0366, pp. 1–12. Bathe, K.J., 1996. Finite Element Procedures. Prentice-Hall, Englewood Cliffs, NJ. Benim, A.C., Nahavandi, A., Assmann, A., Schubert, D., Feindt, P., Suh, S.H., 2011. Simulation of blood flow in human aorta with emphasis on outlet boundary conditions. Applied Mathematical Modelling 35 (7), 3175–3188. Causin, P., Gerbeau, J.-F., Nobile, F., 2005. Added-mass effect in the design of partitioned algorithms for fluid–structure problems. Computer Methods in Applied Mechanics and Engineering 194, 4506–4527. Chesler, N.C., Kamm, R.D., 1998. Performance analysis of a cardiac assist device in counterpulsation. Journal of Biomechanical Engineering 120, 437–445.

204

S.J. Sonntag et al. / Journal of Fluids and Structures 38 (2013) 187–204

Copeland, J.G., Smith, R.G., Arabia, F.A., Nolan, P.E., Sethi, G.K., Tsau, P.H., McClellan, D., Slepian, M.J., 2004. Cardiac replacement with a total artificial heart as a bridge to transplantation. The New England Journal of Medicine 351, 859–867. Degroote, J., Swillens, A., Bruggeman, P., Haelterman, R., Segers, P., Vierendeels, J., 2010. Simulation of fluid–structure interaction with the interface artificial compressibility method. International Journal for Numerical Methods in Biomedical Engineering 26, 276–289. Degroote, J., Bruggeman, P., Haelterman, R., Vierendeels, J., 2008. Stability of a coupling technique for partitioned solvers in FSI applications. Composite Structures 86, 2224–2234. Donea, J., Huerta, A., Ponthot, J.-P., Rodriguez-Ferran, A., 2004. Arbitrary Lagrangian–Eulerian methods, in: Stein, E., de Borst, R., Hughes, T.J.R. (Eds.), Encyclopedia of Computational Mechanics. vol. 1. John Wiley & Sons Ltd, New York, pp. 1–25. Fundamentals. Doyle, M.G., Vergniaud, J.-B., Tavoularis, S., Bourgault, Y., 2008. Numerical simulations of blood flow in artificial and natural hearts with fluid–structure interaction. Artificial Organs 32 (11), 870–879. Duffield, N., Lund, C., Thorup, M., 2003. Estimating flow distributions from sampled flow statistics. Proceedings of the ACM SIGCOMM, 325–336. Felippa, C.A., Park, K.C., Farhat, C., 2001. Partitioned analysis of coupled mechanical systems. Computer Methods in Applied Mechanics and Engineering 190, 3247–3270. Finocchiaro, T., Butschen, T, Kwant, P., Steinseifer, U., Schmitz-Rode, T., Hameyer, K., Leßmann, M., 2008. New linear motor concepts for artificial hearts. IEEE Transactions on Magnetics 44 (6), 678–681. Firouzi, F., Fatouraee, N., Najarian, S., 2005. Simulation of blood flow in the sac-type ventricular assist device using computational fluid dynamic. Iranian Journal of Biomedical Engineering 1 (2), 129–142. ¨ Forster, C., Wall, W.A., Ramm, E., 2007. Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows. Computer Methods in Applied Mechanics and Engineering 196 (7), 1278–1293. Guilmineau, E., Deng, G., Wackers, J., 2011. Numerical simulation with a DES approach for automotive flows. Journal of Fluids and Structures 27, 807–816. Hochareon, P., Manning, K.B., Fontaine, A.A., Deutsch, S., Tarbell, J.M., 2003. Diaphragm motion affects flow patterns in an artificial heart. Artificial Organs 27 (12), 1102–1109. Hou, G., Wang, J., Layton, A., 2012. Numerical methods for fluid–structure interaction—a review. Communications in Computational Physics 12, 337–377. Jansen, K., Shakib, F., Hughes, T., 1992. Fast projection algorithm for unstructured meshes, in: Atluri, S. (Ed.), Computational Nonlinear Mechanics in Aerospace Engineering, American Institute of Aeronautics and Astronautics. American Institute of Aeronautics & Astronautics, Georgia Inst. of Technology, Atlanta. Kaufmann, T.A., Schmitz-Rode, T., Steinseifer, U., 2012. Implementation of cerebral autoregulation into computational fluid dynamics studies of cardiopulmonary bypass. Artificial Organs 36 (8), 754–758. Kiris, C., Kwak, D., Rogers, S., Chang, I.D., 1997. Computational approach for probing the flow through artificial heart devices. Journal of Biomechanical Engineering 19 (4), 452–460. ¨ Konig, C.S., Clark, C., Mokhtarzadeh-Dehghan, M.R., 1999. Investigation of unsteady flow in a model of a ventricular assist device by numerical modeling and comparison with experiment. Medical Engineering and Physics 21 (1), 53–64. ¨ ¨ Kwant, P.B., Finocchiaro, T., Forster, F., Reul, H., Rau, G., Morshuis, M., El Banayosi, A., Korfer, R., Schmitz-Rode, T., Steinseifer, U., 2007. The MiniACcor: constructive redesign of an implantable total artificial heart, initial laboratory testing and further steps. International Journal of Artificial Organs 30 (4), 345–351. Laursen, T.A., 2002. Computational Contact and Impact Mechanics. Springer-Verlag, Berlin. Laursen, T.A., Simo, J.C., 1993. Algorithmic symmetrization of Coulomb frictional problems using augmented Lagrangians. Computer Methods in Applied Mechanics and Engineering 108 (1–2), 133–146. Mendis, S., Puska, P., Norrving, B., 2011. Global Atlas on Cardiovascular Disease Prevention and Control. World Health Organization, Geneva. Menter, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal 32 (8), 1598–1605. Mittal, R., Simmons, S.P., Udaykumar, H.S., 2001. Application of large-Eddy simulation to the study of pulsatile flow in a modeled arterial stenosis. Journal of Biomechanical Engineering 123, 325–332. Mok, D.P., Wall, W.A., 2001. Partitioned analysis schemes for the transient interaction of incompressible flows and nonlinear flexible structures, in: Wall, W.A., Bletzinger, K.-U., Schweitzerhof, K. (Eds.), Proceedings of Trends in Computational Structural Mechanics. CIMNE, Barcelona, Spain, pp. 689–698. Moosavi, M.H., Fatouraee, N., 2007. Flow simulation of a diaphragm-type ventricular assist device with structural interactions. Conference Proceedings IEEE Engineering in Medicine and Biology Society, pp. 1027–1030. Okamoto, E., Hashimoto, T., Mitamura, Y., 2003. Design of a miniature implantable left ventricular assist device using CAD/CAM technology. Artificial Organs 6, 162–167. Organ Procurement and Transplantation Network (OPTN) and Scientific Registry of Transplant Recipients (SRTR), 2011. OPTN/SRTR 2010 Annual Data Report. Rockville, MD: Department of Health and Human Services. Health Resources and Services Administration. Healthcare Systems Bureau, Division of Transplantation. ¨ Raback, P., Ruokolainen, J., Lyly, M., Jarvinen, E., 2001. Fluid–structure interaction boundary conditions by artificial compressibility. ECCOMAS Computational Fluid Dynamics Conference, Swansea, Wales. Sato, K., Orihashi, K., Kurosaki, T., Tokumine, A., Fukunaga, S., Ninomiya, S., 2009. Analysis of flow patterns in a ventricular assist device: a comparative study of particle image velocimetry and computational fluid dynamics. Artificial Organs 33, 352–359. Shim, E.B., Yeo, J.Y., Ko, H.J., Youn, C.H., Lee, Y.R., Park, C.Y., Min, B.G., Sun, K., 2003. Numerical Analysis of the three-dimensional blood flow in the Korean artificial heart. Artificial Organs 27 (1), 49–60. Simo, J.C., Laursen, T.A., 1992. An augmented Lagrangian treatment of contact problems involving friction. Composite Structures 42 (1), 97–116. Visscher, J., Pettersen, B., Andersson, H.I., 2011. Experimental study on the wake behind tapered circular cylinders. Journal of Fluids and Structure 27, 1228–1237. Wall, W.A., Mok, D.P., Ramm, E., 1999. Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures, in: Wunderlich, W. (Ed.), Solids, Structures and Coupled Problems in Engineering. Proceedings of the European Conference on Computer Mechanics, Munich. Wang, S., Ingham, D.B., Ma, L., Pourkashanian, M., Tao, Z., 2012. Turbulence modeling of deep dynamic stall at relatively low Reynolds number. Journal of Fluids and Structures 33, 191–209. Wilson, E.L., 2004. Static and Dynamic Analysis of Structures. Section: Stiffness and Mass Proportional Damping, 4th ed. Computers and Structures, Inc., Berkeley, CA, p. 231. ¨ Wuchner, R., Kupzok, A., Bletzinger, K.U., 2007. A framework for stabilized partitioned analysis of thin membrane–wind interaction. International Journal for Numerical Methods in Fluids 54, 945–963. Yang, X.L., Liu, Y., Yang, J.M., 2007. Unsteady flow and diaphragm motion in total artificial heart. Journal of Mechanical Science and Technology 21, 1869–1875. Zhang, Q., Hisada, T., 2001. Analysis of fluid–structure interaction problems with structural buckling and large domain changes by ALE finite element method. Computer Methods in Applied Mechanics 190, 6341–6357.