<figure><img src="image_1.png" /></figure>

HAL
open science

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

Heng Wei, Faisal Amlani, Niema Pahlevan

► To cite this version:

Heng Wei, Faisal Amlani, Niema Pahlevan. Direct 0D-3D coupling of a lattice Boltzmann methodology for fluid-structure aortic flow simulations. International Journal for Numerical Methods in Biomedical Engineering, 2023, 39 (5), pp.e3683. 10.1002/cnm.3683. hal-03935439

HAL Id: hal-03935439

https://hal.science/hal-03935439

Submitted on 24 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Direct 0D-3D coupling of a lattice Boltzmann methodology for fluid-structure aortic flow simulations

Heng Wei<fnref n="1" /> | Faisal Amlani<fnref n="1" />,<fnref n="2" /> | Niema M. Pahlevan<fnref n="1" />,<fnref n="3" />

**Correspondence**

Niema M. Pahlevan, Department of
Aerospace and Mechanical Engineering,
University of Southern California, Los
Angeles, CA, USA.

Email: pahlevan@usc.edu

**Abstract**

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.

**KEYWORDS**

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.<fnref n="1" /> Computational fluid dynamics (CFD) is a widely-used approach for simulating blood flow in the circulatory system,<fnref n="2" />⁻¹¹ which includes applications to 1D,¹⁰⁻¹² 2D¹³ or 3D⁵⁻⁹,¹⁴,¹⁵ formulations. The lattice Boltzmann (LB) method,¹⁶⁻¹⁹ 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

This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.

© 2023 The Authors. *International Journal for Numerical Methods in Biomedical Engineering* published by John Wiley & Sons Ltd.

<footnote n="1">¹Department of Aerospace and Mechanical Engineering, University of Southern California, Los Angeles, California, USA</footnote>

<footnote n="2">²Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS - Laboratoire de Mécanique Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote n="3">³School of Medicine, University of Southern California, Los Angeles, California, USA</footnote>
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.²⁰,²¹ The accuracy and usefulness of the LB method have been demonstrated in a variety of fluid dynamics problems including turbulence²² and multiphase flow.¹⁸ As highlighted in previous studies,¹⁸,¹⁹,²²⁻²⁴ 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²⁵,²⁶ or more complex hybrid ODE-Dirichlet systems such as time-varying elastance organ models¹⁰ (both known as 0D models due to the absence of spatial dependence). Previous contributions on the 0D-3D coupling for finite element methods²⁷,²⁸ have been implicit and iterative, and for lattice Boltzmann²⁹,³⁰ 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¹⁰).

Additionally, recent work³⁰ 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³⁰ (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¹⁰ (that switches between an ODE and a Dirichlet boundary condition in a “non-stationary” fashion¹⁰) in order to generate physiologically-accurate hemodynamic conditions (instead of simply assigning a given inlet flow or pressure, as is commonly done³⁰⁻³²). Such a coupled model represents the most complicated boundary configuration found in the circulatory system¹⁰: 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¹⁰,²⁵,²⁶ (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¹⁶,³³,³⁴ 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.

## 2 | MATERIALS AND METHODS

### 2.1 | 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 computational domain $\bar{\Omega} = \Omega \cup \partial\Omega_1 \cup \partial\Omega_2 \subset \mathbb{R}^3$ is presented in Figure 1 for the (interior) fluid domain $\Omega$, the solid wall $\partial\Omega_1$ and the 0D-3D coupled domain $\partial\Omega_2$.

### 2.1.1 | 3D lattice Boltzmann equations

For the fluid domain $\Omega$, 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.¹⁶ This distribution function enforces mass and momentum conservation as well as ensuring that the fluid is Galilean invariant and isotropic.³⁵ In the present work, a single-relaxation-time (SRT) incompressible LB method is used to solve the incompressible flow.³⁶ 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

$$f_i(\mathbf{x} + \mathbf{e}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{1}{\tau} [f_i(\mathbf{x}, t) - f_i^{eq}(\mathbf{x}, t)] + \Delta t F_i(\mathbf{x}, t), \quad i = 0, \dots, N_0 - 1, \quad (1)$$

where $f_i(\mathbf{x}, t)$ are distribution functions of the particles in phase space; $\mathbf{e}_i$ are discrete velocities at position $\mathbf{x}$ and time $t$; $\tau$ is a non-dimensional relaxation time; $f_i^{eq}(\mathbf{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 3D-axisymmetric cases, i.e., $N_0 = 9$). The non-dimensional relaxation time $\tau$ is related to fluid viscosity $\mu$ by the expression

$$\mu = \rho\nu = \rho c_s^2 \left( \tau - \frac{1}{2} \right) \Delta t, \quad (2)$$

where $\nu$ is the kinematic viscosity, $\rho$ is the incompressible fluid density (e.g., blood density), and $c_s = \Delta x / (\Delta t \sqrt{3})$ is the lattice sound speed. Uniform discretizations are employed throughout this work for both time ($\Delta t$) and lattice space ($\Delta x$), chosen such that $\Delta x / \Delta t = 1$ (corresponding to $c_s = 1 / \sqrt{3}$).

The equilibrium distribution functions $f_i^{eq}(\mathbf{x}, t)$ for an incompressible Lattice Boltzmann model³⁶ and the forcing terms $F_i(\mathbf{x}, t)^{37}$ are respectively defined as

$$f_i^{eq}(\mathbf{x}, t) = \frac{\omega_i P(\mathbf{x}, t)}{c_s^2} + \omega_i \rho \left[ \frac{\mathbf{e}_i \cdot \mathbf{v}(\mathbf{x}, t)}{c_s^2} + \frac{(\mathbf{e}_i \cdot \mathbf{v}(\mathbf{x}, t))^2}{2c_s^4} - \frac{\mathbf{v}(\mathbf{x}, t)^2}{2c_s^2} \right] \quad \text{and} \quad (3)$$

<figure><img src="image_3.png" /><figcaption>FIGURE 1 A representative illustration of the complete 3D computational domain defined by $\bar{\Omega} = \Omega \cup \partial\Omega_1 \cup \partial\Omega_2$, where $\Omega$ denotes the fluid interior (governed by lattice Boltzmann equations), $\partial\Omega_1$ denotes the compliant solid boundary (governed by elastodynamics PDEs and incorporated by any appropriate fluid-structure interaction algorithm), and $\partial\Omega_2$ denotes the coupled 0D-3D boundary (governed by time-dependent ODEs)</figcaption></figure>
$$F_i(\mathbf{x}, t) = \left(1 - \frac{1}{2\tau}\right) \times \omega_i \times \left[ \left( \frac{\mathbf{e}_i - \mathbf{v}(\mathbf{x}, t)}{c_s^2} + \frac{\mathbf{e}_i \cdot \mathbf{v}(\mathbf{x}, t)}{c_s^4} \mathbf{e}_i \right) \cdot \mathbf{b}(\mathbf{x}, t) \right], \quad (4)$$

where $\mathbf{v}(\mathbf{x}, t)$ is the fluid velocity; $P(\mathbf{x}, t)$ is the macroscopic pressure; $\omega_i$ are weighting factors (whose values are adopted from previous study³⁶); and $\mathbf{b}(\mathbf{x}, t)$ is the force density in Eulerian coordinates. Pressure $P(\mathbf{x}, t)$ and velocity $\mathbf{v}(\mathbf{x}, t)$ can be calculated from the distribution functions $f_i(\mathbf{x}, t)$ via the expressions

$$P(\mathbf{x}, t) = c_s^2 \sum_i f_i(\mathbf{x}, t), \quad (5)$$

$$\mathbf{v}(\mathbf{x}, t) = \frac{1}{\rho} \sum_i \mathbf{e}_i f_i(\mathbf{x}, t) + \frac{1}{2\rho} \mathbf{b}(\mathbf{x}, t) \Delta t. \quad (6)$$

The numerical implementations and algorithmic pseudocode are provided in Section 2.3.

## 2.1.2 | 0D fluid boundary equations: LV-elastance and lumped parameter models

For the fluid inlet boundary condition on the coupled 0D-3D interface $\partial\Omega_{inlet} \subset \partial\Omega_2$, when the heart valve is open, the corresponding time-dependent boundary condition for $\mathbf{x} \in \Omega_{inlet}$ that governs the dynamics of the left ventricle is given by the ODE¹⁰,³⁸

$$\frac{\partial P_v}{\partial t} = -\frac{1}{C_v(t)} \left[ \frac{\partial C_v(t)}{\partial t} P_v(t) + Q(\mathbf{x}, t) \right], \quad (7)$$

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¹⁰,³⁸). Hence once the pressure $P_v(t)$ in the ventricle is greater than that of the inlet fluid domain boundary $P(\mathbf{x} \in \partial\Omega_{inlet}, t)$, the valve opens and $P(\mathbf{x}, t) = P_v(t)$ with the corresponding flow condition $Q(t)$ given naturally by the fluid solver via the area integral given by

$$Q(t) = \iint_{\partial\Omega_2} \frac{\sum \mathbf{e}_i f_i(\mathbf{x}, t)}{\sum f_i(\mathbf{x}, t)} dA. \quad (8)$$

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¹⁰). 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 $\partial\Omega_{outlet} \subset \partial\Omega_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¹⁰ (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(\mathbf{x}, t)$ at the outlet boundary $\mathbf{x} \in \partial\Omega_{outlet}$ through an ODE given by

$$\frac{\partial P_w}{\partial t}(t) = \frac{1}{C_w Z_w} P(\mathbf{x}, t) - \frac{R_w + Z_w}{C_w R_w Z_w} P_w(t). \quad (9)$$

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(\mathbf{x}, t)$ at the outlet $\mathbf{x} \in \partial\Omega_{outlet}$ is given by
$$Q(\mathbf{x},t) = \frac{1}{Z_w} (P(\mathbf{x},t) - P_w(t)). \qquad (10)$$

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.¹⁰

### 2.1.3 | 3D elastic wall equations

In order to account for fluid-structure interactions in a vessel, the solid boundary $\partial\Omega_1$ is assumed to be a thin wall that can be described by deformation of a compliant (elastic) wall in a Lagrangian coordinate system,³⁹ that is,

$$\rho_s h \frac{\partial^2 \mathbf{X}}{\partial t^2} = \sum_{i,j=1}^{2} \frac{\partial}{\partial s_i} \left[ Eh \varphi_{ij} \left( \sqrt{\frac{\partial \mathbf{X}}{\partial s_i} \cdot \frac{\partial \mathbf{X}}{\partial s_j}} - \delta_{ij} \right) \frac{\partial \mathbf{X}}{\partial s_j} - \frac{\partial}{\partial s_j} \left( EI \gamma_{ij} \frac{\partial^2 \mathbf{X}}{\partial s_i \partial s_j} \right) \right] + \mathbf{B}(\mathbf{s}, t), \qquad (11)$$

where $\rho_s$ is the density of the solid wall; $h$ is a constant wall thickness; $\delta_{ij}$ is the Kronecker delta; $\mathbf{X}(\mathbf{s}, t) \in \mathbb{R}^3$ is the position of the solid wall; $\mathbf{s} = (s_1, s_2) \in \mathbb{R}^2$ are the Lagrangian coordinates along the solid wall; $\mathbf{B}(\mathbf{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 $(\varphi)_{ij}$ and $(\gamma)_{ij}$ represent in- and out-of-plane effects and, for a Poisson's ratio $\hat{\nu}$, are respectively given by

$$(\varphi)_{ij} = \begin{pmatrix} 1 & \frac{1}{2(1 + \hat{\nu})} \\ \frac{1}{2(1 + \hat{\nu})} & 1 \end{pmatrix} \quad \text{and} \quad (\gamma)_{ij} = \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}. \qquad (12)$$

The boundary condition of the solid wall (as a shell) at a simply-supported fixed end is given by

$$\mathbf{X} = \mathbf{X}_0, \quad \frac{\partial^2 \mathbf{X}}{\partial s_i^2} = (0,0,0)^T, \quad i=1,2, \qquad (13)$$

where $\mathbf{X}_0$ denotes the displacement coordinates of the fixed boundaries.

### 2.1.4 | Non-dimensionalization

All the above formulations for both the fluid and the solid can be non-dimensionalized via the reference quantities $\rho, U_\infty = Q/A$ (A is the inlet area of $\partial\Omega_2$) and effective length $L$ (which is equal to a diameter $D$ if the inlet is cylindrical). The corresponding non-dimensional parameters considered in our simulations are given by Reynolds number $Re = \rho U_\infty L/\mu$, Womersley number $W_o = \sqrt{\rho\omega L^2/4\mu}$, bending coefficient $K = EI/\rho U_\infty^2 L^3$, tension coefficient $S = Eh/\rho U_\infty^2 L$, and mass ratio of the solid wall to the fluid given by $M = \rho_s h/\rho L$.

## 2.2 | 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,¹⁶,³³ which has been introduced as an alternative to typical “bounce-back” methods⁴⁰ 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.³³ In particular, such a method is ideally suited for curved fluid boundaries³⁴,⁴¹ (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¹⁶,³³,³⁴). In order to update the distribution functions from a timestep $t_n$ to a timestep $t_{n+1} = t_n + \Delta t$ of a fluid point $\mathbf{x}_f \in \Omega$ that is streamed from the corresponding 0D-3D interface (boundary) node $\mathbf{x} \in \partial\Omega_2$ of the lattice (i.e., $\mathbf{x}_f = \mathbf{x} + \mathbf{e}_i\Delta t$), such a formulation can be derived by separating collision and streaming operations of Equation (1) into

$$f_i^*(\mathbf{x}) = f_i(\mathbf{x}, t_n) - \frac{1}{\tau} [f_i(\mathbf{x}, t_n) - f_i^{eq}(\mathbf{x}, t_n)] + \Delta t F_i(\mathbf{x}, t_n), \quad \mathbf{x} \in \partial\Omega_2, \quad (14)$$

and

$$f_i(\mathbf{x}_f, t_{n+1}) = f_i(\mathbf{x} + \mathbf{e}_i \Delta t, t + \Delta t) = f_i^*(\mathbf{x}), \quad \mathbf{x} \in \partial\Omega_2, \quad \mathbf{x}_f \in \Omega, \quad (15)$$

respectively. The distribution function for $\mathbf{x} \in \partial\Omega_2$ in Equation (14) can be decomposed into an equilibrium and non-equilibrium part, i.e.,

$$f_i(\mathbf{x}, t_n) = f_i^{eq}(\mathbf{x}, t_n) + f_i^{neq}(\mathbf{x}, t_n), \quad (16)$$

which, upon substitution into Equation (14) and re-ordering terms, gives the expression for the collision as

$$f_i^*(\mathbf{x}) = f_i^{eq}(\mathbf{x}, t_n) + \left(1 - \frac{1}{\tau}\right) f_i^{neq}(\mathbf{x}, t_n) + \Delta t F_i(\mathbf{x}, t_n). \quad (17)$$

The non-equilibrium component on the boundary can be approximated via a Chapman-Enskog asymptotic expansion³³ at a distance $0 < \epsilon \ll 1$ from the neighboring fluid node $\mathbf{x}_f$ as

$$f_i^{neq}(\mathbf{x}, t_n) = f_i(\mathbf{x}_f, t_n) - f_i^{eq}(\mathbf{x}_f, t_n) + \mathcal{O}(\epsilon^2), \quad (18)$$

whose second-order accuracy in $\epsilon$ is consistent with the order of the LB method. In such an asymptotic analysis, $\epsilon$ is a “smallness parameter” that label's each term's order in the dimensionless Knudsen number $Kn$,⁴² that is, $\epsilon^m = Kn^m$. Inserting the above expression into Equation (16) gives a second-order approximation of the distribution function on the boundary node $\mathbf{x}$ as

$$f_i(\mathbf{x}, t_n) = f_i^{eq}(\mathbf{x}, t_n) + f_i(\mathbf{x}_f, t_n) - f_i^{eq}(\mathbf{x}_f, t_n). \quad (19)$$

From the definition of the equilibrium distribution function given in Equation (3), the above expression yields a “post-stream” state given by

$$f_i(\mathbf{x}, t_n) = f_i^{eq}(P^*(\mathbf{x}), \mathbf{v}^*(\mathbf{x})) + f_i(\mathbf{x}_f, t_n) - f_i^{eq}(P(\mathbf{x}_f, t_n), \mathbf{v}(\mathbf{x}_f, t_n)), \quad (20)$$
where $P^*(\mathbf{x})$, $\mathbf{v}^*(\mathbf{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

$$f_i^*(\mathbf{x}) = f_i^{eq}(P^*(\mathbf{x}), \mathbf{v}^*(\mathbf{x})) + \left(1 - \frac{1}{\tau}\right) (f_i(\mathbf{x}_f, t_n)) - f_i^{eq}(P(\mathbf{x}_f, t_n), \mathbf{v}(\mathbf{x}_f, t_n)) + \Delta t F_i(\mathbf{x}_f, t_n), \quad (21)$$

and hence, from Equation (15), gives a complete update of the fluid point $\mathbf{x}_f$ from the neighboring boundary value $\mathbf{x}$ as

$$f_i(\mathbf{x}_f, t_{n+1}) = f_i^{eq}(P^*(\mathbf{x}), \mathbf{v}^*(\mathbf{x})) + \left(1 - \frac{1}{\tau}\right) (f_i(\mathbf{x}_f, t_n)) - f_i^{eq}(P(\mathbf{x}_f, t_n), \mathbf{v}(\mathbf{x}_f, t_n)) + \Delta t F_i(\mathbf{x}_f, t_n). \quad (22)$$

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 $\mathbf{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

$$(P^*(\mathbf{x}), \mathbf{v}^*(\mathbf{x})) = (P_{0D}(t_{n+1}), \mathbf{v}(\mathbf{x}_f, t_n)), \quad (23)$$

where $P^*$ comes from the evolved 0D boundary condition (e.g., an ODE), and $\mathbf{v}^*$ is approximated by the corresponding value of the streamed fluid node $\mathbf{x}_f$ at the current timestep. Similarly, if the 0D model generates a velocity $\mathbf{v}_{0D}(t)$ at the boundary, then

$$(P^*(\mathbf{x}), \mathbf{v}^*(\mathbf{x})) = (P(\mathbf{x}_f, t_n), \mathbf{v}_{0D}(t_{n+1})). \quad (24)$$

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¹⁰). 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 LV-elastance 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).

### 2.2.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.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 $\partial\Omega_2$. That is, for a point $\mathbf{x}_f \in \Omega$ neighboring a boundary node $\mathbf{x} \in \partial\Omega_2$ with a velocity solution during the systolic phase ($t \in [0, T_d]$) defined as $\mathbf{v}_0(t) = \mathbf{v}(\mathbf{x}_f, t)$, the velocity over a complete cardiac cycle of length $T$ can be expressed as

$$\mathbf{v}(t) = \begin{cases} \mathbf{v}_0(t), & t \in [0, T_d] \text{ (valve open)}, \\ 0, & t \in (T_d, T] \text{ (valve closed)}, \end{cases} \quad (25)$$
<figure><img src="image_4.png" /><figcaption>FIGURE 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 $\partial\Omega_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)</figcaption></figure>

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^\infty$ smoothing-to-zero function $S(t)$ and can be derived from an exponential or erf-based partition-of-unity as

$$S(t) = S(t, t_0, L_s) = \begin{cases} 1, & t < t_0, \\ \exp\left(\frac{2e^{-1/u(t)}}{u(t)-1}\right), & u(t) = \frac{t-t_0}{L_s}, \quad t_0 \le t \le t_0 + L_s, \\ 0, & t > t_0 + L_s, \end{cases} \quad (26)$$

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 $\Delta t$ (corresponding to an interval $L_s = N_s \Delta t$) by multiplying the velocity by the smoothing function as

$$v_{\text{new}}(t) = \begin{cases} v_0(t), & t \in [0, T_d], \\ v_0(T_d)S(t, T_d, N_s\Delta t), & t \in (T_d, T]. \end{cases} \quad (27)$$

An illustrative example of the smoothing function $S$ for closure time $t_0 = T_d = 0.5$ s, $\Delta 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

$$\lim_{t \to T_d} v_{\text{new}}(t) = v_0(T_d), \quad (28)$$

and, from the right, as
$$ 
\begin{aligned}
\lim_{t \to T_d^+} \mathbf{v}_{\text{new}}(t) &= \lim_{t \to T_d^+} \mathbf{v}_0(T_d)S(t) \\
&= \mathbf{v}_0(T_d) \lim_{t \to T_d^+} S(t) \qquad (29) \\
&= \mathbf{v}_0(T_d),
\end{aligned}
$$

where we have taken the notational license $S(t) = S(t, T_d, N_S \Delta t)$. Hence $\mathbf{v}_{\text{new}}(t)$ is continuous everywhere including at $t = T_d$, that is, $\mathbf{v}_{\text{new}}(t) \subset \mathbb{R} \to \mathbb{R}^3$ is of class $C^0([0, T])$. A flowchart summarizing the complete implementation and smooth switching of the 0D LV-elastance model with the $C^\infty$ smoothing employed in this work is presented in Figure 4.

## 2.3 | 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.

### 2.3.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:

1. Compute the macroscopic moments $P(\mathbf{x}, t)$ and $\mathbf{v}(\mathbf{x}, t)$ from $f_i(\mathbf{x}, t)$ via Equations (5) and (6).

2. Obtain the equilibrium distribution $f_i^{\text{eq}}(\mathbf{x}, t)$ from Equation (3).

3. Perform collision (relaxation) and streaming (propagation) to update $f_i(\mathbf{x}, t)$ via Equation (1).

Further details about the implementation of the 3D LB method can be found elsewhere.16,23,24,36

<figure><img src="image_5.png" /><figcaption>FIGURE 3 An illustration of the exponential-based $C^\infty$-function $S(t, t_d, W)$ given by (26) which is employed to smoothly reduce lattice Boltzmann velocity amplitudes as the valve closes in the 0D LV-elastance model (i.e., the ODE-Dirichlet switch)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 4 A flowchart describing the implementation of the particular hybrid ODE-Dirichlet 0D LV-elastance model that is of interest</figcaption></figure>

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

$$f_i(\mathbf{x} + \mathbf{e}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{1}{\tau} [f_i(\mathbf{x}, t) - f_i^{eq}(\mathbf{x}, t)] + \Delta t F_i(\mathbf{x}, t) + H_i(\mathbf{x}, t), \quad (30)$$

where the a source term $H_i(\mathbf{x}, t) = \Delta t h_i^{(1)}(\mathbf{x}, t) + \Delta t^2 h_i^{(2)}(\mathbf{x}, t)$ is incorporated for microscopic evaluation as$^{43}$

$$h_i^{(1)} = -\frac{\omega_i \rho u_r}{r} \quad \text{and} \quad h_i^{(2)} = -\omega_i \frac{3\nu}{r} [\partial_y P + \rho \partial_x u_x u_r + \rho \partial_r u_r u_r + \rho (\partial_r u_x - \partial_x u_r) e_{ix}]. \quad (31)$$

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-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).

## 2.3.2 | 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
study,⁴⁹ 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 fluid-structure interaction (FSI) problems incorporating elastic structures.³⁹,⁵⁰⁻⁵² 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).⁵³ 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.⁴⁹

For the 3D-axisymmetric performance studies included in Section 3.1, a self-implemented staggered grid finite difference (SGFD) methodology⁵⁴ is employed in the Lagrangian coordinate system (where $\mathbf{s}=s \in \mathbb{R}$ is the arc length), where only the tension force given by

$$ \mathbf{T} = Eh \left( \sqrt{\frac{\partial \mathbf{X}}{\partial s} \cdot \frac{\partial \mathbf{X}}{\partial s}} - 1 \right) \qquad (32) $$

is defined on the interface (the displacement variable $\mathbf{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.⁵⁴⁻⁵⁶ That is, for an arbitrary variable, the central, downwind and upwind difference approximations to the first-order derivatives, are given by

$$ \left( \begin{array}{l} D_s^0 \mathbf{X} = (\mathbf{X}(s+\Delta s/2) - \mathbf{X}(s-\Delta s/2))/\Delta s, \\ D_s^+ \mathbf{X} = (\mathbf{X}(s+\Delta s) - \mathbf{X}(s))/\Delta s, \\ D_s^- \mathbf{X} = (\mathbf{X}(s) - \mathbf{X}(s-\Delta s))/\Delta s, \end{array} \right) \qquad (33) $$

such that the corresponding second-order central difference approximation can be defined as

$$ D_s^+ D_s^- \mathbf{X} = (\mathbf{X}(s+\Delta s) - 2\mathbf{X}(s) + \mathbf{X}(s-\Delta s))/\Delta s^2, \qquad (34) $$

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

$$ D_s(\mathbf{T}D_s\mathbf{X}) = D_s^0(\mathbf{T}D_s^0\mathbf{X}) = \frac{\mathbf{T}D_s^0\mathbf{X}_{s+\Delta s/2} - \mathbf{T}D_s^0\mathbf{X}_{s-\Delta s/2}}{\Delta s}. \qquad (35) $$

Similarly, the bending force term in Equation (11) can be approximated as

$$ -D_{ss}(EID_{ss}\mathbf{X}) = -D_s^+ D_s^- (EID_s^+ D_s^- \mathbf{X}) = -EI \frac{D_s^+ D_s^- \mathbf{X}_{s+\Delta s} - 2D_s^+ D_s^- \mathbf{X}_s + D_s^+ D_s^- \mathbf{X}_{s-\Delta s}}{\Delta s^2}. \qquad (36) $$

For coupling the fluid and solid systems, any suitable FSI coupling strategy can be employed. For this particular work, the immersed boundary (IB) method⁵⁷ is used to couple the LB method of the fluid with the 3D FEM (or 3D-axisymmetric FDM) of the solid.³⁹,⁵⁴ This method has been extensively used to simulate FSI problems in cardiovascular biomechanics.⁴⁵,⁵⁸⁻⁶⁰ The body force term $\mathbf{b}(\mathbf{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, $\mathbf{B}(\mathbf{s}, t)$, is then calculated by a penalty scheme⁵⁷ using the expression given by

$$ \mathbf{B}(\mathbf{s}, t) = \alpha \left( \int_{0}^{t} \left( \mathbf{V}_{f}(\mathbf{s}, t') - \mathbf{V}_{s}(\mathbf{s}, t') \right) dt' \right) + \beta (\mathbf{V}_{f}(\mathbf{s}, t) - \mathbf{V}_{s}(\mathbf{s}, t)), \qquad (37) $$

where $\alpha \in \mathbb{R}$ and $\beta \in \mathbb{R}$ are negative penalty parameters (adopted in this work from previous studies³⁹,⁵¹,⁵²,⁶¹,⁶²); $\mathbf{V}_s(\mathbf{s}, t) = \partial\mathbf{X}/\partial t$ is the velocity of Lagrangian material point of the solid wall; and $\mathbf{V}_f(\mathbf{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

$$V_f(\mathbf{s}, t) = \int \mathbf{v}(\mathbf{x}, t) \delta(\mathbf{x} - \mathbf{X}(\mathbf{s}, t)) d\mathbf{x}, \qquad (38)$$

where $\delta(\mathbf{x} - \mathbf{X}(\mathbf{s}, t))$ is a Dirac delta function. The body force $\mathbf{b}(\mathbf{x}, t)$ in Eulerian coordinates can be calculated from the corresponding Lagrangian body force via the expression

$$\mathbf{b}(\mathbf{x}, t) = - \int \mathbf{B}(\mathbf{s}, t) \delta(\mathbf{x} - \mathbf{X}(\mathbf{s}, t)) d\mathbf{s}. \qquad (39)$$

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,³⁹,⁵¹,⁶¹,⁶² including those governed by the dynamics of fluid flow over a circular flexible plate³⁹ and an inverted flexible plate.⁵¹ 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 \neq 0$, respectively.

# 3 | RESULTS AND DISCUSSION

## 3.1 | 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).

### 3.1.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

$$X_c = (-10D, A_m \cos(\omega_e t))^T, \qquad (40)$$

where *D* is the cylinder diameter; and $A_m \in \mathbb{R}$, $\omega_e \in \mathbb{R}$ are the constant oscillation amplitude and frequency, respectively. In the cases considered in this section, the Reynolds number is defined as $Re = \nu D/\nu = 185$ for velocity $\nu$ (m/s), fluid viscosity $\nu$ and $D=1$ m. Defining $\omega_0 = 0.4\pi$ 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 $\omega_e/\omega_0 = 0.9$ and $\omega_e/\omega_0 = 1.1$ (note that $\omega_e/\omega_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.⁶³ 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 $\Delta x = 1/32D$, $\Delta 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 $\omega_e/\omega_0 = 0.9$ (left) and $\omega_e/\omega_0 = 1.1$ (right). As a further validation, particularly at a higher Reynolds number that better corresponds to the
<figure><img src="image_7.png" /><figcaption>FIGURE 5 A fluid-structure interaction (FSI) procedure facilitated by the immersed boundary (IB) method</figcaption></figure>

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 $\Delta x = 1/64D$, $\Delta 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.⁵⁴,⁶³⁻⁷¹

### 3.1.2 | 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.⁵⁴ Adopting the same arbitrary units, the filament is initially held stationary at an angle $k = 0.1\pi$ radians from the vertical, where the physical parameters correspond to a Froude number $Fr = v/\sqrt{gL} = 10.0$ (g is the gravitation constant) and a bending rigidity $EI = 0.01$ Pa m³. 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 $\Delta 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.⁵⁴ 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.⁵⁴

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
**Algorithm 1 Summary of a complete 3D IB-LBM solver with 0D coupling**

**Input** fluid material parameters Re, W₀, μ

**Input** solid material parameters Eh, EI

**Input** characteristic domain parameters D, Lₓᵢ

**Input** LV-elastance parameters C_v(t), ∂C_v/∂t(t)

**Input** numerical discretizations Δxᵢ, Δsᵢ, Δ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** // via Equations (5) and (6)

4: Obtain macroscopic fields P and v from distribution functions fᵢ // via Equation (37)

5: Compute Lagrangian interaction force **B**(**s**, t + Δt) on ∂Ω₁ via IB // via Equation (39)

6: Compute the corresponding body force **b**(**x**, t + Δt) in Ω via IB // via Equation (39)

7: Evolve ODE (i.e., 0D-model) to t + Δt. // via Figure 2

8: Couple ODE with LBM to evolve fᵢ, fₐg to t + Δt on ∂Ω₂ // via Equation (1)

9: Perform collision & streaming (LBM) to evolve fᵢ, fₐg to t + Δt in Ω // via Equation (11)

10: Evolve solid wall position X to t + Δt on ∂Ω₁ via FDM/FEM // via Equation (11)

11: Update the macroscopic variables P, v to t + Δt in Ω // via Equation (11)

12: t → t + Δt

13: **end while**

**Output** the numerical solutions for v(**x**, t), P(**x**, t), **X**(**s**, t) on Ω̄ × [0, T_max]

<figure><img src="image_8.png" /><figcaption>**FIGURE 6** Numerical validation of the steady flow (non-moving cylinder) test case. (Left) Diagram of the 2D computational domain. (Right) The corresponding temporal evolution of drag (solid blue) and lift (dotted red) coefficients at Re = 185, demonstrating excellent agreement with previous studies<sup>63,64</sup></figcaption></figure>

used for validating other hemodynamics solvers.<sup>10,72,73</sup> 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
<figure><img src="image_9.png" /><figcaption>FIGURE 7 Numerical validation of the unsteady flow (oscillating cylinder) test case. (Left) The temporal evolution of drag (solid blue) and lift (dotted red) coefficients at $Re = 185, A_m/D = 0.2$ and $f_e/f_0 = 0.9$. (Right) The temporal evolution of drag and lift coefficients at $Re = 185, A_m/D = 0.2, f_e/f_0 = 1.1$. Both simulated cases are in excellent agreement with results provided by previous studies$^{54,63,64}$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 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–69}$</figcaption></figure>

$$X(s,t) = (s, A \sin(\pi s) \cos(\pi t))^T, \quad (41)$$

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 $\Delta t = 1 \cdot 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
<table><caption>TABLE 1 Quantitative comparison between the presented fluid solver and those of literature<sup>54,63–71</sup> 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_{\ell,rms}$, and maximum lift coefficient $C_{\ell,max}$</caption><thead><tr><th></th><th>Re = 185</th><th>C<sub>d</sub></th><th>C<sub>ℓ,rms</sub></th></tr></thead><tbody><tr><td rowspan="4">Non-oscillatory</td><td>Present</td><td>1.40</td><td>0.464</td></tr><tr><td>Sun, et al.<sup>70</sup></td><td>1.37</td><td>(not reported)</td></tr><tr><td>Lu and Dalton<sup>71</sup></td><td>1.31</td><td>0.422</td></tr><tr><td>Guilmineau and Queutey<sup>63</sup></td><td>1.29</td><td>0.443</td></tr><tr><td rowspan="4">f<sub>e</sub>/f<sub>o</sub> = 0.9</td><td>Present</td><td>1.42</td><td>0.17</td></tr><tr><td>Huang, et al.<sup>54</sup></td><td>1.35</td><td>0.15</td></tr><tr><td>Guilmineau and Queutey<sup>63</sup></td><td>1.33</td><td>0.19</td></tr><tr><td>Kim and Choi<sup>64</sup></td><td>1.37</td><td>0.17</td></tr><tr><td rowspan="4">f<sub>e</sub>/f<sub>o</sub> = 1.1</td><td>Present</td><td>1.44</td><td>0.92</td></tr><tr><td>Huang, et al.<sup>54</sup></td><td>1.41</td><td>0.90</td></tr><tr><td>Guilmineau and Queutey<sup>63</sup></td><td>1.36</td><td>0.87</td></tr><tr><td>Kim and Choi<sup>64</sup></td><td>1.40</td><td>0.87</td></tr><tr><td>Re = 1000</td><td></td><td>C<sub>d</sub></td><td>C<sub>ℓ,max</sub></td></tr><tr><td rowspan="6">Non-oscillatory</td><td>Present</td><td>1.58</td><td>1.40</td></tr><tr><td>Cai, Shang-Gui, et al.<sup>65</sup></td><td>1.55</td><td>1.46</td></tr><tr><td>Mittal and Kumar<sup>66</sup></td><td>1.48</td><td>1.65</td></tr><tr><td>Apte, et al.<sup>67</sup></td><td>1.50</td><td>(not reported)</td></tr><tr><td>Mittal, et al.<sup>68</sup></td><td>1.48</td><td>(not reported)</td></tr><tr><td>Mimeau, et al.<sup>69</sup></td><td>1.51</td><td>1.54</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>FIGURE 9 Numerical validation of the solid solver using an elastic filament. Superposition of simulated filament positions (from left to right) at successive times, demonstrating excellent agreement with positions provided by Huang et al<sup>54</sup></figcaption></figure>

(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).
<figure><img src="image_12.png" /><figcaption>FIGURE 10 (Left) Snapshots in time of a manufactured solution (Equation (41)) of both the exact values (dashed lines) and those produced by our numerical simulations (solid lines). (Right) The corresponding $L^\infty$ errors between the numerical and analytical solutions, which demonstrate the expected second-order accuracy of the solid solver</figcaption></figure>

### 3.1.3 | Fluid-solid solver: Flapping flag

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.⁵⁶ 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 11 present the transverse displacements at the lowest corner of the trailing edge (demarcated point A) and the mid-point of the trailing edge (demarcated point B). As can be observed, results are in very good agreement with those of Huang and Sung⁵⁶ (which are overlaid in red dots). For further validations on other configurations of the complete immersed boundary FSI solver employed in this work, the reader is referred elsewhere.⁴⁵,⁵²,⁷⁴

### 3.1.4 | 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 25D) presented in Figure 12. Adopting parameters of the LV-elastance hybrid ODE-dirichlet model from previous work¹⁰ 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 $\Delta 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), $\mu = 3.5$ cP and $\rho = 1000$ kg/m³ 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^\infty$ 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^\infty$ errors (relative to the finest discretization of $\Delta 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.
<figure><img src="image_13.png" /><figcaption>FIGURE 11 (Left) Problem configuration for benchmarking the coupled FSI solver. (Right) Transverse displacements (solid lines) corresponding to locations A and B (respectively demarcated in the left figure), demonstrate excellent agreement with those values (presented as dots) obtained by Huang and Sung⁵⁶</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 12 Diagram of the simplified straight aorta test case coupled to the LV-elastance model</figcaption></figure>

## 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.³⁸,⁷⁵,⁷⁶ 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 near-wall shear stress (WSS). Indeed, mechanical experiments⁷⁷ and both in vivo/in vitro studies⁷⁸⁻⁸¹ 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.⁷⁷

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
<figure><img src="image_15.png" /><figcaption>FIGURE 13 (Left) Physiological pressure profiles at the inlet for successively-refined discretizations of a straight aorta as produced by the LV-elastance model. (Right) The corresponding $L^\infty$ errors (relative to the finest solution)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 14 (Left) Flow profiles at the inlet for successively-refined discretizations of a straight aorta as produced by the LV-elastance model. (Right) The corresponding $L^\infty$ errors (relative to the finest solution)</figcaption></figure>

by finite elements and the fluid by lattice Boltzmann as described in Section 2.3.2. For an effective aortic diameter $D$ (taken to be unity in the non-dimensionalized configuration), such a domain corresponds to Cartesian coordinates given by $\mathbf{x} \in [0,18D] \times [-1.5D,1.5D] \times [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 $\text{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 $\Delta x = 1/32D$, $\Delta 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
<figure><img src="image_17.png" /><figcaption>FIGURE 15 (Left) Pressure profiles at the midpoint for successively-refined discretizations of a straight aorta as produced by the LV-elasticity model. (Right) The corresponding $L^\infty$ errors (relative to the finest solution)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 16 (Left) Diagram of a physiologically-relevant 3D aortic domain (with carotid and renal branches) coupled to an LV-elasticity model at the aortic inlet and a lumped-parameter Windkessel model at the aortic outlet. (Right) A temporal snapshot of the normalized flow velocity magnitude produced by the solver</figcaption></figure>

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,⁸² 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<sup>77</sup>), 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.<sup>77</sup> 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><img src="image_19.png" /><figcaption>FIGURE 17 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)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 18 (Left) Pressure profiles at the inlet, midpoint and outlet of the 3D aorta model, demonstrating the expected amplification as flow propagates downstream. (Right) Corresponding flow profiles simulated at the inlet, midpoint and outlet, demonstrating the expected decrease in flow amplitude as the wave propagates downstream</figcaption></figure>

BPM<sup>77</sup>) is fixed. 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<sup>77</sup> via the expression

$$ \text{WSS}(\mathbf{x}, t) = \mu \frac{\partial v_1}{\partial x_3}(\mathbf{x}, t), \qquad (42) $$

where $\mu$ 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,<sup>77</sup> 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.<sup>77</sup>
<figure><img src="image_21.png" /><figcaption>FIGURE 19 (Left) Experimental data from an in vitro LV-aortic simulator.⁸² (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</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIGURE 20 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</figcaption></figure>

4 | 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
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.

## ACKNOWLEDGMENTS

The authors would like to acknowledge the helpful discussions with Arian Aghilinejad at the University of Southern California.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

1. Taylor CA, Figueroa C. Patient-specific modeling of cardiovascular mechanics. *Annu Rev Biomed Eng.* 2009;11:109-134.

2. Taylor CA, Steinman DA. Image-based modeling of blood flow and vessel wall dynamics: applications, methods and future directions. *Ann Biomed Eng.* 2010;38(3):1188-1203.

3. Pennati G, Migliavacca F, Dubini G, Bove EL. Modeling of systemic-to-pulmonary shunts in newborns with a univentricular circulation: state of the art and future directions. *Prog Pediatr Cardiol.* 2010;30(1-2):23-29.

4. Vignon-Clementel IE, Marsden AL, Feinstein JA. A primer on computational simulation in congenital heart disease for the clinician. *Prog Pediatr Cardiol.* 2010;30(1-2):3-13.

5. Mittal R, Simmons S, Udaykumar H. Application of large-eddy simulation to the study of pulsatile flow in amodeled arterial stenosis. *J Biomech Eng.* 2001;123(4):325-332.

6. Asgharzadeh H, Borazjani I. A non-dimensional parameter for classification of the flow in intracranial aneurysms. I simplified geometries. *Phys Fluids.* 2019;31(3):031904.

7. Asgharzadeh H, Asadi H, Meng H, Borazjani I. A non-dimensional parameter for classification of the flow in intracranial aneurysms. II patient-specific geometries. *Phys Fluids.* 2019;31(3):031905.

8. Seo JH, Bakhshaee H, Garreau G, et al. A method for the computational modeling of the physics of heart murmurs. *J Comput Phys.* 2017;336:546-568.

9. Wu J, Shadden SC. Coupled simulation of hemodynamics and vascular growth and remodeling in a subject-specific geometry. *Ann Biomed Eng.* 2015;43(7):1543-1554.

10. Amlani F, Pahlevan NM. A stable high-order FC-based methodology for hemodynamic wave propagation. *J Comput Phys.* 2020;405:109130.

11. Aghilinejad A, Amlani F, King KS, Pahlevan NM. Dynamic effects of aortic arch stiffening on pulsatile energy transmission to cerebral vasculature as a determinant of brain-heart coupling. *Sci Rep.* 2020;10(1):1-12.

12. Shi Y, Lawford P, Hose R. Review of zero-D and 1-D models of blood flow in the cardiovascular system. *Biomed Eng Online.* 2011;10(1):1-38.

13. Tasciyan TA, Banerjee R, Cho YI, Kim R. Two-dimensional pulsatile hemodynamic analysis in the magnetic resonance angiography interpretation of a stenosed carotid arterial bifurcation. *Med Phys.* 1993;20(4):1059-1070.

14. Steinman DA, Taylor CA. Flow imaging and computing: large artery hemodynamics. *Ann Biomed Eng.* 2005;33(12):1704-1709.

15. Steinman DA. Image-based computational fluid dynamics modeling in realistic arterial geometries. *Ann Biomed Eng.* 2002;30(4):483-497.

16. Krueger T, Kusumaatmaja H, Kuzmin A, Shardt O, Silva G, Viggen EM. The lattice Boltzmann method: principles and practice. *Graduate Texts in Physics*. Springer; 2016.

17. Benzi R, Succi S, Vergassola M. The lattice Boltzmann equation: theory and applications. *Phys Rep.* 1992;222(3):145-197.

18. Shan X, Chen H. Lattice Boltzmann model for simulating flows with multiple phases and components. *Phys Rev E.* 1993;47(3):1815-1819.

19. Fang H, Wang Z, Lin Z, Liu M. Lattice Boltzmann method for simulating the viscous flow in large distensible blood vessels. *Phys Rev E.* 2002;65(5):051925.

20. Sriram K, Intaglietta M, Tartakovsky DM. Non-Newtonian flow of blood in arterioles: consequences for wall shear stress measurements. *Microcirculation.* 2014;21(7):628-639.

21. Sankar D, Hemalatha K. A non-Newtonian fluid flow model for blood flow through a catheterized artery—steady flow. *Appl Math Model.* 2007;31(9):1847-1864.

22. Cosgrove J, Buick J, Tonge S, Munro C, Greated C, Campbell D. Application of the lattice Boltzmann method to transition in oscillatory channel flow. *J Phys A Math Gen.* 2003;36(10):2609-2620.

23. Boyd J, Buick JM, Green S. Analysis of the Casson and Carreau-Yasuda non-Newtonian blood models in steady and oscillatory flows using the lattice Boltzmann method. *Phys Fluids.* 2007;19(9):093103.

24. Wei H, Cheng AL, Pahlevan NM. On the significance of blood flow shear-rate-dependency in modeling of Fontan hemodynamics. *Eur J Mech B Fluids.* 2020;84:1-14.

25. Westerhof N, Bosman F, De Vries CJ, Noordergraaf A. Analog studies of the human systemic arterial tree. *J Biomech.* 1969;2(2):121-143.
26. Vignon IE, Taylor CA. Outflow boundary conditions for one-dimensional finite element modeling of blood flow and pressure waves in arteries. *Wave Motion*. 2004;39(4):361-374.

27. Moghadam ME, Vignon-Clementel IE, Figliola R, et al. A modular numerical method for implicit 0D/3D coupling in cardiovascular finite element simulations. *J Comput Phys*. 2013;244:63-79.

28. Kim HJ, Vignon-Clementel IE, Figueroa CA, et al. On coupling a lumped parameter heart model and a three-dimensional finite element aorta model. *Ann Biomed Eng*. 2009;37(11):2153-2169.

29. Moore S, Halupka K, Zhuk S. *Towards RealTime 3D Coronary Hemodynamics Simulations During Cardiac Catheterisation* 45. IEEE; 2018:1-4.

30. Sadeghi R, Khodaei S, Ganame J, Keshavarz-Motamed Z. Towards non-invasive computational-mechanics and imaging based diagnostic framework for personalized cardiology for coarctation. *Sci Rep*. 2020;10(1):1-19.

31. Vignon-Clementel IE, Figueroa CA, Jansen KE, Taylor CA. Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. *Comput Methods Appl Mech Eng*. 2006;195(29-32):3776-3796.

32. Vignon-Clementel I, Figueroa C, Marsden A, Feinstein J, Jansen K, Taylor C. Outflow boundary conditions for three dimensional simulations of non-periodic blood flow and pressure fields in deformable arteries. *J Biomech*. 2006;39:S431.

33. Zhao-Li G, Chu-Guang Z, Bao-Chang S. Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. *Chin Phys*. 2002;11(4):366-374.

34. Kang X, Liao Q, Zhu X, Yang Y. Non-equilibrium extrapolation method in the lattice Boltzmann simulations of flows with curved boundaries (non-equilibrium extrapolation of LBM). *Appl Therm Eng*. 2010;30(13):1790-1796.

35. Wolf-Gladrow DA. *Lattice-Gas Cellular Automata and Lattice Boltzmann Models: an Introduction*. Springer; 2004.

36. He X, Luo LS. Lattice Boltzmann model for the incompressible Navier-Stokes equation. *J Stat Phys*. 1997;88(3-4):927-944.

37. He X, Shan X, Doolen GD. Discrete Boltzmann equation model for nonideal gases. *Phys Rev E*. 1998;57(1):R13-R16.

38. Berger DS, Li JK, Noordergraaf A. Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study. *Am J Phys Heart Circ Phys*. 1994;266(4):H1626-H1642.

39. Hua RN, Zhu L, Lu XY. Dynamics of fluid flow over a circular flexible plate. *J Fluid Mech*. 2014;759:56-72.

40. Cornubert R, d'Humières D, Levermore D. A Knudsen layer theory for lattice gases. *Phys D: Nonlinear Phenom*. 1991;47(1-2):241-259.

41. Guo Z, Zheng C, Shi B. An extrapolation method for boundary conditions in lattice Boltzmann method. *Phys Fluids*. 2002;14(6):2007-2010.

42. Krüger T, Kusumaatmaja H, Kuzmin A, Shardt O, Silva G, Viggen EM. *The Lattice Boltzmann Method: Principles and Practice*. Springer; 2016.

43. Lee TS, Huang H, Shu C. An axisymmetric incompressible lattice BGK model for simulation of the pulsatile flow in a circular pipe. *Int J Numer Methods Fluids*. 2005;49(1):99-116.

44. Pahlevan NM, Amlani F, Gorji MH, Hussain F, Gharib M. A physiologically relevant, simple outflow boundary model for truncated vasculature. *Ann Biomed Eng*. 2011;39(5):1470-1481.

45. Wei H, Herrington CS, Cleveland JD, Starnes VA, Pahlevan NM. Hemodynamically efficient artificial right atrium design for ventricular heart patients. *Phys Rev Fluids*. 2021;6:123103.

46. Aghilinejad A, Alavi R, Rogers B, Amlani F, Pahlevan NM. Effects of vessel wall mechanics on non-invasive evaluation of cardiovascular intrinsic frequencies. *J Biomech*. 2021;129:110852.

47. Pahlevan NM, Tavallali P, Rinderknecht DG, et al. Intrinsic frequency for a systems approach to haemodynamic waveform analysis with clinical applications. *J R Soc Interface*. 2014;11(98):20140617.

48. Kang J, Aghilinejad A, Pahlevan NM. On the accuracy of displacement-based wave intensity analysis: effect of vessel wall viscoelasticity and nonlinearity. *PLoS One*. 2019;14(11):e0224390.

49. Doyle JF. *Nonlinear Analysis of Thin-Walled Structures: Statics, Dynamics, and Stability*. Springer Science & BusinessMedia; 2013.

50. Dai H, Luo H, Doyle JF. Dynamic pitching of an elastic rectangular wing in hovering motion. *J Fluid Mech*. 2012;693:473-499.

51. Tang C, Liu NS, Lu XY. Dynamics of an inverted flexible plate in a uniform flow. *Phys Fluids*. 2015;27(7):073601.

52. Huang H, Wei H, Lu XY. Coupling performance of tandem flexible inverted flags in a uniform flow. *J Fluid Mech*. 2018;837:461-476.

53. Batoz JL, Bathe KJ, Ho LW. A study of three-node triangular plate bending elements. *Int J Numer Methods Eng*. 1980;15(12):1771-1812.

54. Huang WX, Shin SJ, Sung HJ. Simulation of flexible filaments in a uniform flow by the immersed boundary method. *J Comput Phys*. 2007;226(2):2206-2228.

55. Connell BS, Yue DK. Flapping dynamics of a flag in a uniform stream. *J Fluid Mech*. 2007;581:33-67.

56. Huang WX, Sung HJ. Three-dimensional simulation of a flapping flag in a uniform flow. *J Fluid Mech*. 2010;653:301-336.

57. Goldstein D, Handler R, Sirovich L. Modeling a no-slip flow boundary with an external force field. *J Comput Phys*. 1993;105(2):354-366.

58. Peskin CS. The immersed boundary method. *Acta Numer*. 2002;11:479-517.

59. Mittal R, Iaccarino G. Immersed boundary methods. *Annu Rev Fluid Mech*. 2005;37:239-261.

60. Seo JH, Vedula V, Abraham T, Mittal R. Multiphysics computational models for cardiac flow and virtual cardiography. *Int J Numer Method Biomed Eng*. 2013;29(8):850-869.

61. Ye H, Wei H, Huang H, Xy L. Two tandem flexible loops in a viscous flow. *Phys Fluids*. 2017;29(2):021902.

62. Zhang C, Huang H, Lu XY. Free locomotion of a flexible plate near the ground. *Phys Fluids*. 2017;29(4):041903.

63. Guilmineau E, Queutey P. A numerical simulation of vortex shedding from an oscillating circular cylinder. *J Fluids Struct*. 2002;16(6):773-794.
64. Kim D, Choi H. Immersed boundary method for flow around an arbitrarily moving body. *J Comput Phys*. 2006;212(2):662-680.

65. Cai SG, Ouahsine A FJ, Hoarau Y. Moving immersed boundary method. *Int J Numer Methods Fluids*. 2017;85(5):288-323.

66. Mittal S, Kumar V. Flow-induced vibrations of a light circular cylinder at Reynolds numbers 103 to 104. *J Sound Vib*. 2001;245(5):923-946.

67. Apte SV, Martin M, Patankar NA. A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows. *J Comput Phys*. 2009;228(8):2712-2738.

68. Mittal R, Dong H, Bozkurttas M, Najjar F, Vargas A, Von Loebbecke A. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. *J Comput Phys*. 2008;227(10):4825-4852.

69. Mimeau C, Gallizio F, Cottet GH, Mortazavi I. Vortex penalization method for bluff body flows. *Int J Numer Methods Fluids*. 2015;79(2):55-83.

70. Sun X, Zhang JZ, Ren XL. Characteristic-based split (CBS) finite element method for incompressible viscous flow with moving boundaries. *Eng Appl Comput Fluid Mech*. 2012;6(3):461-474.

71. Lu XY, Dalton C. Calculation of the timing of vortex formation from an oscillating cylinder. *J Fluids Struct*. 1996;10(5):527-541.

72. Raghu R, Taylor C. Verification of a one-dimensional finite element method for modeling blood flow in the cardiovascular system incorporating a viscoelastic wall model. *Finite Elem Anal Des*. 2011;47(6):586-592.

73. Raghu R, Vignon-Clementel IE, Figueroa CA, Taylor CA. Comparative study of viscoelastic arterial wall models in nonlinear one-dimensional finite element simulations of blood flow. *J Biomech Eng*. 2011;133(8):081003.

74. Aghilinejad A, Wei H, Magee GA, Pahlevan NM. Model-based fluid-structure interaction approach for evaluation of thoracic endovascular aortic repair endograft length in type B aortic dissection. *Front Bioeng Biotechnol*. 2022;10:825015.

75. Berger DS, Robinson KA, Shroff SG. Wave propagation in coupled left ventricle-arterial system. *Hypertension*. 1996;27(5):1079-1089.

76. Campbell KB, Lee LC, Frasch HF, Noordergraaf A. Pulse reflection sites and effective length of the arterial system. *Am J Phys Heart Circ Phys*. 1989;256(6):H1684-H1689.

77. Gharib M, Beizaie M. Correlation between negative near-wall shear stress in human aorta and various stages of congestive heart failure. *Ann Biomed Eng*. 2003;31(6):678-685.

78. Moore JE Jr, Xu C, Glagov S, Zarins CK, Ku DN. Fluid wall shear stress measurements in a model of the human abdominal aorta: oscillatory behavior and relationship to atherosclerosis. *Atherosclerosis*. 1994;110(2):225-240.

79. Pedersen EM, Sung HW, Burlson AC, Yoganathan AP. Two-dimensional velocity measurements in a pulsatile flow model of the normal abdominal aorta simulating different hemodynamic conditions. *J Biomech*. 1993;26(10):1237-1247.

80. Oyre S, Pedersen E, Ringgaard S, Boesiger P, Paaske W. In vivo wall shear stress measured by magnetic resonance velocity mapping in the normal human abdominal aorta. *Eur J Vasc Endovasc Surg*. 1997;13(3):263-271.

81. Pedersen E, Kozerke S, Ringgaard S, Scheidegger M, Boesiger P. Quantitative abdominal aortic flow measurements at controlled levels of ergometer exercise. *Magn Reson Imaging*. 1999;17(4):489-494.

82. Alavi R, Aghilinejad A, Wei H, Niroumandi S, Wieman S, Pahlevan NM. A coupled atrioventricular-aortic setup for in vitro hemodynamic study of the systemic circulation: design, fabrication, and physiological relevancy. *PLoS One*. 2022;17(11):e0267765.

**How to cite this article:** Wei H, Amlani F, Pahlevan NM. Direct 0D-3D coupling of a lattice Boltzmann methodology for fluid-structure aortic flow simulations. *Int J Numer Methods Biomed Engng*. 2023;e3683. doi:10.1002/cnm.3683