A well-balanced finite volume scheme for 1D hemodynamic simulations

with A 0 = k √ A 0 and where A(x, t) is the cross-section area (A = πR 2 with R the radius of the arteria), Q(x, t) = A(x, t)u(x, t) the flow rate or the discharge, u(t, x) the mean flow velocity, ρ the blood density, A 0 (x) the cross section at rest and k(x) the stiffness of the artery. System (1) is into the form of the Saint-Venant problem with variable pressure presented in [3]. We have to mention that arterial pulse wavelengths are long enough to justify the use of a 1D model rather than a 3D model when a global simulation of blood flow in the cardiovascular system is needed.

Since [2,8], it is well known (in the shallow water community) that the scheme should be well-balanced for good source term treatment, i.e. the scheme should preserve at least some steady states. For system (1), we should preserve at least the “man at eternal rest” or “dead man equilibrium” [6] (without artifacts such as [10]), it writes

this means that steady states at rest are preserved (this is the analogous of the “lake at rest” equilibrium). Thus we use the scheme proposed in [3, p.93-94] for that kind of model. This is a finite volume scheme with a modification of the hydrostatic reconstruction (introduced in [1,3] for the shallow water model).

For the homogeneous system

which is (1) with:

an explicit first order conservative scheme writes

where

i refers to the cell

and n to time t n with t n+1 -t n = ∆t. The two points numerical flux

, is an approximation of the flux function F (U, Z) at the cell interface i + 1/2. This numerical flux will be detailled in subsection 1.3.

In system (1), the terms A(∂ x A 0 -2 √ A∂ x k/3)/( √ πρ) are involved in steady states preservation, they need a well-balanced treatment: the variables are reconstructed locally thanks to a variant of the hydrostatic reconstruction [3, p.93-94]

with

). For consistency, the scheme ( 4) is modified as follows

where

and P(A, k) = kA 3/2 /(3ρ √ π). Thus the variation of the radius and the varying elasticity are treated under a well-balanced way. In system (1), the friction term -C f Q/A is treated semi-implicitly. This treatment is classical in shallow water simulations [4,11] and had shown to be efficient in blood flow simulation as well [6]. This treatment does not break the “dead man” equilibrium. It consists in using first (6) as a prediction step without friction, i.e.:

then we apply a semi-implicit friction correction on the predicted values (U * i ):

Thus we get the corrected velocity u n+1 i and we have

As presented in [6], several numerical fluxes might be used (Rusanov, HLL, VFRoe-ncv and kinetic fluxes). In this work we will use the HLL flux (Harten Lax and van Leer [9]) because it is the best compromise between accuracy and CPU time consuming (see [5, 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 blow up of the numerical values, we impose the following CFL (Courant, Friedrichs, Levy) condition

,

2 Some numerical results

In this test, we consider a configuration with no flow and with a change of artery elasticity k(x), this is the case for a dead man with a stented artery (see Figure 1 left). The section of the artery is constant R 0 (x) = 4.0 10 -3 m and the velocity is u(x, t) = 0m/s. We use the following numerical values: J = 50 cells, C f = 0, ρ = 1060m 3 , L = 0.14m, T end = 5s. As initial conditions, we take a fluid at rest The steady state at rest is perfectly preserved in time, we do not notice any spurious oscillation (see Figure 1 right).

We now observe the reflexion and transmission of a pulse through a sudden change of artery elasticity (from k R to k L with k L > k R ) in an elastic tube of constant radius (see Figure 2 left). We take the following numerical values: J = 1500 cells, C f = 0, k L = 1.6 10 8 Pa/m, k R = 1. 10 8 Pa/m, ∆k = 6. 10 7 Pa/m, ρ = 1060m 3 , R 0 = 4.0 10 -3 m, L = 0.16m, T end = 8.0 10 -3 s, c L = k L R/(2ρ) ≃ 17.37m/s and c R = k R R/(2ρ) ≃ 13.74m/s. We take a decreasing elasticity on a rather small scale:

k R else , with x 1 = 19L/40 and x 2 = L/2. As initial conditions, we consider a fluid at rest Q(x, 0) = 0m 3 /s and the following perturbation of radius:

, with ǫ = 1.0 10 -2 . The expression for the pressure is

where p 0 is the external pressure.

The numerical results perfectly match with the predictions for a linearized flow. We get the predicted amplitudes both for the transmitted and the reflected waves (see Figure 2 right).

In this case, the elasticity is constant in space. We consider the viscous term in the linearized momentum equation. A periodic signal is imposed at inflow as a perturbation of a st