Low-Shapiro hydrostatic reconstruction technique for blood flow simulation in large arteries with varying geometrical and mechanical properties

TL;DR

This paper introduces a low-Shapiro hydrostatic reconstruction technique for simulating blood flow in large arteries, offering a balance of simplicity, efficiency, and accuracy in capturing wave dynamics with varying arterial properties.

Contribution

The paper develops and compares two hydrostatic reconstruction methods tailored for blood flow, highlighting HR-LS's efficiency and accuracy in large arterial network simulations.

Findings

HR-S is the most accurate but computationally intensive.

HR-LS effectively captures low-Shapiro steady states and wave dynamics.

HR is inadequate for large arteries with property variations.

Abstract

The purpose of this work is to construct a simple, efficient and accurate well-balanced numerical scheme for one-dimensional (1D) blood flow in large arteries with varying geometrical and mechanical properties. As the steady states at rest are not relevant for blood flow, we construct two well-balanced hydrostatic reconstruction techniques designed to preserve low-Shapiro number steady states that may occur in large network simulations. The Shapiro number S h = u/c is the equivalent of the Froude number for shallow water equations and the Mach number for compressible Euler equations. The first is the low-Shapiro hydrostatic reconstruction (HR-LS), which is a simple and efficient method, inspired from the hydrostatic reconstruction technique (HR). The second is the subsonic hydrostatic reconstruction (HR-S), adapted here to blood flow and designed to exactly preserve all subcritical steady states. We systematically compare HR, HR-LS and HR-S in a series of single artery and arterial network numerical tests designed to evaluate their well-balanced and wave-capturing properties. The results indicate that HR is not adapted to compute blood flow in large arteries as it is unable to capture wave reflections and transmissions when large variations of the arteries' geometrical and mechanical properties are considered. On the contrary, HR-S is exactly well-balanced and is the most accurate hydrostatic reconstruction technique. However, HR-LS is able to compute low-Shapiro number steady states as well as wave reflections and transmissions with satisfying accuracy and is simpler and computationally less expensive than HR-S. We therefore recommend using HR-LS for 1D blood flow simulations in large arterial network simulations.

Figures

Tables

Table 1. Table 1: Parameters describing the artery used in the different test-cases, given in "cgs": the density ρ 𝜌 \rho , the length L 𝐿 L , the inlet radius R i ​ n subscript 𝑅 𝑖 𝑛 R_{in} and the inlet rigidity K i ​ n subscript 𝐾 𝑖 𝑛 K_{in} .

$\rho$ [$g.c{m}^{-3}$]	$L$ [$cm$]	$R_{in}$ [$cm$]	$K_{in}$ [$g.c{m}^{-2}.s^{-2}$]
1	10	0.5	$1\times {10}^5$

Table 2. Table 2: Values of the inlet Shapiro number S h , i ​ n subscript 𝑆 ℎ 𝑖 𝑛 S_{h,in} and the wall deformation parameter Δ ​ 𝒢 Δ 𝒢 \Delta\mathcal{G} used in the single artery test-cases. These values are chosen to test the well-balanced methods in the limits of the low-Shapiro number flow regime.

$S_{h,in}$	$0$	$1\times {10}^{-3}$	$1\times {10}^{-2}$	$1\times {10}^{-1}$

Table 3. Table 3: Steady solutions: Relative errors L 1 ​ [ Q ] superscript 𝐿 1 delimited-[] 𝑄 L^{1}\left[Q\right] and L 1 ​ [ E ] superscript 𝐿 1 delimited-[] 𝐸 L^{1}\left[E\right] computed in the stenosis ( 83 ) and the step ( 84 ) for N = 50 𝑁 50 N=50 cells for all combinations of values of the inlet Shapiro number S h , i ​ n subscript 𝑆 ℎ 𝑖 𝑛 S_{h,in} and the wall deformation parameter Δ ​ 𝒢 Δ 𝒢 \Delta\mathcal{G} taken for table 2 . Only HR-S is exactly well-balanced, but HR-LS is more accurate than HR.

		Stenosis			Step
$S_{h,in}$		0
$\Delta 𝒢$		$1\%$	$10\%$	$30\%$	$1\%$	$10\%$	$30\%$
	HR	$0$	$0$	$0$	$0$	$0$	$0$
$L^1\left[Q\right]$	HR-LS	$0$	$0$	$0$	$0$	$0$	$0$
	HR-S	$0$	$0$	$0$	$0$	$0$	$0$
	HR	$0$	$0$	$0$	$0$	$0$	$0$
$L^1\left[E\right]$	HR-LS	$0$	$0$	$0$	$0$	$0$	$0$
	HR-S	$0$	$0$	$0$	$0$	$0$	$0$
$S_{h,in}$		$1\times {10}^{-3}$
$\Delta 𝒢$		$1\%$	$10\%$	$30\%$	$1\%$	$10\%$	$30\%$
	HR	$4.0\times {10}^{-4}$	$4.2\times {10}^{-3}$	$1.4\times {10}^{-2}$	$2.2\times {10}^{-4}$	$2.3\times {10}^{-3}$	$7.4\times {10}^{-3}$
$L^1\left[Q\right]$	HR-LS	$3.6\times {10}^{-7}$	$4.1\times {10}^{-6}$	$1.9\times {10}^{-5}$	$1.8\times {10}^{-7}$	$2.1\times {10}^{-6}$	$9.4\times {10}^{-6}$
	HR-S	$5.4\times {10}^{-13}$	$5.3\times {10}^{-13}$	$4.2\times {10}^{-14}$	$2.9\times {10}^{-13}$	$3.1\times {10}^{-13}$	$5.8\times {10}^{-13}$
	HR	$3.0\times {10}^{-4}$	$5.1\times {10}^{-3}$	$4.2\times {10}^{-2}$	$2.1\times {10}^{-4}$	$9.4\times {10}^{-3}$	$1.3\times {10}^{-1}$
$L^1\left[E\right]$	HR-LS	$2.1\times {10}^{-7}$	$2.6\times {10}^{-6}$	$1.5\times {10}^{-5}$	$1.1\times {10}^{-7}$	$1.4\times {10}^{-6}$	$1.0\times {10}^{-5}$
	HR-S	$4.6\times {10}^{-13}$	$4.9\times {10}^{-13}$	$6.1\times {10}^{-13}$	$6.7\times {10}^{-13}$	$6.5\times {10}^{-13}$	$1.4\times {10}^{-12}$
$S_{h,in}$		$1\times {10}^{-2}$
$\Delta 𝒢$		$1\%$	$10\%$	$30\%$	$1\%$	$10\%$	$30\%$
	HR	$4.0\times {10}^{-4}$	$4.2\times {10}^{-3}$	$1.4\times {10}^{-2}$	$2.3\times {10}^{-4}$	$2.3\times {10}^{-3}$	$7.4\times {10}^{-3}$
$L^1\left[Q\right]$	HR-LS	$3.6\times {10}^{-6}$	$4.1\times {10}^{-5}$	$1.9\times {10}^{-4}$	$1.8\times {10}^{-6}$	$2.1\times {10}^{-5}$	$9.4\times {10}^{-5}$
	HR-S	$2.6\times {10}^{-13}$	$2.7\times {10}^{-13}$	$9.6\times {10}^{-14}$	$2.4\times {10}^{-13}$	$9.5\times {10}^{-14}$	$1.8\times {10}^{-13}$
	HR	$3.0\times {10}^{-4}$	$5.1\times {10}^{-3}$	$4.2\times {10}^{-2}$	$2.1\times {10}^{-4}$	$9.4\times {10}^{-3}$	$1.2\times {10}^{-1}$
$L^1\left[E\right]$	HR-LS	$2.1\times {10}^{-6}$	$2.6\times {10}^{-5}$	$1.5\times {10}^{-4}$	$1.1\times {10}^{-6}$	$1.4\times {10}^{-5}$	$8.1\times {10}^{-5}$
	HR-S	$2.7\times {10}^{-13}$	$2.7\times {10}^{-13}$	$3.3\times {10}^{-13}$	$2.7\times {10}^{-13}$	$3.4\times {10}^{-13}$	$5.9\times {10}^{-13}$
$S_{h,in}$		$1\times {10}^{-1}$
$\Delta 𝒢$		$1\%$	$10\%$	$30\%$	$1\%$	$10\%$	$30\%$
	HR	$4.0\times {10}^{-4}$	$4.2\times {10}^{-3}$	$1.4\times {10}^{-2}$	$2.3\times {10}^{-4}$	$2.3\times {10}^{-3}$	$7.5\times {10}^{-3}$
$L^1\left[Q\right]$	HR-LS	$3.6\times {10}^{-5}$	$4.1\times {10}^{-4}$	$1.8\times {10}^{-3}$	$1.8\times {10}^{-5}$	$2.1\times {10}^{-4}$	$9.0\times {10}^{-4}$
	HR-S	$2.6\times {10}^{-13}$	$3.4\times {10}^{-13}$	$2.0\times {10}^{-13}$	$2.8\times {10}^{-13}$	$2.4\times {10}^{-13}$	$1.4\times {10}^{-13}$
	HR	$3.2\times {10}^{-4}$	$5.4\times {10}^{-3}$	$4.4\times {10}^{-2}$	$2.2\times {10}^{-4}$	$9.9\times {10}^{-3}$	$1.2\times {10}^{-1}$
$L^1\left[E\right]$	HR-LS	$2.2\times {10}^{-5}$	$2.8\times {10}^{-4}$	$1.8\times {10}^{-3}$	$1.2\times {10}^{-5}$	$2.0\times {10}^{-4}$	$2.2\times {10}^{-3}$
	HR-S	$2.3\times {10}^{-13}$	$2.4\times {10}^{-13}$	$2.9\times {10}^{-13}$	$2.3\times {10}^{-13}$	$2.9\times {10}^{-13}$	$3.4\times {10}^{-13}$

Table 4. Table 4: Steady solutions: Relative errors L 1 ​ [ Q ] superscript 𝐿 1 delimited-[] 𝑄 L^{1}\left[Q\right] and L 1 ​ [ E ] superscript 𝐿 1 delimited-[] 𝐸 L^{1}\left[E\right] computed in the stenosis ( 83 ) and the step ( 84 ) for S h , i ​ n = 1 × 10 − 2 subscript 𝑆 ℎ 𝑖 𝑛 1 superscript 10 2 S_{h,in}=1\times 10^{-2} and Δ ​ 𝒢 = 10 % Δ 𝒢 percent 10 \Delta\mathcal{G}=10\% obtained for N ∈ { 50 , 100 , 200 , 400 } 𝑁 50 100 200 400 N\in\left\{50,100,200,400\right\} . HR and HR-LS converge with order 1 whereas HR-S is exactly well-balanced up to machine precision.

	Stenosis				Step
	HR
$N$	$L^1\left[Q\right]$	Order	$L^1\left[E\right]$	Order	$L^1\left[Q\right]$	Order	$L^1\left[E\right]$	Order
$50$	$4.22\times {10}^{-3}$	$-$	$5.09\times {10}^{-3}$	$-$	$2.34\times {10}^{-3}$	$-$	$9.41\times {10}^{-3}$	$-$
$100$	$2.11\times {10}^{-3}$	$-1.01$	$2.56\times {10}^{-3}$	$-1.01$	$1.17\times {10}^{-3}$	$-1.01$	$8.64\times {10}^{-3}$	$-0.12$
$200$	$1.05\times {10}^{-3}$	$-1.01$	$1.28\times {10}^{-3}$	$-1.01$	$5.86\times {10}^{-4}$	$-1.01$	$8.26\times {10}^{-3}$	$-0.07$
$400$	$5.26\times {10}^{-4}$	$-1.00$	$6.38\times {10}^{-4}$	$-1.00$	$2.93\times {10}^{-4}$	$-1.00$	$8.07\times {10}^{-3}$	$-0.03$
	HR-LS
$N$	$L^1\left[Q\right]$	Order	$L^1\left[E\right]$	Order	$L^1\left[Q\right]$	Order	$L^1\left[E\right]$	Order
$50$	$4.14\times {10}^{-5}$	$-$	$2.61\times {10}^{-5}$	$-$	$2.08\times {10}^{-5}$	$-$	$1.39\times {10}^{-5}$	$-$
$100$	$2.07\times {10}^{-5}$	$-1.01$	$1.31\times {10}^{-5}$	$-1.01$	$1.04\times {10}^{-5}$	$-1.01$	$7.24\times {10}^{-6}$	$-0.96$
$200$	$1.04\times {10}^{-5}$	$-1.01$	$6.58\times {10}^{-6}$	$-1.00$	$5.19\times {10}^{-6}$	$-1.01$	$3.91\times {10}^{-6}$	$-0.90$
$400$	$5.19\times {10}^{-6}$	$-1.00$	$3.30\times {10}^{-6}$	$-1.00$	$2.59\times {10}^{-6}$	$-1.00$	$2.24\times {10}^{-6}$	$-0.80$
	HR-S
$N$	$L^1\left[Q\right]$	Order	$L^1\left[E\right]$	Order	$L^1\left[Q\right]$	Order	$L^1\left[E\right]$	Order
$50$	$2.68\times {10}^{-13}$	$-$	$2.73\times {10}^{-13}$	$-$	$9.53\times {10}^{-14}$	$-$	$3.43\times {10}^{-13}$	$-$
$100$	$1.40\times {10}^{-15}$	$-$	$3.39\times {10}^{-13}$	$-$	$9.20\times {10}^{-14}$	$-$	$3.97\times {10}^{-13}$	$-$
$200$	$1.94\times {10}^{-12}$	$-$	$7.30\times {10}^{-13}$	$-$	$2.26\times {10}^{-12}$	$-$	$8.44\times {10}^{-13}$	$-$
$400$	$8.83\times {10}^{-12}$	$-$	$1.45\times {10}^{-12}$	$-$	$1.01\times {10}^{-11}$	$-$	$1.63\times {10}^{-12}$	$-$

Table 5. Table 5: Wave propagation: Relative error L 1 ​ [ Q ] superscript 𝐿 1 delimited-[] 𝑄 L^{1}\left[Q\right] computed in the step ( 84 ) for values of S h , i ​ n subscript 𝑆 ℎ 𝑖 𝑛 S_{h,in} and Δ ​ 𝒢 Δ 𝒢 \Delta\mathcal{G} taken from table 2 obtained for N = 1600 𝑁 1600 N=1600 . HR, HR-LS and HR-S present similar results except for Δ ​ 𝒢 = 30 % Δ 𝒢 percent 30 \Delta\mathcal{G}=30\% .

$S_{h,in}$		$1\times {10}^{-3}$
$\Delta 𝒢$		$1\%$	$10\%$	$30\%$
	HR	$2.3\times {10}^{-2}$	$5.5\times {10}^{-2}$	$5.5\times {10}^{-1}$
$L^1\left[Q\right]$	HR-LS	$2.3\times {10}^{-2}$	$2.8\times {10}^{-2}$	$6.6\times {10}^{-2}$
	HR-S	$2.3\times {10}^{-2}$	$2.8\times {10}^{-2}$	$6.6\times {10}^{-2}$
$S_{h,in}$		$1\times {10}^{-2}$
$\Delta 𝒢$		$1\%$	$10\%$	$30\%$
	HR	$2.3\times {10}^{-2}$	$5.5\times {10}^{-2}$	$5.5\times {10}^{-1}$
$L^1\left[Q\right]$	HR-LS	$2.3\times {10}^{-2}$	$2.8\times {10}^{-2}$	$6.6\times {10}^{-2}$
	HR-S	$2.3\times {10}^{-2}$	$2.8\times {10}^{-2}$	$6.6\times {10}^{-2}$
$S_{h,in}$		$1\times {10}^{-1}$
$\Delta 𝒢$		$1\%$	$10\%$	$30\%$
	HR	$2.9\times {10}^{-2}$	$6.1\times {10}^{-2}$	$5.1\times {10}^{-1}$
$L^1\left[Q\right]$	HR-LS	$2.9\times {10}^{-2}$	$3.5\times {10}^{-2}$	$7.6\times {10}^{-2}$
	HR-S	$2.9\times {10}^{-2}$	$3.5\times {10}^{-2}$	$7.5\times {10}^{-2}$

Equations

$$\{\begin{array}{rl}\frac{\partial A}{\partial t}& +\frac{\partial Q}{\partial x}=0\\ \frac{\partial Q}{\partial t}& +\frac{\partial }{\partial x}\left(\frac{Q^2}{A}\right)=-\frac{A}{\rho }\frac{\partial P}{\partial x}.\end{array}$$

$$P\left(x,t\right)=P_0+K\left(x\right)\left(\sqrt{A\left(x,t\right)}-\sqrt{A_0\left(x\right)}\right),$$

$$\{\begin{array}{rl}& \frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=0\\ & \frac{\partial Q}{\partial t}+\frac{\partial F}{\partial x}=S_T,\end{array}$$

$$F=\frac{Q^2}{A}+\frac{K}{3\rho }{A}^{\frac{3}{2}},$$

$$S_T=\frac{A}{\rho }\left(\frac{\partial }{\partial x}\left(K\sqrt{A_0}\right)-\frac{2}{3}\sqrt{A}\frac{\partial K}{\partial x}\right).$$

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial }{\partial x}\left[\mathbf{F}\left(\mathbf{U},K\right)\right]=\mathbf{S}\left(\mathbf{U},K\right)\frac{\partial \mathbf{\sigma }}{\partial x}.$$

$$\mathbf{U}=\left[\begin{array}{c}A\\ Q\end{array}\right],\mathbf{F}\left(\mathbf{U},K\right)=\left[\begin{array}{c}Q\\ F\end{array}\right],$$

$$\mathbf{\sigma }=\left[\begin{array}{c}K\\ Z\end{array}\right]=\left[\begin{array}{c}K\\ K\sqrt{A_0}\end{array}\right],\mathbf{S}\left(\mathbf{U}\right)=\left[\begin{array}{cc}0& 0\\ -\frac{2}{3}\frac{A^{\frac{3}{2}}}{\rho }& \frac{A}{\rho }\end{array}\right].$$

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial }{\partial x}\left[\mathbf{F}\left(\mathbf{U},K\right)\right]=0.$$

$$\mathbf{J}\left(\mathbf{U},K\right)=\frac{\partial \mathbf{F}}{\partial \mathbf{U}}=\left[\begin{array}{cc}0& 1\\ \frac{K\sqrt{A}}{2\rho }-\frac{Q^2}{A^2}& \frac{2Q}{A}\end{array}\right].$$

$${\lambda }_1=\frac{Q}{A}-c,{\lambda }_2=\frac{Q}{A}+c,{\mathbf{R}}_1=\left[\begin{array}{c}1\\ {\lambda }_1\end{array}\right],{\mathbf{R}}_2=\left[\begin{array}{c}1\\ {\lambda }_2\end{array}\right].$$

$$c=\sqrt{\frac{K}{2\rho }\sqrt{A}}.$$

$$S_h=\frac{u}{c}=\frac{1}{c}\frac{Q}{A}.$$

$$\{\begin{array}{rl}{\mathbb{U}}_{sub}=& \left\{\frac{Q}{A}\in \mathbb{R}|A>0,K>0,Z>0,S_h<1\right\}\\ {\mathbb{U}}_{sup}=& \left\{\frac{Q}{A}\in \mathbb{R}|A>0,K>0,Z>0,S_h>1\right\}.\end{array}$$

$$\frac{\partial \mathbf{W}}{\partial t}+\mathbf{\Delta }\frac{\partial \mathbf{W}}{\partial x}=0.$$

$$\mathbf{W}=\left[\begin{array}{c}{W}_1\\ W_2\end{array}\right]=\left[\begin{array}{c}\frac{Q}{A}-4c\\ \\ \frac{Q}{A}+4c\end{array}\right].$$

$$\{\begin{array}{rl}& A={\left(\frac{2\rho }{K}\right)}^2{\left(\frac{W_2-W_1}{8}\right)}^4\\ & Q=A\frac{W_1+W_2}{2}.\end{array}$$

$$\frac{\partial \eta }{\partial t}+\frac{\partial G}{\partial x}\le 0,$$

$$\{\begin{array}{rl}& \eta \left(\mathbf{U},K\right)=\frac{Q^2}{2A}+\frac{2}{3}\frac{K}{\rho }{A}^{\frac{3}{2}}\\ & G\left(\mathbf{U},K\right)=\left(\frac{Q^2}{2A}+\frac{K}{\rho }{A}^{\frac{3}{2}}\right)\frac{Q}{A}.\end{array}$$

$$\{\begin{array}{rl}& \tilde{\eta }\left(\mathbf{U},\mathbf{\sigma }\right)=\eta \left(\mathbf{U},K\right)-\frac{Z}{\rho }A\\ & \tilde{G}\left(\mathbf{U},\mathbf{\sigma }\right)=G\left(\mathbf{U},K\right)-\frac{Z}{\rho }Q.\end{array}$$

$$\{\begin{array}{rl}& Q=C_1\\ & \frac{1}{2}\frac{Q^2}{A^2}+\frac{1}{\rho }\left(K\sqrt{A}-Z\right)=C_2,\end{array}$$

$$\{\begin{array}{rl}& Q=0\\ & K\sqrt{A}-Z=C_2.\end{array}$$

$$C_i=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]=\left[x_i-\frac{\Delta x}{2},x_i+\frac{\Delta x}{2}\right],i\in \left[1,N\right],$$

$$t^n=n\Delta t,n\in \mathbb{N}.$$

$$\begin{array}{rl}& {\int }_{C_i}\left[\mathbf{U}\left(x,t^{n+1}\right)-\mathbf{U}\left(x,t^n\right)\right]\mathrm{d}x+\\ & {\int }_{t^n}^{t^{n+1}}\left[\mathbf{F}\left(\mathbf{U}\left(x_{i+\frac{1}{2}},t\right),K\left(x_{i+\frac{1}{2}}\right)\right)-\mathbf{F}\left(\mathbf{U}\left(x_{i-\frac{1}{2}},t\right),K\left(x_{i-\frac{1}{2}}\right)\right)\right]\mathrm{d}t=0.\end{array}$$

$${\mathbf{U}}_{\mathbf{i}}^{\mathbf{n}}\approx \frac{1}{\Delta x}{\int }_{C_i}\mathbf{U}\left(x,t^n\right)\text{d}x,$$

$${\mathbf{F}}_{\mathbf{i}+\frac{1}{2}}^{\mathbf{n}}\approx \frac{1}{\Delta t}{\int }_{t^n}^{t^{n+1}}\mathbf{F}\left(\mathbf{U}\left(x_{i+\frac{1}{2}},t\right),K\left(x_{i+\frac{1}{2}}\right)\right)\mathrm{d}t.$$

$${\mathbf{U}}_{\mathbf{i}}^{\mathbf{n}+1}={\mathbf{U}}_{\mathbf{i}}^{\mathbf{n}}-\frac{\Delta t}{\Delta x}\left[{\mathbf{F}}_{\mathbf{i}+\frac{1}{2}}^{\mathbf{n}}-{\mathbf{F}}_{\mathbf{i}-\frac{1}{2}}^{\mathbf{n}}\right].$$

$${\mathbf{F}}_{\mathbf{i}+\frac{1}{2}}^{\mathbf{n}}=\mathcal{F}\left({\mathbf{U}}_{\mathbf{L}},{\mathbf{U}}_{\mathbf{R}}\right)=\left[\begin{array}{c}{\mathcal{F}}_{\mathcal{A}}\left({\mathbf{U}}_{\mathbf{L}},{\mathbf{U}}_{\mathbf{R}}\right)\\ {\mathcal{F}}_{\mathcal{Q}}\left({\mathbf{U}}_{\mathbf{L}},{\mathbf{U}}_{\mathbf{R}}\right)\end{array}\right].$$

$${\mathbf{U}}_{\mathbf{L}}={\mathbf{U}}_{\mathbf{i}}^{\mathbf{n}},{\mathbf{U}}_{\mathbf{R}}={\mathbf{U}}_{\mathbf{i}+1}^{\mathbf{n}}.$$

Full text

1. 1 Introduction

2. 2 Mathematical model

3. 2.1 Model derivation

4. 2.2 Hyperbolic system

5. 3 Numerical scheme for the homogeneous conservative system

6. 3.1 Finite volume numerical scheme

7. 3.2 Kinetic numerical flux

8. 3.3 Initial condition

9. 3.4 Subcritical boundary condition

10. 3.4.1 Inlet boundary condition: imposed flow rate $Q_{in}$

11. 3.4.2 Unsteady outlet boundary condition: reflection of the outgoing characteristic

12. 3.4.3 Steady outlet boundary condition: imposed cross-sectional area $A_{out}$

13. 3.4.4 Junction boundary condition: conservation of mass and continuity of total pressure

14. 4 Hydrostatic reconstruction

15. 4.1 The hydrostatic reconstruction: HR

16. 4.2 The low-Shapiro hydrostatic reconstruction: HR-LS

17. 4.3 The subsonic hydrostatic reconstruction: HR-S

18. 4.3.1 Well-balanced subsonic positivity-preserving reconstruction procedure for the cross-sectional area $A$

19. 4.3.2 Well-balanced subsonic first-order numerical scheme

20. 5 Physiological examples in a single artery

21. 5.1 Steady solutions

22. 5.2 Single wave propagation

23. 5.2.1 The step configuration

24. 5.2.2 The stenosis configuration

25. 6 A realistic example: stenosis of the iliac artery in a 55 arteries network

26. 7 Conclusion

27. 8 Acknowledgments

1.   Introduction

Since the early work of Euler [24], one-dimensional (1D) models have been successfully used to describe the flow of blood in the large arteries of the systemic network [47, 3, 57, 42, 10]. They have proved to be valuable and efficient tools to capture pulse wave propagation in large network simulations and obtain satisfactory evaluations of average quantities such as the cross-sectional area ($A$), the flow rate ($Q$) or the pressure ($P$) [2, 49]. In recent works, 1D models have also been used to compute inverse problems to obtain patient specific parameters [34, 38, 23]. Due to their simplicity, efficiency and the reduced number of parameters they require, we hope that in the near future these 1D models will be intensively used by medical practitioners to make pre-and-post operative diagnosis and perform patient specific simulations of surgeries.

In physiological situations, the mechanical and geometrical properties of the arterial wall can vary locally. These variations can be caused by tapering (figure 1 left), pathologies such as stenoses (figure 1 center) or aneurysms (figure 1 right) and endovascular prosthesis (stent). Mathematically, they result in a source term in the momentum conservation equation that prevents from writing the system in a conservation-law form. A naive discretization of this nonconservative source term can lead to spurious oscillations of the numerical solution and the failure of the numerical method, especially close to steady states [21]. This problem was originally pointed out by Roe [51] for the scalar equation with source terms and reflects a truncation error between the discretization of the conservative flux gradient and the nonconservative source term that does not vanish close to steady states. Since the works of Bermúdez and Vázquez [9] and LeRoux [28, 29] in the context of shallow-water equations, numerical schemes that preserve some steady states at a discrete level are called well-balanced.

The aim of this study is to propose a simple, robust and efficient well-balanced numerical method for blood flow in an artery with variations of its mechanical and geometrical properties. As blood flow equations are mathematically similar to shallow water equations, several well-balanced numerical schemes have been derived for 1D blood flow equations with varying geometrical and mechanical properties. A popular approach consists in expressing the system in terms of primitive variables, namely the cross-sectional area ($A$) and the flow velocity $u$. The resulting system can be written in a conservation-law form, even in the presence of varying geometrical and mechanical properties. However, it has been proved for shallow water equations that this formulation is not mass-conservative and can lead to erroneous estimations of the wave celerity [55]. This analysis is also valid for blood flow equations and the numerical solutions obtained with a nonconservative system will be incorrect in the presence of elastic jumps. Indeed, the Rankine-Hugoniot jump relation of the nonconservative form is different from the one of the conservative form. C̆anić [15] and Sherwin [53] were among the first to address the issue of the nonconservative source term for blood flow simulation. C̆anić proposed to treat the nonconservative product in this source term through jump conditions, while Sherwin used a two-rarefaction Riemann solver when the material properties varied abruptly. More recently, Toro and Siviglia [56] reformulated the 1D conservative system with varying geometrical and mechanical properties as a homogeneous quasi-linear system and solved the associated Riemann problem. To do so, they introduced an auxiliary steady variable containing the geometrical and mechanical properties of the artery, and also included variations of the external pressure. In the framework of path-conservative methods [45], Müller and Toro [41] used this augmented quasi-linear system to propose an exactly well-balanced numerical scheme for all steady states (subcritical, transcritical and supercritical). Murillo and García-Navarro [43] derived an energy balanced numerical scheme in the framework of augmented solvers for arteries with varying mechanical and geometrical properties, and also variations of the external pressure. In [21], Delestre and Lagrée successfully applied the hydrostatic reconstruction (HR), proposed in [6] for shallow water equations, to compute blood flow in arteries with varying cross-sectional area. In more recent work [20], Delestre extended the hydrostatic reconstruction (HR) to arteries with varying cross-sectional area and arterial wall rigidity.

The hydrostatic reconstruction (HR) meets the simplicity and efficiency requirements for 1D blood flow simulation and will be the reference well-balanced method used in this study. The hydrostatic reconstruction (HR) can be used with any finite-volume numerical flux for a conservative problem and guarantees the following natural properties of shallow water flows:

• •

well-balanced for the steady states at rest, or hydrostatic equilibria;

• •

the conservation of mass;

• •

the non-negativity of the water-height $h$;

• •

the ability to compute dry states and transcritical flows;

• •

a discrete or semi-discrete entropy inequality, which enables to compute the entropic solution in the presence of a discontinuity.

Unfortunately, the steady states at rest preserved by the hydrostatic reconstruction (HR) are not relevant for blood flow as they only occur in "dead men" [21]. We propose two extensions of the hydrostatic reconstruction adapted to blood flow simulation in large arteries.

By relaxing some of the properties of the hydrostatic reconstruction (HR) such as the ability to compute dry states, we derive an extension of the hydrostatic reconstruction, that we refer to as the "low-Shapiro" hydrostatic reconstruction (HR-LS). The low-Shapiro hydrostatic reconstruction (HR-LS) accurately preserves low-Shapiro number steady states that may occur in large network simulations. The Shapiro number $S=u/c$ is the equivalent of the Froude number for shallow water equations and the Mach number for compressible Euler equations. We also adapt the subsonic hydrostatic reconstruction (HR-S), proposed by Bouchut [13], to blood flow equations with variable geometrical and mechanical properties. The subsonic hydrostatic reconstruction (HR-S) exactly preserves all subcritical steady states, including low-Shapiro number steady states. By construction, both the low-Shapiro hydrostatic reconstruction (HR-LS) and the subsonic hydrostatic reconstruction (HR-S) are able to accurately compute wave reflections and transmissions. The different numerical methods are then tested and compared in a series of steady and unsteady physiological flow configurations, where both the geometrical and mechanical wall properties vary.

This work is organized as follows. In section 2 we derive the hyperbolic system of equations that describes the flow of blood in large arteries and recall its main mathematical properties. In section 3, we present a kinetic numerical scheme for the homogeneous problem and the boundary conditions used in the examples presented in this study. In section 4, we introduce the low-Shapiro hydrostatic reconstruction (HR-LS) and the subsonic hydrostatic reconstruction (HR-S) for blood flow in arteries with varying mechanical and geometrical wall properties. In sections 5 and 6, we present a series a steady and unsteady test cases for a single artery and a 55 arteries network, in which we evaluate the performances of the different hydrostatic reconstruction techniques.

2.   Mathematical model

2.1.   Model derivation

The 1D models for blood flow are derived by averaging over the cross-sectional area of each artery a simplified Navier-Stokes system of equations. These simplified equations are obtained using the long wave approximation (${D/\lambda}\ll 1$, ratio between the averaged diameter of an artery $D$ and the average wave length of the pulse wave $\lambda$) and supposing the axial symmetry of blood flow ($\partial_{\theta}=0$). We recall that in arteries the ratio ${D/\lambda}$ is of the order of $1\times 10^{-2}$, therefore the long wave hypothesis is asymptotically valid. Because blood and wall viscosities will damp the effects we want to highlight, namely pulse wave propagation, we neglect them in the rest of this study. We use the inviscid system of equations describing the flow of blood in an elastic artery at the axial position $x$ and time $t$

[TABLE]

The variables $Q$, $A$ and $P$ are respectively the flow rate, the cross-sectional area and the blood pressure. We also introduce the flow velocity $u=\frac{Q}{A}$. The parameter $\rho$ is the density of blood and is supposed constant. For a description of the derivation of system (1) we refer the reader to [35, 8, 31]. To close the system of equations, the variation of pressure is linked to the deformation of the artery. Assuming that the arterial wall is a homogeneous, incompressible Hookean solid and that the artery is represented by a thin-cylinder whose sections move independently of one another, the following wall law is obtained, describing the elastic, or spring-like, behavior of the arterial wall

[TABLE]

where $A_{0}$ is the cross-sectional area at rest of the artery and $K$ is the arterial wall rigidity. Both quantities can vary with the axial position $x$. More complex and general pressure laws can be used (for example in veins [46]), yet equation (2) contains sufficient information to describe the main features of blood flow in large arteries [57, 49]. Combining both system (1) and equation (2) we obtain the final 1D nonconservative system of equations

[TABLE]

where $F$ is the momentum flux

[TABLE]

and $S_{T}$ is a source term taking into account the possible variations of the geometrical and mechanical properties of the arterial wall

[TABLE]

2.2.   Hyperbolic system

System (3) can be written as a system of balance laws

[TABLE]

$\mathbf{U}$ and $\mathbf{F}$ are respectively the vector of conservative variables and the vector of mass and momentum flux

[TABLE]

and the vector $\bm{\sigma}$ and the matrix $\mathbf{S}$ are defined as:

[TABLE]

The main difficulty of system (6) lies in the presence of the nonconservative source term $\mathbf{S}\frac{\partial\bm{\sigma}}{\partial x}$. This nonconservative term vanishes when the cross-sectional area at rest $A_{0}$ and the arterial wall rigidity $K$ are constant, and system (6) is reduced to the following system of conservation laws

[TABLE]

The conservative system (9) has been thoroughly studied by many authors and we only briefly recall its properties. Additional details can be found in [25]. To analyze the mathematical properties of the system (9), we compute the Jacobian matrix of the flux vector $\mathbf{F}$

[TABLE]

$\mathbf{J}\left(\mathbf{U,K}\right)$ has two real eigenvalues $\lambda_{1}$ and $\lambda_{2}$, respectively associated to two right eigenvectors $\mathbf{R_{1}}$ and $\mathbf{R_{2}}$

[TABLE]

The variable $c$ is the Moens-Korteweg wave speed [39, 33] and corresponds to the speed of pulse waves in an artery

[TABLE]

The hyperbolicity of the system is characterized by the Shapiro number $S_{h}$, introduced by Shapiro in [52]

[TABLE]

$S_{h}$ as the analogue of the Froude number $F_{r}$ for the shallow-water equations or of the Mach number $M_{a}$ for compressible flows. Depending on the value of $S_{h}$, we distinguish two flow regimes, represented respectively by the subcritical velocity domain $\mathbb{U}_{sub}$ and the supercritical velocity domain $\mathbb{U}_{sup}$

[TABLE]

In both regions $\mathbb{U}_{sub}$ and $\mathbb{U}_{sup}$, system (9) is strictly hyperbolic as $\lambda_{1}\neq\lambda_{2}$ and the right eigenvectors $\mathbf{R_{1}}$ and $\mathbf{R_{2}}$ are linearly independent. However, when $S_{h}=1$ the flow is critical and the system looses its strict hyperbolicity. In this case resonance phenomena can occur, leading to a possible loss of uniqueness of the solution [37, 32, 36, 30].

In physiological conditions, blood flow is almost always subcritical. Nevertheless, very specific pathologies may lead to supercritical flows but will not be the subject of this study. Only subcritical solutions of system (9) and more generally of system (6) in $\mathbb{U}_{sub}$ will be considered here.

For solutions of system (9) in $\mathbb{U}_{sub}$, linear algebra shows that the Jacobian matrix $\mathbf{J}$ is diagonalizable in the form $\mathbf{J}=\mathbf{R}\mathbf{\Delta}\mathbf{R}^{-1}$, where the columns of $\mathbf{R}$ are the right eigenvectors $\mathbf{R_{1}}$ and $\mathbf{R_{2}}$ and $\mathbf{\Delta}$ is a diagonal matrix containing the eigenvalues of $\mathbf{J}$. Introducing a new vector $\mathbf{W}=\left[W_{1},W_{2}\right]^{T}$ such that $\partial_{\mathbf{U}}\mathbf{W}=\mathbf{R}^{-1}$, system (9) can be written as:

[TABLE]

Finally, by integrating the equation $\partial_{\mathbf{U}}\mathbf{W}=\mathbf{R}^{-1}$, the following expression for $\mathbf{W}$ is obtained

[TABLE]

The vector $\mathbf{W}$ is often referred to as the Riemann invariant vector and is linked to the conservative variables

[TABLE]

The relations (17) are useful to define the boundary conditions at the inlet and outlet of the computational domain.

The vector $\mathbf{U}$, solution of system (9), satisfies an entropy inequality linked to the entropy pair $\left(\eta,G\right)$

[TABLE]

where $\eta$ is the entropy and $G$ is the entropy flux

[TABLE]

This entropy inequality is extended to solutions of system (6) through a new entropy pair $\left(\tilde{\eta},\tilde{G}\right)$ taking into account the vector $\bm{\sigma}$

[TABLE]

This entropy inequality is closely linked to the variation of the physical energy of the system. The existence of such an inequality is essential in order to select the correct physical solution across discontinuities [27].

System (6) admits non-trivial steady solutions, verifying the following steady state system of equations

[TABLE]

where $C_{1}$ and $C_{2}$ are two constants. In the following, we note $E=\frac{1}{2}\frac{Q^{2}}{A^{2}}+\frac{1}{\rho}\left(K\sqrt{A}-Z\right)$ the energy discharge. A particular family of steady states are the steady states at rest, or "man at eternal rest" equilibria, defined by

[TABLE]

For shallow water flows, steady states mainly occur in lakes and verify the "man at eternal rest" equilibria (22). In arteries, steady or quasi-steady flow regimes are observed in small segments when the frequency of the pulse wave is greatly reduced due to a high resistance of the flow, for example after severe stenoses or in smaller arteries. In these cases, the relevant equilibria are no longer the steady states at rest but the non-zero flow steady states described by system (21).

3.   Numerical scheme for the homogeneous conservative system

In this section we describe the finite volume numerical scheme used to solve the homogeneous conservative system (9). The spatial domain is discretized in a series of cells $C_{i}$ defined as

[TABLE]

where $\Delta x$ is the cell size, supposed constant for simplicity. The time domain is also discretized using a constant time step $\Delta t$ and the discrete times are defined as

[TABLE]

3.1.   Finite volume numerical scheme

We first derive the integral form of the conservative system (9) by integrating it with respect to $t$ and $x$ over $\left]t^{n},\>t^{n+1}\right[\times C_{i}$ [36]

[TABLE]

We then approximate the integrals in (25) using the discrete variable $\mathbf{U_{i}^{n}}$ and the numerical flux $\mathbf{F^{n}_{i+\frac{1}{2}}}$, corresponding respectively to an approximation of the space average of the exact solution $\mathbf{U}$ over the cell $C_{i}$ at time $t^{n}$

[TABLE]

and to an approximation of the time average of $\mathbf{F}$ at the cell interface $I_{i+\frac{1}{2}}$

[TABLE]

Using these definitions, we obtain the following explicit finite volume numerical scheme

[TABLE]

We define $\mathbf{F^{n}_{i+\frac{1}{2}}}$ as a two-points numerical flux vector, namely

[TABLE]

As we focus only on first-order finite volume numerical schemes, the vectors $\mathbf{U_{L}}$ and $\mathbf{U_{R}}$ at the cell interface $I_{i+\frac{1}{2}}$ at time $t^{n}$ are defined as

[TABLE]

The choice of the function $\mathcal{F}$ defines the numerical flux and thus the finite volume scheme. Several possibilities exist, and a review of the most common ones applied to blood flow equations can be found in [21, 41, 40, 57, 43, 5].

3.2.   Kinetic numerical flux

We choose to compute the function $\mathcal{F}$ using a kinetic numerical flux, and a review of this method applied to different systems of equations can be found in [11]. The kinetic method was first introduced for shallow water equations in [48], combined with the hydrostatic reconstruction (HR) in [7] and adapted to the blood flow in [21, 5]. The principal motivations for choosing a kinetic numerical flux are that it preserves the positivity of the cross-sectional area and its numerical diffusion is better suited to compute resonant solutions [12, 4]. We briefly recall the classical kinetic approach.

Following [48, 7], we introduce the real, positive, even and compactly supported function $\chi\left(w\right)$, verifying the following properties

[TABLE]

We choose the following expression for the function $\chi\left(w\right)$

[TABLE]

Using this function, we define the kinetic Maxwellian, or so-called Gibbs equilibrium, which represents the density of microscopic particles moving at the velocity $\xi\in\mathbb{R}$

[TABLE]

where

[TABLE]

Noticing that the integral and the first and second moments on $\mathbb{R}$ of $M$ respectively allow to recover $A$, $Q$ and $F$, it can be proved [48] that $\mathbf{U}$ is solution of system (9) if and only if $M$ satisfies the following linear kinetic equation

[TABLE]

where $\mathcal{Q}\left(x,t,\xi\right)$ is a collision term that satisfies

[TABLE]

As the equation (35) is linear, it can be approximated by a simple upwind scheme. The flux function $\mathcal{F}$ is then obtained using the integral and the first moment of the upwind numerical flux used to solve the linear kinetic equation (35), and writes

[TABLE]

where $\mathbf{U_{L}}$ and $\mathbf{U_{R}}$ are defined as in (30). The fluxes $\mathcal{F}^{+}\left(\mathbf{U}\right)$ and $\mathcal{F}^{-}\left(\mathbf{U}\right)$ are defined as

[TABLE]

After some computation, we find that

[TABLE]

with

[TABLE]

The stability of the scheme is ensured if at each time $t^{n}$, the time step $\Delta t$ verifies the following CFL (Courant, Friedrichs and Lewy) [19] condition

[TABLE]

3.3.   Initial condition

All numerical simulations presented in this study are initialized by the following solution of the steady state at rest system (22)

[TABLE]

and the initial vector of conservative variable in the cell $C_{i}$ is then

[TABLE]

3.4.   Subcritical boundary condition

In each artery at time $t^{n}$, boundary conditions are imposed in inlet and outlet ghost cells, respectively noted $C_{in}$ and $C_{out}$, by setting the value of their associated vector of conservative variable $\mathbf{U_{in}^{n}}$ and $\mathbf{U_{out}^{n}}$. As we compute subcritical solutions of system (6) in $\mathbb{U}_{sub}$, one boundary condition is imposed in the inlet ghost cell $C_{in}$ and one boundary condition is imposed in the outlet ghost cell $C_{out}$, respectively allowing to determine one component of $\mathbf{U_{in}^{n}}$ and one component of $\mathbf{U_{out}^{n}}$. To compute the remaining unknown components of $\mathbf{U_{in}^{n}}$ and $\mathbf{U_{out}^{n}}$, we follow the methodology proposed by Bristeau and Coussin [14] and Alastruey [1]. In the following, we assume that in each cell $C_{i}$ at time $t^{n}$, the discrete vector of conservative variables $\mathbf{U_{i}^{n}}$ is known.

3.4.1.   Inlet boundary condition: imposed flow rate $Q_{in}$

We describe here a methodology to impose the flow rate $Q_{in}\left(t^{n}\right)=Q_{in}^{n}$ at the interface between the first cell of the computational domain $C_{1}$ and the inlet ghost cell $C_{in}$, namely

[TABLE]

Taking advantage of the fact that the kinetic flux function $\mathcal{F}$ can be split in two, equation (44) can be expressed as

[TABLE]

To ensure the stability of the scheme, this condition is imposed in an upwind manner. Following [14], we define the quantity

[TABLE]

Two possible cases exist:

• •

If $a_{1}\leq 0$, the dominant part of the information is coming from inside the computational domain. As we are performing an upwind evaluation of the inlet boundary condition, we impose

[TABLE]

• •

If $a_{1}>0$, the dominant part of the information is coming from outside the computational domain. In this case, we impose

[TABLE]

An additional equation is required to completely determine $\mathbf{U_{in}^{n}}$. We take advantage of the characteristic structure of the problem: as we are in the subcritical case, there exists an outgoing characteristic on which the Riemann invariant $W_{1}$ is constant. Using this property, we assume that a correct estimation of the cell average value of the outgoing Riemann invariant $W_{1}\left(\mathbf{U_{in}^{n}}\right)$ is $W_{1}\left(\mathbf{U_{1}^{n}}\right)$. Finally, we impose

[TABLE]

$\mathbf{U_{in}^{n}}$ is obtained by solving either system (47) or system (48). This can be done using a classic Newton’s method in a limited number of iterations ($\sim 5$).

3.4.2.   Unsteady outlet boundary condition: reflection of the outgoing characteristic

We propose here a methodology to characterize the incoming information at the outlet of the computational domain. Indeed, as we are in the subcritical regime, there exists an outgoing characteristic on which the Riemann invariant $W_{2}$ is constant and an incoming characteristic on which propagates the Riemann invariant $W_{1}$. As in the previous case, the cell average value of the outgoing Riemann invariant $W_{2}\left(\mathbf{U_{out}^{n}}\right)$ can be estimated by $W_{2}\left(\mathbf{U_{N}^{n}}\right)$ and we impose

[TABLE]

The value of the incoming Riemann invariant $W_{1}\left(\mathbf{U_{out}^{n}}\right)$ is unknown as it propagates on a characteristic coming from outside the computational domain. In large artery simulations, it is common to estimate the incoming Riemann invariant $W_{1}\left(\mathbf{U_{out}^{n}}\right)$ as a fraction of the outgoing Riemann invariant $W_{2}\left(\mathbf{U_{out}^{n}}\right)$ [57, 1, 2, 43]. This fraction is quantified by a reflection coefficient $R_{t}$ such that

[TABLE]

where $W_{1}\left(\mathbf{U_{out}^{0}}\right)$ and $W_{2}\left(\mathbf{U_{out}^{0}}\right)$ are the initial Riemann invariants of the ghost cell $C_{out}$. The reflection coefficient $R_{t}$, whose value ranges between [math] and $1$, models the reflective and resistive behavior of the network that is not taken into account in the numerical simulation and lies distal (anatomically located far from the point of reference) to the outlet of the computational domain. Finally, using the relations (17), we solve the following system of equations to obtain $\mathbf{U_{out}^{n}}$

[TABLE]

When we wish to remove any incoming information, or equivalently any distal reflection, we set $R_{t}=0$.

3.4.3.   Steady outlet boundary condition: imposed cross-sectional area $A_{out}$

We describe here a methodology to impose the cross-sectional area $A_{out}$ at the outlet of the computational domain. Indeed, when the flow rate is imposed at the inlet of the computational domain and the vector $\bm{\sigma}$ is known, the steady states verifying system (21) are completely determined if we impose the value of the cross-sectional area at the outlet of the computational domain. We set in the outlet ghost cell $C_{out}$ a constant cross-sectional area $A_{out}$

[TABLE]

We need only to compute $Q_{out}^{n}$ to completely determine the outlet vector of conservative variables $\mathbf{U_{out}^{n}}$. To do so, we estimate as in the previous section $W_{2}\left(\mathbf{U_{out}^{n}}\right)$ by

[TABLE]

and using the relations (17), we compute $W_{1}\left(\mathbf{U_{out}^{n}}\right)$ and then $Q_{out}^{n}$

[TABLE]

3.4.4.   Junction boundary condition: conservation of mass and continuity of total pressure

In large network simulations, boundary conditions must also be provided at every junction point, where arteries that are neither the inlet segment nor the terminal segments are connected together. At these junctions points, the outlets of the parent arteries are linked to the inlets of the daughter arteries through the conservation of mass and the continuity of total pressure, or equivalently the energy discharge [54]. We consider here a junction point where a single parent artery $A_{P}$ is connected to $N_{D}$ daughter arteries $\left(A_{Di}\right)_{i=1}^{N_{D}}$. The values of $\mathbf{U_{out}^{n}}|_{A_{P}}$ and of $\left(\mathbf{U_{in}^{n}}|_{A_{Di}}\right)_{i=1}^{N_{D}}$ must be computed, with a total of $2\left(N_{D}+1\right)$ unknowns. $N_{D}+1$ equations are obtained by estimating the outgoing Riemann invariant of the parent and daughter arteries as

[TABLE]

The missing equations are provided by the conservation of mass and total pressure, or equivalently the energy discharge, at each junction point [50, 53]

[TABLE]

In practice, since in physiological conditions the flow is always subcritical, we can simplify the problem and impose only the continuity of pressure [2], neglecting the advection terms in the second equation of system (54)

[TABLE]

This set of equations allows to accurately compute wave reflections and transmissions if a change of impedance occurs between the parent and the daughter arteries. Accurately computing reflected waves is crucial to obtain physiological wave forms in the simulated network, as the observed pressure and flow waves are the sum of the incoming and the multiple reflected waves [2, 3, 49].

4.   Hydrostatic reconstruction

In many physiological configurations, the geometrical and mechanical properties of an artery vary significantly with its length. In the scope of this paper, these geometrical and mechanical gradients are limited to variations of the cross-sectional area at rest $A_{0}$ and the arterial wall rigidity $K$. To prevent spurious oscillations of the numerical solution of system (6) close to steady states, a well-balanced numerical scheme is required to properly balance the source term $S_{T}$ and the flux gradient $\partial_{x}\mathbf{F}$.

In order to make an explicit analogy with the well-balanced methods derived for shallow water equations, we introduce the following notations

[TABLE]

With these notations, we have $\frac{\partial\mathcal{E}}{\partial A}=\frac{\mathcal{P}}{A^{2}}$ and the flux vector $\mathbf{F}$ can be expressed as

[TABLE]

Moreover, the steady state systems (21) and (22) can respectively be written as

[TABLE]

and

[TABLE]

4.1.   The hydrostatic reconstruction: HR

The hydrostatic reconstruction (HR) was introduced by Audusse [6] for shallow water equations and applied to blood flow equations by Delestre [21, 20]. Through a reconstruction of the conservative variables, HR allows to obtain a simple and efficient well-balanced numerical scheme given any finite volume numerical flux for the homogeneous problem (9). It is simple to implement and can easily be adapted to different pressure laws with multiple varying parameters, which is useful when considering veins, collapsible tubes and external pressure variations [46, 18, 42]. This technique allows to preserve at a discrete level the steady states at rest (59) and guarantees that the scheme verifies some natural properties of the shallow water equations (listed as bullets in the introduction), such as the positivity of the water height (equivalent of the cross-sectional area $A$), the ability to compute dry states and transcritical flows and a semi-discrete entropy inequality. This last property is necessary to select the admissible entropy solution across a discontinuity, as explained in [27].

On both sides of each cell interface $I_{i+\frac{1}{2}}$, reconstructed conservative variables are defined to preserve the following system of equations, which coincides with the steady states at rest system (59) when the flow rate $Q$ or the velocity $u$ are zero

[TABLE]

Details on the derivation of HR for blood flow in an artery with variable cross-sectional area $A$ and variable arterial wall rigidity $K$ can be found in [20].

In large arteries, the steady states at rest preserved by HR only occur for "dead men" or distal to an obliterated segment and are of little interest when simulating blood flow in the systemic network. However, in regions of large flow resistance such as small arteries, arterioles or arteries severely constricted by a stenosis, the flow looses its pulsatility and reaches steady or near-steady states with a non-zero flow rate. These quasi-steady flow configurations can occur in large network simulations when the level of arterial precision extends to small arteries and arterioles or in the presence of a very severe stenosis. They are described by the steady state system (58). Therefore, a modification of HR is necessary to capture the relevant steady states for blood flow in large arteries, described by system (58).

4.2.   The low-Shapiro hydrostatic reconstruction: HR-LS

System (58) is nonlinear and difficult to solve in practice. However, in physiological conditions, blood flow is subcritical with a Shapiro number of the order of $S_{h}\approx 1\times 10^{-2}$. Therefore, the nonlinear advection term $\frac{1}{2}\frac{Q^{2}}{A^{2}}$ in system (58) can be neglected at first order with respect to the term $\mathcal{E}(A,K)+\frac{\mathcal{P}(A,K)}{A}-\frac{Z}{\rho}$ that scales as $c^{2}$. Doing so, we obtain the following simplified low-Shapiro number steady state system of equations

[TABLE]

System (61) coincides with the steady state at rest system (59) when $Q$ or $u$ are zero and is an asymptotically correct approximation of the steady state system (58) in low-Shapiro number flow regimes. It also contains the correct conservation properties to obtain low-Shapiro number wave reflections if a change of impedance occurs at the interface between two cells of the computational domain. Indeed, the conservation properties of system (61) are identical to those of system (55), which have proved to be adequate to compute wave reflections and transmissions at junction points [2, 57]. System (61) is the basis for the derivation of the modification of HR we propose in this study, referred to as the low-Shapiro hydrostatic reconstruction (HR-LS) and better suited to compute blood flow in physiological conditions.

HR-LS aims at preserving low-Shapiro number steady states system (61) in an artery with a varying cross-sectional area at rest $A_{0}$ and arterial wall rigidity $K$. Similarly to HR, the well-balanced property is enforced by defining reconstructed variables on both sides of each cell interface $I_{i+\frac{1}{2}}$ according to the reconstruction procedure (61). In the following, variables noted with $"^{*}"$ will refer to the reconstructed variables. Given the vectors of conservative variables $\mathbf{U_{L}}$ and $\mathbf{U_{R}}$ and the vectors $\bm{\sigma_{L}}$ and $\bm{\sigma_{R}}$ at the left and right of the interface $I_{i+\frac{1}{2}}$ between cells $C_{i}$ and $C_{i+1}$, the discrete analogue of system (61) writes

[TABLE]

By solving system (62) and preserving the positivity of $H$, we obtain the following reconstructed variables

[TABLE]

The reconstructed variable $Z^{*}$ is chosen considering nonlinear stability arguments that require that

[TABLE]

to preserve the positivity of $H$. A simple choice is the downwind value

[TABLE]

In order to obtain the reconstructed values $A_{L}^{*}$ and $A_{R}^{*}$, we must select a reconstruction for $K^{*}$. Following [12, 20] we choose

[TABLE]

Therefore, we directly have

[TABLE]

Finally, at each side of the interface $I_{i+\frac{1}{2}}$, we obtain the reconstructed conservative vectors

[TABLE]

that will be used to compute the numerical flux $\mathcal{F}\left(\mathbf{U_{L}^{*}},\mathbf{U_{R}}^{*}\right)$.

A conservative formulation for the source term $S_{T}$ is obtained by integrating over the cell $C_{i}$ the steady flux gradient in which the nonlinear advection term is neglected. This approximation is valid in low-Shapiro number flow regimes, and therefore particularly appropriate for blood flow in large arteries. The following conservative expression for $S_{T}$ is obtained, expressed in terms of the reconstructed conservative vector $U^{*}$

[TABLE]

where $\left({A_{L,i+\frac{1}{2}}^{*}},{A_{R,i-\frac{1}{2}}^{*}}\right)$ are the reconstructed cross-sectional areas at the left of the cell interface $I_{i+\frac{1}{2}}$ and at the right of the cell interface $I_{i-\frac{1}{2}}$ respectively and $\left({K_{i+\frac{1}{2}}^{*}},{K_{i-\frac{1}{2}}^{*}}\right)$ are the reconstructed arterial wall rigidities at the cell interfaces $I_{i+\frac{1}{2}}$ and $I_{i-\frac{1}{2}}$ respectively. For consistency reasons, we modify the previous expression and write

[TABLE]

To simplify the expression, we introduce the notation

[TABLE]

With these notations, the first order well-balanced finite-volume numerical scheme for system (6) is simply

[TABLE]

with

[TABLE]

It is straightforward to see that HR-LS is well-balanced for the steady states at rest system (59) and provides a good evaluation of low-Shapiro number steady states system (61). It also guarantees the following natural properties of blood flow equations:

• •

the conservation of mass;

• •

the non-negativity of the cross-sectional area $A$;

• •

correct reflection and transmission conditions when variations of vessel impedance occur.

In physiological conditions, the arteries never completely collapse, therefore the numerical scheme no longer needs to be able to compute dry states. Furthermore, as the flow is subcritical and the heart input signal is not discontinuous, transcritical or supercritical regimes and discontinuities of the conservative variables do no occur. Hence the semi-discrete entropy inequality as well as the ability to compute transcritical flows are no longer crucial requirements of the numerical scheme. Finally, the viscosity of the blood and of the arterial wall, that are not taken into account in the theoretical part of this study, are of great importance in arteries and have diffusive and dissipative effects that remove high frequency components and therefore any discontinuity in the conservative variables. This last point will be addressed in the last example section 6.

4.3.   The subsonic hydrostatic reconstruction: HR-S

In [13], Bouchut proposed an extension of HR, referred to as the subsonic hydrostatic reconstruction (HR-S), ideal for blood flow simulations in large arteries. HR-S is well-balanced for all subcritical steady states (58) and also preserves the good properties of HR (listed as bullets in the introduction), that is the positivity of the water height (equivalent of the cross-sectional area $A$), the ability to compute dry states and transcritical flows and a semi-discrete entropy inequality. HR-S is also able to correctly capture wave reflections and transmissions in regions where the impedance of the arterial wall changes. Indeed, the subcritical steady states system (58) coincides with the junction conservation properties (54). However, HR-S requires the resolution of the nonlinear steady state system (58) at each time step at every cell interface presenting a gradient of the artery’s geometrical or mechanical properties. This increases the computational cost compared to HR and HR-LS, especially if the region requiring a well-balanced treatment is not limited to a few cells.

In this section, we present the derivation of HR-S adapted to blood flow in an artery where both variations of cross-sectional area at rest $A_{0}$ and variations of the arterial wall rigidity $K$ are taken into account. HR-S will serve as a reference exactly well-balanced method to be compared to HR and HR-LS. In particular, HR-S will allow us to assess if relaxing the dry-state property and the semi-discrete entropy inequality in HR-LS impacts solutions of blood flow in physiological conditions. With the notations (56) and (58), we are in the framework introduced by Bouchut [13]. Therefore, we will only briefly recall the main steps of the derivation of HR-S and additional details can be found in the cited publication.

4.3.1.   Well-balanced subsonic positivity-preserving reconstruction procedure for the cross-sectional area $A$

Similarly to HR and HR-LS, the well-balanced property is enforced by defining reconstructed variables on both sides of each cell interface $I_{i+\frac{1}{2}}$ according to the reconstruction procedure (58). Variables noted with $"^{*}"$ will refer to the reconstructed variables. Following [13], we introduce the function $f$

[TABLE]

and given the vectors of conservative variables $\mathbf{U_{L}}$ and $\mathbf{U_{R}}$ and the vectors $\bm{\sigma_{L}}$ and $\bm{\sigma_{R}}$ at the left and right of the interface $I_{i+\frac{1}{2}}$ between cells $C_{i}$ and $C_{i+1}$, the discrete analogue of system (58) writes

[TABLE]

with

[TABLE]

Similarly to HR-LS, the reconstruction of the flow rate $Q^{*}$ is straightforward. However, contrary to HR and HR-LS, system (74) is nonlinear in $A^{*}$ and is difficult to solve analytically. To help the resolution of system (74), we recall the following properties (see [13] for details).

For fixed values of $Q$ and $K$, function $f$ admits a minimum in $A_{s}\left(Q,K\right)$ and $m_{s}\left(Q,K\right)$ is the minimum value of $f$

[TABLE]

For fixed values of $Q$ and $K$ and since the function $f$ is convex, system (74) admits a subcritical and a supercritical solution for the cross-sectional area $A$ if $f\left(Q,A,K\right)>m_{s}\left(Q,K\right)$. Furthermore, if $A>A_{s}\left(Q,K\right)$ the flow is subcritical with $U\in\mathbb{U}_{sub}$ and inversely if $A<A_{s}\left(Q,K\right)$ the flow is supercritical with $U\in\mathbb{U}_{sup}$ (see figure 2).

Using these properties, Bouchut [13] proposed a reconstruction procedure for the cross-sectional area $A^{*}$. The first step is to select reconstructions of the variables $Z^{*}$ and $K^{*}$ that preserve the positivity of $A$ and select the subcritical solution of system (74). The following inequalities must be verified to respectively preserve the positivity of $A$ and select the subcritical solution of system (74)

[TABLE]

and

[TABLE]

The inequalities (78) are naturally verified as we consider only subcritical flow configurations. On the contrary, the inequalities (77) are verified if inequalities (78) are true and if $Z^{*}$ and $K^{*}$ are chosen such that $\delta_{L,R}\leq 0$. A simple choice for $Z^{*}$ and $K^{*}$ is

[TABLE]

Given the expressions (79) for $Z^{*}$ and $K^{*}$, we adapted the reconstruction procedure for the cross-sectional area $A^{*}$ proposed by Bouchut [13] to blood flow in arteries with variable cross-sectional area $A_{0}$ and variable arterial wall rigidity $K$. It is summarized in figure 2 and is presented in the algorithm 1. The algorithm 1 describes the steps that need to be followed to obtain the reconstructed cross-sectional area $A_{L}^{*}$, solution of system (74). The same algorithm can be applied to reconstruct $A_{R}^{*}$.

4.3.2.   Well-balanced subsonic first-order numerical scheme

Similarly to HR and HR-LS, a conservative formulation for the source term $S_{T}$ is obtained by integrating over the cell $C_{i}$ the steady flux gradient. However, the nonlinear advection term is no longer neglected and an additional flux term is introduced to take it into account

[TABLE]

where $\left({A_{L,i+\frac{1}{2}}^{*}},{A_{R,i-\frac{1}{2}}^{*}}\right)$ are the reconstructed cross-sectional areas at the left the cell interface $I_{i+\frac{1}{2}}$ and at the right of the cell interface $I_{i-\frac{1}{2}}$ respectively and $\left({K_{i+\frac{1}{2}}^{*}},{K_{i-\frac{1}{2}}^{*}}\right)$ are the reconstructed arterial wall rigidities at the cell interfaces $I_{i+\frac{1}{2}}$ and $I_{i-\frac{1}{2}}$ respectively. The additional fluxes $\mathcal{T}_{L}$ and $\mathcal{T}_{R}$ are chosen such that the numerical scheme satisfies an entropy inequality by interface (see [13] for details). The computation of $\mathcal{T}_{L}$ and $\mathcal{T}_{R}$ is presented in the algorithm 2. Only the steps that need to be followed to obtain $\mathcal{T}_{L}$ are detailed in 2 but similar results are obtained for $\mathcal{T}_{R}$.

Finally, the first-order well-balanced finite-volume numerical scheme for system (6) is still

[TABLE]

with

[TABLE]

In the following section, we present a series of numerical test-cases were we systematically compare HR, HR-LS and HR-S.

5.   Physiological examples in a single artery

In this section we present a series of numerical computations designed to evaluate the performances in physiological conditions of the low-Shapiro hydrostatic reconstruction (HR-LS) in comparison with the hydrostatic reconstruction (HR) and the subsonic hydrostatic reconstruction (HR-S). All quantities are represented in centimeters, grams and seconds, or equivalently "cgs", which are the natural units to describe blood flow. Indeed, the density of blood is close to 1 in "cgs".

The following numerical simulations are performed in a single artery representative of a large artery such as the aorta. Table 1 summarizes the values of the characteristic properties of blood and of the artery, namely the blood density $\rho$, the length $L$ of the artery and the inlet radius at rest and arterial wall rigidity $R_{in}$ and $K_{in}$, all written in "cgs".

We study two geometrical configurations in which both the cross-sectional area at rest $A_{0}$ and the arterial wall rigidity $K$ vary. Both are idealized representations of variations of arteries’ geometrical and mechanical properties encountered in arterial networks. The first configuration is a smooth stenosis and corresponds to a local reduction of the cross-sectional area at rest $A_{0}$. It is a classical arterial pathology caused by the formation of plaque that deposits on the arterial wall and slowly obliterates the vessel. The stenosis is represented in figure 3 and is defined by the following radius at rest $R_{0}$ and arterial wall rigidity $K$

[TABLE]

We choose $x_{s}=\frac{3L}{10}$ and $x_{f}=\frac{7L}{10}$. The second configuration we investigate is a decreasing step, or decreasing discontinuity. It is an idealized representation of a pointwize transition between a parent artery and a smaller daughter artery and is useful to evaluate the reflecting behavior of a numerical method. The decreasing step is represented in figure 3 and is defined by the following radius at rest $R_{0}$ and arterial wall rigidity $K$

[TABLE]

We choose $x_{m}=\frac{L}{2}$. In both configuration, the amplitude of the geometrical and mechanical variations depends on the wall deformation parameter $\Delta\mathcal{G}$. The values of $\Delta\mathcal{G}$ used in the following simulations are taken from table 2 and are chosen to test the limits of the well-balanced methods while staying in the subcritical flow regime. From a well-balanced point of view, each of these two configurations has a different behavior with respect to the cell size $\Delta x$. Indeed, the step configuration is a discontinuity of the cross-sectional area at rest $A_{0}$ and of the arterial wall rigidity $K$, and therefore the amplitude of the variation of the geometrical and mechanical properties of the artery, proportional to $\Delta\mathcal{G}$, is independent of $\Delta x$. On the contrary, the stenosis configuration is a smooth variation of $A_{0}$ and $K$, and therefore the local variation of the artery’s geometrical and mechanical properties at each cell interface will decrease with the cell size $\Delta x$.

We now provide the values of the conservative variables at the inlet and outlet of the computational domain, based on methods detailed in section 3.4. We impose the flow rate $Q_{in}$ at the inlet of the computational domain, in $x=0$. In reality, to control the flow regime, we give the value of the inlet Shapiro number $S_{h,in}$ and compute the inlet flow rate $Q_{in}$ as a function of $S_{h,in}$

[TABLE]

$A_{in}$ and $c_{in}$ are respectively the inlet cross-sectional area and Moens-Korteweg wave speed (12) and are unknown. However, a dimensional analysis of system (9) allows us to show that the inlet Shapiro number $S_{h,in}$ scales as the ratio of the perturbation of the wall’s radius $\Delta R=R-R_{0}$ over the radius at rest $R_{0}$. With this scaling law, we can estimate a value of the inlet cross-sectional area $A_{in}$ consistent with the inlet Shapiro number $S_{h,in}$ and we compute $A_{in}$ as

[TABLE]

At the outlet of the computational domain, in $x=L$, we either impose the reflection coefficient $R_{t}=0$ or the cross-sectional area $A_{out}$, depending on the test case. Similarly to the inlet cross-sectional area $A_{in}$, we compute the outlet cross-sectional area as a function of $S_{h,in}$

[TABLE]

The values of the inlet Shapiro number $S_{h,in}$ and the wall deformation parameter $\Delta\mathcal{G}$ used in the following simulations are presented in table 2. They cover a wide range of physiological configurations, allowing us to assess the behavior of the three numerical schemes under consideration in the limit of the low-Shapiro number flow regime. We recall that in arteries the average Shapiro number is of the order of $S_{h}=1\times 10^{-2}$.

5.1.   Steady solutions

We evaluate the well-balanced properties of HR, HR-LS and HR-S by computing steady solutions of system (6) in the smooth stenosis (83) and the decreasing step (84). Steady flow configurations in arterial geometries similar to the stenosis (83) have been studied by Müller [41], where only variations of the wall rigidity $K$ are taken into account. In [43], the authors computed steady solutions in tapered tubes. In the context of the shallow water equations, steady flow solutions over a bump (analogue of the stenosis) or a step have been studied by many authors [16, 44, 17, 22].

The steady numerical solutions are obtained for $t=200$ $s$. The time step $\Delta t$ is constant and chosen such that the CFL condition (41) is always satisfied. We impose the flow rate $Q_{in}$ (85) at the inlet and the cross-sectional area $A_{out}$ (87) at the outlet. We therefore select a specific steady state characterized by its associated flow rate $Q_{st}$ and energy discharge $E_{st}$. These values can be computed analytically and provide exact solutions to compare with our numerical results

[TABLE]

In both configurations (83) and (84), we perform a series of 12 numerical computations for all combinations of the inlet Shapiro number $S_{h,in}$ and the wall deformation parameter $\Delta\mathcal{G}$ taken from table 2. Table 3 shows $L^{1}$ relative errors between the analytic solutions and the results obtained with HR, HR-LS and HR-S for a fixed number of cells $N=50$.

In both the stenosis (83) and the step (84) configurations, the results are similar and indicate that, as expected, each numerical method is exactly well-balanced for the steady states at rest ($S_{h,in}=0$). Only HR-S is exactly well-balanced for all considered subcritical steady states. For the low-Shapiro number steady states ($S_{h,in}=10^{-3},10^{-2},10^{-1}$), HR-LS is more accurate than HR. However, the accuracy of HR-LS diminishes when the values of $S_{h,in}$ and $\Delta\mathcal{G}$ increase, and for $S_{h,in}=1\times 10^{-1}$ and $\Delta\mathcal{G}=30\%$, in the limit of the low-Shapiro number flow regime, HR-LS is only one order of magnitude more accurate than HR. Interestingly, the errors obtained with HR are independent of the inlet Shapiro number $S_{h,in}$, but increase significantly with the wall deformation parameter $\Delta\mathcal{G}$.

To test the consistency and the order of convergence of the different methods, we perform a convergence study for the average low-Shapiro steady configuration $S_{h,in}=1\times 10^{-2}$ and $\Delta\mathcal{G}=10\%$ in both the stenosis and the step configurations. $L^{1}$ relative errors with analytic solutions are presented in table 4 for the following number of cells $N\in\left\{50,100,200,400\right\}$.

In the stenosis configuration (83), both HR and HR-LS converge with order 1, whereas in the step configuration (84), they do not achieve order 1 convergence. Indeed, in the stenosis configuration, the variations of the artery’s geometrical and mechanical properties at each cell interface decrease when the number of cells $N$ increases, enabling the convergence of both methods. On the contrary, the geometrical and mechanical variations remain unchanged in the step configuration when the number of cells $N$ increases. These observations are illustrated by figures 4 and 5, where we respectively plot the spatial evolution of the flow rate $Q$ and the energy discharge $E$ with the number of cells in the stenosis and step configurations.

In both configurations, the values of the errors obtained in table 4 with HR-S are of the order of machine precision, indicating that HR-S is exactly well-balanced for the considered low-Shapiro steady state. However, the errors increase slightly with the number of cells. Similar behaviors are observed in convergence studies presented in [17] for an exactly well-balanced method. In our case, this phenomenon is due to a small error between the computed boundary conditions and those required to obtain the desired steady state, and is not caused by HR-S.

The results indicate that among the three well-balanced methods considered, HR is the least accurate when computing low-Shapiro number steady solutions in an artery presenting smooth and discontinuous variations of its cross-sectional area at rest $A_{0}$ and of its arterial wall rigidity $K$. On the contrary, HR-S is the only exactly well-balanced method for the considered low-Shapiro number steady states. Finally, even though HR-LS is not exactly well-balanced for the considered low-Shapiro number steady states, it allows to compute with satisfying accuracy steady solutions for smooth and discontinuous variations of the artery’s geometrical and mechanical properties. These results show that the system (61) (used by HR-LS) is a better approximation than system (60) (used by HR) of the steady state system (58) (used by HR-S) in low-Shapiro flow configurations.

5.2.   Single wave propagation

The wave-capturing properties of HR, HR-LS and HR-S are now evaluated. We simulate the propagation of a single wave in the smooth stenosis (83) and the decreasing step (84). The step configuration was studied in [21, 20, 58] for an artery with only variations of its cross-sectional area at rest $A_{0}$.

The results are obtained for $t=0.045$ s. The time step $\Delta t$ is constant and chosen such that the CFL condition (41) is always satisfied. We impose a single pulse of flow at the inlet of the computational domain and the unsteady inlet flow rate $Q_{in}\left(t\right)$ is defined as

[TABLE]

We choose $T_{pulse}=0.04$ s to artificially reduce the wave length of the pulse for visualization purposes and the value of $Q_{pulse}$ is a function of the inlet Shapiro number $S_{h,in}$ and is defined as in equation (85). Figure 6 represents the function $Q_{in}$ for $S_{h,in}=1\times 10^{-2}$. At the outlet of the computational domain, we set the reflection coefficient $R_{t}=0$ to remove any terminal reflection.

5.2.1.   The step configuration

We focus on the decreasing step configuration (84). Given the inlet condition (89), the pulse wave propagates in the artery starting from the left-hand side of the domain until it reaches the step. The change of impedance of the vessel creates reflected and transmitted waves that need to be captured by the numerical scheme. A linear analytic solution was proposed in [50] and validated in [21, 20, 57], and gives the expression of the reflection coefficient $R_{t}$ and the transmission coefficient $T_{t}$, based on the conservation properties (55)

[TABLE]

where $Y=A/\left(\rho c\right)$ is the vessel admittance. Subscripts $L$ and $R$ respectively refer to the values at the left and right of the step. As the coefficients $R_{t}$ and $T_{t}$ do not depend on the frequency of the incoming wave, we can analytically predict the position, shape and amplitude of the linear reflected and transmitted waves. However, as the inlet Shapiro number $S_{h,in}$ is non-zero, the flow is nonlinear and the linear analytic solution (90) is only valid in the asymptotic limit $S_{h,in}\rightarrow 0$. To evaluate the quality of the results obtained with HR, HR-LS and HR-S, we compute reference solutions, obtained with HR-S for $N=25600$ and values of $S_{h,in}$ and $\Delta\mathcal{G}$ taken from table 2. To assess the validity of these reference solutions, we compare them to the linear analytic solutions (90) in figure 7. We observe that for low values of the inlet Shapiro number $S_{h,in}$ (figure 7 left), for which the linear approximation is valid, the analytic and reference solutions match. As expected, for higher values of $S_{h,in}$, the flow is no longer linear and the propagation speed as well as the amplitude of the reflected and transmitted waves change (figure 7 center and right).

We present results only for the flow rate $Q$ to reduce the number of variables and simplify the analysis of the results. Similar conclusions to those presented hereafter would have been drawn if we had considered the pressure $P$ or the wall perturbation $R-R_{0}$.

We perform a series of 9 numerical computations with different combinations of the non-zero inlet Shapiro number $S_{h,in}$ and the wall deformation parameter $\Delta\mathcal{G}$ taken from table 2. Table 5 shows $L^{1}\left[Q\right]$ relative errors between the reference solutions and the results obtained with HR, HR-LS and HR-S for a fixed number of cells $N=1600$. We choose a high value of $N$ to reduce the numerical dissipation and highlight the effects of the well-balanced methods.

The results obtained with HR, HR-LS and HR-S are almost identical and indicate that each method is able to correctly compute the expected reflected and transmitted waves. For each method, the error $L^{1}\left[Q\right]$ is independent of the inlet Shapiro number $S_{h,in}$ but increases with the wall deformation parameter $\Delta\mathcal{G}$. However, the error obtained with HR increases faster with $\Delta\mathcal{G}$ than with the other methods. In particular, for $\Delta\mathcal{G}=30\%$, the value of $L^{1}\left[Q\right]$ obtained with HR is one order of magnitude higher than the one obtained with HR-LS or HR-S.

This last point is corroborated by figures 8, 9 and 10, where we represent the spatial evolution of the flow rate $Q$ at $t=0.045$ s, obtained using $N=100$ (left) and $N=1600$ (right) for $S_{h,in}=1\times 10^{-2}$ and $\Delta\mathcal{G}=\left\{10\%,30\%,60\%\right\}$ respectively. In each figure, we compare the results obtained using HR, HR-LS and HR-S to the corresponding reference solution and observe if increasing the number of cells allows the numerical solution to converge towards the reference solution. In figure 8, the results obtained for $\Delta\mathcal{G}=10\%$ with HR, HR-LS and HR-S are similar and indicate that each numerical solution converges towards the reference solution. On the contrary, in figure 9 for $\Delta\mathcal{G}=30\%$ and in figure 10 for $\Delta\mathcal{G}=60\%$, only the solutions obtained with HR-LS and HR-S converge towards the reference solution. HR is unable to compute the expected amplitude of the reflected and transmitted waves and overestimates the amplitude of the reflected wave and underestimates the amplitude of the transmitted wave.

The results indicate that HR-LS and HR-S are able to compute wave reflections and transmissions in an artery presenting an arbitrary large discontinuous variation of its cross-sectional area at rest $A_{0}$ and arterial wall rigidity $K$. On the contrary, HR is unable to compute the correct amplitude of the reflected and transmitted waves when the discontinuous variation of the artery’s geometrical and mechanical properties is too large, independently of the number of cells $N$. Moreover, these results show that the system (61) (used by HR-LS) has the appropriate conservation properties to compute wave reflections for arbitrary large discontinuous geometrical and mechanical variations in low-Shapiro number flow regimes. On the contrary, HR, using the system (60), is only able to compute wave reflections for small discontinuous variations of the artery’s properties ($\Delta\mathcal{G}=10\%$, see figure 8). This last point can be problematic as large variations of the artery’s geometrical and mechanical properties can be encountered when modeling arterial pathologies such as stenoses.

5.2.2.   The stenosis configuration

In this subsection we focus on the stenosis configuration (83). To evaluate the quality of the results obtained with HR, HR-LS and HR-S, we compute reference solutions, obtained with HR-S for $N=25600$ and values of $S_{h,in}$ and $\Delta\mathcal{G}$ taken from table 2. As the variation of geometrical and mechanical properties of the artery is smooth, the observed flow rate is constituted of a continuum of reflected and transmitted waves that are created at each cell interface, where the artery’s geometrical and mechanical properties are discontinuous.

Similar results to those of subsection 5.2.1 are obtained, and therefore we do not completely repeat the previous analysis. In figures 11, 12 and 13, we present the spatial evolution of the flow rate $Q$ at $t=0.045$ s, obtained using $N=100$ (left) and $N=1600$ (right) for $S_{h,in}=1\times 10^{-2}$ and $\Delta\mathcal{G}=\left\{10\%,30\%,60\%\right\}$ respectively. In each figure, we compare the results obtained using HR, HR-LS and HR-S to the corresponding reference solution and observe if increasing the number of cells allows the numerical solution to converge towards the reference solution. Contrary to the step configuration studied in subsection 5.2.1, the results obtained with HR, HR-LS and HR-S are similar and indicate that each numerical solution converges towards the reference solution. However, for $\Delta\mathcal{G}=\left\{30\%,60\%\right\}$ and $N=100$, HR is less accurate than HR-LS and HR-S.

These results are coherent with those of subsection 5.2.1. Indeed, when studying the step configuration, we showed that contrary to HR-LS and HR-S, HR overestimates the amplitude of the reflected wave and underestimates the amplitude of the transmitted wave when a large discontinuous variation of the artery’s geometrical and mechanical properties is considered ($\Delta\mathcal{G}=\left\{30\%,60\%\right\}$). As the stenosis is a smooth variation of the cross-sectional area at rest $A_{0}$ and of the arterial wall rigidity $K$, discontinuous variations of the arterial wall’s geometrical and mechanical properties occur at each cell interface and the amplitude of these variations decreases with the number of cells $N$. Hence, for $\Delta\mathcal{G}=\left\{30\%,60\%\right\}$ and $N=100$, the local discontinuous variations of the artery’s properties are large enough for HR to be inaccurate. On the contrary, for $N=1600$, the local discontinuous variation of the artery’s geometrical and mechanical properties are sufficiently small for HR to be as accurate as HR-LS and HR-S.

We have studied the wave capturing behavior of HR, HR-LS and HR-S. We showed that for arbitrary large smooth or discontinuous variations of the artery’s cross-sectional area at rest $A_{0}$ and arterial wall rigidity $K$, both HR-LS and HR-S are able to compute the expected reflected and transmitted waves. On the contrary, HR is unable to correctly compute reflected and transmitted waves when large discontinuous variations of the artery’s properties are considered. In particular, HR overestimates the reflected wave and underestimates the transmitted wave. Therefore, HR-LS and HR-S are good choices to compute wave reflections and transmissions in low-Shapiro flow regimes. In the following subsection, we will analyze the behavior of the different well-balanced methods in large network computations, where multiple effects come into play.

6.   A realistic example: stenosis of the iliac artery in a 55 arteries network

We study the response at the systemic level of a model network to the presence of a pathology. Indeed, the observed pressure and flow waveforms in the systemic network are the superposition of multiple reflected and transmitted waves, generated at each arterial bifurcation and dampened and diffused by the viscosity of the blood and the arterial wall. The presence of a pathology creates additional reflected and transmitted waves that change the reflection pattern and therefore the shape and amplitude of the observed waveforms. When such pathologies are modeled in the network, a well-balanced method is required to take into account the geometrical and mechanical source term induced by the local variations of the cross-sectional area and arterial wall rigidity representing the pathology.

In the purpose of performing large network blood flow simulations, we use the arterial network proposed by Sherwin in [54] which was adapted from Westerhof [59], describing 55 of the great arteries of the systemic network (human arterial network). This model was more recently used by Wang [57] to perform viscoelastic blood flow simulations using different numerical methods. The network under consideration is represented in figure 14. The parameters of the model were obtained using physiological data and in each artery the geometrical and mechanical parameters do not vary with the axial position $x$. Therefore, in the absence of an arterial pathology, a well-balanced method is not required to compute blood flow in the considered network. The details of the parameters of the model are not recalled here and we refer the reader to the cited publications.

The pathology considered is a stenosis of the right iliac artery (artery 49 in figure 14). We consider two possible shapes for the stenosis. The first corresponds to a succession of an increasing and a decreasing step and will be referred to as the square stenosis. It is defined by the following radius at rest $R_{0}$ and arterial wall rigidity $K$

[TABLE]

We choose $x_{s}=6.25$ cm and $x_{f}=8.25$ cm. The second geometry is the stenosis (83) presented in subsections 5.1 and 5.2.2. Its radius at rest $R_{0}$ and arterial wall rigidity $K$ vary as

[TABLE]

We choose $x_{s}=5.5$ cm and $x_{f}=9.5$ cm. We will refer to this stenosis as the cos stenosis. The cos stenosis is twice as long as the square stenosis to match the deformation area of the square stenosis. However, the maximal amplitudes of both configuration are the same and are proportional to $\Delta\mathcal{G}$.

In [26], the authors studied a similar pathological network and showed that the presence of the stenosis has a noticeable impact on the global hemodynamics for large values of $\Delta\mathcal{G}$. To that effect, we choose $\Delta\mathcal{G}=65\%$.

The results presented in this section are obtained using a time step $\Delta t=5\times 10^{-5}$ and mesh size $\Delta x=0.2$ cm and are compared to results obtained with the 55 arteries network without the pathology. This network will be referred to as the "Sane" network and does not require the use of any well-balanced method. We focus on four measurement points corresponding to typical measurement points used by medical practitioners during surgery. These points are situated in the middle of the following arteries, and the numbers indicated correspond to the numbering of the arteries in figure 14: the Left Subclavian II (20), Left Femoral (45), Right Femoral (51) and Right External Iliac (49), before the stenosis. Furthermore, as the pressure $P$ is the most common and simple variable to measure in vivo, we present only pressure waveforms in the following.

In figure 15, we compare the pressure waveforms obtained using HR, HR-LS and HR-S for the cos stenosis. The results obtained using HR-LS and HR-S are almost identical in the arteries 20 and 45, suited far from the stenosis, but also in the artery 49, located before the stenosis. In artery 51, situated after the stenosis, small differences exist between the results obtained with HR-LS and HR-S, especially during diastole ($t=7.5$ s and $t=8$ s). Moreover, the results obtained with HR-LS and HR-S in arteries 20, 45 and 49 are very close to those obtained with the Sane network, indicating that in these arteries, the resistive behavior of the stenosis is negligible compared to the global resistance of the network. However, in artery 51, the results obtained with HR-LS and HR-S slightly differ from those obtained with the Sane network, indicating that the presence of the stenosis has a local effect, especially during diastole ($t=7.5$ s to $t=8$ s). HR produces significantly different results from HR-LS and HR-S in each artery considered. In arteries 20, 45 and 49, HR overestimates the amplitude of the pressure waveform, whereas in artery 51 it underestimates it. These results are in good accord with the observations made in figures 10 and 13 in subsections 5.2.1 and 5.2.2.

In figure 16, we compare the pressure waveforms obtained using HR-LS for the cos and the square stenosis. The results indicate that in arteries 20, 45 and 49, there are no significant differences when using either the cos or the square stenosis. Only in artery 51 is the influence of the shape noticeable.

In figure 15, the effects of the flow viscosity and the wall viscoelasticity are neglected. However, they play an important role in the global hemodynamics and need to be taken into account to obtain an accurate description of pressure and flow waves in a network. In figure 17 , we present similar results to those obtained in figure 15, but now viscous and viscoelastic effects are taken into account. For the implementation of the viscous and viscoelastic terms, we refer the reader to [57]. The results indicate that viscosity and viscoelasticity have a dissipative and diffusive effect and in the presence of such effects, the results obtained with HR-LS and HR-S overlap, even in artery 51. However, HR still overestimates the amplitude of the pressure waves in arteries 20, 45 and 49 and underestimates it in artery 51.

The results presented in this section indicate that the pressure waveform is sensitive to the choice of the well-balanced method. Even though HR-LS, HR-S have comparable behaviors, the results obtained with these methods are very different from those obtained with HR. On the contrary, changing the shape of the pathology has little effect on the shape and amplitude of the pressure waveforms. The small differences between the results are erased when blood viscosity and wall viscoelasticity are taken into account, due to the damping and diffusion behavior of the fluid and wall viscosities. Overall, HR-LS and HR-S behave similarly and produce satisfying results.

7.   Conclusion

We introduced two well-balanced hydrostatic reconstruction techniques for blood flow in large arteries with varying geometrical and mechanical properties. The low-Shapiro hydrostatic reconstruction (HR-LS) is a simple and efficient well-balanced reconstruction technique, inspired from the hydrostatic reconstruction technique (HR) proposed in [6, 21]. It accurately preserves low-Shapiro number (equivalent of the Froude number for shallow water equations and the Mach number for compressible Euler equations) steady states that may occur in large network simulations and are the appropriate conserved properties at discontinuities of the geometrical and mechanical properties of the artery. The subsonic hydrostatic reconstruction (HR-S), introduced in [13] and adapted here to blood flow, exactly preserves all subcritical steady states. We performed a series a numerical computations to compare the properties of HR, HR-LS and HR-S. In all numerical computations, HR was the least accurate method and was unable to correctly compute wave reflection and transmission when large variations of the artery’s geometrical and mechanical properties were considered. HR-S proved to be exactly well-balanced for all low-Shapiro number steady states and the most accurate reconstruction technique. We showed that HR-LS is well-balanced only for steady states at rest, but provides satisfactory approximations of low-Shapiro steady states. HR-LS is also able to capture wave reflections and transmissions for arbitrary large variations of the artery’s geometrical and mechanical properties, which is an essential property to compute realistic flow and pressure waveforms. We have also evaluated the sensitivity of the model to well-balanced methods and to the shape of the pathology in an 55 arteries network simulation. We showed that the model is not sensitive to the geometry of the pathology. However, important differences were observed between HR and the other well-balanced methods, namely HR-LS and HR-S, due to the fact that HR is unable to capture wave reflection and transmission. Finally, we observed that the small differences between HR-LS and HR-S are erased when adding viscous and viscoelastic effects, which are required to obtain realistic pressure and flow waveforms. This analysis allows us to conclude that both HR-LS and HR-S are adequate well-balanced methods to compute blood flow in large arteries with varying cross-sectional area at rest and arterial wall rigidity. However, in large networks where many arteries present variations of their geometrical and mechanical properties, the extra iterations required by HR-S increase the computational cost compared to HR-LS. We therefore recommended using HR-LS in this case, as it is a good compromise between simplicity, numerical accuracy and efficiency. In future works, we will investigate further the properties of HR-LS and propose an extension of the method to higher order.

8.   Acknowledgments

The authors are grateful to thank F. Bouchut and E. Audusse for their helpful remarks and comments.