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

0.15 MB·English

by  Olivier Delestre

A well-balanced finite volume scheme for 1D
hemodynamic simulations
∗
Olivier Delestre†,†‡and Pierre-Yves Lagr´ee§
2
1 January 11, 2012
0
2

n
a Abstract
J
English version: We are interested in simulating blood flow in arter-
0
ies with variable elasticity with a one dimensional model. We present a
1
well-balanced finite volume scheme based on the recent developments in

] shallowwaterequationscontext. Wethusgetamassconservativescheme
A whichalsopreservesequilibriaofQ=0. Thisnumericalmethodistested
N on analytical tests.
VersionFranc¸aise: Nousnousint´eressons`alasimulationd’´ecoulements
.
h sanguins dans des art`eres dont les parois sont `a ´elasticit´e variable. Ceci
t est mod´elis´e `a l’aide d’un mod`ele unidimensionnel. Nous pr´esentons un
a
m sch´ema ”volume fini ´equilibr´e” bas´e sur les d´eveloppements r´ecents ef-
fectu´espourlar´esolutiondusyst`emedeSaint-Venant. Ainsi,nousobtenons
[
unsch´emaquipr´eservelevolumedefluideainsiqueles´equilibresaurepos:

2 Q=0. Le sch´ema introduit est test´e surdes solutions analytiques.
v
0
1 Introduction
2
6
Weconsiderthefollowingsystemofmassandmomentumconservationwithnon
.
8 dimensionless parametersandvariables,which is the 1D model ofblood flow in
0 anarteryoravesselwithnonuniformelasticity(itisrewritteninaconservative
1
form compared to what we usually find in litterature)
1
:
v ∂ A+∂ Q=0
t x
Xi  ∂ Q+∂ Q2 + 1 kA3/2 = A ∂ 2√A∂ k C Q , (1)
 t x x 0 x f
(cid:20) A 3√πρ (cid:21) √πρ(cid:18) A − 3 (cid:19)− A
r
a

with = k√A and where A(x,t) is the cross-section area (A = πR2 with R
0 0
A
the radiusofthe arteria), Q(x,t)=A(x,t)u(x,t) the flowrateorthe discharge,
u(t,x) the mean flow velocity, ρ the blood density, A (x) the cross section at
0
∗WethankENDOCOMANRforfinancialsupport,wethankFranc¸oisBouchutforfruitful
discussions.
†CNRS& UPMC Universit´eParis 06, UMR 7190, 4 place Jussieu, Institut Jean le Rond
d’Alembert,Boˆıte162,F-75005Paris,France;[email protected]
‡Presentlyat: LaboratoireJ.A.Dieudonn´e&EPUNice-Sophia,Universit´edeNiceSophia
Antipolis,France;[email protected]
§CNRS& UPMC Universit´eParis 06, UMR 7190, 4 place Jussieu, Institut Jean le Rond
d’Alembert,Boˆıte162,F-75005Paris,France;[email protected]
1
restandk(x)thestiffnessoftheartery. System(1)isintotheformoftheSaint-
Venantproblemwithvariablepressurepresentedin[3]. Wehavetomentionthat
arterialpulsewavelengthsarelongenoughtojustifytheuseofa1Dmodelrather
than a 3D model when a global simulation of blood flow in the cardiovascular
system is needed.
1 Numerical method
Since [2, 8], it is well known(in the shallow water community) that the scheme
shouldbe well-balancedforgoodsourcetermtreatment, i.e. the schemeshould
preserveatleastsomesteadystates. Forsystem(1),weshouldpreserveatleast
the ”manateternalrest”or”deadmanequilibrium”[6](withoutartifactssuch
as [10]), it writes
u=0
 1 A , (2)
A3/2k =Cst
 √πρ − √πρA0
this means that steadystates at rest are preserved(this is the analogousof the
”lakeatrest”equilibrium). Thusweuse theschemeproposedin[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).
1.1 Convective step
For the homogeneous system
∂ U +∂ F(U,Z)=0, (3)
t x
which is (1) with:
A Q
U = , Z = A0 andF(U,k)= ,
(cid:18) Q (cid:19) (cid:18) k (cid:19) (cid:18) Q2/A+kA3/2/(3√πρ) (cid:19)
an explicit first order conservative scheme writes
Un+1 Un Fn Fn
i − i + i+1/2− i−1/2 =0, (4)
∆t ∆x
where Un is an approximation of U
i
1 xi+1/2
Un U(x,t )dx.
i ≃ ∆xZ n
xi−1/2
i refers to the cell C = (x ,x ) = (x ,x +∆x) and n to time
i i−1/2 i+1/2 i−1/2 i−1/2
t with t t =∆t.
n n+1 n
−
The two points numerical flux
Fn = (Un,Un ,k∗ ),
i+1/2 F i i+1 i+1/2
with k∗ = max(k ,k ), is an approximation of the flux function F(U,Z)
i+1/2 i i+1
at the cell interface i+1/2. This numerical flux will be detailled in subsection
1.3.
2
1.2 Source terms treatment
In system (1), the terms A(∂ 2√A∂ k/3)/(√πρ) are involved in steady
x 0 x
A −
states preservation, they need a well-balanced treatment: the variables are re-
constructed locally thanks to a variant of the hydrostatic reconstruction [3,
p.93-94]
A =max(k √A +min(∆ ,0),0)/k∗
i+1/2L i i A0i+1/2 i+1/2
 pUi+1/2L =(Ai+1/2L,Ai+1/2L.ui)t , (5)
 Ai+1/2R =max(ki+1 Ai+1−max(∆A0i+1/2,0),0)/ki∗+1/2
Up =(A ,Ap .u )t
i+1/2R i+1/2R i+1/2R i+1

with∆ = =k A k √A andk∗ =max(k ,k ).
A0i+1/2 A0i+1−A0i i+1 0i+1− i 0i i+1/2 i i+1
For consistency, the scheme (4) is mopdified as follows
∆t
Un+1 =Un Fn Fn , (6)
i i − ∆x(cid:16) i+1/2L− i+1/2R(cid:17)
where
Fn =Fn +S
i+1/2L i+1/2 i+1/2L ,
Fn =Fn +S
i−1/2R i−1/2 i−1/2R
with
Fn = U ,U ,k∗
i+1/2 F(cid:16) i+1/2L i+1/2R i+1/2(cid:17)
0
S =
i+1/2L (cid:18) P(Ani,ki)−P(Ani+1/2L,ki∗+1/2) (cid:19)
0
S =
i−1/2R (cid:18) P(Ani,ki)−P(Ani−1/2R,ki∗−1/2) (cid:19)
and (A,k)=kA3/2/(3ρ√π). Thus the variationofthe radiusandthe varying
P
elasticityaretreatedunderawell-balancedway. Insystem(1),thefrictionterm
C Q/Ais treatedsemi-implicitly. This treatmentis classicalin shallowwater
f
−
simulations[4,11]andhadshowntobeefficientinbloodflowsimulationaswell
[6]. This treatment does not break the ”dead man” equilibrium. It consists in
using first (6) as a prediction step without friction, i.e.:
∆t
U∗ =Un Fn Fn ,
i i − ∆x(cid:16) i+1/2L− i−1/2R(cid:17)
then we apply a semi-implicit friction correction on the predicted values (U∗):
i
un+1 u∗
A∗ i − i = C un+1.
i (cid:18) ∆t (cid:19) − f i
Thus we get the corrected velocity un+1 and we have An+1 =A∗.
i i i
1.3 HLL numerical flux
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
3
1.7e+08
the elasticity k 1e-10
1.6e+08 velocity, t=0s
velocity, t=5s
1.5e+08
5e-11
1.4e+08
m]
k [Pa/ 11..23ee++0088 u [m/s] 0
1.1e+08
-5e-11
1e+08
9e+07
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 -1e-10
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
x [m]
x [m]
Figure1: Nospuriousflows(right)aregeneratedbyachangeofelasticity(left).
Lax and van Leer [9]) because it is the best compromise between accuracy and
CPU time consuming (see [5, chapter 2]). It writes:
F(U ,k∗) if0 c
L 1
F(UL,UR,k∗)= c2F(UL,kc∗2)−cc11F(UR,k∗) + c2c1c2c1(UR−UL) ifc1≤<0<c2 ,
F(U ,k∗) − − ifc 0
R 2
≤

with
c = inf ( inf λ (U,k∗))andc = sup ( sup λ (U,k∗)),
1 j 2 j
U=UL,UR j∈{1,2} U=UL,UR j∈{1,2}
where λ (U,k∗) and λ (U,k∗) are the eigenvalues of the system and k∗ =
1 2
max(k ,k ).
L R
To prevent blow up of the numerical values, we impose the following CFL
(Courant, Friedrichs, Levy) condition
∆x
∆t n ,
CFL
≤ max(u +c )
i i
i | |
where c = k √A /(2ρ√π) and n =1.
i i i CFL
p
2 Some numerical results
2.1 ”The stented man at eternal rest”
Inthistest,weconsideraconfigurationwithnoflowandwithachangeofartery
elasticity k(x), this is the case fora deadman witha stented artery(see Figure
1left). ThesectionofthearteryisconstantR (x)=4.010−3mandthevelocity
0
is u(x,t)=0m/s. We use the following numerical values: J =50 cells, C =0,
f
ρ=1060m3, L=0.14m, T =5s. As initialconditions,wetakeafluidatrest
end
4
Q(x,0)=0m3/s and
k ifx [0:x ] [x :L]
0 1 4
∈ ∪
 ∆k x x1
k + sin − π π/2 +1 ifx ]x :x [
k(x)= k00+∆2k (cid:18) (cid:18)x2−x1 − (cid:19) (cid:19) ifx∈[x12 :x23] ,
∈
∆k x x
3
k + cos − π +1 ifx ]x :x [
 0 2 (cid:18) (cid:18)x4−x3 (cid:19) (cid:19) ∈ 3 4
with k =1.0108Pa/m, ∆k =6.0107Pa/m, x =1.010−2m, x =3.0510−2m,
0 1 2
x =4.9510−2m and x =7.010−2m.
3 4
The steady state at rest is perfectly preservedin time, we do notnotice any
spurious oscillation (see Figure 1 right).
2.2 Wave reflection-transmission in a stented artery
We now observe the reflexion and transmission of a pulse through a sudden
change of artery elasticity (from k to k with k > k ) in an elastic tube of
R L L R
constantradius(seeFigure2left). Wetakethefollowingnumericalvalues: J =
1500cells,C =0,k =1.6108Pa/m,k =1.108Pa/m,∆k =6.107Pa/m,ρ=
f L R
1060m3, R = 4.0 10−3m, L = 0.16m, T = 8.0 10−3s, c = k R/(2ρ)
0 end L L
≃
17.37m/s and c = k R/(2ρ) 13.74m/s. We take a decreapsing elasticity
R R
≃
on a rather small scaple:
k +∆k ifx [0:x ]
R 1
∈
 ∆k x x1
k(x)= k + 1+cos − π ifx ]x :x ] ,
 R 2 (cid:20) (cid:18)x2 x1 (cid:19)(cid:21) ∈ 1 2
−
k else
R

with x = 19L/40 and x = L/2. As initial conditions, we consider a fluid at
1 2
rest Q(x,0)=0m3/s and the following perturbation of radius:
100 65L
R (x) 1+ǫsin π x ifx [65L/100:85L/100]
R(x,0)= 0 (cid:20) (cid:18)20L (cid:18) − 100(cid:19)(cid:19)(cid:21) ∈ ,
 R (x) else
0

with ǫ=1.010−2. The expression for the pressure is
p(x,t)=p (x)+k(x)(R(x,t) R (x)),
0 0
−
where p is the external pressure.
0
The numerical results perfectly match with the predictions for a linearized
flow. Wegetthepredictedamplitudesbothforthetransmittedandthereflected
waves (see Figure 2 right).
2.3 Wave ”damping”
In this case,the elasticity is constantin space. We consider the viscousterm in
the linearized momentum equation. A periodic signal is imposed at inflow as a
perturbation of a steady state (R = Cst, u = 0) with a constant section at
0 0
rest. WetakeR=R +ǫR andu=0+ǫu ,where(R ,u )istheperturbationof
0 1 1 1 1
5
1.7e+08
the elasticity k 2000
11..56ee++0088 1500 t=t=3TTeenntdd=//044
k [Pa/m] 111...234eee+++000888 p-p [Pa]0 1 500000 (Tpr-(pp0-p(t0=(0t=))0/)2)/2
1.1e+08
1e+08 0
Re(p-p0(t=0))/2
9e+07
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 -500
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16
x [m]
x [m]
Figure 2: Across a discontinuity of k (left), an initial pulse evolves (right) ; it
is transmitted and reflected.
the steadystate. Lookingforprogressivewaves(i.e. underthe formei(ωt−Kx)),
we take for the imposed incoming discharge
Q (t)=Q(t,x=0)=Q sin(ωt)m3/s.
b amp
Thus, we have a damping wave in the domain
0 ifk x>ωt
Q(t,x)= r ,
(cid:26) Qampsin(ωt krx)ekix ifkrx ωt
− ≤
withω =2π/T , T thetime lengthofapulseandK =k +ik thewave
pulse pulse r i
vector.
We use the following numerical values: J =750 cells, ρ=1060m3, L=3m,
k = 1.108Pa/m, R = 4.10−3m and T = 5s. We consider both C = 0 and
0 end f
C = 0.005053. As initial conditions, we take a fluid at rest Q(x,0) = 0m3/s
f
and as input boundary condition
Q (t)=Q sin(ωt),
b amp
with ω = 2π/T = 2π/0.5s and Q = 3.45 10−7m3/s. The output is an
pulse amp
outgoing wave.
The results are closed to the analytical solution (see Figure 3). We notice a
small numerical diffusion for C =0.005053.
f
3 Conclusion
Inthiswork,wehaveconsideredthe1Dmodelofflowinanarterywithvarying
elasticityandconstantsection. Wehavepresentedawell-balancedfinitevolume
scheme. Thus we get a mass conservative scheme. Moreover,the well-balanced
property allows to have a good treatment of the source, i.e. we do not get
numerical artifacts. This numerical method gave good results on numerical
tests. In future works, we will have to add some extra source terms in order
to get a more realistic model. These extra terms will require to develop a low
diffusive high order scheme in the spirit of [7]. Moreover, this will improve the
accuracy of the scheme. And we will also have to test more complex cases such
as bifurcations and networks.
6
0
Exact, Cf=0 2e-08
HLL1-SI
-5e-08
0
-1e-07 -2e-08
3Q [m/s] -1-.25ee--0077 3Q [m/s] --64ee--0088
-8e-08
-2.5e-07
-1e-07
-3e-07
-3.5e-07 -1.2e-07 Exact, Cf=0H.0L0L510-5S3I
0 0.5 1 1.5 2 2.5 3 -1.4e-07
0 0.5 1 1.5 2 2.5 3
x [m]
x [m]
Figure 3: Response to an oscillating input condition, left with no viscosity and
right with viscosity.
References
[1] E. Audusse, F. Bouchut, M.-O. Bristeau, R. Klein, and B. Perthame. A
fast and stable well-balanced scheme with hydrostatic reconstruction for
shallow water flows. SIAM J. Sci. Comput., 25(6):2050–2065,2004.
[2] Alfredo Bermudez and M. Elena Vazquez. Upwind methods for hyperbolic
conservation laws with source terms. Computers & Fluids, 23(8):1049 –
1071,1994.
[3] F. Bouchut. Nonlinear stability of finite volume methods for hyperbolic
conservation laws, and well-balanced schemes for sources, volume 2/2004.
Birkh¨auser Basel, 2004.
[4] M.-O. Bristeau and Benoˆıt Coussin. Boundary conditions for the shallow
waterequationssolvedby kinetic schemes. TechnicalReport4282,INRIA,
October 2001.
[5] OlivierDelestre.Simulationduruissellementd’eaudepluiesurdessurfaces
agricoles/ rain water overland flow on agricultural fields simulation. PhD
thesis, Universit´e d’Orl´eans (in French), available from TEL: tel.archives-
ouvertes.fr/INSMI/tel-00531377/fr,July 2010.
[6] Olivier Delestre and Pierre-Yves Lagr´ee. A ”well balanced” finite volume
scheme for blood flow simulation. submitted.
[7] Olivier Delestre and Fabien Marche. A numerical scheme for a viscous
shallow water model with friction. J. Sci. Comput., DOI 10.1007/s10915-
010-9393-y,2010.
[8] J. M. Greenberg and A.-Y. LeRoux. A well-balanced scheme for the nu-
merical processing of source terms in hyperbolic equation. SIAM Journal
on Numerical Analysis, 33:1–16,1996.
[9] Amiram Harten, Peter D. Lax, and Bram van Leer. On upstream differ-
encing and godunov-type schemes for hyperbolic conservation laws. SIAM
Review, 25(1):35–61,January 1983.
7
[10] Robert Kirkman, Tony Moore, and Charlie Adlard. The Walking Dead.
Image Comics, 2003.
[11] Qiuhua Liang and Fabien Marche. Numerical resolution of well-balanced
shallow water equations with complex source terms. Advances in Water
Resources, 32(6):873 – 884, 2009.
8