The application of the Lattice Boltzmann method to the one-dimensional blood flow dynamics in elastic vessels

TL;DR

This paper demonstrates that the Lattice Boltzmann method can effectively model one-dimensional blood flow in elastic vessels, capturing nonlinear wave propagation and vessel stiffening with high accuracy.

Contribution

It introduces a Lattice Boltzmann model for blood flow in elastic vessels, linking kinetic theory with cardiovascular dynamics, and validates it through various test problems.

Findings

Accurate modeling of nonlinear wave propagation in elastic vessels

Effective simulation of artery stiffening effects

Good agreement with previous results in test cases

Abstract

The one-dimensional nonlinear equations for the blood flow motion in distensible vessels are considered using the kinetic approach. It is shown that the Lattice Boltzmann (LB) model for non-ideal gas is asymptotically equivalent to the blood flow equations for compliant vessels at the limit of low Knudsen numbers. The equations of state for non-ideal gas are transformed to the pressure-luminal area response. This property allows to model arbitrary pressure-luminal area relations. Several test problems are considered: the propagation of a sole nonlinear wave in an elastic vessel, the propagation of a pulse wave in a vessel with varying mechanical properties (artery stiffening) and in an artery bifurcation, in the last problem Resistor-Capacitor-Resistor (RCR) boundary conditions are considered. The comparison with the previous results show a good precision.

Figures

Equations

$$\frac{\partial A}{\partial t}+\frac{\partial Au}{\partial x}=0,\frac{\partial u}{\partial t}+\frac{\partial u^2/2}{\partial x}=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial x}+\frac{\nu }{A}\frac{{\partial }^{2}Au}{\partial x^2}+f_{dr},$$

$$p=f(A),$$

$$c_{pulse}^2(A)=\frac{A}{{\rho }_0}\frac{\partial p}{\partial A}.$$

$$\frac{\partial A}{\partial t}+\frac{\partial Au}{\partial x}=0,\frac{\partial u}{\partial t}+\frac{\partial u^2/2}{\partial x}=-\frac{c_{pulse}^2(A)}{A}\frac{\partial A}{\partial x}+\frac{\nu }{A}\frac{{\partial }^{2}Au}{\partial x^2}+f_{dr}.$$

$$f_{-1}(t+\Delta t,x-c\Delta t)-f_{-1}(t,x)=\frac{\Delta t}{\tau +\frac{\Delta t}{2}}\left(f_{-1}^{eq}(t,x)-f_{-1}(t,x)\right)+\frac{\tau \Delta t}{\tau +\frac{\Delta t}{2}}{F}_{-1}(t,x),$$

$$f_0(t+\Delta t,x)-f_0(t,x)=\frac{\Delta t}{\tau +\frac{\Delta t}{2}}\left(f_0^{eq}(t,x)-f_0(t,x)\right)+\frac{\tau \Delta t}{\tau +\frac{\Delta t}{2}}{F}_0(t,x),$$

$$f_1(t+\Delta t,x+c\Delta t)-f_1(t,x)=\frac{\Delta t}{\tau +\frac{\Delta t}{2}}\left(f_1^{eq}(t,x)-f_1(t,x)\right)+\frac{\tau \Delta t}{\tau +\frac{\Delta t}{2}}{F}_1(t,x),$$

$$f_{-1}^{eq}(t,x)=w_{-1}\rho (t,x)\left(1-3\frac{u(t,x)}{c}+3\frac{u(t,x{)}^2}{c^2}\right),$$

$$f_0^{eq}(t,x)=w_0\rho (t,x)\left(1-3\frac{u(t,x{)}^2}{2c^2}\right),$$

$$f_1^{eq}(t,x)=w_1\rho (t,x)\left(1+3\frac{u(t,x)}{c}+3\frac{u(t,x{)}^2}{c^2}\right),$$

$$\rho (t,x)=f_{-1}(t,x)+f_0(t,x)+f_{+1}(t,x),$$

$$\rho (t,x)u(t,x)=(-f_{-1}(t,x)c+f_{+1}(t,x)c)+\frac{\Delta t}{2}\sum _i{F}_i(t,x)c_i,$$

$$\rho (t,x)(u(t,x{)}^2+c_s^2)=(f_{-1}(t,x)+f_{+1}(t,x))c$$

$$c_s=\sqrt{\frac{1}{3}}c.$$

$$F_i=w_i(c_i/c_s^2)a,i=\pm 1,0,$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho u}{\partial x}=0,\frac{\partial u}{\partial t}+\frac{\partial u^2/2}{\partial x}=-\frac{c_s^2}{\rho }\frac{\partial \rho }{\partial x}+\frac{2\nu }{\rho }\frac{\partial }{\partial x}\left(\rho \frac{\partial u}{\partial x}\right)+\frac{a}{\rho }+O(u^3)+O(\Delta t^2),$$

$$F=-w_i(c_i/c_s^2)\frac{\partial U}{\partial x},U=h(\rho )-\rho c_s^2,$$

$$\frac{\partial h(A)}{\partial A}=c_{pulse}^2(A)=\frac{A}{{\rho }_0}\frac{\partial p(A)}{\partial A}$$

$$\frac{\partial U}{\partial x}=-2\Phi \frac{\partial \Phi }{\partial x}\approx -2\Phi (x)\frac{\Phi (x+\Delta x)-\Phi (x-\Delta x)}{2\Delta x},$$

$$\frac{\partial h(A)}{\partial A}\le c^2+c_s^2=\frac{4}{3}{c}^2.$$

$$f_i(t+\Delta t,x+c_i\Delta t)=CFL\left(f_i^{\ast }(t,x)+\frac{(1-CFL)}{2}\delta f_i^{\ast }(t,x)\right)+$$

$$+(1-CFL)\left(f_i^{\ast }(t,x+c_i\Delta t)-\frac{CFL}{2}\delta f_i^{\ast }(t,x+c_i\Delta t)\right),$$

$$f_{-1}(t+\Delta t,x-c\Delta t)-f_{-1}(t,x)=\frac{\Delta t}{\tau +\frac{\Delta t}{2}}\left(f_{-1}^{eq}(t,x)-f_{-1}(t,x)\right)-$$

$$-\frac{\tau}{2c^{2}\left(\tau+\frac{\Delta t}{2}\right)}\Phi(x)\bigl{(}\Phi(x+c\Delta t)-\Phi(x-c\Delta t)\bigr{)},$$

$$f_0(t+\Delta t,x)-f_0(t,x)=\frac{\Delta t}{\tau +\frac{\Delta t}{2}}\left(f_0^{eq}(t,x)-f_0(t,x)\right),$$

$$f_1(t+\Delta t,x+c\Delta t)-f_{+1}(t,x)=\frac{\Delta t}{\tau +\frac{\Delta t}{2}}\left(f_1^{eq}(t,x)-f_1(t,x)\right)+$$

$$+\frac{\tau}{2c^{2}\left(\tau+\frac{\Delta t}{2}\right)}\Phi(x)\bigl{(}\Phi(x+c\Delta t)-\Phi(x-c\Delta t)\bigr{)},$$

$$\frac{\partial h(A)}{\partial A}=c_{pulse}^2(A)=\frac{A}{{\rho }_0}\frac{\partial p(A)}{\partial A},$$

$$f_{\pm 1}^{eq}=\frac{A}{6}\left(1\pm 3\frac{u}{c}+3\frac{u^2}{c^2}\right),f_0^{eq}=\frac{4A}{6}\left(1-3\frac{u^2}{2c^2}\right),$$

$$A(t,x)=(f_{-1}+f_0+f_1)(t,x),$$

Full text

1 Introduction

The models of blood motion in cardiovascular system vary from the [math]D lumped models, $1$D pulse propagation equations to $3$D viscous flow equations [1]. In many cases the $3$D approach based on the solution of the Navier-Stokes equations is too detailed while [math]D lumped models are oversimplified and applicable only for the distal vasculature. In the case of the $1$D models it is assumed that the radial velocity is negligible [2]. Then, integrating the Navier-Stokes equations over the radial variable the $1$D nonlinear system of equations (depending only on one spatial variable, axial coordinate) for luminal area change (or pressure) and axial blood velocity is derived [3, 2]. Several numerical approaches for solving of these equations have been proposed including variations of Galerkin discretization [2, 4, 5], McCormack and Lax-Wendroff schemes [6, 7, 8, 9], trapezoidal schemes [10] and a relatively novel well-balance discretization method [11, 12]. The comparison results of the numerical methods can be found in [13, 14]. Among the open-source numerical solvers one can mention [math]D framework pyNS [15], $1$D finite-volume OpenBF solver [16], Lax-Wendroff method based $1$D solver VaMpy [17], finite element method based Artery.FE framework [18].

The hydrodynamic equations can be modeled using the kinetic methods like the Lattice Boltzmann approach [19, 20]. This method describes the motion of particles on a Cartesian spatial lattice (advection part), the collision of the particles is modeled by assuming that the velocity distribution of the particles tends to a local equilibrium state. The Lattice Bolzmann (LB) method correctly reproduces low-Mach incompressible flows like blood motion and can be used for the modeling of the flow in cardiovascular network.

The Lattice Boltzmann simulations of the blood flow dynamics in 2D and 3D vessel geometry have gained some popularity recently [21, 22, 23, 24, 25, 26]. The boundary conditions for LB models applied to blood flow problems in elastic tubes have been proposed in [21]. Several results are devoted to the estimation of blood endothelial (wall) shear stress in complicated vessel networks using LB approach, or using LB models coupled with Monte-Carlo simulations [22, 23, 24, 25, 27].

In the present paper an alternative way to model 1D blood motion equations based on the Lattice Boltzmann method is proposed. It is shown that there exists an analogy between the one-dimensional LB model D1Q3 with a virtual force (such approach is popular in modeling of non-ideal gas and multiphase flows using LB approach) at the limit of small Knudsen numbers (small times between collisions) and the blood flow equations in elastic vessels. More precisely, the LB hydrodynamics for D1Q3 model is equivalent to the 1D blood motion equations if one changes the sound velocity in the LB model by the pulse propagation velocity and the density by the vessel luminal area. The addition of the fictitious force allows to model arbitrary area-pressure vessel responses. The presented method is relatively simple in implementation.

Several test problems are considered. The method correctly describes the change in shape of the initial pulse wave while propagating along a vessel, the analytical solution profile is very close to the LB modeling results. The second test problem concerns the propagation of pulse waves with various amplitudes and duration in a vessel with a prosthesis [4, 28], in the next problem an artery network (5 vessels, 2 bifurcations) is considered [2]. Finally, the inclusion of the RCR boundary conditions for the LB model is discussed and LB simulation results are compared with McCormack difference scheme for the vessel branching supplemented with RCR boundary conditions at the distal ends of the daughter vessels.

2 Lattice Boltzmann approach for non-ideal gas and 1D blood flow equations in elastic vessels

2.1 The mass and momentum conservation equations for inviscid 1D blood flow in elastic vessels

Consider a long elastic vessel filled with incompressible fluid (blood). The conservation equations of mass and momentum in the elastic vessel are obtained from the Navier–Stokes equations by the integration over the cross-sectional spatial coordinate (radial variable). Assume that the radial component of the blood velocity is close to zero, and the velocity profile is flat. Then the following equations are obtained [2]

[TABLE]

where $\rho_{0}$ is the constant blood density (incompressible fluid), $A,u$ is the luminal area, blood velocity respectively and $f_{dr}$ is the viscous drag force. Here $x$ is the axial coordinate and $p$ is the blood pressure.

It is obvious that one has three unknowns $u,A,p$ and only the two governing equations. To close the system (1) one needs to introduce a pressure-area relation

[TABLE]

where $f$ is a some function, its form should be defined from the elastic properties of the considered vessel. For realistic vessels this dependence can be complicated and saturation effects are observed [29]. In practice the Laplace law is popular [1, 2].

From the pressure-area relation one can find the pulse-wave velocity using the formula

[TABLE]

Using the pulse wave velocity definition Eqs (1)-(2) can be rewritten in the following equivalent form

[TABLE]

In the present paper the drag force is not considered, therefore $f_{dr}=0$.

2.2 The Lattice Boltzmann method and hydrodynamics

It will be sufficient to consider the simplest three-velocity model $D1Q3$ [19]. This model describes three populations of particles traveling on a 1D lattice with equal spacing $\Delta x=c\Delta t$ between the lattice nodes. Here $c$ is the lattice velocity, $\Delta t$ is the lattice time step. During the discrete time step $\Delta t$ the particles at a node $x$ can hop to adjacent nodes $x\pm c\Delta t$ or stay at the node $x$. The concentrations of the particles traveling with the velocities $\pm c,0$ are given by $f_{\pm 1}(t,x)$ and $f_{0}(t,x)$ respectively. The dynamics of $f_{\pm 1},f_{0}$ obeys the following discrete equations in space and time

[TABLE]

[TABLE]

[TABLE]

where $\tau$ is the relaxation time which can be considered as a free parameter, $F_{i}$ are the components of an external force and $f^{eq}_{\pm 1},f^{eq}_{0}$ are the equilibrium states (analogs of the Maxwell distribution)

[TABLE]

[TABLE]

[TABLE]

where $w_{\pm 1}=1/6,w_{0}=4/6$ are the lattice weights and $\rho,u$ are the density and flow velocity

[TABLE]

[TABLE]

moreover, one can define the full energy of the lattice gas in every node by the formula

[TABLE]

and here $c_{s}$ is the sound velocity of the lattice gas defined by

[TABLE]

Here the external force is taken in the linearized form [19]

[TABLE]

where $a$ is the force magnitude.

At the limit of small values of $\tau$ the lattice gas can be considered as continuous media since the time between collisions tends to zero ( $\tau$ is proportional to the expected time between particle collisions). In the continuous limit the considered LB model describes some hydrodynamics which can be obtained using the Chapman-Enskog expansion on a small parameter, the detailed analysis of this procedure for LB method in 1D case can be found in [30]. One has

[TABLE]

where $O(u^{3})$ is the spurious term which appears due to the quadratic form of the local equilibrium on the bulk velocity $u$, $O(\Delta t^{2})$ is space-time discretization error ( LB scheme has second order accuracy) and $\nu=c_{s}^{2}\tau$ is the gas viscosity. If the flow is slow $u<<c_{s}$ (i.e the Mach number $Ma=u/c_{s}$ is small) and the time step is small then one can assert that the LB model correctly describes 1D hydrodynamics. It should be mentioned that the Navier-Stokes equations without spurious $u^{3}$ term can be obtained if high-order LB models are considered [31, 32]. In the present case the term proportional $u^{3}$ is unimportant since the blood flow is significantly subsonic.

2.3 The Lattice Boltzmann equation as a model of hemodynamics

One can see that the luminal area $A$ in the equations (4) is very similar to the density $\rho$ in the equations (5). If one changes the density $\rho$ by the luminal area $A$ and the sound velocity $c_{s}$ by the pulse wave velocity $c_{pulse}$ then Eqs (5) turn into Eqs (4) if the flow is slow, discretization error is negligible. These conditions are realistic due to a) for a typical blood flow the Mach number is over $0.1$, $u<<c_{pulse}$ [2] b) the discretization error can be controlled by the choice of $\Delta t$.

The lattice Boltzmann approach gives slightly different longitudinal viscosity term $\frac{\partial}{\partial x}\left(A\frac{\partial u}{\partial x}\right)$ than $\frac{\partial}{\partial x}\left(\frac{\partial Au}{\partial x}\right)$ in Eqs (1)-(2) . This difference is unimportant since the terms are very close (same in the linear theory), moreover stream-wise diffusion effects play small role in hemodynamics and usually omitted from the consideration. In the present case the diffusion is an inherent property of the LB method.

One should emphasize that the force free case in the equations (5) describes only the particular case of constant pulse velocity $c_{pulse}$ since the sound velocity $c_{s}$ is constant, the blood flow equations are derived for the pressure-area relation in the form $p\sim\log(A)$ (this results from the definition of the pulse wave velocity (3)). The logarithmic dependence is only qualitatively correct, a more popular Laplace law states that $p\sim\sqrt{A}$, moreover the dependence of the pressure on the luminal area in real arteries is saturating [29], hence a generalization of the model is desirable.

The idea of the generalization of the presented approach is to choose the force term in such a way that the case of the non-constant pulse velocity is captured. This approach is similar to the method of modeling of non-ideal gases and phase transitions using LB schemes [33, 34, 35, 36, 37, 38]. Following the articles [35, 36, 38] the force is introduced using a pseudo-potential $U$

[TABLE]

where $h(\rho)$ is a some function. Using the expressions (6) the terms $-\frac{c_{s}^{2}}{\rho}\frac{\partial\rho}{\partial x}+\frac{a}{\rho}$ become $-\frac{1}{\rho}\frac{\partial h}{\partial\rho}\frac{\partial\rho}{\partial x}$ in the hydrodynamic equations (5), next changing $\rho$ by $A$ and selecting $h(A)$ in a such way that

[TABLE]

where $p(A)$ is the desired pressure-area law one obtains the equations (4).

Following [36] the pseudo-potential $U$ can be written in the form $U=-\Phi^{2}$, $\Phi=\sqrt{Ac_{s}^{2}-h(A)}$ and for the LB model the central discretization of the spatial derivatives is adopted

[TABLE]

where $\Delta x=c\Delta t$ is the lattice spatial step. For the LB D1Q3 model with the force discretization in the form (7) the stability conditions are known [36]

[TABLE]

From the relation (8) one can see that stability can be controlled by adjusting the lattice velocity $c=\Delta x/\Delta t$.

The model implementation at the each time step consists of the two parts. Firstly, the collision term (right hand side of LB model) is computed in all spatial nodes and then the post-collision distribution functions $f_{i}^{*}(t,x)$ are obtained. Next, the streaming step is performed $f(t+\Delta t,x+c_{i}\Delta t)=f_{i}^{*}(t,x)$.

2.4 Stability improvement

For the conventional LB method Courant-Friedrichs-Lewy (CFL) number (defined as $c\Delta t/\Delta x$ ) equals to unity. This means that LB method is only marginally stable. This feature is well-known and several efforts have been paid to improve the stability properties, the possible solution is to rewrite the streaming step in such a form that LB model becomes a finite-difference scheme with controllable CFL number [39, 40, 41, 42, 43, 44].

The simplest way to decrease CFL number is the implementation of the fractional time step. In this method [40] the streaming part $f_{i}(t+\Delta t,x+c_{i}\Delta t)=f_{i}^{*}(t,x)$, where $c_{i}=(-c,0,c)$ and $f_{i}^{*}$ is the post-collision distribution function is replaced by

[TABLE]

[TABLE]

where $\delta f_{i}^{*}(t,x)=f_{i}^{*}(t,x+c_{i}\Delta t)-f_{i}^{*}(t,x)$. For this realization of the streaming component the particles travel a distance $CFLc\Delta t$ during a time step $\Delta t$ and the parameter $0<CFL\leq 0$ is Courant-Friedrichs-Lewy number. The form (9) guarantees second order accuracy in space and time. Finally, one should mention that for the LB model with the streaming step (9) the viscosity is governed by the relation $\nu=CFLc_{s}^{2}\tau$ [40].

In the present paper for the test problems in Paragraphs 3.1-3.3 the conventional streaming part ($CFL=1$) was applied, for the test problem in Paragraph 3.4 the fractional time step with $CFL=0.5$ was implemented.

3 The model application.

From the results of the previous paragraph one derives the following LB model

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\Phi(x)=\sqrt{Ac_{s}^{2}-h(A)}$ and $c_{s}^{2}=(c^{2}/3)$, $c=\Delta x/\Delta t$ . The function $h(A)$ is obtained from the equation

[TABLE]

here $c_{pulse}(A),p(A)$ are the target pressure pulse velocity and pressure-area relation. The local equilibrium states $f^{eq}_{\pm 1},f^{eq}_{0}$ are defined as

[TABLE]

where

[TABLE]

[TABLE]

The relaxation time $\tau$ should be taken small enough if the inviscid flow is modeled. Note that wall friction forces are not taken in account but potentially they can be implemented as an additional external force. The pulse pressure waves presented in figures are recovered with use of the pressure-area relation $p(A)$.

The case of LB model without external force corresponds to the following pressure-area relation

[TABLE]

For the hemodynamic modeling the following tube law is popular

[TABLE]

where $D_{0},p_{0},\rho_{0}$ are the vessel distensibility, blood pressure at rest, blood density and $n>0$. The case $n=1/2$ corresponds to the Laplace law. The value of $p_{0}$ can be taken arbitrary since only $p-p_{0}$ is computed. It will be more convenient to use a linearized pressure pulse velocity $c_{pulse}(A_{0})$ instead of the distensibility $D_{0}$

[TABLE]

3.1 A sole forward pulse wave propagation.

The evolution of the nonlinear wave which has triangle shaped form in the initial moment of time is considered in a semi-infinite vessel $x>0$ with constant elastic properties.

For a forward traveling wave one has the relation between the wave velocity $u$ and the luminal area $A$ [45]

[TABLE]

then the equations (1)-(2) reduce to a differential equation for $A$

[TABLE]

The equation (16) is supplemented with the symmetric triangle-shaped initial-boundary condition at $x=0$ for $t\in[0,T_{0}]$, where $T_{0}$ is one heartbeat

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

The pressure-area relation (14),$n=1/2$, which corresponds to the Laplace law is considered. The analytical solutions to this problem are obtained in [45]. For the sake of brevity this solution is not given here. The LB equations (10)-(12) are applied with the initial-boundary conditions (17)-(20) and the modeling results are compared with the analytical solutions at $t=0.35\,s$. The modeling constants are given in Fig. 1 caption. A very good agreement between the solutions is observed ( Fig. 1).

3.2 A pulse wave propagation in a vessel with prosthesis.

A semi-infinite vessel with an interior stiff region (stent or prosthesis) is considered. Similar test problem was investigated using Galerkin and Taylor-Galerkin finite difference methods in the papers [28, 4].

The mass flow and full energy are continuous across the junctions of the stiff and elastic regions [46, 2]. Then

[TABLE]

where $L$ and $R$ denote left and right side of the junction interface respectively and $\rho_{0}$ is the blood density. Since the blood flow is slow then $p>>\frac{\rho_{0}u^{2}}{2}$ and one can consider the pressure continuity property $p_{L}=p_{R}$ instead of the full energy. Consider the Laplace area-pressure relation (14), $n=1/2$. Using the linearized pressure pulse wave velocity $c_{pulse}(A_{0})=\sqrt{A_{0}^{1/2}/\rho_{0}D_{0}}$ the Laplace law reads

[TABLE]

The continuity condition at the junction is considered for the pressure-area relation linearized near $A_{0}$, i.e. $p=p_{0}+\rho_{0}c_{pulse}^{2}(A_{0})\left(\frac{A}{A_{0}}-1\right)$. In a similar way one can linearize the logarithmic tube law (13).

In terms of LB variables the continuity of the mass flux and the linearized pressure take the following form

[TABLE]

[TABLE]

[TABLE]

where the junction is placed between the lattice nodes $x_{n},x_{n+1}$; $c_{L}=\frac{\Delta x_{L}}{\Delta t},c_{R}=\frac{\Delta x_{R}}{\Delta t}$ are the LB velocities on the left and right sides from the junction; $\Delta x_{L},\Delta x_{R}$ are the LB spatial steps and $\Delta t$ is the LB time step; $a_{L},a_{R}$ are the force magnitudes on the left and right sides of the junction interface; $c_{pulse,L}^{2}(A_{0})$ and $c_{pulse,R}^{2}(A_{0})$ are the linearized pulse wave velocities. For the present problem the lattice velocities are taken as $c_{L}^{2}=3c_{pulse,L}^{2}(A_{0}),c_{R}^{2}=3c_{pulse,R}^{2}(A_{0})$. Such a calibration is convenient and guarantees that the lattice sound velocity equals to the linearized pulse wave velocity.

The distribution functions $f_{-1}(t,x_{n}),f_{1}(t,x_{n+1})$ are unknowns, they are obtained from the continuity conditions

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

For the test problems considered here it is assumed that the undisturbed vessel cross-sectional area $A_{0}$ equals $10^{-4}\,m^{2}$, the blood density $\rho_{0}$ is $10^{3}\,kg/m^{3}$.

The first vessel geometry is depicted in Fig. 2. The dynamics of a half-sinusoidal pressure pulse wave with the maximal amplitude of $0.2\,kPa$ and the duration of $0.165\,s$ starting at the point $x=0$ is studied. The pressure-area relation (13) is adopted.

The distensibility of the prosthesis is assumed to be $10$ times smaller then in the other part of the vessel, therefore the linearized pulse wave velocity is $\sqrt{10}$ is larger in the stiffening than in the remaining part of the vessel ( $4\,m/s$ and $12.6\,m/s$), Fig. 2.

The simulation results show that the pressure maximal magnitude is amplified (in the comparison with the wave at $x=0\,m$) near the left junction with a subsequent pressure drop below the undisturbed level (Fig. 3). This behavior is in a good agreement with the results from [28, 4]. The pressure increase over the initial pressure values happens due to the reflection of the wave entering the stiff region while the pressure drop is observed due to the second reflection of the wave leaving the stiff region.

The second considered vessel geometry is depicted in Fig. 4. A relatively short vessel segment (Fig. 4) and a small pulse wave amplitude $0.055\,kPa$ are taken. Such a choice guarantees that the nonlinear effects do not distort initial wave shapes and shock waves do not appear. The distensiblity of the stiff region is $100$ smaller than the elastic parts of the vessel, this means that the pulse wave velocity in the stiff region is $10$ times larger than in the the elastic segments. The Laplace pressure-area relation is adopted.

The modeling results are presented in Fig. 5. Multiple reflections propagating from the junctions are clearly visible. The positive amplitudes correspond to the reflections occurred at the left junction while the waves with negative pressure amplitudes correspond to the reflections at the right junction. Again, the wave form behavior is very similar to the results from [28, 4].

3.3 Bifurcations of arteries

Consider a bifurcation of $N+1$ vessels: an undisturbed luminal area and linearized pulse wave velocity for a parent vessel are defined as $A_{0}^{(0)},c_{pulse}(A_{0}^{(0)})$ and the same quantities for daughter vessels are defined as $A_{0}^{(i)},c_{pulse}(A_{0}^{(i)}),i=1\ldots N$. At the junction point it is assumed that the continuity of the pressure and the conservation of the flow are satisfied [2]

[TABLE]

[TABLE]

where $u^{(i)},\Delta p^{(i)}=p-p_{0}$ are the blood velocity and pressure alteration from the diastolic value $p_{0}$ in the vessels $i=0\ldots N$.

For simplicity, the pressure continuity condition is considered in the linearized form for the pressure-area relation (14), then one has $\Delta p(A)\equiv p(A)-p_{0}=\rho_{0}c_{pulse}^{2}(A_{0})(A/A_{0}-1).$ In terms of the LB variables the equations (21)-(22) read as

[TABLE]

[TABLE]

where $i=1\ldots N$ and

[TABLE]

where $c^{(i)}$ is the lattice step in $i$-s vessel, moreover

[TABLE]

where $a^{(i)},i=0\ldots N$ is the amplitude of the external force in $i$-s vessel (the definition of the external virtual force is given in (6)). The time step $\Delta t$ is constant in all vessels, while the lattice velocities $c^{(i)}=\Delta x^{(i)}/\Delta t$ can be different, here $\Delta x^{(i)}$ is the lattice spatial step. The lattice velocities $c^{(i)}$ are taken in a such way that

[TABLE]

the latter property guarantees that the pulse wave propagation speed will be different if $c_{pulse}(A_{0}^{(i)})$ are non-equal in the vessels.

At the junction point $x$ the value of $f_{-1}^{(0)}(t,x)$ is unknown for the parent vessel and the values of $f_{1}^{(i)}(t,x)$ are unknown for the daughter vessels $i1\ldots N$, they should be found from the equations (23)-(24). One has

[TABLE]

where $f^{(i)}_{1},i=1\ldots N$ are calculated from the linear system

[TABLE]

[TABLE]

here $i=1\ldots N$, and the following notations are introduced

[TABLE]

It is assumed that the wave is absorbed by the vessel boundaries (the vessels are well-matched with the distal vasculature). In order to absorb the incident wave the impedance boundary conditions should be stated [47].

The modeling artery geometry is depicted in Fig. 6. The input velocity wave is prescribed to the left inlet of the vessel $1$ and has the following form $0.25*sin(\pi t/0.33)\,\,m/s$ for $0\leq t\leq 0.33\,s$. The linearized pulse wave velocities are taken same in all vessels: in the first case the value $3\,m/s$ is used, in the second case this velocity equals $6\,m/s$. The case of the increased pulse velocity corresponds to the more stiff vessels which results to the smaller amplitudes of luminal area changes (Fig. 7-8). The Laplace pressure-area relation is adopted for the present problem. For the first (left) bifurcation the linear reflection coefficient [46] is

[TABLE]

and equals [math] for the forward wave traveling to the junction from the the vessel $1$, on the other hand, for the backward wave traveling to this junction from the vessel $2$ or $3$ this coefficient equals $-0.5$. The linear reflection coefficient at the second (right) junction equals $0.5$ for the forward wave in the vessel $3$. The backward waves at the distal ends of the vessels $4$ and $5$ do not exist since the outlets of these vessels absorb incident waves. Thus, the wave is partially ”trapped” in the vessel 3.

The wave dynamics and multiple wave reflections are measured at the middle of the vessel $3$, see Fig. 7-8. In Fig. 8 the modeling results are presented for the geometry in Fig. 6 where the vessels 4 and 5 have increased radius ($D/\sqrt{3}$). This results in smaller reflection coefficient $0.2$ and therefore smaller amplitudes of the reflected waves in comparison with the profiles in Fig. 7.

The problem with the geometry in Fig. 6 was also considered in the paper [2]. The presented area and velocity profiles are very similar to the data from [2].

3.4 RCR boundary conditions

Resistor-Capacitor-Resistor (RCR) boundary conditions [48, 49, 50] are very popular in hemodynamical modeling. These boundary conditions are analogous to electric circuit with two resistances $R_{1},R_{2}$ and one capacitor $C$ (Fig. 9). The first resistance $R_{1}$ accounts for the distal part of the truncated vessel and is usually selected in such a way that there are no wave reflection in a junction between a vessel and RCR model, while the second resistor $R_{2}$ accounts for the drag forces in the distal arterioles, the capacitor $C$ term is responsible for the compliance of the arterioles. One has for the pressures and velocities $p,q$ at the distal ends of the vessels 2,3 the following equation

[TABLE]

the first order discretization of this equation is following

[TABLE]

where $R=R_{1}+R_{2}$ and $\Delta t$ is the time step for the difference scheme. It is assumed that the values $p(t),q(t)$ are known.

Any traveling wave in a vessel can be separated into the forward and backward components by the rule

[TABLE]

Using (27) one has

[TABLE]

[TABLE]

then

[TABLE]

Next, we mention that

[TABLE]

since $p_{f}$ can be considered as the forward traveling wave with the velocity $c_{pulse}$, where $x=L$ is the distal point of the considered vessel. Therefore, we have expressed $p_{f}(t+\Delta t,L),p_{b}(t+\Delta t,L)$ via known functions $p_{f}(t,L),p_{b}(t,L),p_{f}(t,x-c_{pulse}\Delta t)$. Then

[TABLE]

[TABLE]

Now one expresses $p,q$ via the lattice Boltzmann variables

[TABLE]

[TABLE]

where $a$ is the amplitude of the external force, $C_{0}=\frac{A^{(0)}}{\rho c_{pulse}^{2}}$, where $c_{pulse},A^{(0)}$ are the linearized pulse velocity and undisturbed luminal area, $a(L)$ is approximated as $2a(L-1)-a(L-2)$. The values $f_{0}(t+\Delta t,L),f_{-1}(t+\Delta t,L)$ are unknowns at the distal end of the vessel, they are found from the equations

[TABLE]

[TABLE]

The relations (32)-(33) define RCR boundary conditions for the LB model. Finally one can see that the unknown values of LB distribution functions at the right end of the vessel (distal end) $x=L$ are defined by $p(t+\Delta t,L),q(t+\Delta t,L)$, they can be represented as a sum of the forward and backward components (30)-(31) which in turn depend on $p_{f}(t,L),p_{b}(t,L),p_{f}(t,x-c_{pulse}\Delta t)$ (the formulas (28), (29)).

As a test problem the propagation of a pulse wave in a vessel network consisting of three vessels (one parent vessel and two daughter vessels) and one bifurcation connected with RCR models at the distal outlets of the daughter vessels is considered. The parent vessel has the length $0.1\,m$ , the luminal radius equals $0.01\,m$, the linearized pulse wave velocity equals $4\,m/s$. For the daughter vessels the values of the vessel length, luminal radius, linearized pulse wave velocity are taken as $0.05\,m,0.009\,m,3\,m/s$ respectively. The Laplace law is adopted for the pressure-area relation. The initial pressure profile has the duration of $0.05\,s$ and the maximal amplitude corresponding to $10\%$ luminal area increase in the vessel 1. The value of $R_{1}$ equals $\rho_{0}c_{pulse}(A^{(0)})/A^{(0)}$, where $c_{pulse}(A^{(0)}),A^{(0)}$ are the linearized pulse velocity and undisturbed luminal area respectively. The full resistance $R_{1}+R_{2}$ and compliance $C$ are taken as $10^{8}Pa\cdot s/m^{3},10^{-9}\,m^{3}/Pa$ respectively, then the diastolic decay time $\tau\equiv R_{2}C\approx 0.09\,s$.

The choice of the modeling parameters results in relatively complicated behavior of the pressure wave: since the linear reflection coefficient is negative for the vessel 1 at the branching point ($\approx-0.367$) the negative pressure wave component appears with subsequent positive backward wave caused by the microcirculation (Fig. 10, upper figure); the pressure component in the vessel 2 contains the forward pressure wave passed through the bifuracation and the backward pressure wave caused by the microcirculation (Fig. 10, lower figure).

As a benchmark solution the modeling results for the second order McCormack difference scheme [7] are adopted. For the McCormack scheme at the bifurcation point the continuity equations for the flux and pressure are considered jointly with the requirement that the characteristic variables are constant [2]. The realization of the boundary conditions at the vessel junction for the McCormack scheme differs from the LB approach (25)-(26) in which the characteristics are not considered.

The pressure and velocity profiles measured at the middle of the vessels $1,2$ are depicted in Fig.10 and Fig.11 for LB model and McCormack difference scheme. Both methods give very similar profiles, small differences are observed at some points, they can be attributed to the different approximation methods of the boundary conditions at the vessel bifurcations.

4 Conclusion

In the present paper the method for modeling of the 1D blood flow equations in elastic vessels using kinetic LB approach is proposed. The inclusion of the virtual external force in LB model allows to mimic arbitrary vessel luminal area response to the exerted blood pressure. Several test problems are considered.

In comparison to the conventional numerical methods of 1D blood flow modeling LB approach has the following distinctive features.

1. In LB approach the viscosity of the fluid is the inherent property and its value is linearly dependent on $\tau$. This property can be valuable for the modeling of lengthy vessel networks in which the diffusion effects are non-negligible.
2. The implementation of the boundary conditions for junctions presented in the paper differs from the conventional realization. In the presented method the characteristics are not considered, the conditions for the pressure and the velocity in the junction are written explicitly in terms of LB variables. It is also possible to implement these boundary conditions in a standard way by expressing the characteristics in terms of LB variables (in a similar way like RCR conditions are introduced). The modeling results of LB approach are very similar (not exactly the same) to the pressure and velocity profiles obtained with use of the MacCormack scheme with the application of the conventional boundary conditions at the junctions.
3. The abrupt changes in materials properties are allowed in the presented approach. Note that for some discrete methods the discontinuities in elastic properties are inadmissible [4]. On the other hand the inclusion of the vessel parts in which the pulse wave velocity changes smoothly over a some spatial region is not trivial for the presented method. 4. In LB approach the advection is linear (streaming part is linear) while for the conventional methods the discretization of the nonlinear terms $u\frac{\partial\rho}{\partial x},u\frac{\partial u}{\partial x}$ should be considered. The linearity of the advection part is an attractive property for the parallel multi-CPU implementation of LB method.

The implementation of viscoelastic vessel response was not considered in the present paper. Potentially the viscoelastic effect can be inserted in the virtual force term. This term will contain an approximation of the vessel area time derivative (Kelvin-Voigt model), therefore an additional theoretical question about the stability of the scheme should be considered. This question is leaved for the future study.