A Steady State Preserving Numerical Scheme for One-Dimensional Blood Flow Model

1. Introduction

Blood flow in arteries can be described mathematically by the Navier-Stokes equations in three dimensions (3D). However, numerical approximations for three-dimensional models are computationally expensive. Therefore, the simpler one-dimensional (1D) models are considered, which provide a realistic description of certain parts of the cardiovascular system and also provide boundary conditions for the more complex three-dimensional models. A simplified 1D model that describes the blood flow in arteries is given by the following system of partial differential equations (see e.g. [1])

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial A}{\partial t}}+{\displaystyle \frac{\partial \left(AU\right)}{\partial x}}=0\hfill \\ {\displaystyle \frac{\partial U}{\partial t}}+U{\displaystyle \frac{\partial U}{\partial x}}+{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial P}{\partial x}}={\displaystyle \frac{\mathcal{F}}{\rho A}},\hfill \end{array}\right.$$

(1)

where $A(x,t)=\pi R^2(x,t)$ is the cross-sectional area of the vessel at spatial position x and time t with $R(x,t)$ being the radius, $U(x,t)$ is the mean blood velocity in the axial direction, $\rho$ is the blood density, assumed to be constant and equal to $1060\mathrm{kg}/{\mathrm{m}}^3$, $P=P\left(A\right)$ is the internal pressure and $\mathcal{F}(x,t)$ is the friction force per unit length. The present work is restricted to the inviscid limit where $\mathcal{F}=0$.

Several expressions can be found in the literature to describe the pressure. In the present work we use the algebraic expression proposed in [2]:

$$P=P_{ext}+\beta (\sqrt{A}-\sqrt{A_0}),$$

(2)

where the external pressure $P_{\mathrm{ext}}$, the coefficient $\beta$ stands for the arterial stiffness and is supposedly constant and $A_0=A_0\left(x\right)$ is the vessel cross-sectional area at rest, which may be variable in the case of some pathologies.

We refer to equations (1) as the system $(A,U)$. Such equations are an example of a system of balance laws

$$\frac{\partial \mathbf{w}}{\partial t}+\frac{\partial \mathbf{f}\left(\mathbf{w}\right)}{\partial x}=\mathbf{S}\left(\mathbf{w}\right)$$

(3)

where

$$\begin{array}{cc}\hfill & \mathbf{w}=\left[\begin{array}{c}A\\ U\end{array}\right],\mathbf{f}\left(\mathbf{w}\right)=\left[\begin{array}{c}AU\\ {\displaystyle \frac{U^2}{2}}+{\displaystyle \frac{\beta \sqrt{A}}{\rho }}\end{array}\right],\mathbf{S}\left(\mathbf{w}\right)=\left[\begin{array}{c}0\\ {\displaystyle \frac{\beta }{\rho }}{\displaystyle \frac{\partial \sqrt{A_0}}{\partial x}}\end{array}\right]\hfill \end{array}$$

correspond to the vector of unknowns, flux function and source term, respectively. In order to analyse the characteristic information it is convenient to write the system in the non-conservative form

$${\displaystyle \frac{\partial \mathbf{w}}{\partial t}}+\mathbf{J}\left(\mathbf{w}\right){\displaystyle \frac{\partial \mathbf{w}}{\partial x}}=\mathbf{S}(x,\mathbf{w}),$$

where

$$\mathbf{J}\left(\mathbf{w}\right)=\left[\begin{array}{cc}U& A\\ c^2/A& U\end{array}\right]$$

is the Jacobian matrix of the flux function with eigenvalues

$${\lambda }_1=U-c\mathrm{y}{\lambda }_2=U+c$$

and the corresponding right eigenvectors are given by

$${\mathbf{r}}_1={[-1,c/A]}^{\top }\mathrm{y}{\mathbf{r}}_2={[1,c/A]}^{\top }.$$

(4)

Here, $c=\sqrt{\frac{\beta \sqrt{A}}{2\rho }}$ is the Moens-Korteweg wave speed.

It is well known that solutions of balance laws may contain discontinuities, even for smooth initial data. Consequently, the solutions of such systems must be considered in weak sense. Moreover, the weak solutions may not be unique, then it is necessary to impose an additional admissibility criterion or entropy conditions to select among the possible solutions, the one that is physically relevant [3]. In particular, the system $(A,U)$ admits the entropy inequality

$${\displaystyle \frac{\partial \widehat{\eta }\left(\mathbf{w}\right)}{\partial t}}+{\displaystyle \frac{\partial \widehat{G}\left(\mathbf{w}\right)}{\partial x}}\le 0,$$

(5)

where $(\widehat{\eta },\widehat{G})$ is the extended entropy pair introduced in [4]:

$$\widehat{\eta }=\eta -{\displaystyle \frac{\beta \sqrt{A_0}}{\rho }}A,\widehat{G}=G-{\displaystyle \frac{\beta \sqrt{A_0}}{\rho }}AU,$$

(6)

and $(\eta ,G)$ is the entropy pair $(\eta ,G)$ given by

$$\eta =\frac{1}{2}A{U}^2+{\displaystyle \frac{2\beta A^{3/2}}{3\rho }},G=\frac{1}{2}A{U}^3+{\displaystyle \frac{\beta U{A}^{3/2}}{\rho }}$$

(7)

Using the entropy pair we may obtain the entropy variables

$$\mathbf{v}\left(\mathbf{w}\right)=\nabla \eta ={\left[\frac{U^2}{2}+{\displaystyle \frac{\beta \sqrt{A}}{\rho }},AU\right]}^{\top }$$

(8)

with the corresponding entropy potential

$$\begin{array}{cc}\hfill \psi \left(\mathbf{v}\right)& :={\mathbf{v}}^{\top }\mathbf{f}\left(\mathbf{w}\left(\mathbf{v}\right)\right)-G\left(\mathbf{w}\left(\mathbf{v}\right)\right)={\displaystyle \frac{A{U}^3}{2}}+{\displaystyle \frac{\beta U{A}^{3/2}}{\rho }}.\hfill \end{array}$$

(9)

In the extended case, we have

$$\widehat{\mathbf{v}}\left(\mathbf{w}\right)=\nabla \widehat{\eta }=\mathbf{v}-{[\beta \sqrt{A_0}/\rho ,0]}^{\top }.$$

(10)

Systems as (3) are of special interest because of the delicate balance that must exist between the convective term and the source term during the time evolution. Specifically, the steady-states solutions for the system (3) are the solutions $\mathbf{w}\left(x\right)$ which are time-independent. Thus, steady-state solutions satisfy the system

$$\frac{\partial \mathbf{f}\left(\mathbf{w}\right)}{\partial x}=\mathbf{S}\left(\mathbf{w}\right).$$

(11)

In particular, the system $(A,U)$ admits the following steady-state solution, known as the (non-zero pressure) man-at-eternal-rest steady state or dead-man equilibrium [5]:

$$\begin{array}{c}\hfill U=0\mathrm{y}\sqrt{A}-\sqrt{A_0}=C=\mathrm{constant}.\end{array}$$

(12)

When $C=0$, the (zero pressure) man-at-eternal-rest steady state is obtained, namely

$$\begin{array}{c}\hfill U=0\mathrm{and}A=A_0.\end{array}$$

(13)

It is crucial that numerical methods to approximate steady-state solutions also satisfy a discrete version of (11). Such methods are termed well-balanced schemes.

2. Numerical method

Since, in general, it is not possible to find an analytical solution of the system (3), it is imperative to approximate the solution numerically. To this end, we consider a semi-discrete finite volume scheme for (3) on a uniform spatial mesh with nodes $x_j=j\Delta x$, $j\in \mathbb{Z}$:

$$\frac{\mathrm{d}{\mathbf{w}}_j\left(t\right)}{\mathrm{d}t}=-\frac{1}{\Delta x}\left({\mathbf{F}}_{j+1/2}-{\mathbf{F}}_{j-1/2}\right)+{\mathbf{S}}_j,$$

(14)

where ${\mathbf{w}}_j\left(t\right)$ is the cell average on $I_j=[x_{j-1/2},x_{j+1/2})$, ${\mathbf{F}}_{j+1/2}$ is the numerical flux associated with $x_{j+1/2}$ and ${\mathbf{S}}_j$ is a discretization of the source term. It should be noted that the first term on the right-hand side corresponds to the discretisation of the spatial derivative $-\mathbf{f}{\left(\mathbf{w}\right)}_x$. The discretization of the time derivative will be discussed later.

The scheme (14) is entropy stable with respect to the entropy pair $(\widehat{\eta },\widehat{G})$ if it satisfies a discrete version of the entropy inequality (5), that is,

$$\frac{\mathrm{d}}{\mathrm{d}t}\widehat{\eta }\left({\mathbf{w}}_j\left(t\right)\right)+\frac{1}{\Delta x}\left({\widehat{\mathcal{G}}}_{j+1/2}-{\widehat{\mathcal{G}}}_{j-1/2}\right)\le 0$$

(15)

for some numerical entropy flux ${\widehat{\mathcal{G}}}_{j+1/2}$ consistent with the entropy flux $\widehat{G}$. If equality holds in (15), then the scheme (14) is called entropy conservative.

From now on, we use

$${\left[\left[a\right]\right]}_{j+1/2}:=a_{j+1}-a_j,{\left\{\left\{a\right\}\right\}}_{j+1/2}:=\frac{1}{2}(a_{j+1}+a_j).$$

to denote, respectively, the jump of a across the interface $x_{j+1/2}$ and the arithmetic mean of a quantity.

In this work we will used a numerical flux of the form [7]:

$${\mathbf{F}}_{j+1/2}={\tilde{\mathbf{F}}}_{j+1/2}-\frac{1}{2}{R}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2},$$

(16)

where

(i).

${\tilde{\mathbf{F}}}_{j+1/2}={[{\tilde{F}}_{1,j+1/2},{\tilde{F}}_{2,j+1/2}]}^{\top }$ is the second order and entropy conservative flux for the homogeneous case of (3) given by

$${\tilde{F}}_{1,j+1/2}={\left\{\left\{AU\right\}\right\}}_{j+1/2},{\tilde{F}}_{2,j+1/2}=\frac{1}{2}{\left\{\left\{U^2\right\}\right\}}_{j+1/2}+\frac{\beta }{\rho }{\left\{\left\{\sqrt{A}\right\}\right\}}_{j+1/2},$$

(17)

with the corresponding numerical entropy flux

$${{\tilde{G}}_{j+1/2}:={\left\{\left\{{\mathbf{v}}^{\top }\right\}\right\}}_{j+1/2}{\tilde{\mathbf{F}}}_{j+1/2}-\left\{\right\{\psi \left\}\right\}}_{j+1/2}.$$

(18)

The numerical flux (17) was obtained in [8] from the condition (see [6])

$${\left[\left[{\mathbf{v}}^{\top }\right]\right]}_{j+1/2}{\tilde{\mathbf{F}}}_{j+1/2}={\left[\left[\psi \right]\right]}_{j+1/2}.$$

(19)

(ii).

$R_{j+1/2}$ is the matrix of right eigenvectors of the Jacobian matrix $\mathbf{J}\left({\mathbf{w}}_{j+1/2}\right)$ being evaluated at the average state ${\mathbf{w}}_{j+1/2}:=({\mathbf{w}}_j+{\mathbf{w}}_{j+1})/2$, ${\Lambda }_{j+1/2}=diag\left(\right|{\lambda }_1|,|{\lambda }_2\left|\right)$ is a Roe-type matrix and ${\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}={\mathbf{z}}_{j+1/2}^{+}-{\mathbf{z}}_{j+1/2}^{-}$ ${\mathbf{z}}_{j+1/2}^{-}$ y ${\mathbf{z}}_{j+1/2}^{+}$ denoting, respectively, the left and right limiting values of the scaled entropy variables $\mathbf{z}:=R_{j+1/2}^{\top }\widehat{\mathbf{v}}$ at interface $x_{j+1/2}$, obtained by ENO reconstruction. The choice of this ENO method is due to the fact that it satisfies the so-called sign property [9]:

$$\mathrm{sign}\left({\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}\right)=\mathrm{sign}\left({\left[\left[\mathbf{z}\right]\right]}_{j+1/2}\right),$$

(20)

which will be useful in proving entropy stability of the flux (16).

The discretization of the second component of the source term $B\left(x\right):={\partial }_x\sqrt{A_0\left(x\right)}$ is performed using the simple and well-known second-order centered difference approximation

$$B^{\prime }\left(x\right)\approx \frac{B(x+\Delta x)-B(x-\Delta x)}{2\Delta x}.$$

Thus, the discretization of the source term yields

$${\mathbf{S}}_j=\frac{\beta }{\rho }\left[\begin{array}{c}0\\ {\displaystyle \frac{B_{j+1}-B_{j-1}}{2\Delta x}},\end{array}\right]$$

(21)

where $B_j$ denotes the discrete approximation of $B\left(x_j\right)=\sqrt{A_0\left(x_j\right)}$.

3. Theoretical results and numerical experiments

In this section it is shown that the scheme (14) with numerical flux (16) and discretization of the source term (21) is entropy stable and well- balanced. The proof is similar to that developed for the Saint-Venant equations and the relativistic magnetohydrodynamics equations in [10] and [11], respectively. Then, some examples are included to check that the numerical scheme verifies these properties at the discrete level.

Theorem 1.

The numerical scheme (14) with numerical flux (16) and discretization of the source term (21) satisfies the following properties:

(i) 

It is entropy stable, i.e., satisfies the discrete entropy inequality (15), where $\widehat{\eta }$ is the entropy function given by (6) and the corresponding numerical entropy flux is

$${\widehat{\mathcal{G}}}_{j+1/2}={\tilde{\mathcal{G}}}_{j+1/2}-{\tilde{\mathcal{D}}}_{j+1/2}$$

(22)

with

$$\begin{array}{cc}\hfill {\tilde{\mathcal{G}}}_{j+1/2}& ={\tilde{G}}_{j+1/2}-\frac{\beta (B_j{A}_{j+1}{U}_{j+1}+B_{j+1}{A}_j{U}_j)}{2\rho },\hfill \end{array}$$

(23)

$${\tilde{\mathcal{D}}}_{j+1/2}=\frac{1}{2}{\left\{\left\{{\widehat{\mathbf{v}}}^{\mathrm{T}}\right\}\right\}}_{j+1/2}{R}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}.$$

(24)

(ii) 

it preserves the discrete version of the man-at-eternal-rest (12), this means that given the initial data

$$U_j=0,\sqrt{A_j}-B_j\equiv C\forall j,$$

(25)

with C a constant, then the solution obtained with the scheme (14) satisfies

$$\frac{d{\mathbf{w}}_j}{dt}=0\forall j.$$

Proof. 

Let us first prove entropy stability. Multiplying (14) by ${\widehat{\mathbf{v}}}_j^{\top }$ yields

$$\begin{array}{cc}\hfill {\widehat{\mathbf{v}}}_j^{\top }\frac{\mathrm{d}{\mathbf{w}}_j\left(t\right)}{\mathrm{d}t}& =-\frac{1}{\Delta x}{\widehat{\mathbf{v}}}_j^{\top }\left({\mathbf{F}}_{j+1/2}-{\mathbf{F}}_{j-1/2}\right)+{\widehat{\mathbf{v}}}_j^{\top }{\mathbf{S}}_j\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =-\frac{1}{\Delta x}{\widehat{\mathbf{v}}}_j^{\top }\left({\tilde{\mathbf{F}}}_{j+1/2}-{\tilde{\mathbf{F}}}_{j-1/2}\right)\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & -\frac{1}{\Delta x}{\widehat{\mathbf{v}}}_j^{\top }\left(-\frac{1}{2}{R}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}+\frac{1}{2}{R}_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2}\right)\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & +{\widehat{\mathbf{v}}}_j^{\top }{\mathbf{S}}_j\hfill \end{array}$$

(29)

The left-hand side of (26) gives

$${\widehat{\mathbf{v}}}_j^{\top }\frac{\mathrm{d}{\mathbf{w}}_j}{\mathrm{d}t}=\nabla \widehat{\eta }\frac{\mathrm{d}{\mathbf{w}}_j}{\mathrm{d}t}=\frac{\mathrm{d}\widehat{\eta }\left({\mathbf{w}}_j\right)}{\mathrm{d}t}$$

(30)

We now turn to the right-hand side. Using (10), the numerical flux ${\tilde{\mathbf{F}}}_{j+1/2}$ explicitly given by (17), the equation (18) and the definition of ${\tilde{G}}_{j+1/2}$ (18), we get

$$\begin{array}{cc}\hfill {\widehat{\mathbf{v}}}_j^{\top }\left({\tilde{\mathbf{F}}}_{j+1/2}-{\tilde{\mathbf{F}}}_{j-1/2}\right)& ={\mathbf{v}}_j^{\top }\left({\tilde{\mathbf{F}}}_{j+1/2}-{\tilde{\mathbf{F}}}_{j-1/2}\right)-[\beta B_j/\rho ,0]\left({\tilde{\mathbf{F}}}_{j+1/2}-{\tilde{\mathbf{F}}}_{j-1/2}\right)\hfill \\ \hfill & ={\left\{\left\{{\mathbf{v}}^{\top }\right\}\right\}}_{j+1/2}{\tilde{\mathbf{F}}}_{j+1/2}-\frac{1}{2}{\left[\left[{\mathbf{v}}^{\top }\right]\right]}_{j+1/2}{\tilde{\mathbf{F}}}_{j+1/2}\hfill \\ \hfill & -{\left\{\left\{{\mathbf{v}}^{\top }\right\}\right\}}_{j-1/2}{\tilde{\mathbf{F}}}_{j-1/2}-\frac{1}{2}{\left[\left[{\mathbf{v}}^{\top }\right]\right]}_{j-1/2}{\tilde{\mathbf{F}}}_{j-1/2}\hfill \\ \hfill & -\frac{\beta B_j}{\rho }({\left\{\left\{AU\right\}\right\}}_{j+1/2}-{\left\{\left\{AU\right\}\right\}}_{j-1/2})\hfill \\ \hfill & ={\left\{\left\{{\mathbf{v}}^{\top }\right\}\right\}}_{j+1/2}{\tilde{\mathbf{F}}}_{j+1/2}-\frac{1}{2}{\left[\left[\psi \right]\right]}_{j+1/2}-{\left\{\left\{{\mathbf{v}}^{\top }\right\}\right\}}_{j-1/2}{\tilde{\mathbf{F}}}_{j-1/2}-\frac{1}{2}{\left[\left[\psi \right]\right]}_{j-1/2}\hfill \\ \hfill & -\frac{\beta }{2\rho }(B_j{A}_{j+1}{U}_{j+1}-B_j{A}_{j-1}{U}_{j-1})\hfill \\ \hfill & ={\left\{\left\{{\mathbf{v}}^{\top }\right\}\right\}}_{j+1/2}{\tilde{\mathbf{F}}}_{j+1/2}-{\left\{\left\{\psi \right\}\right\}}_{j+1/2}-{\left\{\left\{{\mathbf{v}}^{\top }\right\}\right\}}_{j-1/2}{\tilde{\mathbf{F}}}_{j-1/2}+{\left\{\left\{\psi \right\}\right\}}_{j-1/2}\hfill \\ \hfill & -\frac{\beta }{2\rho }(B_j{A}_{j+1}{U}_{j+1}-B_j{A}_{j-1}{U}_{j-1})\hfill \\ \hfill & =\left({\tilde{G}}_{j+1/2}-{\tilde{G}}_{j-1/2}\right)-\frac{\beta }{2\rho }(B_j{A}_{j+1}{U}_{j+1}-B_j{A}_{j-1}{U}_{j-1}).\hfill \end{array}$$

(31)

To deal with the factor (28) in parenthesis, we take into account the definition of the scaled entropy variables and (24) to obtain

$$\begin{array}{cc}\hfill & -\frac{1}{2}{\widehat{\mathbf{v}}}_j^{\top }(R_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}-R_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2})\hfill \\ \hfill & =-\frac{1}{2}{\left\{\left\{{\widehat{\mathbf{v}}}^{\top }\right\}\right\}}_{j+1/2}{R}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}+\frac{1}{2}{\left\{\left\{{\widehat{\mathbf{v}}}^{\top }\right\}\right\}}_{j-1/2}{R}_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2}\hfill \\ \hfill & +\frac{1}{4}{\left[\left[{\widehat{\mathbf{v}}}^{\top }\right]\right]}_{j+1/2}{R}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}+\frac{1}{4}{\left[\left[{\widehat{\mathbf{v}}}^{\top }\right]\right]}_{j-1/2}{R}_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2}\hfill \\ \hfill & =-({\tilde{\mathcal{D}}}_{j+1/2}-{\tilde{\mathcal{D}}}_{j-1/2})+\frac{1}{4}{\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}+\frac{1}{4}{\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2}\hfill \end{array}$$

(32)

Using the discretization of the source term bive by (21) we can rewrite (29) as

$$\begin{array}{cc}\hfill {\widehat{\mathbf{v}}}_j^{\top }{\mathbf{S}}_j& =({\mathbf{v}}^{\top }-[\beta \sqrt{A_0}/\rho ,0])\left[\begin{array}{c}0\\ {\displaystyle \frac{\beta (B_{j+1}-B_{j-1})}{2\rho \Delta x}}\end{array}\right]\hfill \\ \hfill & =\left[\frac{U_j^2}{2}+\frac{\beta \sqrt{A_j}}{\rho },A_j{U}_j\right]\left[\begin{array}{c}0\\ {\displaystyle \frac{\beta (B_{j+1}-B_{j-1})}{2\rho \Delta x}}\end{array}\right]\hfill \\ \hfill & =\frac{\beta }{2\rho \Delta x}(B_{j+1}{A}_j{U}_j-B_{j-1}{A}_j{U}_j).\hfill \end{array}$$

(33)

Inserting (30), (31), (32) and (33) into (26) an dusing the definition of ${\tilde{\mathcal{G}}}_{j+1/2}$ given by (22) we get

$$\begin{array}{cc}\hfill & \frac{\mathrm{d}\widehat{\eta }\left({\mathbf{w}}_j\right)}{\mathrm{d}t}=\hfill \\ \hfill & -\frac{1}{\Delta x}\left[\left({\tilde{G}}_{j+1/2}-{\displaystyle \frac{\beta \left(B_j{A}_{j+1}{U}_{j+1}+B_{j+1}{A}_j{U}_j\right)}{2\rho }}\right)-\left({\tilde{G}}_{j-1/2}-{\displaystyle \frac{\beta \left(B_{j-1}{A}_j{U}_j+B_j{A}_{j-1}{U}_{j-1}\right)}{2\rho }}\right)\right]\hfill \\ \hfill & +\frac{1}{\Delta x}({\tilde{\mathcal{D}}}_{j+1/2}-{\tilde{\mathcal{D}}}_{j-1/2})-\frac{1}{4\Delta x}{\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}-\frac{1}{4\Delta x}{\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2}\hfill \\ \hfill & =-\frac{1}{\Delta x}\left(({\tilde{\mathcal{G}}}_{j+1/2}-{\tilde{\mathcal{D}}}_{j+1/2})-({\tilde{\mathcal{G}}}_{j-1/2}-{\tilde{\mathcal{D}}}_{j-1/2})\right)\hfill \\ \hfill & -\frac{1}{4\Delta x}\left({\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}\right.+\left.{\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2}\right)\hfill \\ \hfill & =-\frac{1}{\Delta x}({\widehat{\mathcal{G}}}_{j+1/2}-{\widehat{\mathcal{G}}}_{j-1/2})-\frac{1}{4\Delta x}\left({\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}\right.+\left.{\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2}\right)\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill & \frac{\mathrm{d}\widehat{\eta }\left({\mathbf{w}}_j\right)}{\mathrm{d}t}+\frac{1}{\Delta x}({\widehat{\mathcal{G}}}_{j+1/2}-{\widehat{\mathcal{G}}}_{j-1/2})\hfill \\ \hfill & =-\frac{1}{4\Delta x}\left({\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j+1/2}{\Lambda }_{j+1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}\right.+\left.{\left[\left[{\mathbf{z}}^{\top }\right]\right]}_{j-1/2}{\Lambda }_{j-1/2}{\langle \langle \mathbf{z}\rangle \rangle }_{j-1/2}\right)\hfill \\ \hfill & \le 0\hfill \end{array}$$

where the last inequality is obtained by using the sign property. Summing up,

$$\frac{\mathrm{d}\widehat{\eta }\left({\mathbf{w}}_j\right)}{\mathrm{d}t}+\frac{1}{\Delta x}({\widehat{\mathcal{G}}}_{j+1/2}-{\widehat{\mathcal{G}}}_{j-1/2})\le 0,$$

which proves entropy stability.

We next show the second part of the theorem. From (25) it follows easily that

$${\tilde{\mathbf{F}}}_{j+1/2}={\left[0,{\displaystyle \frac{\beta }{\rho }}{\left\{\left\{\sqrt{A}\right\}\right\}}_{j+1/2}\right]}^{\top },{\widehat{\mathbf{v}}}_j={\left[{\displaystyle \frac{\beta C}{\rho }},0\right]}^{\top }.$$

The second equality implies that ${\widehat{\mathbf{v}}}_j$ is constant, thus ${\langle \langle \mathbf{z}\rangle \rangle }_{j+1/2}=\mathbf{0}$, and consequently the diffusion term in (16) vanishes. On account of this remark, we have ${\mathbf{F}}_{j+1/2}={\tilde{\mathbf{F}}}_{j+1/2}$. Hence,

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathbf{w}}_j\left(t\right)}{\mathrm{d}t}& =-\frac{\beta }{\rho \Delta x}\left({\tilde{\mathbf{F}}}_{j+1/2}-{\tilde{\mathbf{F}}}_{j-1/2}\right)+{\mathbf{S}}_j\hfill \\ \hfill & =-\frac{\beta }{\rho \Delta x}\left[\begin{array}{c}0\\ {\left\{\left\{\sqrt{A}\right\}\right\}}_{j+1/2}-{\left\{\left\{\sqrt{A}\right\}\right\}}_{j-1/2}\end{array}\right]+\frac{\beta }{2\rho \Delta x}\left[\begin{array}{c}0\\ B_{j+1}-B_{j-1}\end{array}\right]\hfill \\ \hfill & =-\frac{\beta }{2\rho \Delta x}\left[\begin{array}{c}0\\ ({\sqrt{A}}_{j+1}-B_{j+1})-({\sqrt{A}}_{j-1}-B_{j-1})\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}0\\ 0\end{array}\right],\hfill \end{array}$$

the last equality being a consequence of (25). □

Discretización temporal

To discretize the temporal derivative in (14), the explicit second-order strong stability preserving Runge-Kutta method SSPRK (Heun’s method) will be used. This method is described by the following steps

$$\begin{array}{cc}\hfill {\mathbf{w}}^{\left(1\right)}& ={\mathbf{w}}^n+\Delta t\mathcal{L}\left({\mathbf{w}}^n\right),\hfill \\ \hfill {\mathbf{w}}^{\left(2\right)}& ={\mathbf{w}}^{\left(1\right)}+\Delta t\mathcal{L}\left({\mathbf{w}}^{\left(1\right)}\right),\hfill \\ \hfill {\mathbf{w}}^{n+1}& =\frac{1}{2}\left({\mathbf{w}}^n+{\mathbf{w}}^{\left(2\right)}\right),\hfill \end{array}$$

where

$$\begin{array}{c}\hfill {\left[\mathcal{L}\left(\mathbf{w}\right)\right]}_j:=-\frac{1}{\Delta x}\left({\mathbf{F}}_{j+1/2}-{\mathbf{F}}_{j-1/2}\right)+{\mathbf{S}}_j,\end{array}$$

with ${\mathbf{F}}_{j+1/2}$ and ${\mathbf{S}}_j$ given by (16) and (21), respectively. In order to guarantee that the explicit scheme obtained to be stable, the time and space discretization steps must obey the CFL condition. For this purpose, the value of $\Delta t$ is computed adaptively for each step. More exactly, the solution ${\mathbf{w}}^{n+1}$ at $t_{n+1}=t_n+\Delta t$ is calculated from ${\mathbf{w}}^n$ by using the time step $\Delta t=CFL\ast \Delta x/{\alpha }_{max}^n$, where ${\alpha }_{max}$ is an estimate of the maximal characteristic velocity for $\mathbf{w}$. All numerical experiments are carried out with a CFL number of $0.5$.

3.1. Numerical tests

In this section, some examples are presented to validate that the scheme correctly captures the steady-state solutions of system (1) and to verify that it is entropy stable at discrete level. In all examples transmissive boundary conditions are imposed.

3.1.1. Ejemplo 1.

In the first example taken from [12], we consider three different configurations for the Riemann problem with initial data

$$\mathbf{w}(x,0)=\left\{\begin{array}{cc}{\mathbf{w}}_l={[A_l,U_l]}^{\top }\hfill & \mathrm{si}x\le x_d\hfill \\ {\mathbf{w}}_r={[A_r,U_r]}^{\top }\hfill & \mathrm{si}x>x_d,\hfill \end{array}\right.$$

with the spatial domain of length $L=0.2\mathrm{m}$ and $x_d=0.1\mathrm{m}$ the point where the discontinuity is located. The numerical solutions are computed on a mesh with $M=200$ uniform cells and for comparison purposes reference solutions are computed on a refined mesh with 3200 uniform cells. The other parameters are $A_0=3.14\times {10}^{-4}{\mathrm{m}}^2$, $\beta =3.31\times {10}^6\mathrm{Pa}/\mathrm{m}$ along with different configurations of the initial data to obtain different types of waves. One can observe in Figure 1 the two rarefactions obtained at the final time $t=0.009\mathrm{s}$ for the initial data ${\mathbf{w}}_l={[2A_0,-1]}^{\top }$, ${\mathbf{w}}_r={[2A_0,1]}^{\top }$. In Figure 2 it is observed that the scheme can capture the two shocks waves obtained at the time $t=0.012\mathrm{s}$ for the initial data ${\mathbf{w}}_l={[A_0,1]}^{\top }$, ${\mathbf{w}}_r={[2A_0,-1]}^{\top }$. The third configuration is ${\mathbf{w}}_l={[A_0,0]}^{\top }$, ${\mathbf{w}}_r={[2A_0,0]}^{\top }$ and the results at the time $t=0.012\mathrm{s}$ are displayed in Figure 3. In this case, the solutions develop a left discontinuity follows by a right rarefaction. Due to the presence of two types of waves we also present in Figure 4 the total entropy

$$\begin{array}{c}\hfill \mathcal{E}\left(t_n\right):=\Delta x\sum _{j=1}^M\eta \left({\mathbf{w}}_j\left(t_n\right)\right).\end{array}$$

(34)

As expected, this quantity decreases over time, that is, the scheme is entropy stable at the discrete level.

3.1.2. Ejemplo 2.

The purpose of this test problem taken from [13] is to check that the proposed scheme preserves the zero pressure man-at-eternal-rest steady state (13) at a discrete level. The configuration exhibits with no flow and includes a change in the section of the artery. This is the case for a dead man with an aneurism. The initial conditions are $U(x,0)=0$ and $A(x,0)=A_0\left(x\right)=\pi R_0^2\left(x\right)$, where

$$R_0\left(x\right)=\left\{\begin{array}{c}\tilde{R}\mathrm{si}x\in [0,x_1]\cup [x_4,L],\hfill \\ \tilde{R}+\frac{\Delta R\left[sin\left(\frac{x-x_1}{x_2-x_1}\pi -\frac{\pi }{2}\right)+1\right]}{2}\mathrm{si}x\in ]x_1,x_2[,\hfill \\ \tilde{R}+\Delta R\mathrm{si}x\in [x_2,x_3],\hfill \\ \tilde{R}+\frac{\Delta R\left[cos\left(\frac{x-x_3}{x_4-x_3}\pi \right)+1\right]}{2}\mathrm{si}x\in ]x_3,x_4[,\hfill \end{array}\right.$$

with $\tilde{R}=4\times {10}^{-3}\mathrm{m}$, $\Delta R=1.0\times {10}^{-3}\mathrm{m}$, $x_1=1.0\times {10}^{-2}\mathrm{m}$, $x_2=3.05\times {10}^{-2}\mathrm{m}$, $x_3=4.95\times {10}^{-2}\mathrm{m}$, $x_4=7.0\times {10}^{-2}\mathrm{m}$ and $L=0.14\mathrm{m}$. The computational domain is $[0,L]$ and $\beta ={\pi }^{-1}\times {10}^8\mathrm{Pa}/\mathrm{m}$. The numerical results computed on a mesh with 200 cells at time $t=5\mathrm{s}$ are plotted in Figure 5. In Figure 5 (left) can be observed that the area remains the same as the area at rest as expected. In Figure 5 (right) we compare the velocity obtained by using the well-balanced and entropy stable scheme proposed in this work (hereafter referred as WB-ES) and the velocity obtained by using the Lax-Friedrichs numerical flux combined with discretization of the source term given by (21). Notice that the WB-ES scheme preserves the steady-state exactly while Lax-Friedrichs produces unacceptable spurious oscillations. This kind of anomalies was also reported before in reference [5].

3.1.3. Ejemplo 3

This example (see [13]) corresponds to the case of a dead man with stenosis. Stenosis occurs when the artery narrows and it leads to reduced blood flow from the heart to the rest of the body. The cross-sectional area at rest is $A_0\left(x\right)=\pi R_0^2\left(x\right)$, where

$$R_0\left(x\right)=\left\{\begin{array}{cc}\tilde{R}+\Delta R\hfill & \mathrm{for}x\in [0,x_1]\cup [x_4,L],\hfill \\ \tilde{R}-{\displaystyle \frac{\Delta R}{2}}\left(sin\left({\displaystyle \frac{x-x_1}{x_2-x_1}}\pi -{\displaystyle \frac{\pi }{2}}\right)-1\right]\hfill & \mathrm{for}x\in (x_1,x_2),\hfill \\ \tilde{R}\hfill & \mathrm{for}x\in [x_2,x_3],\hfill \\ \tilde{R}-{\displaystyle \frac{\Delta R}{2}}\left(cos\left({\displaystyle \frac{x-x_3}{x_4-x_3}}\pi \right)-1\right)\hfill & \mathrm{for}x\in (x_3,x_4)\hfill \end{array}\right.$$

with $\tilde{R}=4\times {10}^{-3}\mathrm{m}$, $\Delta R=1.0\times {10}^{-3}\mathrm{m}$, an artery of length $L=0.14\mathrm{m}$, $x_1=9L/40$, $x_2=L/4$, $x_3=3L/40$, $x_4=31L/40$. The initial conditions are

$$A(x,0)={\left(C+{\pi }^{1/2}{R}_0\left(x\right)\right)}^2,U(x,0)=0,$$

where $C={10}^{-3}\mathrm{m}$.

Numerical solutions are computed at $t=1$ s. As observed in Figure 6 the proposed scheme preserves the steady-state (12) at the discrete level, i.e., the quantity $\sqrt{A(x,t)}-\sqrt{A_0\left(x\right)}$ remains constant and the velocity is null.

4. Conclusions

We present a numerical scheme to approximate the solutions of a blood flow model in arteries. It was shown analytically that the scheme is entropy stable and preserves the steady state solutions. These properties were also tested numerically with some standard examples. Although higher order schemes can be found in the literature related to the numerical approximation of blood flows, the value of the second order discretization of the source term proposed in this work lies in the fact that using linear combinations of them combined with high-order entropy stable fluxes it is possible to construct arbitrarily high order entropy stable schemes that preserve the steady states solutions.