Direct 0D‐3D coupling of a lattice Boltzmann methodology for fluid–structure aortic flow simulations

This work introduces a numerical approach and implementation for the direct coupling of arbitrary complex ordinary differential equation‐ (ODE‐)governed zero‐dimensional (0D) boundary conditions to three‐dimensional (3D) lattice Boltzmann‐based fluid–structure systems for hemodynamics studies. In particular, a most complex configuration is treated by considering a dynamic left ventricle‐ (LV‐)elastance heart model which is governed by (and applied as) a nonlinear, non‐stationary hybrid ODE‐Dirichlet system. Other ODE‐based boundary conditions, such as lumped parameter Windkessel models for truncated vasculature, are also considered. Performance studies of the complete 0D‐3D solver, including its treatment of the lattice Boltzmann fluid equations and elastodynamics equations as well as their interactions, is conducted through a variety of benchmark and convergence studies that demonstrate the ability of the coupled 0D‐3D methodology in generating physiological pressure and flow waveforms—ultimately enabling the exploration of various physical and physiological parameters for hemodynamics studies of the coupled LV‐arterial system. The methods proposed in this paper can be easily applied to other ODE‐based boundary conditions as well as to other fluid problems that are modeled by 3D lattice Boltzmann equations and that require direct coupling of dynamic 0D boundary conditions.

a most complex configuration is treated by considering a dynamic left ventricle-(LV-)elastance heart model which is governed by (and applied as) a nonlinear, non-stationary hybrid ODE-Dirichlet system. Other ODE-based boundary conditions, such as lumped parameter Windkessel models for truncated vasculature, are also considered. Performance studies of the complete 0D-3D solver, including its treatment of the lattice Boltzmann fluid equations and elastodynamics equations as well as their interactions, is conducted through a variety of benchmark and convergence studies that demonstrate the ability of the coupled 0D-3D methodology in generating physiological pressure and flow waveforms-ultimately enabling the exploration of various physical and physiological parameters for hemodynamics studies of the coupled LV-arterial system. The methods proposed in this paper can be easily applied to other ODE-based boundary conditions as well as to other fluid problems that are modeled by 3D lattice Boltzmann equations and that require direct coupling of dynamic 0D boundary conditions.

K E Y W O R D S
0D-3D coupling, 3D blood flow, fluid-structure interactions, lattice Boltzmann equations, LV-arterial hemodynamics, mathematical physiology 1 | INTRODUCTION 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. 1 Computational fluid dynamics (CFD) is a widely-used approach for simulating blood flow in the circulatory system, [2][3][4][5][6][7][8][9][10][11] which includes applications to 1D, 10-12 2D 13 or 3D [5][6][7][8][9]14,15 formulations. The lattice Boltzmann (LB) method, [16][17][18][19] 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. 20,21 The accuracy and usefulness of the LB method have been demonstrated in a variety of fluid dynamics problems including turbulence 22 and multiphase flow. 18 As highlighted in previous studies, 18,19,[22][23][24] 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 (and the local nature of LBM collision calculations enables highly-parallelizable implementations).
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 25,26 or more complex hybrid ODE-Dirichlet systems such as time-varying elastance organ models 10 (both known as 0D models due to the absence of spatial dependence). Previous contributions on the 0D-3D coupling for finite element methods 27,28 have been implicit and iterative, and for lattice Boltzmann 29,30 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 nonstationary, i.e., switching, boundary conditions 10 ).
Additionally, recent work 30 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 30 (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 fully fluid-structure 3D pulsatile blood flow solvers based on LB and elastodynamics equations. In particular, the 0D equations considered in this work govern a highly-complex and non-stationary dynamic left ventricle (LV-)elastance heart model 10 (that switches between an ODE and a Dirichlet boundary condition in a "non-stationary" fashion 10 ) in order to generate physiologically-accurate hemodynamic conditions (instead of simply assigning a given inlet flow or pressure, as is commonly done [30][31][32]. Such a coupled model represents the most complicated boundary configuration found in the circulatory system 10 : 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 10,25,26 (also treated in this contribution).
This work presents a numerical approach for directly coupling these 0D LV-elastance and Windkessel boundary conditions (or any simpler ODE-based boundary condition) to a 3D LB-based fluid-structure interaction solver for hemodynamics (where the solid is governed by elastic equations). The methodology introduces, for both ODE-based as well as non-stationary boundaries, a discrete explicit-in-time modification to an LB non-equilibrium extrapolation method 16,33,34 that has been previously proposed for fluid-only problems (i.e., no solid interaction) and only for given (often analytical) Dirichlet-based pressure or velocity boundary conditions (a particular novelty here is in the non-stationary switching between pressure and velocity). 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.1 presents the governing equations for the fluid, the solid and the LV-elastance models. Section 2.2 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 2.3 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 3.1 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 3.2 considers a sample physiological study of oscillatory wall shear stress in a 3D aorta with carotid and renal branches.

| Governing formulations
This section presents the governing equations employed in the numerical solver described in Section 2.3: those for the fluid domain of a vessel (governed by lattice Boltzmann equations, Section 2.1.1); those for the LV-elastance model for the fluid inlet (governed by hybrid ODE-Dirichlet equations, Section 2.1.2); and those for the solid vessel walls (governed by elastodynamics equations, Section 2.1.3). An illustration of the complete coupled fluid-structure computa- Figure 1 for the (interior) fluid domain Ω, the solid wall ∂Ω 1 and the 0D-3D coupled domain ∂Ω 2 .

| 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. 16 This distribution function enforces mass and momentum conservation as well as ensuring that the fluid is Galilean invariant and isotropic. 35 In the present work, a single-relaxation-time (SRT) incompressible LB method is used to solve the incompressible flow. 36 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 f i x, t ð Þ are distribution functions of the particles in phase space; e i are discrete velocities at position x and time t; τ is a non-dimensional relaxation time; f eq i x, t ð Þ are equilibrium distribution functions; and F i are forcing terms. Here, N 0 ¼ 19 since a D3Q19 (19 discrete velocity vectors) stencil is applied (and a D2Q9 stencil is employed for the 3Daxisymmetric cases, i.e., N 0 ¼ 9.). 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 c s ¼ Δx= Δt ffiffi ffi 3 p À Á is the lattice sound speed. Uniform discretizations are employed throughout this work for both time (Δt) and lattice space (Δx), chosen such that Δx=Δt ¼ 1 (corresponding to c s ¼ 1= ffiffi ffi 3 p ). The equilibrium distribution functions f eq i x, t ð Þ for an incompressible Lattice Boltzmann model 36 and the forcing terms F i x,t ð Þ 37 are resepctively defined as where v x, t ð Þ is the fluid velocity; P x, t ð Þ is the macroscopic pressure; ω i are weighting factors (whose values are adopted from previous study 36 ); and b x, t ð Þ is the force density in Eulerian coordinates. Pressure P x, t ð Þ and velocity v x, t ð Þ can be calculated from the distribution functions f i x, t ð Þ via the expressions The numerical implementations and algorithmic pseudocode are provided in Section 2.3.

| 0D fluid boundary equations: LV-elastance and lumped parameter models
For the fluid inlet boundary condition on the coupled 0D-3D interface ∂Ω inlet & ∂Ω 2 , when the heart valve is open, the corresponding time-dependent boundary condition for x Ω inlet that governs the dynamics of the left ventricle is given by the ODE 10,38 for pressure inside the ventricle P v t ð Þ and time-varying compliance C v t ð Þ (the inverse of which is the corresponding time-varying elastance whose maximum is a measure of LV contractility 10,38 ). Hence once the pressure P v t ð Þ in the ventricle is greater than that of the inlet fluid domain boundary P x ∂Ω inlet ,t ð Þ , the valve opens and P x, t ð Þ¼P v t ð Þ with the corresponding flow condition Q t ð Þ given naturally by the fluid solver via the area integral given by Once the inflow reaches zero (or, numerically, the time at which Q t ð Þ < 0), the valve closes and the 0D boundary condition remains Q t ð Þ ¼ 0 (a Dirichlet-type flow condition). Generally, C v t ð Þ 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 10 ). An analysis of this nonlinear, non-stationary LV boundary condition and details on the algorithmic implementation of such a switching configuration are provided in Section 2.2.
At the outlet boundary ∂Ω outlet & ∂Ω 2 of the 3D physiological aorta considered in Section 3.2, a conventional 0D lumped parameter model, a so-called circuit-like Windkessel model, is employed to represent the effects of truncated vasculature 10 (i.e., eliminated periphery). The outlet of the fluid domain is coupled to such a model through a matching characteristic impedance Z w that is related to the fluid inductance and overall outgoing aortic compliance. Together with an effective chamber compliance C w and a total peripheral resistance R w , the pressure P w in the terminal compliance chamber is related to the aortic pressure P x, t ð Þ at the outlet boundary x ∂Ω outlet through an ODE given by The form of Equation (9) can be derived from an electrical analogy; it effectively governs an impedance-matching (at the outgoing vessel) circuit in order to affect steady flow (i.e., pressure corresponds to voltage and flow to current). The corresponding 0D outflow Q x, t ð Þ at the outlet x ∂Ω outlet is given by Further details on the parameters, usage and implementation of both this outgoing Windkessel model, as well as the LV-elastance heart model above, can be found elsewhere. 10

| 3D elastic wall equations
In order to account for fluid-structure interactions in a vessel, the solid boundary ∂Ω 1 is assumed to be a thin wall that can be described by deformation of a compliant (elastic) wall in a Lagrangian coordinate system, 39 that is, where ρ s is the density of the solid wall; h is a constant wall thickness; δ ij is the Kronecker delta; X s, t ð Þ ℝ 3 is the position of the solid wall; s ¼ s 1 , s 2 ð Þ ℝ 2 are the Lagrangian coordinates along the solid wall; B s,t ð Þ is the Lagrangian force exerted on the solid wall by the fluid; and Eh, EI are stretching and bending stiffnesses (respectively). The matrices φ ð Þ ij and γ ð Þ ij represent in-and out-of-plane effects and, for a Poisson's ratio b ν, are respectively given by The boundary condition of the solid wall (as a shell) at a simply-supported fixed end is given by where X 0 denotes the displacement coordinates of the fixed boundaries.

| Non-dimensionalization
All the above formulations for both the fluid and the solid can be non-dimensionalized via the reference quantities ρ, The corresponding non-dimensional parameters considered in our simulations are given by Reynolds number ∞ L, and mass ratio of the solid wall to the fluid given by M ¼ ρ s h=ρL.

| 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 2.2.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, 16,33 which has been introduced as an alternative to typical "bounce-back" methods 40 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. 33 In particular, such a method is ideally suited for curved fluid boundaries 34,41 (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 modification 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 16,33,34 ). In order to update the distribution functions from a timestep t n to a timestep t nþ1 ¼ t n þ Δt of a fluid point x f Ω that is streamed from the corresponding 0D-3D interface (boundary) node x ∂Ω 2 of the lattice (i.e., x f ¼ x þ e i Δt), such a formulation can be derived by separating collision and streaming operations of Equation (1) into and respectively. The distribution function for x ∂Ω 2 in Equation (14) can be decomposed into an equilibrium and nonequilibrium part, i.e., which, upon substitution into Equation (14) 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 33 at a distance 0 < ϵ ( 1 from the neighboring fluid node x f as whose second-order accuracy in ϵ is consistent with the order of the LB method. In such an asymptotic analysis, ϵ is a "smallness parameter" that label's each term's order in the dimensionless Knudsen number Kn, 42 that is, ϵ m ¼ Kn m : Inserting the above expression into Equation (16) gives a second-order approximation of the distribution function on the boundary node x as From the definition of the equilibrium distribution function given in Equation (3), the above expression yields a "poststream" state given by where P Ã x ð Þ,v Ã x ð Þ 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 (14) yields the complete post-collision stream update of the boundary values as and hence, from Equation (15), gives a complete update of the fluid point x f from the neighboring boundary value x as Other non-equilibrium extrapolation implementations have mostly considered purely fluid domains, and hence do not consider a forcing F i 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 P Ã and v Ã . That is, if the 0D model generates a pressure P 0D t ð Þ from the flow solution (e.g., the LV-elastance model of this work), the coupling can be instituted in Equation (22) as where P Ã comes from the evolved 0D boundary condition (e.g., an ODE), and v Ã is approximated by the corresponding value of the streamed fluid node x f at the current timestep. Similarly, if the 0D model generates a velocity v 0D t ð Þ 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 10 ). 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. For the LVelastance model of Section 2.1.2, the flow (i.e., velocity integrated over the inlet area) is input into the 0D model, which returns a singular (internal) pressure before valve closure or a fully-zero flow velocity after valve closure (hence the 0D output is prescribed everywhere at the inlet of the fluid domain with no additional treatment).

| 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.1.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 t ¼ T d 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 ∂Ω 2 . That is, for a point x f Ω neighboring a boundary node x ∂Ω 2 with a velocity solution during the systolic phase (t 0, , the velocity over a complete cardiac cycle of length T 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, that is, at t ¼ T d . 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 C ∞ smoothing-to-zero function S t ð Þ and can be derived from an exponential or erf-based partition-of-unity as where t 0 ¼ T d is the location of the start of the smoothing-to-zero (i.e., the time of valve closure) and L S is the interval over which to smoothly transition to zero (indeed, there is a non-instantaneous valve closure in physical reality). For example, the velocity v 0 t ð Þ in Equation (25) can be smoothly brought to zero over a discrete interval of N S timesteps of size Δt (corresponding to an inteval L S ¼ N S Δt) by multiplying the velocity by the smoothing function as An illustrative example of the smoothing function S for closure time t 0 ¼ T d ¼ 0:5 s, Δt ¼ 0:005 s, and N S ¼ 15 is shown in Figure 3. One can easily verify that the limit from the left can be found as and, from the right, as F I G U R E 2 An illustrative example diagram of the direct (explicit) coupling between a 0D model (e.g., the ODEs corresponding to the LV-heart model and the Windkessel model) and the 3D lattice Boltzmann (LB) model. The flow in ∂Ω 2 is computed in terms of the LB distribution functions at a timestep n which is fed into ODEs governing the 0D model. The corresponding pressure produced by the 0D model is then re-translated into distribution functions on the boundary via the non-equilibrium extrapolation described by Equation (20) lim where we have taken the notational license S t A flowchart summarizing the complete implementation and smooth switching of the 0D LV-elastance model with the C ∞ smoothing employed in this work is presented in Figure 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 2.3.2, can be facilitated by any appropriate fluid-structure algorithm: in this work, an immersed boundary method is used.

| 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: 1. Compute the macroscopic moments P x,t ð Þ and v x, t ð Þ from f i x, t ð Þ via Equations (5) and (6). 2. Obtain the equilibrium distribution f eq i x, t ð Þ from Equation (3). 3. Perform collision (relaxation) and streaming (propagation) to update f i x, t ð Þ via Equation (1). Further details about the implementation of the 3D LB method can be found elsewhere. 16,23,24,36  For some of the performance studies discussed in Section 3.1, 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 x ¼ x, r ð Þ that describe 3D axisymmetric flow, Equation (1) can be transformed into where the a source term H i x, t ð Þ¼Δth With the inflow conditions given by the complex hybrid ODE-Dirichlet system of the LV-elastance model (Section 2.1.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. 44 Such a model has been successfully utilized for hemodynamics studies. [45][46][47][48] Other models can be used, including lumped parameter Windkessel ODEs 10,25,26 (which can be easily implemented using the same direct 0D-3D coupling introduced earlier).

| Fluid-structure interactions
For all 3D simulations of this work (including the physiological example study of Section 3.2), the solid deformation given by Equation (11) is numerically simulated by the nonlinear finite element method (FEM) solver of previous F I G U R E 4 A flowchart describing the implementation of the particular hybrid ODE-Dirichlet 0D LV-elastance model that is of interest study, 49 where 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 fluidstructure interaction (FSI) problems incorporating elastic structures. 39,[50][51][52] 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). 53 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 elsewhere. 49 For the 3D-axisymmetric performance studies included in Section 3.1, a self-implemented staggered grid finite difference (SGFD) methodology 54 is employed in the Lagrangian coordinate system (where s ¼ s ℝ is the arc length), where only the tension force given by is defined on the interface (the displacement variable X s, t ð Þ, for example, is defined on all the nodes). The solid deformation governed by Equation (11) is subsequently solved by such a finite difference methodology in a strong form. [54][55][56] 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 where the same difference approximation is applied for the time derivative. The tension force term of Equation (11) (given by Equation 32) is hence approximated as Similarly. the bending force term in Equation (11) 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 57 is used to couple the LB method of the fluid with the 3D FEM (or 3Daxisymmetric FDM) of the solid. 39,54 This method has been extensively used to simulate FSI problems in cardiovascular biomechanics. 45,[58][59][60] The body force term b x, t ð Þ 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, B s, t ð Þ, is then calculated by a penalty scheme 57 using the expression given by where α ℝ and β ℝ are negative penalty parameters (adopted in this work from previous studies 39,51,52,61,62 ); V s s, t ð Þ¼ ∂X=∂t is the velocity of Lagrangian material point of the solid wall; and V f s, t ð Þ is the fluid velocity at the position X. The latter can be obtained through transforming the Eulerian fluid velocity v x, t ð Þ into Lagrangian coordinates via the integral where δ x À X s, t ð Þ ð Þis a Dirac delta function. The body force b x,t ð Þ in Eulerian coordinates can be calculated from the corresponding Lagrangian body force via the expression The Lagrangian interaction force B 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, 39,51,61,62 including those governed by the dynamics of fluid flow over a circular flexible plate 39 and an inverted flexible plate. 51 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 U 0 ¼ 0, Q 0 ¼ 0 and an arbitrary P 0 ≠ 0, respectively.

| 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 classic case of a hanging elastic filament as well as the method of manufactured solutions (MMS).

| 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 fluidstructure coupling such that the solid position of the center of the cylinder is given by where D is the cylinder diameter; and A m ℝ, ω e ℝ are the constant oscillation amplitude and frequency, respectively. In the cases considered in this section, the Reynolds number is defined as Re ¼ vD=ν ¼ 185 for velocity v (m/s), fluid viscosity ν and D ¼ 1 m. Defining ω 0 ¼ 0:4π radians as the natural shedding frequency for a stationary cylinder, the computation is performed with the parameters A m ¼ 0:2 m, and two oscillating cases ω e =ω 0 ¼ 0:9 and ω e =ω 0 ¼ 1:1 (note that ω e =ω 0 ¼ 0:0 for a non-oscillating cylinder). These test cases and their corresponding parameters are adopted from the benchmark cases proposed by previous study. 63 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 C d , C l of the steady case as simulated by our solver up to a final time of t ¼ 400 s. The numerical discretization corresponds to Δx ¼ 1=32D, Δt ¼ 1=100s (where D ¼ 1 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 ω e =ω 0 ¼ 0:9 (left) and ω e =ω 0 ¼ 1:1 (right). As a further validation, particularly at a higher Reynolds number that better corresponds to the range for realistic blood flow configurations, Figure 8 additionally presents the drag and lift coefficients for a steady case at Re ¼ 1000. The numerical discretization in this case corresponds to Δx ¼ 1=64D, Δt ¼ 1=200s. Table 1 presents a comparison of various drag and lift-based coefficients between the simulations described above and those reported in literature. All cases (the steady and moving at Re ¼ 185 and the steady at Re ¼ 1000) demonstrate excellent agreement with previous studies. 54,63-71

| Solid solver: filament under gravity and manufactured solutions
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).
The first case considers a hanging filament of length L ¼ 1 m without an ambient fluid (in order to independently test the solid solver alone), as proposed by Huang et al. 54 Adopting the same arbitrary units, the filament is initially held stationary at an angle k ¼ 0:1π radians from the vertical, where the physical parameters correspond to a Froude number Fr ¼ v= ffiffiffiffiffi gL p ¼ 10:0 (g is the gravitation constant) and a bending rigidity EI ¼ 0:01 Pa m 3 . Snapshots of the simulated positions of the filament over a time period of 0.8 s are illustrated in Figure 9, where a spatial discretization of N ¼ 100 elements is utilized together with a timestep of Δt ¼ 0:0001 s. Again, the solutions produced by the solver of this work are in excellent agreement with those results presented by Huang et al. 54 In particular, the maximum difference in the free ends (leftmost and rightmost snapshots of Figure 9) is found to be 0:67 versus 0:68. 54 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 F I G U R E 5 A fluid-structure interaction (FSI) procedure facilitated by the immersed boundary (IB) method used for validating other hemodynamics solvers. 10,72,73 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 X s, t ð Þ as Algorithm 1 Summary of a complete 3D IB-LBM solver with 0D coupling Input fluid material parameters Re, W o , μ Input solid material parameters Eh, EI Input characteristic domain parameters D, L x i Input LV-elastance parameters C v t ð Þ, ∂C v =∂t t ð Þ Input numerical discretizations Δx i , Δs i , Δt Input total number of cardiac cycles to simulate () total time T max ) 1: Initialize the pressure P and velocity v for the fluid // initial time t ¼ 0 2: Initialize X for the solid structure position // initial time t ¼ 0 3: while t < T max do 4: Obtain macroscopic fields P and v from distribution functions f i // via Equations (5) and (6) where A ¼ 0:01 m is the maximum amplitude of displacement in the vertical component. The solid structure considered is again an elastic filament of length L ¼ 1 m with elastic material parameters corresponding to a Young's modulus of E ¼ 0:5 MPa and a thickness of h ¼ 1 mm. Employing the same numerical discretization as in the gravity case, Figure 10 (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 (N ¼ 25 elements), the simulation is advanced for 100,000 timesteps in all cases at a fine time-step size of Δt ¼ 1 Á 10 À5 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 (41), overall space and for all timesteps, are presented in Figure 10  F I G U R E 8 Numerical validation at a higher Reynolds number (Re ¼ 1000) for the steady flow (non-moving cylinder) test case, demonstrating excellent agreement with results provided in previous studies [65][66][67][68][69] (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).
T A B L E 1 Quantitative comparison between the presented fluid solver and those of literature 54,[63][64][65][66][67][68][69][70][71] for all the fluid solver validation cases of Section 3.1.1. Results demonstrate good agreement of the proposed LBM implementation for reported mean drag coefficient C d , root-mean-square averaged lift coefficient C ℓ,rms , and maximum lift coefficient C ℓ, max In order to validate coupling the fluid and solid solvers together, one can consider the classical 3D benchmark problem of a (square) flapping flag subjected to a 3D uniform flow. 56 Figure 11 (left) illustrates the problem configuration. For a simulation employing a Reynolds number of 100, a stretching coefficient K s ¼ 1000, and a bending coefficient K b ¼ 0:0001, the right plots of Figure  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 25D) presented in Figure 12. Adopting parameters of the LV-elastance hybrid ODE-dirichlet model from previous work 10 corresponding to an end-systolic LV elastance E es ¼ 2:2 mmHg/ml and a CO ¼ 4:3 L/min, Figure 13 (left) presents the expected physiologically-accurate pressure profiles at the 0D-3D interface as simulated by the complete solver for discretizations corresponding to Δx ¼ 1=20, 1=32,1=64, 1=100 and 1=128, where physical parameters of W o ¼ 16, a non-dimensionalized D ¼ 1 (corresponding to 24 mm), μ ¼ 3:5cP and ρ ¼ 1000kg=m 3 are employed. The timestep is fixed again and is taken small enough so that errors are dominated by the spatial discretization. Figure 13 (right) presents the corresponding L ∞ 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 14 (left) and Figure 15 (left) additionally present the simulated physiological flow profiles at the inlet and the pressure profiles at the midpoint of the vessel, respectively. The corresponding L ∞ errors (relative to the finest discretization of Δx ¼ 1=128) are presented in Figure 14 (right) and Figure 15 (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. 3.2 | 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. 38,75,76 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.1.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 nearwall shear stress (WSS). Indeed, mechanical experiments 77 and both in vivo/in vitro studies [78][79][80][81] 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. 77 As a demonstration of the applicability of the proposed solver toward exploring these parameters for studying pathophysiological conditions in the cardiovascular system (e.g., CHF), a computational model of a simplified 3D aorta that includes carotid and renal branches is considered and illustrated in Figure 16 (left), where the elastic wall is discretized For an effective aortic diameter D (taken to be unity in the non-dimensionalized configuration), such a domain corresponds to Cartesian coordinates given by x 0, 18D ½ ÂÀ1:5D,1:5D ½ Â0, 8D ½ for a non-dimensionalized D ¼ 1 (corresponding to 24 mm). The fluid is considered Newtonian (although the LBM is well-known to also handle non-Newtonian flow) with a Reynolds number of Re ¼ 884. For the compliant wall, a linear elastic material with Young's modulus of E ¼ 0:5 MPa and a wall thickness corresponding to h ¼ 1 mm is considered. The complete fluid-structure (immersed boundary) solver, where a no-slip condition is imposed at the fluid-structure interface, is coupled to the 0D LV-elastance heart model ( Figure 4) at the inlet and a Windkessel ODE at the outlet (see Section 2.1.2). At all peripheral branch outlets, extension tube boundary models 44 are employed. Figure 16 (right) presents the corresponding normalized velocity magnitudes produced by a simulation that employs discretizations of Δx ¼ 1=32D, Δt ¼ 1=50000 s and is advanced up to a time T ¼ 5 s (where 1 s corresponds to the period of a cardiac cycle). The LV parameters (including the compliance function) and the Windkessel lumped parameters are adopted from previous work 10 and correspond to a healthy case with normal contractility. Additionally, Figure 17 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 P v t ð Þ and the aortic pressure P x, t ð Þ at the coupled boundary. Figure 18 (left) further demonstrates that the simulations capture the expected increase in pressure amplitude as the LV-sourced waves propagate downstream. Figure 18 (right) provides the corresponding flow profiles as simulated at the inlet, midpoint and outlet of the 3D aorta, demonstrating the physiologically-expected decrease in amplitude as flow propagates downstream. For an experimental reference, Figure 19. additionally presents experimental, 82 demonstrating similar morphology characteristics and physiological ranges for both pressure and flow waveforms. Differences between the simulated results of Figures 17 and 18 and the experiments of Figure 19 can be attributed primarily to differences in the physiological parameters employed (HR = 60 BPM, CO = 3.5 L/min for the simulations; HR = 75 BPM, CO = 5.0 L/min for the pressure data; and HR = 75 BPM, CO = 3.0 L/min for the flow data).
In investigating WSS (as a relevant hemodynamic biomarker in CHF 77 ), two contractility cases (representing a low flow rate and a high flow rate) can be considered employing the same Womersely number W o ¼ 16. For normal contractility, the end-systolic LV elastance E es 10 is set to 1:76 mmHg/ml, which corresponds to a cardiac output of CO ≈ 3:5 l/min and is in accordance with values employed in experimental studies. 77 For the high contractility scenario, E es = 2:75 mmHg/ml. Again, for both cases, the Womersely number W o ¼ 16 (corresponding to a heart rate of HR ¼ 60  Figure 20 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 77 via the expression where μ is the fluid viscosity (corresponding to 3.5 centipoise), v 1 is the simulated axial velocity, and x 3 is the dimension normal to the wall. As expected, 77 a negative WSS (corresponding to a retrograde flow) is evident for the normal contractility case parameters, and such effects disappear 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. 77 F I G U R E 1 7 Aortic pressure at the inlet (blue) and the corresponding ventricular pressure (dashed red) for healthy patient parameters, simulated by the 3D FSI solver. As expected, aortic inlet pressure is equal to LV pressure during the systolic phase (when the valve is opened)

| CONCLUSION
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 2.2 can F I G U R E 1 9 (Left) Experimental data from an in vitro LV-aortic simulator. 82 (Right) Experimental flow data (of a different run) from the same setup. The overall morphologies and physiological ranges are in agreement with those of the simulations presented in this work F I G U R E 2 0 Simulated wall shear stress (WSS), for both normal and high contractility cases, at a location between the midpoint and outlet of the 3D aorta model 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.