当前位置：文档库 › Large eddy simulation of powered Fontan hemodynamics

Large eddy simulation of powered Fontan hemodynamics

Y.Delorme a,n,K.Anupindi a,A.E.Kerlo a,D.Shetty a,M.Rodefeld b,J.Chen a,S.Frankel a

a School of Mechanical Engineering,Purdue University,Lafayette,IN,United States

b Department of Surgery,Indiana University School of Medicine,Indianapolis,IN,United States

a r t i c l e i n f o

Article history:

Accepted26October2012

Keywords:

High order large eddy simulation Fontan circulation

Viscous impeller pump a b s t r a c t

Children born with univentricular heart disease typically must undergo three open heart surgeries within the?rst2–3years of life to eventually establish the Fontan circulation.In that case the single working ventricle pumps oxygenated blood to the body and blood returns to the lungs?owing passively through the Total Cavopulmonary Connection(TCPC)rather than being actively pumped by a subpulmonary ventricle.The TCPC is a direct surgical connection between the superior and inferior vena cava and the left and right pulmonary arteries.We have postulated that a mechanical pump inserted into this circulation providing a3–5mmHg pressure augmentation will reestablish bi-ventricular physiology serving as a bridge-to-recovery,bridge-to-transplant or destination therapy as a ‘‘biventricular Fontan’’circulation.The Viscous Impeller Pump(VIP)has been proposed by our group as such an assist device.It is situated in the center of the4-way TCPC intersection and spins pulling blood from the vena cavae and pushing it into the pulmonary arteries.We hypothesized that Large Eddy Simulation(LES)using high-order numerical methods are needed to capture unsteady powered and unpowered Fontan hemodynamics.Inclusion of a mechanical pump into the CFD further complicates matters due to the need to account for rotating machinery.In this study,we focus on predictions from an in-house high-order LES code(WenoHemo TM)for unpowered and VIP-powered idealized TCPC hemodynamics with quantitative comparisons to Stereoscopic Particle Imaging Velocimetry(SPIV) measurements.Results are presented for both instantaneous?ow structures and statistical data. Simulations show good qualitative and quantitative agreement with measured data.

1.Introduction

Compared to normal circulation,children born with single ventricle heart disease typically undergo a series of three staged open heart surgeries within the?rst2–3years of life to eventually live with a Fontan circulation(Gewillig,2005).The anatomy of the Fontan circulation involves a direct connection between the Superior and Inferior Vena Cavae(SVC and IVC)and the Right and Left Pulmonary Arteries(RPA and LPA)forming the Total Cavo-pulmonary Connection or TCPC.Effectively,these patients have a single working ventricle which pumps oxygenated blood to the body and deoxygenated blood passively?ows through the TCPC to the lungs in series,but without a subpulmonary power source.It is an inherently inef?cient circulation which is asso-ciated with elevated systemic venous pressures and reduced ventricular?lling(Deleval,1998).The TCPC hemodynamics effec-tively involve two con?ned impinging jets and complex unsteady chaotic vortical?ow patterns resulting in undesirable pressure loss and irregular streaming.For more than the past decade researchers have focused on reducing pressure losses at the TCPC junction by passive geometric means such as introducing offset (Migliavacca et al.,2003),or splitting the SVC and/or IVC (Soerensen et al.,2007;Marsden et al.,2009).Those solutions prevent the two jets coming from the IVC and SVC from colliding. By doing so,the competition between the two jets is avoided,and so are the strong pressure losses.In contrast,our group has focused on a cavopulmonary assist device,which consists in inserting a mechanical pump into the Fontan circulation to provide the estimated3–5mmHg needed to reestablish normal biventricular pressure levels as a means to bridge-to-recovery or bridge-to-transplant(Rodefeld et al.,2003).Previous designs have been proposed,from axial blood pumps placed in the vena cavae(Rodefeld et al.,2003)to a folding propeller design (Throckmorton et al.,2007).Our current design has evolved substantially from prior device concepts,is based on the von Karman viscous pump principle and involves spinning a double-sided cone featuring6mild vanes at the center of the TCPC junction(Rodefeld et al.,2010).The pump design is motivated by the desire for a percutaneous expandable rotary blood pump.This pump could be used to reduce the number of surgeries required for the survival of the infants from three to two or even one.In addition,this would minimize the Fontan-related disease caused by the current staged procedures while older.This pump would also be used in older patients showing signs of heart failure.Placing a right-sided power source in the univentricular Fontan circulation will

Contents lists available at SciVerse ScienceDirect

journal homepage:https://www.wendangku.net/doc/eb6322710.html,/locate/jbiomech

https://www.wendangku.net/doc/eb6322710.html,

Journal of Biomechanics

n Corresponding author.Tel.:t176********.

E-mail addresses:delorme.yann@https://www.wendangku.net/doc/eb6322710.html,,

ydelorme@https://www.wendangku.net/doc/eb6322710.html,(Y.Delorme).

Journal of Biomechanics46(2013)408–422

restore conditions closer to a normal two-ventricle circulation.This will ease the symptoms of Fontan failure and allow the physician to take care of the patient as a‘‘biventricular Fontan’’.Preliminary performance,biocomptability and CFD studies of a rigid prototype-have demonstrated the feasibility of the Viscous Impeller Pump (VIP)(Kennington et al.,2011;Giridharan et al.,in press).

Computational?uid dynamics(CFD)of pathological and medical device hemodynamics has become very popular over the past decade due to the widespread availability of medical images, associated processing software,and open-source and commercial CFD packages.Many of these studies assume laminar?ow or use Reynolds Averaged Navier–Stokes(RANS)based turbulence models for predictions.This is potentially a problem when dealing with pathological and medical device hemodynamics due to the?ow regime being transitional or low Reynolds number turbulence,and featuring?ow curvature and rotation.Traditional RANS-based turbulence models are known to have accuracy issues in dealing with these?ow features.These issues were recently brought out in paper reporting on an FDA-approved computational and experi-mental study of an idealized medical device(Stewart et al.,2012). In addition,blood damage prediction models typically require accurate predictions of instantaneous wall shear stress and?ow residence time,and this information is typically not available in a RANS simulation.Other simulations of Fontan circulation involved the use of low order?nite element methods allowing an easier coupling between?uid and solid domains(Marsden et al.,2009).

In this study,the focus is on Large Eddy Simulation(LES)and its potential to predict pathological and medical device hemody-namics.In LES the large-scale unsteady three-dimensional?ow features are numerically resolved and the effect of the unresolved small-scale eddies on the large-scales is modeled using a Subgrid-Scale(SGS)model(Pope,2000).While it is debatable,there is a general consensus in the literature that LES calls for the use of a high-order numerical method to accurately capture the still rela-tively wide-range of large eddy length scales.The use of high-order numerical methods often makes application to complex geometries, especially those involving patient-speci?c geometries and rotating machinery,dif?cult.Also,accurate SGS models that can handle transitional?ows are still a challenge.Furthermore,validation of LES predictions through quantitative comparisons to measured velocity pro?les,and not just more global measures like pressure drop or pump performance,is critical.

Our approach is based on the use of a?nite-difference method using a structured Cartesian grid combined with a version of the immersed boundary method for handling complex geometries and rotating machinery.In the next section,the computational method is described in detail.This is followed by a brief descrip-tion of the experimental method used to obtain the measured data for validation.More details on this experiments and a full description of the measurements appear in a separate publication (Kerlo et al.,2011).Following,LES results for unpowered and powered Fontan hemodynamics are presented with comparisons to measured data.A summary concludes the paper.

2.Methods

2.1.LES equations

The?ltered incompressible Navier–Stokes equations solved in our in-house code(WenoHemo TM)are

@u i

tu j

@u i

j ?à

1@2u i

j j

@t ij

e1T

where u i is the i th component of the spatially?ltered velocity vector,p is the

spatially?ltered pressure,x i is the i th component of the spatial domain and t ij is

the sub-grid stress tensor which arises from the?ltering of the equations,and is

de?ned as

t ij?u i u jàu i u jà2u j u jàu j u j

àád

e3T

All the parameters(velocities,spatial variables,y)are dimensionless.In the case

of the TCPC simulations,the velocity is non dimensionalized using the inlet

velocity U inlet and the physical dimensions are non-dimensionalized using the inlet

diameter D.

The Vreman SGS turbulence model is employed in this study.The SGS stress

tensor t ij is related to the?ltered strain rate tensor S ij using an eddy viscosity n T,

where

t ij?à2n T S ije4T

n T?C?Meu!T

S ij?

@u i

@x j

@u j

@x i

>:e5T

The model for the kernel Meu

T,originally proposed by Vreman(2004)and tested

in the present code by Shetty et al.(2010),is robust and easy to implement.The

kernel is a function of the components of the velocity gradient tensor.The

constant C arises from the Vreman model and its value can vary as a function of

the case studied.In our case,the value of the constant used is0.07.More details on

the eddy viscosity formula are provided in the two papers cited above(Vreman,

2004;Shetty et al.,2010).Due to the form of the model kernel,the eddy viscosity

goes to zero in laminar regions,as well as in near wall regions.These are

important features of a SGS model for transitional wall-bounded?ows such as

those studied herein.

2.2.Numerical methods

In order to solve Eqs.(1)and(2),a predictor/corrector scheme is used.During

the?rst stage,a predicted value of the velocity?eld u n is obtained by solving

Eq.(1)without the pressure term.The pressure?eld is then obtained by solving

the following Poisson equation:

@2p

i i

@u n

e6T

The velocity?eld is?nally corrected to obtain u nt1

,which now respects the

continuity equation.

u nt1

?u n

àD t

@x i

e7T

The time integration is performed using a third order accurate backward

difference scheme(Shetty et al.,2011).

The convective terms in Eq.(1)are discretized using a?fth-order accurate

Weighted Essentially Non-Oscillatory(WENO)scheme.The WENO scheme uses

adaptive stencils to balance the desire for high-order accuracy with a non-

oscillatory solution near discontinuities.This scheme uses upwinding based on

the sign of the velocity multiplying the spatial derivative and chooses the

smoothest stencil by avoiding the discontinuities in the interpolation process.

The detailed of the scheme can be found in Zhang’s paper(Zhang and Jackson,

2009).The computation of the viscous terms and SGS terms involves the use of

?rst-order and second-order derivatives,which are discretized using fourth-order

central?nite-difference schemes.

2.3.Immersed boundary method

The use of the above high-order numerical method is facilitated by using a

?xed structured staggered Cartesian grid(Shetty et al.,2010).In order to account

for?ow over or through complex shaped rigid or rotating bodies,an Immersed

Boundary Method(IBM)is employed.Speci?cally,a mirroring of the?uid?ow

inside the solid body is used(Mark and van Wachem,2008).The IBM requires a

triangulated surface mesh for the objects we wish to?ow over or through.The

procedure begins with a computation of the surface normal for each surface point.

Next,the structured Cartesian grid points are separated into two categories:inside

or outside of the immersed body(IB)or bodies.For the solid points,ghost points

(GPs)are de?ned as those having a neighbor inside the?uid region.The ghost

points are projected into the?uid region,normal to the surface of the IB(using the

previously computed normals).This results in a boundary point(BP)on the

Y.Delorme et al./Journal of Biomechanics46(2013)408–422409

interpolating the velocity from the neighbor points using the following equations:a ?x i t1àx IP x i t1àx i b ?y j t1ày IP

j t1j g ?z k t1àz IP k t1

8>>>>>>><>>>>>>>:e8Tu !IP ?a b g u !i ,j ,k te1àa Tb g u !

i t1,j ,k

ta e1àb Tg u !i ,j t1,k ta b e1àg Tu !

i ,j ,k t1

te1àa Te1àb Tg u !

i t1,j t1,k

te1àa Tb e1àg Tu !

i t1,j ,k t1

ta e1àb Te1àg Tu !

i ,j t1,k t1

te1àa Te1àb Te1àg Tu !

i t1,j t1,k t1

e9T

In the above equation,x n ,y n and z n correspond to the coordinates of the 8points surrounding the IPs and x IP ,y IP and z IP correspond to the IP coordinates.The value of the velocity at the GP location is then obtained using u !GP ?2u !BP àu !IP

e10T

This method has been shown to be second order accurate.2.4.Grid stretching

In order to cluster as many points as possible in the ?uid region,we are using an analytical grid clustering formula that allows to stretch the Cartesian mesh at one particular location in each direction.The LES computations are then performed in computational space (x ,Z ,z )instead of the physical space (x ,y ,z ).The grid stretching expression used here is N ?1texp ew Tà1àá?

L e11TD ?1texp eàw Tà1àá?

L e12T

g ?12w ln

N D

e13T

x ?x min t1tsinh ew ex àg TTw g

?c y ?y min t1tsinh ew eZ àg TT

8>>>>>>>>

where w controls the strength of the clustering,c controls the location of the stretching,L is the length of the domain in the stretching direction,and x min ,y min and z min are the minimum values of the domain along the x ,y and z directions.The values of the derivatives in the computational space are then computed using the chain rule (only the conversion from x space to x space is presented)@@x ?@@x @x @x e15T

@2?@2

@x 2

t@x @2x

e16T

2.5.Experimental method

Experiments to measure the ?ow ?eld associated with powered Fontan

hemodynamics using the VIP are also carried out and serve as further validation of the LES predictions.Stereoscopic Particle Image Velocimetry (SPIV)is employed in an in vitro mock ?ow loop using an index-matched blood analog ?uid and the same idealized TCPC model as in the LES to measure the three velocity components within two-dimensional planes at different locations.The experimental setup enables the study of the ?ow characteristics in a horizontal plane near and away from the VIP along the LPA or RPA.Hollow glass beads are uniformly mixed to the blood analog ?uid as seeding particles.A dual-head Nd:YAG pulse laser (Quantel Twins BSL140)periodically illuminates the central horizontal plane of the test section (y –z

plane).

Fig.1.Schematic of the Immersed Boundary Method:the shaded region corre-sponds to the solid body.The ‘‘mirroring’’characteristic of the method corre-sponds here to a zero velocity boundary condition at the surface of the

body.Fig.2.Presentation of the idealized TCPC.(a)Schematic of the Idealized TCPC with

Y.Delorme et al./Journal of Biomechanics 46(2013)408–422

410

The laser beam is directed toward the area of interest by forming1mm-thick laser sheets through a combination of optical lenses.The two CCD cameras required for the SPIV measurements point towards the area of interest at tilt angle of a?301.The cameras and the laser are synchronized using a multi-channel pulse generator (Quantum Composers).SPIV data are obtained at the central horizontal plane and for two consecutive downstream locations along the x-axis of the RPA.At each location 8000snapshots of particle images are recorded from each camera,which gives once processed4000instantaneous velocity?eld enabling statistical analysis.

3.Results

3.1.Simulation details

The hybrid high-order?nite-difference LES/IBM solver allows both internal and external?ow simulations.For the VIP-powered Fontan hemodynamics studied here,both types of?ow simulations are present at the same time with?ow through the TCPC and over the rotating VIP.For the idealized TCPC geometry and the VIP pump head, a surface mesh was created using triangular elements in Gambit(AN-SYS).These mesh(es)are then read by the solver to save each nodes location and compute the normal vector for each triangular element.

At the inlet(s)a velocity pro?le is imposed to match the desired?ow rate.In the case of Fontan hemodynamics,a para-bolic velocity pro?le is used to match the laminar?ow imposed at the entrance of the TCPC in the experiments.The mean velocity is computed using the?ow rate and the imposed velocity pro?le is uerT?2u mean1à

e17

TY.Delorme et al./Journal of Biomechanics46(2013)408–422411

where R is the radius of the Vena Cava (VC).At the outlet,a constant value of 10is imposed for the pressure and a zero-gradient boundary condition is imposed for the velocity assuming fully developed ?ow @u

@n outlet ?0e18T

At the inlet(s)and at the wall(s),a zero gradient boundary condition is imposed for the pressure.If the wall is an IB,the values of the GPs automatically enforce the zero gradient boundary condition for the pressure at the wall.If the wall boundary condition is imposed at a ?ow boundary of the domain,the Poisson solver (Adams,1999)imposes the desired boundary condition automatically.

In order to impose the desired boundary condition for the velocity at the edges of the domain,three layers of points are

added outside of the domain in every directions.The values of the extra points are set to satisfy the boundary conditions.Let us consider the upper boundary of the computational domain and N be the number of grid points in a particular direction:in that case the values of the points outside the domain are set as follow:No slip u N ti ?àu N ài i ?1,2or 3Free slip u N ti ?u N ài i ?1,2or 3In?ow u N ti ?U inflow i ?1,2or 3Out?ow u N ti ?u N à1i ?1,2or 3Periodic

u N ti ?u i

i ?1,

2or 3

The code is parallelized using OpenMP (OMP,2012)and simula-tions are typically run using 8processors over a month for the no VIP case and the static VIP case,and two months for the

rotating

Y.Delorme et al./Journal of Biomechanics 46(2013)408–422

412

VIP cases.This is the time needed to obtain converged statistics to compare to experimental measurements.

3.2.Validation

In order to validate the IBM,the classical case of?ow over a sphere is considered.At low Re number based on the sphere diameter,?ow over a sphere results in the formation of a wake with steady vortices.The location of the center of those vortices and the length of the recirculation bubble are directly related to the Re number.Simulations for four different Re number(50,100, 200and200)are considered and compared to the experimental measurements reported by Johnson and Patel(1999).The com-putational domain is25sphere diameters in the?ow direction and8sphere diameters in the other two directions.A Cartesian grid with384points in the?ow direction and128points in the other two directions is used.Results show very good agreement between the simulations and the measurements for the location of the vortices and the length of the recirculation zone.

3.3.LES of Fontan hemodynamics

The current state-of-the-art in simulating Fontan hemody-namics involves the use of patient-speci?c geometries(Zelicourt et al.,2010).In this study,an idealized TCPC geometry,shown in Fig.2-a,is employed in both the experiments and simulations. The two inlet pipes of inner diameter22mm correspond to the SVC and IVC and the two outlet pipes of inner diameter18mm correspond to the RPA and LPA.The computational domain is a rectangular box that is15inlet diameters along the Pulmonary Arteries(PA)or x-direction and5inlet diameters along the VC or z-direction.The Cartesian grid used is a512?64?128with stretching along the x and z directions to cluster the points towards the center of the domain(?uid region).The stretching coef?cient w is chosen as3,resulting in a clustering of65

points Y.Delorme et al./Journal of Biomechanics46(2013)408–422413

Fig.6.Illustration of the swirl switching phenomena.(a)Evolution of the vorticity at the junction of the TCPC as a function of time.The red line represents the time averaged value.(b)Streamlines colored by vorticity magnitude at the plane x ?0,t ?3.5s.(c)Streamlines colored by vorticity magnitude at the plane x ?0,t ?8.0s.(d)Streamlines colored by vorticity magnitude at the plane x ?0,t ?9.0s.(e)Streamlines colored by vorticity magnitude at the plane x ?0,t ?32.0s.(f)Streamlines colored by vorticity magnitude at the plane x ?0,t ?33.0s.(g)Streamlines colored by vorticity magnitude at the plane x ?0,t ?34.0s.(For interpretation of the references to color in this ?gure caption,the reader is referred to the web version of this paper.)

Y.Delorme et al./Journal of Biomechanics 46(2013)408–422

414

Y.Delorme et al./Journal of Biomechanics46(2013)408–422415

Fig.7.Powered Fontan hemodynamics:static VIP.(a)Iso-surface of l2(l2?à25)colored by vorticity magnitude.(b)Velocity vectors and vorticity magnitude,x/D?0.0.

(c)Velocity vectors and vorticity magnitude,x/D?1.5.(d)Velocity vectors and vorticity magnitude,x/D?3.0.(e)Velocity vectors and vorticity magnitude,x/D?4.5.

(f)Velocity vectors and vorticity magnitude,x/D?3.0,t?74s.(g)Schematic highlighting v and w patterns for static VIP case.(For interpretation of the references to color

across the diameter of the VCs and55points across the diameter of the PAs.Due to the particular geometry of the TCPC,45%of the grid points are located outside the?uid region.

At each inlet,a?ow rate of2.2L minà1is imposed,corre-sponding to a mean velocity of0.0964m sà1.The value of the density and viscosity is chosen to match those used in the experiments and are1060kg mà3and3.5cP,respectively,which are close to physiological values.The Re number based on the inlet diameter and?ow rate is about700.At the outlets an out?ow boundary condition for the u-component of the velocity and a free slip boundary condition for the v and w-components of the velocity are imposed(see Section3.1).For the pressure, a constant value of10is imposed at the outlets.Along the y-direction,a no-slip boundary condition is used for the velocity and a homogeneous Neumann boundary condition for the pres-sure is imposed.The time step is D t n?10à3.This value is chosen to ensure numerical stability and gives a Courant–Friedrichs–Lewy(CFL)number of0.15.

In order to understand the complex?ow structures associated with Fontan hemodynamics,a case without the VIP is considered ?rst.A snapshot from the simulation showing vortical structures depicted by the negative second eigenvalue of the velocity gradient tensor(Jeong and Hussain,1994)and colored by vorti-city magnitude,is shown in Fig.3-a.Two long vortical structures are clearly visible throughout the PA.Four PA axial slices showing vorticy magnitude contours with velocity vectors superimposed are shown in Fig.3-b to e and clearly con?rm this?ow feature. Animations from this simulation show that these two vortices are unstable and interact with each other to produce a single vortex with signi?cant increase in vorticity magnitude.Evidence for this can be seen from a similar set of plots shown in Fig.4at a later time in the simulation.When the two vortices interact

and

Y.Delorme et al./Journal of Biomechanics46(2013)408–422

416

eventually merge,they produce a single stronger vortex and many additional small scale vortical structures form during the breakup of the single vortex.The presence of the strong vortex enhances the mixing at the connection,which is crucial to ensure a correct hepatic factor distribution between the two lungs.On the other hand,a too powerful vortex might induce high stresses resulting in high blood damage.

To facilitate comparisons to the SPIV measurements,mean (time-averaged over30?ow through times)velocity pro?les are extracted at?ve different locations in the RPA(see Fig.2-b).Those pro?les are extracted at the center plane y?0and at x?0.5,1,1.5, 2and2.5in the RPA,and comparisons between LES and PIV are sown in Fig.5.Excellent quantitative agreement is obtained for the x-and z-components of the velocity.The average deviation from the experiments is0.015m/s for the x-component(0.056m/ s maximal deviation)and0.016m/s for the z-component (0.027m/s maximal deviation).The y-component of the velocity obtained numerically and which is perpendicular to both the PA’s and the VC’s shows a?at pro?le with a near zero value.It appears to not be in agreement with the measured pro?les.This can be explained as follows:animations from these simulations show that the rotational direction of the vortices in the PA changes sign during the simulation.A series of snapshots clearly depicting this phenomena are shown in Fig.6-b to g.This is similar to the phenomena of swirl switching as it relates to the low frequency oscillations of Dean vortices in turbulent pipe bend?ows(Rutten et al.,2005).This effect can also be seen in Fig.6-a showing the time evolution of the x-component of the vorticity spatially averaged over a cube of non-dimensional length0.25centered at the TCPC junction(x?y?z?0)together with its time aver-aged value.The oscillatory sign change of the vorticity is obvious from this plot as well as the near zero average value.The SPIV data that is used to compute the mean velocity pro?les is obtained from2different SPIV windows in the PA and are extracted over different time periods.Hence,Fig.5shows pro?les with clockwise as well as counter-clockwise rotation in qualita-tive agreement with the phenomena observed in the simulations.

3.4.LES of powered Fontan hemodynamics

In this section,LES results are presented with the VIP centered at the TCPC junction.The VIP concept is inspired by the Von Karman viscous pump principle(Panton,2005)and the current design resembles a double-sided cone with six mild

vanes Y.Delorme et al./Journal of Biomechanics46(2013)408–422417

(Kennington et al.,2011).The height of the pump is17.6mm and its maximum diameter is20.3mm.Fig.2shows the position of the pump in the TCPC.In practice,the VIP will be inserted via a catheter through a femoral vein or jugular vein,positioned at the TCPC junction and deployed by a novel expandable design currently under development.The VIP will be powered by an external electromagnetic motor.Kerlo et al.(2011)show the hydraulic performances(difference of pressure between the VCs and the PAs)of the VIP over a range of?ow rates and rotation rates.The VIP provides the desired pressure rise for the selected rotation rates,with a relatively?at pro?le over the range of?ow rates.This predicts stable performance over a wide range of physiological conditions.The operating conditions for the pump are between2000RPM and5000RPM(up to7000RPM if needed).LES results are?rst presented for the not rotating VIP case(static)in order to simulate failure mode and then for the rotating VIP cases https://www.wendangku.net/doc/eb6322710.html,parisons to mea-surements are also presented.The pro?les for the following case are extracted at the center plane y?0and at x?0.8,1.2,1.6and 2to match the locations of extracted data from the experiments.

3.4.1.Static case

Instantaneous?ow features,similar to that already shown in Fig.3,are shown in Fig.7.Four vortical structures are clearly discerned over about half of the extent of the PA’s.These resemble the Dean’s vortices associated with?ow through a pipe bend,as the TCPC and the static VIP form a curved passage from the VC’s to the PA’s.Unlike the no pump case,these vortices are highly stable and animations only show a slight oscillation of the location of the four vortices.The vorticity magnitude is also much lower compared to the no VIP case.

A comparison between the LES predictions and the SPIV measure-ments for the mean velocity pro?les are shown in Fig.8.Agreement is not as good as for the no pump case but still reasonable both qualitatively and quantitatively.The mean deviation is0.03m/s

for

Y.Delorme et al./Journal of Biomechanics46(2013)408–422

418

all components and the averaged maximal deviation is around 0.08m/s.The pattern of the out of plane velocity component can be explained by the fact that,as we move forward in time,the four vortices in the PA slowly rotates but stay stable(Fig.7-f).Fig.7-g shows a schematic of the four vortices,with a black line representing the location where the data are extracted.From this schematic, we can visualize the positive/negative/positive/negative pattern of the y-component of the velocity,as well as the negative/positive pattern of the z-component of the velocity.Fig.8shows that the LES simulation is able to capture qualitatively the?ow physics;only the magnitude of the velocity components at these particular locations are off.

3.4.2.Rotating cases

LES results for powered Fontan hemodynamics are shown in Fig.9.The simulations are run using D t n?5?10à4for the case at 1000RPM and D t n?2:5?10à4for the case at2000RPM to capture accurately the motion of the pump and maintain a low CFL number for stability(0.3).Two small vortical structures?lling half of the PA on both sides of the TCPC junction are observed. Also,smaller discrete vortical structure emanating from the impel-ler are observed and are thought to be associated with shedding from the VIP blades.Those small structures are being convected downstream of the pump by the mean?ow.The?ow near the pump is highly turbulent with a high vorticity magnitude.Further downstream,free?ow interactions(vortex shedding)appear.A comparison between the LES predictions and the SPIV measure-ments for the mean velocity pro?les are shown in Figs.10and11. Reasonable qualitative and quantitative agreement is observed.The mean deviation is about0.045m/s for both cases,with an averaged maximal deviation of0.15m/s.Both LES and experiments show a peak of the axial velocity at the centerline that can be explained by the geometry of the VIP.At1000RPM,we observe no recirculating ?ow near the pump(no negative axial velocity pro?les).The effect of the pump on the?ow?eld dissipates quicker for the

experiments Y.Delorme et al./Journal of Biomechanics46(2013)408–422419

than for the LES.Indeed,the amplitudes of the oscillations decrease more as we move in the PA for the SPIV measurements than for the LES simulation.On the other hand,at2000RPM,LES is not capable of predicting the recirculation close to the pump that is measured by the experiments.Even though LES and experiments show reasonable agreements,we have to move away from the pump to observe a better match between numerical and experimental results.This may be due to the resolution of the grid that seems to be not high enough to fully capture the geometry and the motion of the VIP.

3.5.Flow and wall shear stress

Flow and wall shear stresses(WSS)associated with medical device hemodynamics are important as they related to blood damage predictions especially hemolysis and thrombosis.High stress levels can damage red blood cells and result in hemolysis and low stress levels are often associated with?ow status and thrombosis.In this section,wall shear stress on the TCPC vessel walls and scalar stress on the VIP surface are presented.The wall shear stresses are computed following the method presented by Mark and van Wachem(2008)and the scalar stress distribution is computed using Eqs.(19)and(20)(Bludszuweit,1995):

s s?

??????????????????????????????????????????????????????????????????????

es iiàs jjTes iiàs jjTt

s ij s ij

e19T

s ij?m@u i

@x j

@u j

@x i

tr u0

e20T

A comparison of WSS for all cases is presented in Fig.12-a to d. For the no VIP case and the static VIP case,the wall shear stresses are low at the junction and the highest(0.7Pa)at the entry of the PAs.When the pump is rotating,the maximum value of the wall shear stress(1.8Pa at1000RPM and 3.5Pa at2000RPM) increases and is located close to the pump,where the distance between the?ns and the blood vessels is the smallest.

Fig.12-e to g shows the scalar stress distribution on the VIP wall.It is generally accepted that a maximum value of450Pa is allowed to prevent the damage of blood cells(Paul et al.,2003). The scalar stress on the static pump are very low(less than0.5Pa) which is a good indicator that in case of failure of the VIP to rotate,the presence of the VIP will not harm the blood.When the VIP is rotating,the maximum value observed is less than 70Pa(respectively300Pa)at1000RPM(respectively2000

RPM).

Y.Delorme et al./Journal of Biomechanics46(2013)408–422

420

For both cases,these values are below the maximum allowed scalar stress.

Fig.13shows the repartition of the scalar stress in the TCPC at the plane z?0for all the studied cases.All four?gures present scalar stresses well below the limit allowed(less than15Pa,4Pa, 100Pa and400Pa for the no VIP case,static VIP case,VIP at 1000RPM and VIP at2000RPM respectively).It is observed that the presence of the static VIP decreases the scalar stresses on the blood by70%.The rotating VIP induces higher scalar stresses localized close the pump walls.For those cases,the stresses drop rapidly downstream of the pump.After1diameter away from the VIP,the stresses reach values of the same order of magnitude as the no VIP case.

3.6.Limitations

In LES,grid independence is related to suf?cient resolution of the large-scale(?ltered)?ow?eld for a given?lter width.While this was demonstrated for the TCPC simulations without the pump,this has not yet been demonstrated for the case with the pump.This could explain the lack of agreement between the LES and the PIV measurements for the case of the rotating VIP.Within the current code framework a signi?cant number of grid points are located outside of the?uid domain in the single block computational domain.Efforts are underway to develop a multi block version of this code to minimize wasted grid and facilitate ?ner resolution near the pump in future studies.

4.Conclusion

The Fontan circulation is currently the only way to palliate children born with one ventricle.Yet,serious challenges remain. We propose here a novel solution to increase the pressure step-up geometry and reduce the pressure loss.Due to the complex?ow features induced by both the TCPC and the VIP,a high order LES code has been developed.Non-assisted Fontan simulations show good agreement with experiments and reveal an interesting swirl switching phenomenon due to the instabilities of the two con-?ned impinging jets.The static VIP case shows that the VIP does not obstruct the Fontan?ow path and instead acts as a?ow diverter decreasing the turbulence and instabilities of the?ow in the PAs.The rotating VIP cases show reasonable qualitative agreement between LES and SPIV.The lack of resolution around the VIP is probably responsible for the discrepancies observed. Ongoing efforts focus on studying patient speci?c geometries for TCPC with more realistic boundary conditions using a lumped parameter model.Other efforts include studying blood damage using particle tracking and hepatic factor distribution in the PAs using scalar mixing.

Con?ict of interest statement

This statement is to declare that we,the authors of manuscript ‘‘Large Eddy Simulation of Powered Fontan Hemodynamics’’do not possess any?nancial relationships that might bias our work. We hereby declare that no con?ict of interest exists in our work. Acknowledgments

The authors will like to acknowledge?nancial support for this work by the National Institute of Health(NIH)(Grant no. HL098353).

References

Fig.13.Fluid scalar stresses inside the TCPC.(a)No pump.(b)Static pump.(c)Pump at1000RPM.(d)Pump at2000RPM.

Y.Delorme et al./Journal of Biomechanics46(2013)408–422421

Deleval,M.,1998.The fontan circulation:what have we learned?what to expect?

Pediatric Cardiology19,316–320.

Gewillig,M.,2005.The Fontan circulation.Congenital Heart Disease91,839–846. Giridharan,G.,Koenig,S.,Kennington,J.,Sobieski,M.,Chen,J.,Frankel,S.,Rodefeld, M.Performance evaluation of a pediatric viscous impeller pump for Fontan cavopulmonary assist.The Journal of Thoracic and Cardiovascular Surgery https://www.wendangku.net/doc/eb6322710.html,/10.1016/j.jtcvs.2012.01.082,in press.

Jeong,J.,Hussain, F.,1994.On the identi?cation of a vortex.Journal of Fluid Mechanics285,69–94.

Johnson,T.,Patel,V.,1999.Flow past a sphere up to a Reynolds number of300.

Journal of Fluid Mechanics378,19–70.

Kennington,J.,Frankel,S.,Chen,J.,Koenig,S.,Sobieski,M.,Giridharan,G.,Rodefeld, M.,2011.Design optimization and performance studies of an adult scale viscous impeller pump for powered Fontan in an idealized total cavopulmon-ary connection.Cardiovascular Engineering and Technology2,237–343. Kerlo,A.-E.,Kennington,J.,Xu,D.,Frankel,S.,Giridharan,G.,Rodefeld,M.,Chen,J., 2011.Experimental study of powered Fontan hemodynamics in idealized total cavopulmonary connection model.In:11th International Conference on Fluid Control,Measurements and Visualization.

Mark,A.,van Wachem,B.,2008.Derivation and validation of a novel implicit second order accurate immersed boundary method.Journal of Computational Physics227,6660–6680.

Marsden,A.,Bernstein,A.,Reddy,M.,Shadden,S.,Spilker,R.,Chan,F.,Taylor,C., Feinstein,J.,2009.Evaluation of a novel y-shaped extracardiac Fontan baf?e using computational?uid dynamics.The Journal of Thoracic and Cardiovas-cular Surgery137,394–403.

Migliavacca, F.,Dubini,G.,Bove, E.,deLeval,M.,https://www.wendangku.net/doc/eb6322710.html,putational?uid dynamics simulations in realistic3d geometries of the total cavopulmonary anastomosis:the in?uence of the inferior caval anastomosis.Journal of Biomechanical Engineering805,805–813.

OMP,2012./https://www.wendangku.net/doc/eb6322710.html,/wp S.

Panton,R.,2005.Incompressible Flow.John Wiley and Sons,Inc.

Paul,R.,Apel,J.,Klaus,S.,Schugner,F.,Schwindke,P.,Reul,H.,2003.Shear stress related blood damage in laminar couette?ow.Arti?cial Organs27,517–529. Pope,S.,2000.Turbulent Flows.Cambridge University Press.Rodefeld,M.,Boyd,J.,Myers,C.,Lalone,B.,Bezruczko,A.,Potter,A.,Brown,J.,2003.

Cavopulmonary assist:circulatory support for the univenricular Fontan circulation.Annals of Thoracic Surgery76,1911–1916.

Rodefeld,M.,Coats,B.,Fisher,T.,Giridharan,G.,Chen,J.,Brown,J.,Frankel,S.,2010.

Cavopulmonary assist for the univentricular Fontan circulation:von Karman viscous impeller pump.The Journal of Thoracic and Cardiovascular Surgery 140.

Rutten,F.,Schroder,W.,Meinke,M.,https://www.wendangku.net/doc/eb6322710.html,rge-eddy simulation of low frequency oscillations of the dean vortices in turbulent pipe bend?ows.Physics of Fluids

17.

Shetty,D.,Fisher,T.,Chunekar,A.,Frankel,S.,2010.High order incompressible large eddy simulation of fully inhomogeneous turbulent?ows.Journal of Computational Physics229,8802–8822.

Shetty,D.,Jie,S.,Chandy,A.,Frankel,S.,2011.A pressure correction scheme for rotational Navier Stokes equations and its application to rotating turbulnet ?https://www.wendangku.net/doc/eb6322710.html,munications in Computational Physics9,740–755.

Soerensen, D.,Pekkan,K.,de Zelicourt, D.,Sharma,S.,Kanter,K.,Fogel,M., Yoganathan,A.,2007.Introduction of a new optimized total cavopulmonary connection.The Annals of Thoracic Surgery83,2182–2190.

Stewart,S.,Patterson,E.,Burgreen,G.,Hariharan,P.,Giarra,M.,Reddy,V.,Day,S., Manning,K.,Deutsch,S.,Berman,M.,Myers,M.,Malinauskas,R.,2012.

Assessment of CFD performance in simulations of an idealized medical device: results of FDA’s?rst computational interlaboratory study.Cardiovascular Engineering and Technology3(2),139–160.

Throckmorton,A.,Ballman,K.,Myers,C.,Litwak,K.,Frankel,S.,Rodefeld,M.,2007.

Mechanical cavopulmonary assist for the univentricular Fontan circulation using a novel folding propeller blood pump.American Society of Arti?cial Internal Organs53,734–741.

Vreman,A.,2004.An eddy-viscosity subgrid-scale model for turbulent shear?ow: algebraic theory and applications.Physics of Fluids16,3670–3681. Zelicourt,D.,Marsden,A.,Fogel,M.,Yogonathan,A.,2010.Imaging and patient speci?c simulations for the Fontan surgery:current methodologies and clinical applications.Progress in Pediatric Cardiology30,31–44.

Zhang,J.,Jackson,T.,2009.A high order incompressible?ow solver with WENO.

Journal of Computational Physics228,2426–2442.

Y.Delorme et al./Journal of Biomechanics46(2013)408–422 422

Reproduced with permission of the copyright owner.Further reproduction prohibited without permission.