Simulation of the dynamic flow field in the left ventricle of the heart during diastolic filling

Fluid–structure interaction (FSI) problems in modeling flexible structures and moving boundaries can be simulated by combining the immersed boundary (IB) method and the Lattice–Boltzmann (LB) method. In this paper, we propose an FSI model of mitral valve leaflets and the left ventricle to simulate the flow field in the left ventricle during diastolic filling based on an IB–LB coupling scheme. Changes in the flow field and mitral valve leaflets in the FSI model are simulated with the IB–LB method, in which two arc-shaped flexible fibers simulate mitral valve leaflets immersed in the flow field. At the same time, a semi-elliptical model is used to simulate the left ventricle, which is simplified as a rigid boundary. The LB method is used to solve the Newtonian flow field, and the IB method is used to simulate the deformation of the flexible fiber interacting with the flow. In this paper, we introduce the basic principles underlying the combination of the LB and the IB methods and elucidate the coupling frame and the left ventricular flow model in detail. Finally, we verify the effectiveness of the coupled models by simulating the effects of diastolic jet flow on the motion of the mitral leaflets in the simplified left ventricular flow model.


I. INTRODUCTION
Cardiology is a bioscience of fluid dynamics, and the primary function of the cardiovascular system is to facilitate the perfusion of blood to the many and varied tissues and organs in the human circulatory system. 1 The anatomy of the heart is complex, functioning as a pump to direct blood flow. The flow field resulting from the pulsatile action of the heart is complicated and varied. Although cardiac blood flow is principally composed of laminar flow, with elements of non-laminar vortex flow seen in certain phases of the cardiac cycle, within every cardiac cycle, certain consistent characteristics of blood flow within a specific phase are stable and reproducible. These characteristics, which include current vortex parameters, flow velocity distribution, and pressure distribution, are unvaried at a specific time in the cycle. However, cardiac pathophysiology alters these cycles and patterns of cardiac flow. To detect the abnormal flow field caused by heart disease, it is very important to understand and characterize the change in and the spatiotemporal distribution of blood flow in the heart. 2 Within either an entire cardiac cycle or certain phases of the cardiac cycle, particular heart diseases will present particular abnormalities in patterns of cardiac flow. Therefore, an in-depth study of the flow field and characteristic changes in flow that correlate with the progression of different diseases of the heart may provide important diagnostic clues, creating a potential tool for early accurate diagnosis of cardiovascular diseases.
Numerical simulation technology is inexpensive and highly efficient and can specifically study the effects of various cardiac geometric deformations, mechanical factors, and fluid-solid coupling factors on hemodynamics. Such simulations can also be used to demonstrate multi-parameter, spatiotemporal distributions, expressing multifactorial results visually based on changing biomechanics that conventional physical experiments cannot. When in vivo or in vitro experimental studies are not feasible, numerical simulation is the only feasible method for observing and studying hemodynamics. Although numerical simulation analyses of hemodynamics alone cannot explain the pathophysiological mechanisms of cardiovascular disease, these in silico results can amply supplement physical experimental research and are helpful in revealing the pathological physiology and progression of cardiovascular disease, providing valuable information for the prevention and treatment of cardiovascular disease.
A principal challenge of numerical simulation is accurately modeling complex solid-fluid boundary conditions and movements. It has been proven that ensuring the accuracy and efficiency of traditional computations of very complicated boundary movements is difficult. In 1972, Peskin proposed a numerical method for addressing problems of the fluid-structure interaction (FSI) called the "immersed boundary" (IB) method; it was applied to the mitral valve motion, simulating the change in the flow field in the heart. 3 The IB method has since been widely used for numerical simulation in the study of biological system-fluid interactions. This technique has been successfully used to simulate the deformability and aggregation of red blood cells, 4-7 the deformation and detachment of biofilms, [8][9][10] insect flight, 11,12 aquatic animal locomotion, [13][14][15] and change in the flow field in a three-dimensional human heart. [16][17][18] In addition, deformations and trajectories of elastic fresh tea leaves in a simple straight channel model have been investigated using a combined IB-Lattice-Boltzmann (LB) method. 19 Chunlin Gong and Zhe Fang introduced a nonlinear finite element method into the IB-LB framework to simulate nonlinear fluid-solid interactions for moving deformable objects in incompressible fluid flow, which not only satisfies the nonslip boundary condition well at the boundary points but also has better accuracy for capturing nonlinearity because of its mature nonlinear finite element method solver. 20 The basic idea behind the IB method is to define the boundary of a complex structure into internal forces in a flow field using a simple Cartesian grid foundation with the boundary of the structure being connected in sequence by movable points. An essential feature is the use of Cartesian grids throughout, rather than a complex of body grids, to calculate the flow field. The modeling of IBs as a force source rather than a boundary condition in the flow field enables the use of Cartesian grids for calculations. Of course, although the IB method has many advantages, there are some shortcomings, 21 including inefficiencies in fluid equation solutions, that remain to be improved. 22 Because the IB method's fluid equations must be solved in a regular Cartesian grid, one advantageous approach would be to take the LB method as a substitute for the finite-difference fluid flow scheme. 23 The advantages of the LB method include the use of relatively simple algorithms, high parallelism, and the ability to handle complex boundary conditions. [24][25][26][27][28] Because the LB method is also a regular lattice-based scheme, it can be effectively combined with the IB method, as initially described by Feng and Michaelides in 2004 as the IB-LB method, 29,30 to deal with FSI problems. Numerous incremental modifications of this approach have been developed. The multirelaxation-time collision scheme proposed by Peng et al., which replaces the single-relaxation-time collision scheme used in the original LB-IB method, reduces computational effort. 31 Niu et al. proposed a momentum exchange-based LB-IB method to improve the calculation of boundary forces on fluid flow. 32 Shu et al. proposed an improved IB velocity correction approach that directly corrects the velocity to better describe the physical boundary conditions. 33 Dupuis et al. investigated two alternatives for coupling boundary force terms with grid nodes, namely, direct forcing and interpolated forcing, validating using these simulations of incompressible flow past an impulsively started cylinder. 34 Zhang et al. developed an LB-IB approach to simulate deformable capsules in flow. 35,36 Furthermore, Hu et al. proposed an efficient LB-IB method for treating three different thermal boundary conditions (Dirichlet, Neumann, and Robin) in thermal flow problems. 37 To date, only a few IB-LB coupling schemes for flexible structure problems are available, and even fewer schemes are available for FSI simulations of intracardiac flow.
In this paper, we present a simplified FSI simulation model of mitral valve leaflets and a rigid left ventricle. The IB-LB method is used to simulate the change in the flow field in the left ventricle and the dynamics of the mitral valve in a rectangular channel, in which two arc-shaped flexible fibers immersed in the flow field simulate the mitral valve. One end of the arc-shaped flexible fiber is fixed, and the other end moves with the blood flow. Simultaneously, a semielliptical model is used to simulate the left ventricle, which is simplified as a rigid boundary. The IB-LB method is used to calculate the deformation of the arc-shaped mitral valve and its interaction with the fluid. The IB-LB method is introduced and described in detail, with the numerical results for the flow field in a two-dimensional left ventricle.

A. LB method
Here, we use a combination of the Bhatnagar-Gross-Krook (BGK) LB method and the IB method to describe the motion of a fluid and the interaction between the fluid and a boundary. 38 The BGK model 39 is an approximation of the collision term of the Boltzmann equation. In the BGK model, the nonlinear collision term of the Boltzmann equation is replaced by a simpler term, which makes the derivation of the transport equation for macroscopic variables much easier. A problem that is easily solved by the BGK model is that of the relaxation of a fluid state to equilibrium. 40 At present, the lattice models commonly used in the LB method are D2Q9, D3Q13, D3Q15, and D3Q19, among others. The D2Q9 model is used to discretize the BGK equation. In the D2Q9 model, D2 represents two-dimensional space (D3 represents threedimensional space in the D3Q13, D3Q15, and D3Q19 models), and Q9 represents the nine different directions in which particles can migrate at each grid node, as shown in Fig. 1. The red dots in Fig. 1 represent LB method grid nodes, and e0-e8, respectively, represent the nine different directions of particle motion. The particle's initial location is represented on the grid node by e0, and e1-e8 represent the direction along which the vector migrates to the neighboring grid node. The vector length is equal to the distance between the current grid node and the adjacent grid node.
Solving a flow field problem using the LB method consists of two key steps: first, the collision and then, the migration of particles. 41 In the collision step, the equilibrium density distribution function of each grid node in all directions can be calculated with the LB equation according to the current particle density distribution. In the subsequent migration step, the particle density distribution function moves in a fixed direction to adjacent grid nodes. In the D2Q9 model, the particles of all grid nodes migrate in nine directions and arrive at neighboring grid nodes at the same time after the migration step; after a recurring collision to achieve balance, the process is repeated. Throughout the process, no particle remains in between any two grid nodes after migration, which is a fundamental assumption of the LB method. In the D2Q9 model, the speed of e0-e8 is defined as where c is the lattice velocity, c = δx/δt, and δx and δt are the grid step and time step, respectively. Herein, δx = δt. The LB method is derived from lattice gas automata, and the evolution equation is shown as 42 where xi is a grid node, eα is as shown in Eq. (1), and α shows the direction. f α(xi, t) is the particle density distribution function on the grid, δx is the grid step, defined as the unit length of the grid, and δt is the unit time step. Ωα(f α) is the impact of the collision on the distribution function f α. Simplifying the collision terms introduces a simple model of a single relaxation factor, Ωα where f eq α is a local equilibrium distribution function, and τ is a dimensionless relaxation time, which is an important parameter related to the physical properties of the fluid in the LB method, and the viscosity of the fluid can be directly determined by τ. Because blood is an incompressible viscous fluid and the mitral valve is a flexible structure in that fluid, in our case, the dimensionless equation combining the LB equation and BGK approximation for the interaction of the fluid and the boundary is shown as where Gα is the external force term. According to Eq. (3), the change in the particle density distribution function at the grid node at position x and time t is equal to the change in the external force per unit time and the particle density distribution function after the collision tends to balance.
The local equilibrium distribution function is shown as where cs is the lattice velocity, and its value is c/ √ 3, ρ is the pressure, u is the velocity, and ωα is the direction weight.
The external force Gα is calculated by The value of ωα is as shown in the following equation: Because the particle density distribution function rather than a macroscopic physical quantity is obtained with the LB method, the macroscopic physical quantity cannot be obtained directly just as the fluid velocity u, the fluid viscosity, the fluid density ρ, and the fluid pressure P. Therefore, the particle density distribution function needs to be converted into a macro-physical quantity, as shown in the following equation: B. IB method The IB method generally uses two sets of independent grid coordinates. First, a rectangular grid is established in the flow field, in which the grid point of the flow field, called an "Euler point," is fixed; the IB is represented by a series of discrete points, "Lagrange points," that move with the deformation of the IB. A schematic of the IB is shown in Fig. 2. Figure 2(a) describes the Lagrangian IB and Eulerian-Cartesian lattice for fluid. In Fig. 2(a), Γ is the enclosed circular boundary, and Ω denotes the flow field. The specific position of Γ can be determined according to the arc length s and time t. X(s, t) represents the Lagrange vector function of the arc length s and time t. The boundary variation at X(s, t) is determined with a singular Lagrangian force density function F(s, t). F(s, t) needs to be converted into the Eulerian variable f to solve the Navier-Stokes (N-S) equation. Finally, the flow velocity and pressure are obtained by solving the N-S equation in the flow field Ω. Because the Lagrangian force F k from the boundary point k spread to surrounding Eulerian fluid nodes, the IB moves at the local fluid flow velocity, as shown by Fig. 2(b).
For two-dimensional incompressible problems, the whole scheme can be expressed with the following series of equations: dX(s, t) dt = U(X(s, t), t) = ∫ Ω u(x, t)δ(x − X(s, t))dx, where lowercase letters indicate Eulerian variables, uppercase letters indicate Lagrangian variables, x(x, y) represents the Eulerian grid coordinates, u represents the velocity vector at the Euler point, U is the velocity vector at the Lagrange point, p is the pressure of the fluid, μ is the viscosity coefficient, δ(t) is the Dirac delta function, and f(x, t) represents the force generated by the boundary. X(s, t) and F(s, t) are the coordinates of Lagrange points and the force generated by the IB, respectively. Equations (8) and (9) are commonly used incompressible viscous N-S equations. Equation (12) and the right side of Eq. (10) are the fluid-boundary interaction equations: the former extends the Lagrangian force to the Eulerian force, and the latter is used to apply velocity on the boundary to obtain the velocity U. Equation (11) shows that the concentrated force on the boundary is determined by the flexibility of the IB at time t, where S f is a boundary force function on the deformation of the boundary that it approximates Hook's law for flexible materials.

C. The modeling of a flexible muscle fiber and force computations
In the present study, to model flow-solid interactions of the mitral valve with blood flow, we define the mitral valve as a flexible muscle fiber and model the change in the flexible force in the interaction process to ensure accurate simulation of the flow field during the FSI.
Assuming that the flexible potential energy ε of the mitral valve immersed in the boundary is generated by the stretching and bending of muscle fibers, the external force F can be divided into two parts, F = Fs + F b , where Fs and F b represent the flexible force generated by stretching and bending, respectively. Similarly, the flexible potential energy ε can also be divided into two parts, ε = εs + ε b , where εs and ε b represent the potential energy generated by stretching and bending, respectively, as shown in the following equations: where Ks represents the intensity coefficient of stretching, K b represents the intensity coefficient of bending, and As and A b are the assumed flexible functions.
Fs and F b , which represent the flexible force generated by stretching and bending, respectively, are as shown in the following equation: where t is the tangential direction of the mitral valve, and n is the normal direction of the mitral valve. A ′ s and A ′ b are derivatives of As and A b , respectively.
During fluid-solid interaction between blood flow and the mitral valve, the Eulerian variable f can be calculated from the Lagrangian force F, as shown in the following equation: where α is the amplitude along the origin of the fiber, and δ(x) is the Dirac function.

D. FSI of blood flow and the mitral valve
The governing equation of the LB method in conjunction with the external force is shown as follows: 41 where f j(x, t) represents the distribution function at the j direction, cs = c/ √ 3 is the speed of sound in the model, and c is the lattice ARTICLE scitation.org/journal/adv velocity in the model; c = 0 when j = 0, c = 1 when j = 1, 2, 3, 4, and c = √ 2 when j = 5, 6, 7, 8. The result obtained by the LB method is the particle density distribution function. The macroscopic physical quantities that can be detected in usual experiments are obtained by converting the particle density distribution. Here, the macroscopic physical quantities ρ(x, t) and (ρu) (x, t) are calculated with the following equations: The stretching potential energy of mitral muscle fibers and their corresponding stretching forces are shown in the following equations: where, l = 1,2, . . ., n f , and n f is the total number of Lagrange points on the mitral muscle fibers. Similarly, the bending potential energy of the mitral muscle fibers and their corresponding bending forces are shown in the following equations: The total Lagrangian force is the sum of the stretching force and bending force, and that is F(α, t) = Fs(α, t) + F b (α, t). After the Lagrangian force spreads to the Euler grid coordinates, the converted Eulerian force is obtained as shown in the following equation: The velocity of the IB can be calculated from the velocity of the surrounding fluid in conjunction with the interpolation as shown in the following equation: The Lagrange coordinates X n+1 of the IB are updated with the following equation: Here, the main computational steps at a time increment n of the FSI between the mitral valve and blood flow based on the LB-IB method are as follows: (1) Calculate the Lagrangian flexible force F n+1 at the boundary points of the mitral valve from X n with Eqs. (21) and (23)

III. SIMULATION OF THE FLOW FIELD IN A TWO-DIMENSIONAL LEFT VENTRICLE AND RESULTS
The mitral valve of the human heart consists of four parts: the anterior and posterior lobes, the chordae, and the papillary muscles. The valve opens and closes about 3 × 10 9 times in a human life and is responsible for ensuring that oxygenated blood from the lungs is pumped through to the aorta and the rest of the body, affecting myriad downstream blood functions. If any one of these components is damaged, it changes the hemodynamics in the heart and triggers valvular prolapse or incomplete closure, contributing to further complications leading to disease. Therefore, the physiological and pathological biomechanical study of the dynamic movement of the mitral valve is of great significance. Here, we present a model and numerical simulation focusing on the FSI of the mitral valve in a rigid left ventricular wall.

A. Geometric model and simulation conditions
In this section, the IB-LB method is used to simulate the change in the flow field in a two-dimensional left ventricle during diastolic filling. The entire simulation is carried out in a regular rectangular channel that conforms well with the LB method; the size of the rectangle is 100 × 70, which is discretized to a 120 × 80 lattice for the Eulerian LB fluid flow. A simplified geometric model of the mitral valve and left ventricle is established, as shown in Fig. 3. This is a symmetrical vertical layout with a horizontal symmetry line in the middle. A no-slip condition is imposed on the moving top and bottom walls. Two arc-shaped flexible fibers immersed in the flow field are used to simulate the mitral valve; one end of the arc-shaped flexible fiber is fixed, and the other end can move with the blood flow. The mitral valve is modeled by a quarter arc, and the initial arc-shaped flexible fibers of the mitral valve are defined as in Eq. (27). At the same time, a semi-elliptical model is used to simulate the left ventricular wall, which is designated a rigid boundary. The semi-elliptical model of the left ventricle, which shapes the wall of the left ventricle, is defined in Eq. (28). The IB-LB method is used to calculate the deformation of the arc-shaped mitral valve and its interaction with blood flow.
In this model, the intersection of the arc and the wall is fixed. According to the properties of the mitral valve, the arc is divided into two parts. 43 One part is the quarter of the arc nearest to the fixed point, which is assumed to be an inflexible rigid segment to imitate the mitral valve ring. The remaining three quarters of the arc is the flexible part with stretching and bending properties that imitate the mechanical characteristics of leaflets. The left ventricular wall simulated by the semi-elliptical model is a rigid, inflexible boundary. We choose the fluid density ρ = 1000 kg/m 3 , fluid viscosity μ = 0.0033 kg/ms, and fluid velocity v = 0.8 m/s. Each mitral valve leaflet is discretized to 80 fiber segments, and the semi-elliptical flexible fiber representing the left ventricle is discretized to 220 fiber segments.

B. Results
Previous studies have shown that blood, bile, and lymph are non-Newtonian fluids, and the shear force and shear strain rate of blood are not linearly related. Here, we assume that blood is an isotropic, incompressible Newtonian fluid with a constant density and viscosity and that the ventricular wall is rigid and impermeable. Our assumption of a Newtonian fluid is meant to simplify calculations. The velocity of non-Newtonian fluid should be slightly lower than that of Newtonian fluid at the same position because of the influence of its surface viscosity and shear coupling, but the overall simulation effect is roughly the same. 44 Therefore, numerical simulation results obtained by simplifying the blood in the left ventricle as a Newtonian fluid and numerical simulation results for non-Newtonian blood fluids are basically the same in the overall outline, especially when it comes to the deformation of the mitral valve and the formation of blood flow vortices.
The cardiac cycle is the activity of the human heart from the beginning of one heartbeat to the beginning of the next. It consists of two periods: one in which the heart muscle relaxes and refills with blood, called "diastole," followed by a period of robust contraction and pumping of blood, called "systole." After emptying, the heart immediately relaxes and expands to receive another influx of blood returning from the lungs and other systems of the body, before again contracting to pump blood to the lungs and those systems. 45 In diastole, the semilunar valves close, and the atrioventricular valves open. Oxygenated blood from the pulmonary veins fills the left atrium. For a typical heart rate of 72 beats/min, each cardiac cycle lasts approximately 0.8 s, with 0.3 s in systole and 0.5 s in diastole. 46 The mitral valve, which is located between the left atrium and left ventricle, prevents oxygenated blood from flowing back into the left atrium. In this work, we dynamically simulate changes in blood flow during diastolic filling with simplified models of the mitral valve and left ventricle. The fluid-solid interaction process from early diastole to the end of diastole in this model is shown in Fig. 4, which is intercepted from the dynamic diastole simulation result. Figure 4 shows the initial grid and six sequential velocity-vector plots at different time intervals of approximately 40-100 ms dictated by crucial changes in the flow field; the entire simulated diastole lasts 480 ms.
When ventricular pressure falls below atrial pressure, the rapidly changing LA-LV pressure gradient leads to the opening of the leaflets. Figure 4(a) is the initial state showing only the structure, and Fig. 4(b) shows t = 1 ms containing the vector velocity of the blood flow. Under the influence of both pressure and blood flow, the leaflets gradually expand and are poised to open at t = 60 ms and t = 100 ms, as shown in Figs. 4(c) and 4(d), with the mitral valve ring unchanged, reflecting its inflexible state. At t = 160 ms, the leaflets have just opened, forming two vortices at the orifice of the leaflets, as shown in Fig. 4(e). At t = 240 ms, the leaflets have opened nearly completely, forming two vortices in the blood flow behind them, as shown in Fig. 4(f). Finally, Fig. 4(g) shows the state after the leaflets have fully opened at t = 300 ms, also indicating two vortices in the blood flow behind the leaflets. In addition, the process took another 180 ms to end after the mitral valve was fully opened. In this process, the flow field in the left ventricle is basically the same as shown in Fig. 4(g). Figure 4 illustrates that the numerical simulation reveals the whole dynamic process of blood flow during diastolic filling. The results of the numerical simulation are very reasonable, accurately reflect the changes in blood flow during left ventricular diastolic filling, 47,48 and are consistent with the results of other simulations. [49][50][51][52][53] The state of the mitral valve from closed to gradually fully open demonstrates the strong FSI process. During the process of simulated diastole, the mitral valve shows good elastic curve characteristics, and the deformation and elastic characteristics of the mitral valve are well simulated. After the injection of blood, first, the mitral valve protrudes outward to form a semicircular shape, and then, the two leaflets detach from each other and gradually expand. Similarly, during the simulation process, with the injection of blood, the change in the flow field in the left ventricle and the formation of the vortex are well confirmed.

IV. CONCLUSIONS
The IB method is an effective scheme for simulating the FSI and has been widely applied to different problems related to the mechanical interactions of biological fluids with solid structures. Because the LB method also uses a rectangular grid, it can be effectively combined with the IB method to deal with FSI problems. Here, we used an IB-LB scheme to simulate change in blood flow during diastolic filling using a simplified geometric model of the mitral valve and left ventricle.
In this study, the left ventricular wall was set as a fixed rigid boundary, and the mitral valve was immersed in a flow field. The construction of the mitral valve and left ventricular wall was based on an elliptical arc. Two arc-shaped flexible fibers immersed in the flow field simulated the mitral valve; one end of the arc-shaped flexible fiber was fixed, and the other end could move with the blood flow. Furthermore, a semi-elliptical model was used to simulate a simplified left ventricle, a rigid boundary. Simultaneously, we examined the flow-solid interaction of the mitral valve and blood flow, considering the characteristics of the leaflets as flexible muscle fibers, where the change in the flexible force in the interaction process was modeled to ensure accurate simulation of the flow field during the FSI. Assuming that the flexible potential energy of the leaflets immersed in the flow field could be generated by the stretching and bending of muscle fibers, the external force F was divided into two parts, F = Fs + F b , where Fs and F b represented the flexible force generated by stretching and bending, respectively. The lattice was divided with the LB method, and the Lagrange coordinates were established with the IB method. The BGK LB method was used to describe the motion of the fluid, and the D2Q9 model was used to discretize the BGK equation. The IB-LB method was used to calculate the deformation of the arc-shaped leaflets and their interaction with blood flow. The results of the numerical simulation showed that blood flow velocity, the movement of the leaflets, and the vortex flow patterns were reasonable and correctly reflected the change in blood flow during diastolic filling in the left ventricle, validating our methods.