A shallow water with variable pressure model for blood flow simulation

We performed numerical simulations of blood flow in arteries with a variable stiffness and cross-section at rest using a finite volume method coupled with a hydrostatic reconstruction of the variables at the interface of each mesh cell. The method was then validated on examples taken from the literature. Asymptotic solutions were computed to highlight the effect of the viscous and viscoelastic source terms. Finally, the blood flow was computed in an artery where the cross-section at rest and the stiffness were varying. In each test case, the hydrostatic reconstruction showed good results where other simpler schemes did not, generating spurious oscillations andnonphysical velocities.


Introduction
In this work we are interested in modeling and simulating blood flow in arteries with varying stiffness and cross-section. The blood flow in the main arteries of the systemic network is governed by the 3D Navier-Stokes equations which can be complicated and time-consuming to solve numerically. Fortunately, using well-known hypothesis valid in the case of blood flow in arteries (long wave approximation D/λ << 1, axial symmetry ∂ θ = 0), this system of equations can be simplified and then integrated over the cross-section of the artery in order to obtain a 1D hyperbolic system of equations, similar to the Saint-Venant system for shallow water flows. Details on the derivation of the model are presented in section (2) and can also be found in [23,38]. Finally, we are left with a set of mass and momentum conservation equations with non dimensionless variables and parameters: with A(x, t) = πR(x, t) 2 the cross-section area (R is the radius of the artery), Q(x, t) = A(x, t)u(x, t) the discharge, u(t, x) the mean flow velocity, ρ the blood density, C f the friction coefficient and A 0 = k √ A 0 with k(x) the stiffness of the artery and A 0 (x) = πR 0 (x) 2 the cross-section at rest.
The vast majority of arteries in the systemic network are tapered, meaning that the cross-section at rest A 0 (x) varies throughout the length of the artery. Similarly, in the presence of arterial pathologies such as aneurysm or stenoses, the stiffness k (x) of the arterial wall can vary locally. As for shallow water equations with topography, the presence of tapper or variable stiffness in an artery modifies the blood flow, and both behaviors are accounted for in (1) through the source term To numerically solve (1), it is necessary, among other things, to discretize this source term. A naive treatment of the topography gradients will most likely generate numerical oscillations, therefore the use of the so-called well-balanced schemes is required to properly balance the fluxes and the source terms. In the following, we will focus on a specific well-balance method, called the hydrostatic reconstruction.
We will first present the derivation of the model and its properties, then the numerical method and in particular the derivation of the well-balanced scheme applied to the case of blood flow in arteries. We will then validate our method on examples taken from the literature and verify asymptotic behaviors of the numerical solution. Finally, we will compute the blood flow in an artery with varying cross-section and stiffness.

Derivation of the 1D blood flow equations
The 1D model for blood flow equations is derived from the conservative form of the Navier-Stokes equations for an incompressible fluid with constant viscosity µ: where u is the velocity vector, ρ the density, supposed constant, p the pressure and τ the stress tensor to be defined. Using the control volume of the Figure  1, we integrate the Navier-Stokes equations over a volume V of cross-section A surrounded by a surface S (V = S ∪ A) and of length dz. We define then the average velocity U and pressure P as From the mass conservation equation (2) we have: We then transform the volume integral using the Green (Divergence) theorem and writing the surface integral as S ∪ A. The surface element is dS = Rdθdz and the two terms are written as We retrieve therefore the first equation of our system For the conservation of momentum equation (3), the temporal term ∂ t ρu becomes and the divergence term In the last two integrals the integration over the surface S is where the term uudS tends to zero. Finally, the integration over the area A gives In terms of the cross-section A and the flow rate Q, we obtained the following system of equations: The viscous effects are contained in f v which is computed by the integration of the shear stress at the wall τ rx over the internal surface dS. Therefore, it depends on the exact flow condition. To close the mathematical problem we need a relation between the pressure P and the cross-section A, P = P (A), called the wall or state law. For f v = C f Q/A and the state law , which corresponds to the elastic response of the artery, we obtain the proposed system of equations. (1).

Conservative hyperbolic system and steady states
Considering an artery with a constant stiffness k and a variable cross-section at rest A 0 (x), (1) reduces to the following system, similar to the shallow water equations with topography: As a reminder, the shallow water system is: with h(x, t) the water height, q(x, t) = h(x, t)u(x, t) the unit discharge, u(x, t) the mean flow velocity, g the constant of gravity, S 0 = −∂ x z the opposite of the slope, z the topography and S f the friction term (which takes the form of Manning's, Stickler's, Chézy's, ... empirical friction law).

Hyperbolic system
The system (5) can be written using the following vectorial form: where U is the vector of the conservative variables, F (U ) is the flux: and S(U ) is the source term, taking into account the shape of the vessel at rest A 0 (x) and the friction term The analogous term for the shallow water equations is the topography source term. The gradient of the flux (8) can be written as the product of the Jacobian matrix J(U ) with the partial derivative of the vector of conservative variables U : When the cross-section A > 0, the Jacobian matrix admits two different real eigenvalues, λ 1 and λ 2 : with c the Moens-Korteweg wave propagation velocity (for the shallow water equations (6), c = √ gh). In this case, the system is said to be strictly hyperbolic, which is a generalization of the advection phenomenon [18,45,30]: a part of the information concerning the flow propagates at the velocity λ 1 and the other part at the velocity λ 2 . For blood flow under physiological conditions, we have λ 1 > 0 and λ 2 < 0, hence the flow is subcritical.

Steady states
Since the works of [3,2] on the shallow water equations, it is well known that if a numerical scheme does not preserve steady states at the discrete level, spurious oscillations and artificial non zero velocities will be generated. The steady states for the system (5) are obtained when considering a stationary flow (i.e. there is no evolution in time) and are governed by the following equations: with b = k/(ρ √ π) constant since we are considering an artery with a constant stiffness k. Neglecting the viscous friction effects (inviscid flow) by setting C f = 0, we obtain the conservation of the discharge and Bernoulli's law for blood flow: In the literature [8,37,44,6], we can find well-balanced numerical methods able to preserve the following steady state: which is the analogous of (13) in the case of the shallow water equations. However, these methods are complicated to handle due to the occurrence of critical points when solving (13) or (14). Therefore we chose to focus on simpler steady states that we call the rest steady states or the "man at eternal rest" equilibrium [13] by analogy with the "lake at rest" (introduced in [1]) or the hydrostatic equilibrium for the shallow water equations: where η is the water level. In this case we have a hydrostatic balance between the hydrostatic pressure and the gravitational acceleration. By analogy, we have the following equilibrium for the blood flow in arteries: Numerical methods able to preserve at least the steady states (16) are said to be "well-balanced" since the work of [19]. A wide panel of well-balanced methods has been developed for shallow water equations. Among others we can mention [29,24,39,27,16,25,1,11,36,17,4,22,5,20]. In [13], we adapted the hydrostatic reconstruction introduced in [1] to the system with constant stiffness (5).
We will now present the hydrostatic reconstruction introduced in [1] adapted to the original system of equations (1) with varying stiffness k(x) and crosssection at rest A 0 (x). By a combination of the mass and momentum equations in (1), under some regularity assumptions, we have: . Considering a stationary flow where the viscous friction is neglected by setting C f = 0, we recover Bernoulli's law (13). The notable difference is that k is now a function of x. In the case of the "man at rest" equilibrium" (without artifacts such as [26,35]) we obtain: The fact that now k is a function of x will influence the way the well-balanced scheme is obtained. In the following section, we will present a well-balanced scheme for system (1), based on the hydrostatic reconstruction for Saint-Venant/shallow water equations with variable pressure [5].

The numerical method 4.1 Numerical context
Several numerical methods have been used to solve the blood flow equations. In [43], they are solved thanks to the Methods of Characteristics (MOC). In [51,50], they use a conservative form of the model with the non-conserved vector (A, u) and equations (19) are solved with a twostep Lax-Wendroff scheme. In [42], a quasi conservative form of the equations (with s(U ) a source term) is solved thanks to a first order explicit in time upwind finite difference scheme.
In [38], they are the first to solve blood flow equations under a conservative form, thanks to a two-step Lax-Wendroff scheme. The solutions of the equations under the form (19) using an upwind Discontinuous Galerkin method (used by [49,48]) and a Taylor Galerkin finite element method (also used in [33,14,34]) have been compared in [41]. A MacCormack finite difference method has been applied in [15] followed by [40]. Finite volume methods seem to be first used to solve these equations in [9,10]. In [13], a well-balanced finite volume method based on the hydrostatic reconstruction (introduced in [1]) is applied on system (5), and this method is compared with a Taylor Galerkin method in [46]. We will present in the following sections the extension of the well-balanced scheme (based on an extension of the hydrostatic reconstruction) we have used to solve the system (1), which can be written under the following vectorial form with and the source terms

Convective step
For the homogeneous system which is (21) without source term, an explicit first order in time conservative scheme can be written as: where i refers to the cell and F i+ 1 2 is an approximation of the flux function F (U, Z) at the cell interface ). This numerical flux will be detailed in subsection 4.4.

Source terms treatment
4.3.1 Topography source term S 1 (U, Z) In the system (21), the term S 1 (U, Z) is involved in the steady state preservation, therefore requires a well-balanced treatment. Following a variant of the hydrostatic reconstruction [5, p.93-94], the variables are reconstructed locally from (18) on both sides of the interface i + 1/2 of the cell C i : . In order to help the understanding of the principle of the hydrostatic reconstruction (26), we present the hydrostatic reconstruction for the shallow water system of equations (6): with z i+1/2 = max(z i , z i+1 ). The water height is reconstructed in a way that allows to have locally the hydrostatic equilibrium h + z = cst on each side of the interface i + 1/2. As mentioned in [1], max(., 0) is there to ensure the positivity of the water height in case of drying and the upwind evaluation of z i+1/2 ensures that 0 ≤ h i+1/2L ≤ h i and 0 ≤ h i+1/2R ≤ h i+1 , which has been proved in [1] to ensures the positivity of the water height. For blood flow equations with a constant stiffness k, the corresponding equilibrium writes A (respectively − √ A 0 ) "plays the role" of h (resp. z), thus in that case the hydrostatic reconstruction writes: As we have − √ A 0 instead of z, we take A 0i+1/2 = min( √ A 0i , A 0i+1 ), thus we have: We can notice that we recover reconstruction (29) if the stiffness k is constant in reconstruction (26). For consistency, the scheme (25) is modified as follows: where . Thus blood flow in a artery with varying cross-section at rest and stiffness is treated in a well-balanced way.

Viscous source term S 2 (U )
In system (21), the friction term −C f Q/A in S 2 (U ) is treated semi-implicitly. This treatment is classical in shallow water simulations [7,31] and has proven efficient in blood flow simulation as well [13]. Furthermore, this treatment preserves the "dead man" equilibrium (18). It consists in using first (30) as a prediction step without friction, i.e.: then applying a semi-implicit friction correction on the predicted values (U * i ): Thus we get the corrected velocity u n+1 i and we have A n+1 i = A * i .

HLL numerical flux
As presented in [13], several numerical fluxes can be used (Rusanov, HLL, VFRoe-ncv and kinetic fluxes) for numerical simulations of blood flow in arteries. Details can be found in [5,12,13]. In this work we will use the HLL flux (Harten Lax and van Leer [21]) because it is the best compromise between accuracy and CPU time consumption (see [12, chapter 2]). It writes: where λ 1 (U, k * ) and λ 2 (U, k * ) are the eigenvalues of the system and k * = max(k L , k R ).
To prevent a blow up of the numerical values, we impose the following CFL (Courant, Friedrichs, Lewy) condition: and n CF L = 1.

Validation of the method
To validate the well-balanced scheme presented in the previous sections for blood flow in arteries with varying stiffness k (x) and cross-section at rest A 0 (x), we applied it to different test cases taken from [13], where arteries with a varying cross-section at rest A 0 (x) and a constant stiffness k were considered. For each of these examples, the rest equilibrium state was: Q = 0 and √ A − √ A 0 = 0 and non-reflecting boundary conditions were set at each end of the computational domain in the form of homogeneous Neumann boundary conditions. The hydrostatic reconstruction scheme as well as a naive centered discretization of the source term were systematically tested to clearly evaluate the benefit of using a well-balanced scheme. According to [13], several Riemann solvers can be used, but we only display results obtained using the HLL flux. In the following, we present the numerical parameters, the analytic solution if it exists and the numerical results. For further details we refer the reader to [13].

"The man at eternal rest"
We considered an artery at its equilibrium state, where there is no flow and the radius of the cross-section at rest R 0 (x) varies throughout the artery, as for example in a dead man with an aneurysm. This equilibrium state is exactly the one well-balanced methods are designed to preserve. If the topography source term is not treated correctly, non-physical velocity may be generated.
We used the following numerical values: L = 0.14 m, J = 50 cells, T end = 5 s, ρ = 1060 kg.m −3 , C f = 0 and k = 4.0 × 10 8 P a.m −1 . We used the equilibrium state as an initial condition, setting Q(x, 0) = 0 and: The results obtained are presented in Figure 2 right. As expected, a naive centered discretization of the topography source term results in nonphysical oscillations of the velocity u (x, t), whereas the well-balanced solution preserves the equilibrium state.

The ideal "Tourniquet"
This test case is the equivalent of the dam break problem for the Shallow Water equations (Stoker's solution in [12]). We considered an artery with a constant radius at rest R 0 , a constant stiffness k and no viscous friction (C f = 0), therefore the governing system of equations was (24). Initially, a tourniquet was applied and then immediately removed. We have a Riemann problem and the method of characteristics allowed us to compute an analytic solution that we compared to the numerical solutions. This Riemann problem has been first introduced in compressible gas dynamic with the Sod tube (for further details we refer the reader to [28,32]) and extended to blood flow in [13].
The results obtained are presented in Figure 3. We can see that the numerical solution obtained with the well balanced scheme is in good agreement with the analytic solution presented in [13]. This is also true for the solution obtained using a centered discretization of the topography source term, which is superposed on the well-balanced solution, since in this case the source term is null.

Wave reflection-transmission of the pulse towards a constriction
In this section we considered the propagation of a pulse towards constriction. This configuration is an idealized representation of a transition between a parent artery and a daughter artery of smaller cross-section. We tested here the ability of the numerical scheme to capture the propagation of a small perturbation of the equilibrium state at the beginning of an artery with a varying radius at rest R 0 (x). In order to accurately compute the numerical solution, the forward and backward traveling waves need to be correctly captured as well as the reflected and transmitted waves generated by the abrupt change in topography at the transition point. To test if these reflections were accurately described, we computed the analytic reflection and transmission coefficients at the transition point and compared them to the amplitude of the numerical reflected waves. For further details we refer the reader to [13].
We considered an artery of length L = 0.16 m and used the following numerical parameters: J = 1500 cells, T end = 8.0 × 10 −3 s, ρ = 1060 kg.m −3 , C f = 0 and k = 1.0 × 10 8 P a.m −1 . The constriction was defined by the following radius of the cross-section at rest: with R R = 4.0 × 10 −3 m, ∆R = 1.0 × 10 −3 m, x 1 = 19 40 L and x 2 = L 2 . We set Q(x, 0) = 0 as an initial condition and we defined the initial perturbation as: with x 3 = 15 100 L < x 1 , x 4 = 35 100 L < x 2 and = 5.0 × 10 −3 a small parameter ensuring that we stayed in the range of small perturbations of the equilibrium state.
The numerical results are plotted in Figure 4. We can see that the propagation of the pulse as well as the wave reflections and transmissions are accurately described using the well balanced scheme (Figure 4 left) whereas spurious waves appear with the centered discretization of the source term (Figure 4 right).

Asymptotic solutions for a uniform vessel
In this section we studied the propagation of a pulse wave in a uniform vessel (k = cst, A 0 = cst) and derived asymptotic solutions of the system of equation (1), following the work of Wang and al. [47]. Small perturbations Q , A 0 + Ã of the base state (Q = 0, A = A 0 ) were considered, resulting in the following linearized system of equations: where c 0 = kR 0 / (2ρ) is the Moens-Korteweg celerity.
In the following numerical examples, we only present results obtained for the hydrostatic reconstruction since we considered a uniform vessel. The numerical parameters were defined as follows: L = 3 m, R 0 = 1.0 × 10 −2 m, J = 1500 cells, T end = 0.5 s, ρ = 1060 kg.m −3 , µ = 3.5 × 10 −3 P a.s and k = 1.0 × 10 7 P a.m −1 . The parameters C f and C v , respectively the viscous coefficient and the viscoelastic coefficient, were set according to the desired test case. Initially, the system was at its equilibrium state Q = 0, A = A 0 = πR 2 0 and an inflow boundary condition was prescribed as Q (x = 0, t) = Q in (t) with: where H(t) is the Heaviside function, T c the period of the sinusoidal wave and Q c the maximum amplitude of the inflow wave. We set Q c = 1.0 × 10 −6 m 3 .s −1 and T c = 0.4 s to insure that only small perturbations from the equilibrium state were considered. The cross-section at the inlet A(x = 0, t) was reconstructed by a matching of the outgoing characteristic, technique that takes advantage of the hyperbolic nature of the problem. A homogeneous Neumann boundary condition was prescribed at the outlet to simplify the computation of the asymptotic solutions and to avoid reflections.

The d'Alembert equation
Following ideas developed in [47], we set C f = 0 in (31) and we obtained the d'Alembert equation, which admits the following pure wave solution c 0Ã0 = Q = Q in (x − c 0 t).
In Figure 5, we can see the propagation of a pulse wave without dissipation or diffusion, as predicted by the analytic solution.

Dissipation due to the viscosity of the blood
We also investigated the effect of the blood viscosity on the propagation of the pulse wave and set C f = 0. Starting from the linearized system of equations (31), we considered the small parameter f = T c C f A0 and performed the change of variables ξ = x − c 0 t and τ = f t to place ourselves in the moving frame at slow times to properly capture the effects of the viscous term. The first order solution obtained in [47] is: t=0.1s t=0.2s t=0.3s t=0.4s t=0.5s where exp − f t 2Tc is the exponential envelop of the pure wave solutionQ 0 (x − c 0 t). To obtain this asymptotic solution numerically, we set C f = 40πν = 4.15 × 10 −4 m 2 .s −1 , therefore f = 0.53.
In Figure 6, we can see the propagation of the pulse with dissipation (or attenuation) of its amplitude due to the viscosity of the blood. The straight doted line represents the exponential envelop exp − f x 2Tcc0 computed previously and is in good agreement with the decrease in amplitude of the pulse wave. One can note that as expected, there is no diffusion, since the wavelength of the pulse does not change while it propagates in the artery.

Diffusion due to the viscoelasticity of the arterial wall
In this section, we set the friction coefficient to zero (C f = 0) and focused on an other important characteristic of the blood flow in the arteries: the viscoelasticity of the arterial wall. We chose here to take into account this time-dependent behavior in our governing system of equations through a very simple lumped model, the Kelvin-Voigt model, resulting in an additional parabolic term in the governing system of equations: where the viscoelastic coefficient C ν is defined as C ν = 2 3 φh ρR0 = 1.57 m 2 .s −1 with φ = 5000 P a.s and h = 5.0 × 10 −3 m. The parabolic term was treated by performing a temporal splitting of the problem. First the purely hyperbolic t=0.1s t=0.2s t=0.3s t=0.4s t=0.5s problem with a non reflecting boundary condition at the outlet was solved, and its solution was then used as an initial condition of the parabolic problem. A Crank-Nicolson scheme coupled with homogeneous Neumann boundary conditions was than used to solve the parabolic problem.
To correctly capture the behavior of this new viscoelastic term, we defined a new small parameter ν = Cv c 2 0 Tc = 8.3 × 10 −2 and applied the same technique as in the previous section. From [47] we have the following first order diffusive analytic solution, which is a solution of the heat equation: The numerical results for several times and the analytic solution at t = 0.4 s are presented in Figure 7. We can see that the viscoelastic term induces a diffusion of the pulse wave, changing its wavelength, and that the numerical solution at t = 0.4 s perfectly matches with the asymptotic solution at t = 0.4 s.

Real artery simulation
In this section, we focused on simulating the propagation of a pulse wave in a tapered artery of length L = 3 m, where the the radius of the cross section at rest R 0 (x) was linearly decreasing from the proximal to the distal end of the artery:  L and x 2 = 16 20 L. Following [47], the stiffness of the arterial wall was defined as k(x) = 4 3 Eh R 2 0 (x) with E the Young's modulus and h the width of the arterial wall. Therefore we were in a configuration where R 0 and k were varying throughout the length of the artery and if the well-balanced scheme was not used, spurious waves might have arisen.
We used the following numerical parameters to mimic the geometrical and mechanical properties of a real artery: J = 1500 cells, T end = 0.5 s, ρ = 1060 kg.m −3 , µ = 3.5 × 10 −3 P a.s, E = 4.0 × 10 5 P a, h = 5.0 × 10 −4 m, C f = 8πν, φ = 5000 P a.s and C v = 2 3 φh ρR0 . We used the same initial inflow condition as for the asymptotic solutions.   t=0.1s t=0.2s t=0.3s t=0.4s t=0.5s The results are presented in figures 8 and 9. We can see that in the absence of friction and viscoelastic effects (figure 8), if the well-balanced scheme is not used (figure 8 left) nonphysical reflections appear. On the contrary, the wellbalanced scheme provides a satisfactory numerical solution, where a continuous reflection phenomena takes place due to the tapering, resulting in a decrease of the amplitude of the backward traveling wave and an increase of the amplitude of the forward traveling wave. Indeed, in the case of a tapered artery, the transmission coefficient T r > 1 and the reflection coefficient R e < 1. When viscous and viscoelastic effects are taken into account (figure 9), all phenomena add up and we recognize the effects of the continuous reflection, the viscous dissipation and the viscoelastic diffusion.

Conclusion and perspectives
In this work we have presented a numerical method based on a well-balanced finite volume scheme for the blood flow equations with variable wall elasticity. This scheme based on an extension of the hydrostatic reconstruction gave very good results on several tests, for which classical methods failed. In further work, we will try to improve the accuracy of the numerical method by raising the order of the numerical method and to apply this method to real network modeling.