Cardiovascular modeling is a challenging fluid-structure interaction problem that involves treatment of complex geometries and boundary conditions in order to effectively capture physiological dynamics Taylor and Figueroa (2009). Computational fluid dynamics (CFD) is a widely-used approach for simulating blood flow in the circulatory system Taylor and Steinman (2010); Pennati et al. (2010); Vignon-Clementel et al. (2010); Mittal et al. (2001); Asgharzadeh and Borazjani (2019); Asgharzadeh et al. (2019); Seo et al. (2017); Wu and Shadden (2015); Amlani and Pahlevan (2020); Aghilinejad et al. (2020), which includes applications to 1D Amlani and Pahlevan (2020); Aghilinejad et al. (2020); Shi et al. (2011), 2D Tasciyan et al. (1993) or 3D Mittal et al. (2001); Asgharzadeh and Borazjani (2019); Asgharzadeh et al. (2019); Seo et al. (2017); Wu and Shadden (2015); Steinman and Taylor (2005); Steinman (2002) formulations. The lattice Boltzmann (LB) method Krueger et al. (2016); Benzi et al. (1992); Shan and Chen (1993); Fang et al. (2002), originating from classical statistical physics, is a powerful alternative to conventional continuum-based CFD methods that use Navier-Stokes equations. The LB method uses simplified kinetic equations combined with a modified molecular-dynamics approach to model both Newtonian and non-Newtonian fluid flow in any complex geometry (the fluid is modeled as particles that stream and collide over a discrete lattice mesh). Indeed, a particular advantage of LB-based hemodynamics solvers is their ability to easily model non-Newtonian effects via its right-hand-side; capturing such effects may be important for small vessels or vessels where the shear rate is low Sriram et al. (2014); Sankar and Hemalatha (2007). The accuracy and usefulness of the LB method have been demonstrated in a variety of fluid dynamics problems including turbulence Cosgrove et al. (2003) and multiphase flow Shan and Chen (1993). As highlighted in previous studies Shan and Chen (1993); Fang et al. (2002); Cosgrove et al. (2003); Boyd et al. (2007); Wei et al. (2020), LB methods have been shown to be particularly suitable for hemodynamics simulations since many flow features of clinical interest may require efficient numerical treatment of fully 3D computational domains.
In order to extend the clinical applicability of fluid-structure blood flow solvers based on LB equations applied to large vessels, this work introduces a direct 0D-3D coupling for the treatment of physiological boundary conditions that are governed by ordinary differential equations (ODEs) such as lumped parameter Windkessel models Westerhof et al. (1969); Vignon and Taylor (2004) or more complex hybrid ODE-Dirichlet systems such as time-varying elastance organ models Amlani and Pahlevan (2020) . Previous contributions on the 0D-3D coupling for finite element methods Moghadam et al. (2013); Kim et al. (2009) have been implicit and iterative, and for lattice Boltzmann Moore et al. (2018); Sadeghi et al. (2020) blood flow models usually only a Dirichlet or Neumann pressure or flow is prescribed during the entirety of a cardiac cycle (precluding the use of more sophisticated and non-stationary, i.e., switching, boundary conditions Amlani and Pahlevan (2020)). Additionally, recent work Sadeghi et al. (2020) on LB-based hemodynamics solvers have assumed only rigid walls, and have applied 0D lumped parameter models externally through an iterative procedure where the heart model is evolved and precomputed entirely independently Sadeghi et al. (2020) (such that the resultant pressure profile is applied on a 3D LB domain simply as a Dirichlet condition, i.e., not a true mathematical coupling).
This work, on the other hand, presents a first direct 0D-3D coupling for 3D pulsatile blood flow solvers based on LB equations. In particular, the 0D equations considered in this work govern a highly-complex and non-stationary dynamic left ventricle (LV-)elastance heart model Amlani and Pahlevan (2020) (that switches between an ODE and a Dirichlet boundary condition in a “non-stationary” fashion Amlani and Pahlevan (2020)) in order to generate physiologically-accurate hemodynamic conditions (instead of simply assigning a given inlet flow or pressure, as is commonly done Sadeghi et al. (2020); Vignon-Clementel et al. (2006b, a)). Such a coupled model represents the most complicated boundary configuration found in the circulatory system Amlani and Pahlevan (2020): a hybrid ODE-Dirichlet boundary condition representing the left ventricle, where the time at which the ODE-governed condition transitions to a Dirichlet condition is itself determined by the corresponding solution of the governing fluid-structure LB system. Hence the methodology introduced in this work can be trivially extended to the application of non-switching ODE-based boundary conditions such as lumped parameter models based on Windkessels Westerhof et al. (1969); Vignon and Taylor (2004).
This work presents a numerical approach for directly coupling this 0D LV-elastance boundary condition (or any simpler ODE-based boundary condition) to a 3D fluid-structure interaction LB solver. The methodology introduces (for both ODE-based as well as non-stationary boundaries) a discrete explicit-in-time extension to an LB non-equilibrium extrapolation method Zhao-Li et al. (2002); Krueger et al. (2016); Kang et al. (2010) for fluid-only problems that has been previously proposed only for given (often analytical) Dirichlet-based pressure/velocity boundary conditions. The ultimate aim is to enable accurate physiological LV-aortic coupling conditions for cardiovascular studies of, for example, the effects of left ventricle contractility on pulsatile hemodynamics in the aorta. Section 2 presents the governing equations for the fluid, the solid and the LV-elastance models. Section 3 details the direct 0D-3D LV-coupling strategy that is introduced in this work, including a discussion of the loss of mathematical regularity of such a model from a discontinuity in the velocity upon valve closure (and the proposition of a smoothing operator in order to ensure a continuous transition). Section 4 provides algorithmic details of the complete solver, including the numerical methods employed for the solid as well as the fluid-structure interactions (both of which can be provided by any number of suitable schemes). Section 5 presents a variety of performance studies attesting to the valid implementation and the accuracy of the fluid and solid solvers presented in this paper. Finally, Section 6 considers a sample physiological study of oscillatory wall shear stress in a 3D aorta.
2 Governing formulations
This section presents the governing equations employed in the numerical solver described in Section 4: those for the fluid domain of a vessel (governed by lattice Boltzmann equations, Section 2.1); those for the LV-elastance model for the fluid inlet (governed by hybrid ODE-Dirichlet equations, Section 2.2); and those for the solid vessel walls (governed by elastodynamics equations, Section 2.3). An illustration of the complete coupled fluid-structure computational domain is presented in Figure 1 for the (interior) fluid domain , the solid wall and the 0D-3D coupled domain .
2.1 3D lattice Boltzmann equations
For the fluid domain , the lattice Boltzmann (LB) equations are employed, where the synchronous motions of fluid particles on a regular lattice are enforced through a particle distribution function Krueger et al. (2016). This distribution function enforces mass and momentum conservation as well as ensuring that the fluid is Galilean invariant and isotropic Wolf-Gladrow (2004). In the present work, a single-relaxation-time (SRT) incompressible LB method is used to solve the incompressible flow He and Luo (1997). The evolution of the distribution functions on the lattice is governed by the discrete Boltzmann equation with the Bhatnagar-Gross-Krook (BGK) collision model, given by
where are distribution functions of the particles in phase space; are discrete velocities at position and time ; is a non-dimensional relaxation time; are equilibrium distribution functions; and are forcing terms. Here,
since a D3Q19 (19 discrete velocity vectors) stencil is applied (and a D2Q9 stencil is employed for the 3D-axisymmetric cases, i.e.,.). The non-dimensional relaxation time is related to fluid viscosity by the expression
where is the kinematic viscosity, is the incompressible fluid density (e.g., blood density), and is the lattice sound speed. Uniform discretizations are employed throughout this work for both time () and lattice space (), chosen such that (corresponding to ).
where is the fluid velocity; is the macroscopic pressure; are weighting factors (whose values are adopted from He and Luo (1997)); and is the force density in Eulerian coordinates. Pressure and velocity can be calculated from the distribution functions via the expressions
The numerical implementations and algorithmic pseudocode are provided in Section 4.
2.2 0D fluid boundary equations
For the fluid inlet boundary condition on the coupled 0D-3D interface , when the heart valve is open, the corresponding time-dependent boundary condition for that governs the dynamics of the left ventricle is given by the ODE Berger et al. (1994); Amlani and Pahlevan (2020)
for pressure inside the ventricle and time-varying compliance (the inverse of which is the corresponding time-varying elastance whose maximum is a measure of LV contractility Berger et al. (1994); Amlani and Pahlevan (2020)). Hence once the pressure in the ventricle is greater than that of the 3D fluid domain boundary , the valve opens and with the corresponding flow condition given naturally by the fluid solver via the area integral given by
Once the inflow reaches zero (or, numerically, the time at which ), the valve closes and the 0D boundary condition remains (a Dirichlet-type flow condition). Generally, is given by clinical parameters either through a look-up table or through a closed-form approximation (adopted throughout this paper from the compliance curve presented in Amlani and Pahlevan (2020)). An analysis of this LV boundary condition and details on the algorithmic implementation of such a switching configuration are provided in Section 3.
2.3 3D elastic wall equations
In order to account for fluid-structure interactions in a vessel, the solid boundary is assumed to be a thin wall that can be described by deformation of a compliant (elastic) wall in a Lagrangian coordinate system Hua et al. (2014), i.e.,
where is the density of the solid wall; is a constant wall thickness; is the Kronecker delta; is the position of the solid wall; are the Lagrangian coordinates along the solid wall; is the Lagrangian force exerted on the solid wall by the fluid; and , are stretching and bending stiffnesses (respectively). The matrices and represent in- and out-of-plane effects and, for a Poisson’s ratio , are respectively given by
The boundary condition of the solid wall (as a shell) at a simply-supported fixed end is given by
where denotes the displacement coordinates of the fixed boundaries.
All the above formulations for both the fluid and the solid can be non-dimensionalized via the reference quantities , ( is the inlet area of ) and effective length (which is equal to a diameter if the inlet is cylindrical). The corresponding non-dimensional parameters considered in our simulations are given by Reynolds number , Womersley number , bending coefficient , tension coefficient , and mass ratio of the solid wall to the fluid given by .
3 A direct 0D-3D coupling for ODE-based boundary equations and lattice Boltzmann solvers
This section proposes a discrete, direct methodology for the coupling of 3D lattice Boltzmann equations with dynamic (ODE-based) 0D models which, as described before, are often found in inflow and outflow conditions for cardiovascular configurations. The hybrid ODE-Dirichlet system governed by Equation (7) represents a most complex form of the myriad such formulations found in cardiovascular modeling, and its particular treatment is discussed in Section 3.1 (although the method proposed in what follows is straightforwardly applicable to other ODE-based 0D boundary equations such as Windkessel models). The general strategy for such multidimensional coupling of the solver presented in this work is based on the non-equilibrium extrapolation method Krueger et al. (2016); Zhao-Li et al. (2002), which has been introduced as an alternative to typical “bounce-back” methods Cornubert et al. (1991) used to implement (given) pressure and velocity boundary conditions for LB methods in order to preserve a consistent order-of-accuracy between the boundaries and the order-of-accuracy inherent to LB formulations Zhao-Li et al. (2002). In particular, such a method is ideally suited for curved fluid boundaries Kang et al. (2010); Guo et al. (2002) (such as those provided by fluid-structure interfaces of interest in this work), where the physical boundary need not coincide with the regular fluid lattice.
In this contribution, we present a discrete explicit-in-time LB extension of the non-equilibrium extrapolation method for the use of ODE-based boundary conditions (previous usages of non-equilibrium extrapolation have been mostly confined to given Dirichlet-based pressure or velocity boundary conditions of the fluid system Krueger et al. (2016); Zhao-Li et al. (2002); Kang et al. (2010)). In order to update the distribution functions from a timestep to a timestep of a fluid point that is streamed from the corresponding 0D-3D interface (boundary) node of the lattice (i.e., ), such a formulation can be derived by separating collision and streaming operations of Equation (1) into
respectively. The distribution function for in Equation (12) can be decomposed into an equilibrium and non-equilibrium part, i.e.,
which, upon substitution into Equation (12) and re-ordering terms, gives the expression for the collision as
The non-equilibrium component on the boundary can be approximated via a Chapman-Enskog asymptotic expansion Zhao-Li et al. (2002) at a distance from the neighboring fluid node as
whose second-order accuracy in is consistent with the order of the LB method. Inserting this expression into Equation (14) gives a second-order approximation of the distribution function on the boundary node as
From the definition of the equilibrium distribution function given in Equation (3), the above expression yields a “post-stream” state given by
where are the (unknown) effective pressure or velocity at the boundary point (given by the boundary condition, as described shortly). Substitution into the collision governed by Equation (12) yields the complete post-collision stream update of the boundary values as
and hence, from Equation (13), gives complete update of the fluid point from the neighboring boundary value as
Other non-equilibrium extrapolation implementations have mostly considered purely fluid domains, and hence do not consider a forcing term in the derivation. However, for the complete solver of this work, such forces (although small in an elastic regime) can be nonzero close to a fluid-structure interface.
A 0D model (correspondingly evolved in time, if an ODE) can then be directly and explicitly coupled with this formulation at each timestep through and . That is, if the 0D model generates a pressure from the flow solution (e.g., the LV-elastance model of this work), the coupling can be instituted in Equation (20) as
where comes from the evolved 0D boundary condition (e.g., an ODE), and is approximated by the corresponding value of the streamed fluid node at the current timestep. Similarly, if the 0D model generates a velocity at the boundary, then
Such a direct 0D-3D coupling enables the use of any sort of mathematical or numerical model for the 0D boundary. For example, an ODE-based boundary condition can be resolved and marched forward-in-time by any suitable integration technique (e.g., forward Euler or higher-order Amlani and Pahlevan (2020)). As an example 0D model that inputs a flow from the 3D fluid solver (in terms of the corresponding equilibrium function solutions) and generates a corresponding pressure (via, e.g., an ODE), an illustrative schematic diagram for such coupling as detailed above is presented in Figure 2.
3.1 On the particularities of the specific hybrid ODE-Dirichlet LV model
The non-stationary switching condition of the specific 0D LV-elastance equations presented in Section 2.2 leads to a loss in regularity of the solution near boundaries governed by the hybrid-ODE Dirichlet system. Indeed, the corresponding velocity at the time of valve closure (also known as a dicrotic notch on the corresponding pressure waveform) is generally non-zero, and hence the switch to a Dirichlet condition for diastole leads to a discontinuity in the velocity solution close to the 0D-3D boundary . That is, for a point neighboring a boundary node with a velocity solution during the systolic phase () defined as , the velocity over a complete cardiac cycle of length can be expressed as
which, again, can be discontinuous. Such a loss in the smoothness of the velocity derived from the switching solutions to the LV-elastance boundary equations may lead to spurious reflections (including in the form of artificial backflow) at the point of closure of the valve, i.e., at In order to ensure there is no spurious backflow or artificial oscillations resulting from immediately setting a zero velocity, we introduce a smooth transition function to the velocity profile upon valve closure. Such a function can be defined as a continuously differentiable smoothing-to-zero function and can be derived from an exponential or erf-based partition-of-unity as
where is the location of the start of the smoothing-to-zero (i.e., the time of valve closure) and is the interval over which to smoothly transition to zero. For example, the velocity in Equation (23) can be smoothly brought to zero over a discrete interval of timesteps of size (corresponding to an inteval ) by multiplying the velocity by the smoothing function as
An illustrative example of the smoothing function for closure time s, s, and is shown in Figure 3. One can easily verify that the limit from the left can be found as
and, from the right, as
where we have taken the notational license Hence is continuous everywhere including at , i.e., is of class . A flowchart summarizing the complete implementation and smooth switching of the 0D LV-elastance model with the smoothing employed in this work is presented in Figure 4.
4 Algorithmic details
The following details the numerical algorithms/methods employed for the LB fluid solver and the solid (elastic) solver (which utilizes finite differences for 3D-axisymmetric and finite elements for 3D). Their interactions, described in Section 4.2, can be facilitated by any appropriate fluid-structure algorithm: in this work, an immersed boundary method is used.
4.1 Lattice Boltzmann method
The core algorithm for solving the 3D LB equations consists of a cyclic sequence of sub-steps (where each cycle corresponds to one overall timestep) that is prescribed as:
2. Obtain the equilibrium distribution from Equation (3).
3. Perform collision (relaxation) and streaming (propagation) to update via Equation (1).
Further details about the implementation of the 3D LB method can be found elsewhere Wei et al. (2020); Boyd et al. (2007); He and Luo (1997); Krueger et al. (2016).
For some of the performance studies discussed in Section 5, a 3D-axisymmetric LB model is implemented, where an incompressible D2Q9 BGK model is used to derive an axisymmetric configuration. In such a formulation, for the pseudo-Cartesian coordinates that describe 3D axisymmetric flow, Equation (1) can be transformed into
where the a source term is incorporated for microscropic evaluation as Lee et al. (2006)
With the inflow conditions given by the complex hybrid ODE-Dirichlet system of the LV-elastance model (Section 2.2), the corresponding outflow conditions can be physiologically modeled by any suitable boundary condition that accounts for downstream physiological effects including the effective compliance, resistance and wave reflection of the truncated vasculature (to approximate the effect of the eliminated peripheral vessels). This is reasonably captured in this work using the extension outflow boundary tube model consisting of an elastic tube terminated in a rigid contraction Pahlevan et al. (2011). Such a model has been successfully utilized for hemodynamics studies Wei et al. (2021); Aghilinejad et al. (2021); Pahlevan et al. (2014); Kang et al. (2019). Other models can be used, including lumped parameter Windkessel ODEs Amlani and Pahlevan (2020); Westerhof et al. (1969); Vignon and Taylor (2004) (which can be easily implemented using the same direct 0D-3D coupling introduced earlier).
4.2 Fluid-structure interactions (FSI)
For all 3D simulations of this work (including the physiological example study of Section 6), the solid deformation given by Equation (9) is numerically simulated by the nonlinear finite element method (FEM) solver of Doyle (2013), where the the large-displacement and small-strain deformation problems are handled by co-rotational schemes. The numerical strategy has been successfully implemented in previous works for resolving a wide range of fluid-structure interaction (FSI) problems incorporating elastic structures Hua et al. (2014); Dai et al. (2012); Tang et al. (2015); Huang et al. (2018)
. Briefly, such a method uses three-node triangular elements to describe the deformation using six degrees of freedom (three displacement components and three angles of rotation)Batoz et al. (1980). An iterative strategy is then used for the time integration of the subsequent nonlinear systems of algebraic equations in order to ensure second-order accuracy. A further detailed description of the particular finite element method employed in this work can be found in Doyle (2013).
For the 3D-axisymmetric performance studies included in Section 5, a self-implemented staggered grid finite difference (SGFD) methodology Huang et al. (2007) is employed in the Lagrangian coordinate system (where is the arc length), where only the tension force given by
is defined on the interface (the displacement variable , for example, is defined on all the nodes). The solid deformation governed by Equation (9) is subsequently solved by such a finite difference methodology in a strong form Huang et al. (2007); Connell and Yue (2007); Huang and Sung (2010). That is, for an arbitrary variable, the central, downwind and upwind difference approximations to the first-order derivatives, are given by
such that the corresponding second-order central difference approximation can be defined as
Similarly. the bending force term in Equation (9) can be approximated as
For coupling the fluid and solid systems, any suitable FSI coupling strategy can be employed. For this particular work, the immersed boundary (IB) method Goldstein et al. (1993) is used to couple the LB method of the fluid with the 3D FEM (or 3D-axisymmetric FDM) of the solid Huang et al. (2007); Hua et al. (2014). This method has been extensively used to simulate FSI problems in cardiovascular biomechanics Wei et al. (2021); Peskin (2002); Mittal and Iaccarino (2005); Seo et al. (2013). The body force term in Equation (4) is used as an interaction force between the fluid and the boundary in order to enforce the no-slip velocity boundary condition at the FSI interface. The Lagrangian force between the fluid and structure, , is then calculated by a penalty scheme Goldstein et al. (1993) using the expression given by
where and are negative penalty parameters (adopted in this work from Hua et al. (2014); Tang et al. (2015); Ye et al. (2017); Zhang et al. (2017); Huang et al. (2018)); is the velocity of Lagrangian material point of the solid wall; and is the fluid velocity at the position . The latter can be obtained through transforming the Eulerian fluid velocity into Lagrangian coordinates via the integral
where is a Dirac delta function. The body force in Eulerian coordinates can be calculated from the corresponding Lagrangian body force via the expression
The Lagrangian interaction force can be explicitly obtained by the penalty IB strategy described above. Such a formulation of the IB numerical strategy has been successfully applied to a wide range of FSI problems Hua et al. (2014); Tang et al. (2015); Ye et al. (2017); Zhang et al. (2017), including those governed by the dynamics of fluid flow over a circular flexible plate Hua et al. (2014) and an inverted flexible plate Tang et al. (2015). The overall numerical algorithm for solving the complete FSI system is summarized in the block-diagram of Figure 5, and its corresponding pseudocode implementation is provided in Algorithm 1. At the start of each numerical simulation, the 0D and 3D domains can be initialized with , and an arbitrary , respectively.
5 Performance evaluation
This section presents a series of performance studies, based on benchmark cases and manufactured solutions, that validate the implementation and convergence of the solid, fluid and interaction components of the complete solver (including the LV-elastance model). Firstly, the fluid solver is validated by non-oscillatory and oscillatory cases of a flow around a cylinder. The solid solver is then validated through a classical case of a hanging elastic filament as well as the method of manufactured solutions (MMS).
5.1 Fluid solver: steady and oscillating cylinders
In order to validate the present LBM solver and its coupling with the IB method for solid wall boundaries, two cases are considered: a uniform flow passing a non-oscillating cylinder and uniform flow passing a transversely oscillating cylinder. In those cases the motion of the solid wall is prescribed (rigidly), and hence the problem becomes a one-way fluid-structure coupling such that the solid position of the center of the cylinder is given by
where is the cylinder diameter; and , are the constant oscillation amplitude and frequency, respectively. In the cases considered in this section, the Reynolds number is defined as for velocity (m/s), fluid viscosity and m. Defining radians as the natural shedding frequency for a stationary cylinder, the computation is performed with the parameters m, and two oscillating cases and (note that for a non-oscillating cylinder). These test cases and their corresponding parameters are adopted from the benchmark cases proposed by Guilmineau and Queutey (2002). Figures 6 (left) illustrates the computational domain for both steady and moving cases.
Figure 6 (right) presents the corresponding time evolution of the drag and lift coefficients , of the steady case as simulated by our solver up to a final time of s. The numerical discretization corresponds to (where m). Using the same discretization and computational domain, Figure 7 additionally presents the time evolution of the drag and lift coefficients of the moving (oscillating) cylinder case for (left) and (right). All of these cases (steady and moving) are in excellent agreement with those results presented in Guilmineau and Queutey (2002); Kim and Choi (2006).
5.2 Solid solver
In order to validate the elastic (solid) component, two case are considered: an elastic filament under gravity and a manufactured solution (i.e., a given right-hand-side).
5.2.1 A filament under gravity
The first case considers a hanging filament of length m without an ambient fluid (in order to independently test the solid solver alone), as proposed by Huang et al. (2007). Adopting the same arbitrary units, the filament is initially held stationary at an angle radians from the vertical, where the physical parameters correspond to a Froude number ( is the gravitation constant) and a bending rigidity Pa m. Snapshots of the simulated positions of the filament over a time period of 0.8 s are illustrated in Figure 8, where a spatial discretization of elements is utilized together with a timestep of s. Again, the solutions produced by the solver of this work are in excellent agreement with those results presented by Huang et al. (2007).
5.2.2 An arbitrarily moving filament
In order to further validate the implementation of the solid solver (in particular, it’s numerical accuracy), the method of manufactured solutions (MMS) is additionally employed. Such a verification procedure has been extensively used for validating other hemodynamics solvers Amlani and Pahlevan (2020); Raghu and Taylor (2011); Raghu et al. (2011). In MMS, one proposes a closed-form smooth solution (i.e., arbitrary movement) and subsequently derives (algebraically) the corresponding right-hand forcing terms and boundary conditions in order to render the postulated function to be an exact solution of the solid equations. Here, we postulate a given displacement function as
where m is the maximum amplitude of displacement in the vertical component. The solid structure considered is again an elastic filament of length m with elastic material parameters corresponding to a Young’s modulus of MPa and a thickness of mm. Employing the same numerical discretization as in the gravity case, Figure 9 (left) presents snapshots in time of the simulated filament positions for both numerical (solid lines) and analytical values (dashed lines).
Using successive discretization sizes that are integer multiples of the coarsest one used ( elements), the simulation is advanced for timesteps in all cases at a fine time-step size of s (in order to ensure that errors are dominated by the spatial discretization). The maximum absolute errors between simulated displacements and the exact manufactured solution of Equation (39), over all space and for all timesteps, are presented in Figure 9 (right). The overlaid slopes in the plots illustrate the expected second-order of accuracy for the elastic (solid) solver employed in this paper (hence verifying its implementation).
5.3 The complete FSI solver coupled to the 0D LV-elastance heart model
In order to verify the complete FSI solver coupled to a 0D hybrid ODE-Dirichlet heart model, one can consider the axisymmetric straight aorta configuration (of length ) presented in Figure 10. Adopting parameters of the LV-elastance hybrid ODE-dirichlet model from Amlani and Pahlevan (2020) corresponding to an end-systolic LV elastance mmHg/ml and a L/min, Figure 11 (left) presents the expected physiologically-accurate pressure profiles at the 0D-3D interface as simulated by the complete solver for discretizations corresponding to and , where physical parameters of , a non-dimensionalized (corresponding to 24 mm), and are employed. The timestep is fixed again and is taken small enough so that errors are dominated by the spatial discretization. Figure 11 (right) presents the corresponding errors (relative to the finest solution), where a convergence between first and second order can be observed (and is expected from the second-order nature of the lattice Boltzmann solver and the first-order discretizations of the LV-elastance ODEs). Figure 12 (left) and Figure 13 (left) additionally present the simulated physiological flow profiles at the inlet and the pressure profiles at the midpoint of the vessel, respectively. The corresponding errors (relative to the finest discretization of ) are presented in Figure 12 (right) and Figure 13 (right), respectively. As before, one can appreciate the convergence and accuracy as expected from the second-order fluid discretization and the first order LV-elastance ODE time integration.
6 An example physiological case: wall shear stress in the aorta
A predominant effect of advanced congestive heart failure (CHF) is reduced blood flow in the aorta that results from a reduction in cardiac output (CO) and a low ejection fraction (in almost half of the patients). Many factors can influence the heart’s pumping ability, including those related to the direct coupling between the LV and the arterial system Berger et al. (1994, 1996); Campbell et al. (1989). The hybrid ODE-Dirichlet boundary condition considered in this paper has been chosen for its ability to model the non-stationary and nonlinear effects of such complex coupling (which is expressed as an alternating boundary condition between systole—an ODE—and diastole—a Dirichlet condition—as described in Section 2.2). The 0D elastance- (compliance-)based LV model enables the generation of physiological pressure and flow waveforms that can account for different contractilities and cardiac outputs, and its corresponding coupling to 3D lattice Boltzmann equations can enable investigation of the heart’s influence on corresponding 3D fluid-structure effects such as those related to near-wall shear stress (WSS). Indeed, mechanical experiments Gharib and Beizaie (2003) and both in-vivo/in-vitro studies Moore Jr et al. (1994); Pedersen et al. (1993); Oyre et al. (1997); Pedersen et al. (1999) have shown there can exist negative WSS corresponding to a retrograde flow during a substantial interval within a cardiac cycle, and such effects have been strongly correlated with the state of CHF Gharib and Beizaie (2003).
As a demonstration of the applicability of the proposed solver towards exploring these parameters for studying pathophysiological conditions in the cardiovascular system (e.g., CHF), a computational model of a simplified 3D aorta (which excludes the branches) is considered and illustrated in Figure 14 (left), where the elastic wall is discretized by finite elements and the fluid by lattice Boltzmann as described in Section 4.2. For an effective aortic diameter (taken to be unity in the non-dimensionalized configuration), such a domain corresponds to Cartesian coordinates given by for a non-dimensionalized (corresponding to mm). For the compliant wall, a linear elastic material with Young’s modulus of MPa and a wall thickness corresponding to mm is considered. The complete fluid-structure (immersed boundary) solver is coupled to the 0D LV-elastance heart model (Figure 4), where a no-slip condition is imposed at the fluid-structure interface. Using an extension tube boundary model Pahlevan et al. (2011) at the outlet (with a constant outlet pressure of mmHg), Figure 14 (right) presents the corresponding normalized velocity magnitudes produced by a simulation that emplys discretizations s and is advanced up to a time s (where s corresponds to the period of a cardiac cycle). The LV parameters and compliance function are adopted from Amlani and Pahlevan (2020) and correspond to a healthy case with normal contractility. Figure 15 (left) additionally illustrates the expected physiological characteristics of the LV and aortic pressures, particularly the equality during systole (in the absence of a diseased valve condition) between ventricular pressure and the aortic pressure at the coupled boundary. Figure 15 (right) provides the corresponding physiological flow profiles as simulated at the inlet, midpoint and outlet of the simplified 3D aorta.
In investigating WSS (as a relevant hemodynamic biomarker in CHF Gharib and Beizaie (2003)), two contractility cases (representing a low flow rate and a high flow rate) can be considered employing the same Womersely number . For normal contractility, the end-systolic LV elastance Amlani and Pahlevan (2020) is set to mmHg/ml, which corresponds to a cardiac output of l/min and is in accordance with values employed in experimental studies Gharib and Beizaie (2003). For the high contractility scenario, = mmHg/ml. Again, for both cases, the Womersely number (corresponding to a heart rate of BPM Gharib and Beizaie (2003)) is fixed. Figure 16 presents the corresponding wall shear stress (WSS) for both normal and high contractility cases, as calculated through the fluid points next to the solid wall Gharib and Beizaie (2003) via the expression
where is the fluid viscosity (corresponding to 3.5 centipoise), is the simulated axial velocity, and is the dimension normal to the wall. As expected Gharib and Beizaie (2003), a negative WSS (corresponding to a retrograde flow) is evident for the normal contractility case parameters, and such effects dissappear in the high contractility case since the corresponding flow rate is very high. These results are in agreement with the experimental results presented in Gharib and Beizaie (2003).
This work presents a direct 0D-3D coupling for dynamic (ODE-based) boundary conditions applied to lattice Boltzmann solvers for hemodynamic flow. Benchmark performance studies and a physiological case of wall shear stress in a simplified 3D aorta are treated in order to validate the proposed methodology and its implementation. In particular, this work treats a most complicated configuration of such coupling conditions: a hybrid non-stationary ODE-Dirichlet boundary condition. Such a methodology produces a physiologically-accurate hemodynamics solver (with a heart model) for studying wave propagation and pulsatile blood flow in arterial vessels. The methodology introduced in this paper can be easily extended to non-switching ODE conditions such as 0D lumped parameter models (e.g., Windkessels for truncating vasculature at vessel outlets), as well as to other methods for treating fluid-structure interactions (facilitated here by an immersed boundary method). Such a direct 0D-3D coupling with the proposed regularization of Section 3 can also be applied to any other fluid problem that is governed by lattice Boltzmann equations and that requires direct time-dependent ODEs as boundary conditions.
- Effects of vessel wall mechanics on non-invasive evaluation of cardiovascular intrinsic frequencies. Journal of Biomechanics, pp. 110852. Cited by: §4.1.
- Dynamic effects of aortic arch stiffening on pulsatile energy transmission to cerebral vasculature as a determinant of brain-heart coupling. Scientific Reports 10 (1), pp. 1–12. Cited by: §1.
- A stable high-order FC-based methodology for hemodynamic wave propagation. Journal of Computational Physics 405, pp. 109130. Cited by: §1, §1, §1, §2.2, §2.2, §3, §4.1, §5.2.2, §5.3, §6, §6.
- A non-dimensional parameter for classification of the flow in intracranial aneurysms. ii. patient-specific geometries. Physics of Fluids 31 (3), pp. 031905. Cited by: §1.
- A non-dimensional parameter for classification of the flow in intracranial aneurysms. i. simplified geometries. Physics of Fluids 31 (3), pp. 031904. Cited by: §1.
- A study of three-node triangular plate bending elements. International Journal for Numerical Methods in Engineering 15 (12), pp. 1771–1812. Cited by: §4.2.
- The lattice boltzmann equation: theory and applications. Physics Reports 222 (3), pp. 145–197. Cited by: §1.
- Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study. American Journal of Physiology-Heart and Circulatory Physiology 266 (4), pp. H1626 – H1642. Cited by: §2.2, §6.
- Wave propagation in coupled left ventricle-arterial system. Hypertension 27 (5), pp. 1079 – 1089. Cited by: §6.
- Analysis of the casson and carreau-yasuda non-newtonian blood models in steady and oscillatory flows using the lattice Boltzmann method. Physics of Fluids 19 (9), pp. 093103. Cited by: §1, §4.1.
- Pulse reflection sites and effective length of the arterial system. American Journal of Physiology-Heart and Circulatory Physiology 256 (6), pp. H1684 – H1689. Cited by: §6.
- Flapping dynamics of a flag in a uniform stream. Journal of fluid mechanics 581, pp. 33–67. Cited by: §4.2.
- A Knudsen layer theory for lattice gases. Physica D: Nonlinear Phenomena 47 (1-2), pp. 241–259. Cited by: §3.
- Application of the lattice Boltzmann method to transition in oscillatory channel flow. Journal of Physics A: Mathematical and General 36 (10), pp. 2609. Cited by: §1.
- Dynamic pitching of an elastic rectangular wing in hovering motion. Journal of Fluid Mechanics 693, pp. 473–499. Cited by: §4.2.
- Nonlinear analysis of thin-walled structures: statics, dynamics, and stability. Springer Science & Business Media. Cited by: §4.2.
- Lattice Boltzmann method for simulating the viscous flow in large distensible blood vessels. Physical Review E 65 (5), pp. 051925. Cited by: §1.
- Correlation between negative near-wall shear stress in human aorta and various stages of congestive heart failure. Annals of Biomedical Engineering 31 (6), pp. 678–685. Cited by: §6, §6.
- Modeling a no-slip flow boundary with an external force field. Journal of Computational Physics 105 (2), pp. 354–366. Cited by: §4.2.
- A numerical simulation of vortex shedding from an oscillating circular cylinder. Journal of Fluids and Structures 16 (6), pp. 773–794. Cited by: Figure 6, Figure 7, §5.1, §5.1.
- An extrapolation method for boundary conditions in lattice Boltzmann method. Physics of Fluids 14 (6), pp. 2007–2010. Cited by: §3.
- Lattice Boltzmann model for the incompressible Navier–Stokes equation. Journal of Statistical Physics 88 (3-4), pp. 927–944. Cited by: §2.1, §2.1, §4.1.
- Discrete Boltzmann equation model for nonideal gases. Physical Review E 57 (1), pp. R13. Cited by: §2.1.
- Dynamics of fluid flow over a circular flexible plate. Journal of Fluid Mechanics 759, pp. 56–72. Cited by: §2.3, §4.2, §4.2.
- Coupling performance of tandem flexible inverted flags in a uniform flow. Journal of Fluid Mechanics 837, pp. 461–476. Cited by: §4.2, §4.2.
- Simulation of flexible filaments in a uniform flow by the immersed boundary method. Journal of Computational Physics 226 (2), pp. 2206–2228. Cited by: §4.2, §4.2, Figure 8, §5.2.1.
- Three-dimensional simulation of a flapping flag in a uniform flow. Journal of Fluid Mechanics 653, pp. 301–336. Cited by: §4.2.
- On the accuracy of displacement-based wave intensity analysis: effect of vessel wall viscoelasticity and nonlinearity. PloS one 14 (11). Cited by: §4.1.
- Non-equilibrium extrapolation method in the lattice Boltzmann simulations of flows with curved boundaries (non-equilibrium extrapolation of LBM). Applied Thermal Engineering 30 (13), pp. 1790–1796. Cited by: §1, §3, §3.
- Immersed boundary method for flow around an arbitrarily moving body. Journal of Computational Physics 212 (2), pp. 662–680. Cited by: Figure 6, Figure 7, §5.1.
- On coupling a lumped parameter heart model and a three-dimensional finite element aorta model. Annals of Biomedical Engineering 37 (11), pp. 2153–2169. Cited by: §1.
- The Lattice Boltzmann Method: Principles and Practice. Graduate Texts in Physics, Springer. Cited by: §1, §1, §2.1, §3, §3, §4.1.
- An axisymmetric incompressible lattice Boltzmann model for pipe flow. International Journal of Modern Physics C 17 (05), pp. 645–661. Cited by: §4.1.
- Application of large-eddy simulation to the study of pulsatile flow in a modeled arterial stenosis. Journal of Biomechanical Engineering 123 (4), pp. 325–332. Cited by: §1.
- Immersed boundary methods. Annual Review of Fluid Mechanics 37, pp. 239–261. Cited by: §4.2.
- A modular numerical method for implicit 0d/3d coupling in cardiovascular finite element simulations. Journal of Computational Physics 244, pp. 63–79. Cited by: §1.
- Fluid wall shear stress measurements in a model of the human abdominal aorta: oscillatory behavior and relationship to atherosclerosis. Atherosclerosis 110 (2), pp. 225–240. Cited by: §6.
- Towards realtime 3d coronary hemodynamics simulations during cardiac catheterisation. In 2018 Computing in Cardiology Conference (CinC), Vol. 45, pp. 1–4. Cited by: §1.
- In vivo wall shear stress measured by magnetic resonance velocity mapping in the normal human abdominal aorta. European Journal of Vascular and Endovascular Surgery 13 (3), pp. 263–271. Cited by: §6.
- A physiologically relevant, simple outflow boundary model for truncated vasculature. Annals of Biomedical Engineering 39 (5), pp. 1470–1481. Cited by: §4.1, §6.
- Intrinsic frequency for a systems approach to haemodynamic waveform analysis with clinical applications. Journal of The Royal Society Interface 11 (98), pp. 20140617. Cited by: §4.1.
- Quantitative abdominal aortic flow measurements at controlled levels of ergometer exercise. Magnetic Resonance Imaging 17 (4), pp. 489–494. Cited by: §6.
- Two-dimensional velocity measurements in a pulsatile flow model of the normal abdominal aorta simulating different hemodynamic conditions. Journal of Biomechanics 26 (10), pp. 1237–1247. Cited by: §6.
- Modeling of systemic-to-pulmonary shunts in newborns with a univentricular circulation: state of the art and future directions. Progress in Pediatric Cardiology 30 (1-2), pp. 23–29. Cited by: §1.
- The immersed boundary method. Acta Numerica 11, pp. 479–517. Cited by: §4.2.
- Verification of a one-dimensional finite element method for modeling blood flow in the cardiovascular system incorporating a viscoelastic wall model. Finite Elements in Analysis and Design 47 (6), pp. 586–592. Cited by: §5.2.2.
- Comparative study of viscoelastic arterial wall models in nonlinear one-dimensional finite element simulations of blood flow. Journal of Biomechanical Engineering 133 (8). Cited by: §5.2.2.
- Towards non-invasive computational-mechanics and imaging-based diagnostic framework for personalized cardiology for coarctation. Scientific Reports 10 (1), pp. 1–19. Cited by: §1, §1.
- A non-Newtonian fluid flow model for blood flow through a catheterized artery—steady flow. Applied Mathematical Modelling 31 (9), pp. 1847–1864. Cited by: §1.
- A method for the computational modeling of the physics of heart murmurs. Journal of Computational Physics 336, pp. 546–568. Cited by: §1.
- Multiphysics computational models for cardiac flow and virtual cardiography. International journal for numerical methods in Biomedical Engineering 29 (8), pp. 850–869. Cited by: §4.2.
- Lattice Boltzmann model for simulating flows with multiple phases and components. Physical Review E 47 (3), pp. 1815. Cited by: §1.
- Review of zero-D and 1-D models of blood flow in the cardiovascular system. Biomedical Engineering online 10 (1), pp. 1–38. Cited by: §1.
- Non-newtonian flow of blood in arterioles: consequences for wall shear stress measurements. Microcirculation 21 (7), pp. 628–639. Cited by: §1.
- Flow imaging and computing: large artery hemodynamics. Annals of Biomedical Engineering 33 (12), pp. 1704–1709. Cited by: §1.
- Image-based computational fluid dynamics modeling in realistic arterial geometries. Annals of Biomedical Engineering 30 (4), pp. 483–497. Cited by: §1.
- Dynamics of an inverted flexible plate in a uniform flow. Physics of Fluids 27 (7), pp. 073601. Cited by: §4.2, §4.2.
- Two-dimensional pulsatile hemodynamic analysis in the magnetic resonance angiography interpretation of a stenosed carotid arterial bifurcation. Medical Physics 20 (4), pp. 1059–1070. Cited by: §1.
- Patient-specific modeling of cardiovascular mechanics. Annual Review of Biomedical Engineering 11, pp. 109–134. Cited by: §1.
- Image-based modeling of blood flow and vessel wall dynamics: applications, methods and future directions. Annals of Biomedical Engineering 38 (3), pp. 1188–1203. Cited by: §1.
- Outflow boundary conditions for one-dimensional finite element modeling of blood flow and pressure waves in arteries. Wave Motion 39 (4), pp. 361–374. Cited by: §1, §1, §4.1.
- Outflow boundary conditions for three-dimensional simulations of non-periodic blood flow and pressure fields in deformable arteries. Journal of Biomechanics (39), pp. S431. Cited by: §1.
- Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Computer methods in applied mechanics and engineering 195 (29-32), pp. 3776–3796. Cited by: §1.
- A primer on computational simulation in congenital heart disease for the clinician. Progress in Pediatric Cardiology 30 (1-2), pp. 3–13. Cited by: §1.
- On the significance of blood flow shear-rate-dependency in modeling of Fontan hemodynamics. European Journal of Mechanics-B/Fluids. Cited by: §1, §4.1.
- Hemodynamically efficient artificial right atrium design for univentricular heart patients. Physical Review Fluids (in press). Cited by: §4.1, §4.2.
- Analog studies of the human systemic arterial tree. Journal of Biomechanics 2 (2), pp. 121–143. Cited by: §1, §1, §4.1.
- Lattice-gas cellular automata and lattice Boltzmann models: an introduction. Springer. Cited by: §2.1.
- Coupled simulation of hemodynamics and vascular growth and remodeling in a subject-specific geometry. Annals of Biomedical Engineering 43 (7), pp. 1543–1554. Cited by: §1.
- Two tandem flexible loops in a viscous flow. Physics of Fluids 29 (2), pp. 021902. Cited by: §4.2.
- Free locomotion of a flexible plate near the ground. Physics of Fluids 29 (4), pp. 041903. Cited by: §4.2.
- Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. Chinese Physics 11 (4), pp. 366. Cited by: §1, §3, §3.