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

Preprint

Article

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

Altmetrics

Downloads

Views

Comments

A peer-reviewed article of this preprint also exists.

This version is not peer-reviewed

Submitted:

23 October 2023

Posted:

25 October 2023

You are already at the latest version

Alerts

Abstract

In this work, an entropy-stable and well-balanced numerical scheme for a one-dimensional blood flow model is presented. Such scheme is obtained from an explicit entropy conservative flux along with a second order discretization of the source term by using centered finite differences. We prove that the scheme is entropy-stable and preserves steady-states solutions. In addition, some numerical examples are included to test the performance of the proposed scheme.

Keywords:

Subject: Computer Science and Mathematics - Applied Mathematics

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])

\{\begin{matrix} \frac{\partial A}{\partial t} + \frac{\partial (A U)}{\partial x} = 0 \\ \frac{\partial U}{\partial t} + U \frac{\partial U}{\partial x} + \frac{1}{ρ} \frac{\partial P}{\partial x} = \frac{F}{ρ A}, \end{matrix}

(1)

where

A (x, t) = π 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,

ρ

is the blood density, assumed to be constant and equal to

1060 kg / m^{3}

P = P (A)

is the internal pressure and

F (x, t)

is the friction force per unit length. The present work is restricted to the inviscid limit where

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_{e x t} + β (\sqrt{A} - \sqrt{A_{0}}),

(2)

where the external pressure

P_{ext}

, the coefficient

β

stands for the arterial stiffness and is supposedly constant and

A_{0} = A_{0} (x)

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 w}{\partial t} + \frac{\partial f (w)}{\partial x} = S (w)

(3)

where

\begin{matrix} w = [\begin{matrix} A \\ U \end{matrix}], f (w) = [\begin{matrix} A U \\ \frac{U^{2}}{2} + \frac{β \sqrt{A}}{ρ} \end{matrix}], S (w) = [\begin{matrix} 0 \\ \frac{β}{ρ} \frac{\partial \sqrt{A_{0}}}{\partial x} \end{matrix}] \end{matrix}

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

\frac{\partial w}{\partial t} + J (w) \frac{\partial w}{\partial x} = S (x, w),

where

J (w) = [\begin{matrix} U & A \\ c^{2} / A & U \end{matrix}]

is the Jacobian matrix of the flux function with eigenvalues

λ_{1} = U - c y λ_{2} = U + c

and the corresponding right eigenvectors are given by

r_{1} = {[- 1, c / A]}^{⊤} y r_{2} = {[1, c / A]}^{⊤} .

(4)

Here,

c = \sqrt{\frac{β \sqrt{A}}{2 ρ}}

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

\frac{\partial \hat{η} (w)}{\partial t} + \frac{\partial \hat{G} (w)}{\partial x} \leq 0,

(5)

where

(\hat{η}, \hat{G})

is the extended entropy pair introduced in [4]:

\hat{η} = η - \frac{β \sqrt{A_{0}}}{ρ} A, \hat{G} = G - \frac{β \sqrt{A_{0}}}{ρ} A U,

(6)

and

(η, G)

is the entropy pair

(η, G)

given by

η = \frac{1}{2} A U^{2} + \frac{2 β A^{3 / 2}}{3 ρ}, G = \frac{1}{2} A U^{3} + \frac{β U A^{3 / 2}}{ρ}

(7)

Using the entropy pair we may obtain the entropy variables

v (w) = \nabla η = {[\frac{U^{2}}{2} + \frac{β \sqrt{A}}{ρ}, A U]}^{⊤}

(8)

with the corresponding entropy potential

\begin{matrix} ψ (v) & : = v^{⊤} f (w (v)) - G (w (v)) = \frac{A U^{3}}{2} + \frac{β U A^{3 / 2}}{ρ} . \end{matrix}

(9)

In the extended case, we have

\hat{v} (w) = \nabla \hat{η} = v - {[β \sqrt{A_{0}} / ρ, 0]}^{⊤} .

(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

w (x)

which are time-independent. Thus, steady-state solutions satisfy the system

\frac{\partial f (w)}{\partial x} = S (w) .

(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{matrix} U = 0 y \sqrt{A} - \sqrt{A_{0}} = C = constant . \end{matrix}

(12)

When

C = 0

, the (zero pressure) man-at-eternal-rest steady state is obtained, namely

\begin{matrix} U = 0 and A = A_{0} . \end{matrix}

(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 Δ x

j \in Z

\frac{d w_{j} (t)}{d t} = - \frac{1}{Δ x} (F_{j + 1 / 2} - F_{j - 1 / 2}) + S_{j},

(14)

where

w_{j} (t)

is the cell average on

I_{j} = [x_{j - 1 / 2}, x_{j + 1 / 2})

F_{j + 1 / 2}

is the numerical flux associated with

x_{j + 1 / 2}

and

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

- f {(w)}_{x}

. The discretization of the time derivative will be discussed later.

The scheme (14) is entropy stable with respect to the entropy pair

(\hat{η}, \hat{G})

if it satisfies a discrete version of the entropy inequality (5), that is,

\frac{d}{d t} \hat{η} (w_{j} (t)) + \frac{1}{Δ x} ({\hat{G}}_{j + 1 / 2} - {\hat{G}}_{j - 1 / 2}) \leq 0

(15)

for some numerical entropy flux

{\hat{G}}_{j + 1 / 2}

consistent with the entropy flux

\hat{G}

. If equality holds in (15), then the scheme (14) is called entropy conservative.

From now on, we use

{[[a]]}_{j + 1 / 2} : = a_{j + 1} - a_{j}, {{a}}_{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]:

F_{j + 1 / 2} = {\tilde{F}}_{j + 1 / 2} - \frac{1}{2} R_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2},

(16)

where

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

${\tilde{F}}_{1, j + 1 / 2} = {{A U}}_{j + 1 / 2}, {\tilde{F}}_{2, j + 1 / 2} = \frac{1}{2} {{U^{2}}}_{j + 1 / 2} + \frac{β}{ρ} {{\sqrt{A}}}_{j + 1 / 2},$

(17)

with the corresponding numerical entropy flux

${\tilde{G}}_{j + 1 / 2} : = {{v^{⊤}}}_{j + 1 / 2} {\tilde{F}}_{j + 1 / 2} - {{ψ}}_{j + 1 / 2} .$

(18)

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

${[[v^{⊤}]]}_{j + 1 / 2} {\tilde{F}}_{j + 1 / 2} = {[[ψ]]}_{j + 1 / 2} .$

(19)
(ii).: $R_{j + 1 / 2}$ is the matrix of right eigenvectors of the Jacobian matrix $J (w_{j + 1 / 2})$ being evaluated at the average state $w_{j + 1 / 2} : = (w_{j} + w_{j + 1}) / 2$ , $Λ_{j + 1 / 2} = d i a g (| λ_{1} |, | λ_{2} |)$ is a Roe-type matrix and ${〈 〈 z 〉 〉}_{j + 1 / 2} = z_{j + 1 / 2}^{+} - z_{j + 1 / 2}^{-}$ $z_{j + 1 / 2}^{-}$ y $z_{j + 1 / 2}^{+}$ denoting, respectively, the left and right limiting values of the scaled entropy variables $z : = R_{j + 1 / 2}^{⊤} \hat{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]:

$sign ({〈 〈 z 〉 〉}_{j + 1 / 2}) = sign ({[[z]]}_{j + 1 / 2}),$

(20)

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

The discretization of the second component of the source term

B (x) : = \partial_{x} \sqrt{A_{0} (x)}

is performed using the simple and well-known second-order centered difference approximation

B^{'} (x) \approx \frac{B (x + Δ x) - B (x - Δ x)}{2 Δ x} .

Thus, the discretization of the source term yields

S_{j} = \frac{β}{ρ} [\begin{matrix} 0 \\ \frac{B_{j + 1} - B_{j - 1}}{2 Δ x}, \end{matrix}]

(21)

where

B_{j}

denotes the discrete approximation of

B (x_{j}) = \sqrt{A_{0} (x_{j})}

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 $\hat{η}$ is the entropy function given by (6) and the corresponding numerical entropy flux is

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

(22)

with

$\begin{matrix} {\tilde{G}}_{j + 1 / 2} & = {\tilde{G}}_{j + 1 / 2} - \frac{β (B_{j} A_{j + 1} U_{j + 1} + B_{j + 1} A_{j} U_{j})}{2 ρ}, \end{matrix}$

(23)

${\tilde{D}}_{j + 1 / 2} = \frac{1}{2} {{{\hat{v}}^{T}}}_{j + 1 / 2} R_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{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 w_{j}}{d t} = 0 \forall j .$

Proof.

Let us first prove entropy stability. Multiplying (14) by

{\hat{v}}_{j}^{⊤}

yields

\begin{matrix} {\hat{v}}_{j}^{⊤} \frac{d w_{j} (t)}{d t} & = - \frac{1}{Δ x} {\hat{v}}_{j}^{⊤} (F_{j + 1 / 2} - F_{j - 1 / 2}) + {\hat{v}}_{j}^{⊤} S_{j} \end{matrix}

(26)

\begin{matrix} = - \frac{1}{Δ x} {\hat{v}}_{j}^{⊤} ({\tilde{F}}_{j + 1 / 2} - {\tilde{F}}_{j - 1 / 2}) \end{matrix}

(27)

\begin{matrix} - \frac{1}{Δ x} {\hat{v}}_{j}^{⊤} (- \frac{1}{2} R_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} + \frac{1}{2} R_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2}) \end{matrix}

(28)

\begin{matrix} + {\hat{v}}_{j}^{⊤} S_{j} \end{matrix}

(29)

The left-hand side of (26) gives

{\hat{v}}_{j}^{⊤} \frac{d w_{j}}{d t} = \nabla \hat{η} \frac{d w_{j}}{d t} = \frac{d \hat{η} (w_{j})}{d t}

(30)

We now turn to the right-hand side. Using (10), the numerical flux

{\tilde{F}}_{j + 1 / 2}

explicitly given by (17), the equation (18) and the definition of

{\tilde{G}}_{j + 1 / 2}

(18), we get

\begin{matrix} {\hat{v}}_{j}^{⊤} ({\tilde{F}}_{j + 1 / 2} - {\tilde{F}}_{j - 1 / 2}) & = v_{j}^{⊤} ({\tilde{F}}_{j + 1 / 2} - {\tilde{F}}_{j - 1 / 2}) - [β B_{j} / ρ, 0] ({\tilde{F}}_{j + 1 / 2} - {\tilde{F}}_{j - 1 / 2}) \\ = {{v^{⊤}}}_{j + 1 / 2} {\tilde{F}}_{j + 1 / 2} - \frac{1}{2} {[[v^{⊤}]]}_{j + 1 / 2} {\tilde{F}}_{j + 1 / 2} \\ - {{v^{⊤}}}_{j - 1 / 2} {\tilde{F}}_{j - 1 / 2} - \frac{1}{2} {[[v^{⊤}]]}_{j - 1 / 2} {\tilde{F}}_{j - 1 / 2} \\ - \frac{β B_{j}}{ρ} ({{A U}}_{j + 1 / 2} - {{A U}}_{j - 1 / 2}) \\ = {{v^{⊤}}}_{j + 1 / 2} {\tilde{F}}_{j + 1 / 2} - \frac{1}{2} {[[ψ]]}_{j + 1 / 2} - {{v^{⊤}}}_{j - 1 / 2} {\tilde{F}}_{j - 1 / 2} - \frac{1}{2} {[[ψ]]}_{j - 1 / 2} \\ - \frac{β}{2 ρ} (B_{j} A_{j + 1} U_{j + 1} - B_{j} A_{j - 1} U_{j - 1}) \\ = {{v^{⊤}}}_{j + 1 / 2} {\tilde{F}}_{j + 1 / 2} - {{ψ}}_{j + 1 / 2} - {{v^{⊤}}}_{j - 1 / 2} {\tilde{F}}_{j - 1 / 2} + {{ψ}}_{j - 1 / 2} \\ - \frac{β}{2 ρ} (B_{j} A_{j + 1} U_{j + 1} - B_{j} A_{j - 1} U_{j - 1}) \\ = ({\tilde{G}}_{j + 1 / 2} - {\tilde{G}}_{j - 1 / 2}) - \frac{β}{2 ρ} (B_{j} A_{j + 1} U_{j + 1} - B_{j} A_{j - 1} U_{j - 1}) . \end{matrix}

(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{matrix} - \frac{1}{2} {\hat{v}}_{j}^{⊤} (R_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} - R_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2}) \\ = - \frac{1}{2} {{{\hat{v}}^{⊤}}}_{j + 1 / 2} R_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} + \frac{1}{2} {{{\hat{v}}^{⊤}}}_{j - 1 / 2} R_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2} \\ + \frac{1}{4} {[[{\hat{v}}^{⊤}]]}_{j + 1 / 2} R_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} + \frac{1}{4} {[[{\hat{v}}^{⊤}]]}_{j - 1 / 2} R_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2} \\ = - ({\tilde{D}}_{j + 1 / 2} - {\tilde{D}}_{j - 1 / 2}) + \frac{1}{4} {[[z^{⊤}]]}_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} + \frac{1}{4} {[[z^{⊤}]]}_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2} \end{matrix}

(32)

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

\begin{matrix} {\hat{v}}_{j}^{⊤} S_{j} & = (v^{⊤} - [β \sqrt{A_{0}} / ρ, 0]) [\begin{matrix} 0 \\ \frac{β (B_{j + 1} - B_{j - 1})}{2 ρ Δ x} \end{matrix}] \\ = [\frac{U_{j}^{2}}{2} + \frac{β \sqrt{A_{j}}}{ρ}, A_{j} U_{j}] [\begin{matrix} 0 \\ \frac{β (B_{j + 1} - B_{j - 1})}{2 ρ Δ x} \end{matrix}] \\ = \frac{β}{2 ρ Δ x} (B_{j + 1} A_{j} U_{j} - B_{j - 1} A_{j} U_{j}) . \end{matrix}

(33)

Inserting (30), (31), (32) and (33) into (26) an dusing the definition of

{\tilde{G}}_{j + 1 / 2}

given by (22) we get

\begin{matrix} \frac{d \hat{η} (w_{j})}{d t} = \\ - \frac{1}{Δ x} [({\tilde{G}}_{j + 1 / 2} - \frac{β (B_{j} A_{j + 1} U_{j + 1} + B_{j + 1} A_{j} U_{j})}{2 ρ}) - ({\tilde{G}}_{j - 1 / 2} - \frac{β (B_{j - 1} A_{j} U_{j} + B_{j} A_{j - 1} U_{j - 1})}{2 ρ})] \\ + \frac{1}{Δ x} ({\tilde{D}}_{j + 1 / 2} - {\tilde{D}}_{j - 1 / 2}) - \frac{1}{4 Δ x} {[[z^{⊤}]]}_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} - \frac{1}{4 Δ x} {[[z^{⊤}]]}_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2} \\ = - \frac{1}{Δ x} (({\tilde{G}}_{j + 1 / 2} - {\tilde{D}}_{j + 1 / 2}) - ({\tilde{G}}_{j - 1 / 2} - {\tilde{D}}_{j - 1 / 2})) \\ - \frac{1}{4 Δ x} ({[[z^{⊤}]]}_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} + {[[z^{⊤}]]}_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2}) \\ = - \frac{1}{Δ x} ({\hat{G}}_{j + 1 / 2} - {\hat{G}}_{j - 1 / 2}) - \frac{1}{4 Δ x} ({[[z^{⊤}]]}_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} + {[[z^{⊤}]]}_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2}) \end{matrix}

Therefore,

\begin{matrix} \frac{d \hat{η} (w_{j})}{d t} + \frac{1}{Δ x} ({\hat{G}}_{j + 1 / 2} - {\hat{G}}_{j - 1 / 2}) \\ = - \frac{1}{4 Δ x} ({[[z^{⊤}]]}_{j + 1 / 2} Λ_{j + 1 / 2} {〈 〈 z 〉 〉}_{j + 1 / 2} + {[[z^{⊤}]]}_{j - 1 / 2} Λ_{j - 1 / 2} {〈 〈 z 〉 〉}_{j - 1 / 2}) \\ \leq 0 \end{matrix}

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

\frac{d \hat{η} (w_{j})}{d t} + \frac{1}{Δ x} ({\hat{G}}_{j + 1 / 2} - {\hat{G}}_{j - 1 / 2}) \leq 0,

which proves entropy stability.

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

{\tilde{F}}_{j + 1 / 2} = {[0, \frac{β}{ρ} {{\sqrt{A}}}_{j + 1 / 2}]}^{⊤}, {\hat{v}}_{j} = {[\frac{β C}{ρ}, 0]}^{⊤} .

The second equality implies that

{\hat{v}}_{j}

is constant, thus

{〈 〈 z 〉 〉}_{j + 1 / 2} = 0

, and consequently the diffusion term in (16) vanishes. On account of this remark, we have

F_{j + 1 / 2} = {\tilde{F}}_{j + 1 / 2}

. Hence,

\begin{matrix} \frac{d w_{j} (t)}{d t} & = - \frac{β}{ρ Δ x} ({\tilde{F}}_{j + 1 / 2} - {\tilde{F}}_{j - 1 / 2}) + S_{j} \\ = - \frac{β}{ρ Δ x} [\begin{matrix} 0 \\ {{\sqrt{A}}}_{j + 1 / 2} - {{\sqrt{A}}}_{j - 1 / 2} \end{matrix}] + \frac{β}{2 ρ Δ x} [\begin{matrix} 0 \\ B_{j + 1} - B_{j - 1} \end{matrix}] \\ = - \frac{β}{2 ρ Δ x} [\begin{matrix} 0 \\ ({\sqrt{A}}_{j + 1} - B_{j + 1}) - ({\sqrt{A}}_{j - 1} - B_{j - 1}) \end{matrix}] \\ = [\begin{matrix} 0 \\ 0 \end{matrix}], \end{matrix}

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{matrix} w^{(1)} & = w^{n} + Δ t L (w^{n}), \\ w^{(2)} & = w^{(1)} + Δ t L (w^{(1)}), \\ w^{n + 1} & = \frac{1}{2} (w^{n} + w^{(2)}), \end{matrix}

where

\begin{matrix} {[L (w)]}_{j} : = - \frac{1}{Δ x} (F_{j + 1 / 2} - F_{j - 1 / 2}) + S_{j}, \end{matrix}

with

F_{j + 1 / 2}

and

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

Δ t

is computed adaptively for each step. More exactly, the solution

w^{n + 1}

t_{n + 1} = t_{n} + Δ t

is calculated from

w^{n}

by using the time step

Δ t = C F L * Δ x / α_{max}^{n}

, where

α_{max}

is an estimate of the maximal characteristic velocity for

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

w (x, 0) = \{\begin{matrix} w_{l} = {[A_{l}, U_{l}]}^{⊤} & si x \leq x_{d} \\ w_{r} = {[A_{r}, U_{r}]}^{⊤} & si x > x_{d}, \end{matrix}

with the spatial domain of length

L = 0.2 m

and

x_{d} = 0.1 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} m^{2}

β = 3.31 \times 10^{6} Pa / 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 s

for the initial data

w_{l} = {[2 A_{0}, - 1]}^{⊤}

w_{r} = {[2 A_{0}, 1]}^{⊤}

. In Figure 2 it is observed that the scheme can capture the two shocks waves obtained at the time

t = 0.012 s

for the initial data

w_{l} = {[A_{0}, 1]}^{⊤}

w_{r} = {[2 A_{0}, - 1]}^{⊤}

. The third configuration is

w_{l} = {[A_{0}, 0]}^{⊤}

w_{r} = {[2 A_{0}, 0]}^{⊤}

and the results at the time

t = 0.012 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{matrix} E (t_{n}) : = Δ x \sum_{j = 1}^{M} η (w_{j} (t_{n})) . \end{matrix}

(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} (x) = π R_{0}^{2} (x)

, where

R_{0} (x) = \{\begin{matrix} \tilde{R} si x \in [0, x_{1}] \cup [x_{4}, L], \\ \tilde{R} + \frac{Δ R [sin (\frac{x - x_{1}}{x_{2} - x_{1}} π - \frac{π}{2}) + 1]}{2} si x \in] x_{1}, x_{2} [, \\ \tilde{R} + Δ R si x \in [x_{2}, x_{3}], \\ \tilde{R} + \frac{Δ R [cos (\frac{x - x_{3}}{x_{4} - x_{3}} π) + 1]}{2} si x \in] x_{3}, x_{4} [, \end{matrix}

with

\tilde{R} = 4 \times 10^{- 3} m

Δ R = 1.0 \times 10^{- 3} m

x_{1} = 1.0 \times 10^{- 2} m

x_{2} = 3.05 \times 10^{- 2} m

x_{3} = 4.95 \times 10^{- 2} m

x_{4} = 7.0 \times 10^{- 2} m

and

L = 0.14 m

. The computational domain is

[0, L]

and

β = π^{- 1} \times 10^{8} Pa / m

. The numerical results computed on a mesh with 200 cells at time

t = 5 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} (x) = π R_{0}^{2} (x)

, where

R_{0} (x) = \{\begin{matrix} \tilde{R} + Δ R & for x \in [0, x_{1}] \cup [x_{4}, L], \\ \tilde{R} - \frac{Δ R}{2} (sin (\frac{x - x_{1}}{x_{2} - x_{1}} π - \frac{π}{2}) - 1] & for x \in (x_{1}, x_{2}), \\ \tilde{R} & for x \in [x_{2}, x_{3}], \\ \tilde{R} - \frac{Δ R}{2} (cos (\frac{x - x_{3}}{x_{4} - x_{3}} π) - 1) & for x \in (x_{3}, x_{4}) \end{matrix}

with

\tilde{R} = 4 \times 10^{- 3} m

Δ R = 1.0 \times 10^{- 3} m

, an artery of length

L = 0.14 m

x_{1} = 9 L / 40

x_{2} = L / 4

x_{3} = 3 L / 40

x_{4} = 31 L / 40

. The initial conditions are

A (x, 0) = {(C + π^{1 / 2} R_{0} (x))}^{2}, U (x, 0) = 0,

where

C = 10^{- 3} 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} (x)}

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.

Author Contributions

Conceptualization, S.V., C.V. and J.B.; Formal Analysis, C.V., Investigation, S.V. and C.V; Software, Validation & Supervision, C.V. and J.B.; Writing-original draft, C.V.; Writing-Review & Editing, S.V. and C.V. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

Sherwin, S. J.; Franke, V.; Peiró, J.; Parker, K. One-dimensional modelling of a vascular network in space–time variables. J. Eng. Math. 2003, 47(3), 217–250. [Google Scholar] [CrossRef]
Formaggia, L.; Nobile, F.; Quarteroni, A.; Veneziani, A. Multiscale modelling of the circulatory system: a preliminar analysis. Comput. Vis. Sci. 1999, 2, 75–83. [Google Scholar] [CrossRef]
Dafermos, C. Hyperbolic conservation laws in continuum physics; Springer: Berlin, 2000. [Google Scholar]
Guigo, A. R.; Delestre, O.; Fullana, J.M.; Lagrée, P.Y. Low-Shapiro hydrostatic reconstruction tecnique for blood flow simulation in large arteries with varying geometrical and mechanical properties. J. Comput. Phys. 2017, 331, 108–136. [Google Scholar] [CrossRef]
Delestre, O.; Lagrée, P. Y. A ’well-balanced’ finite volume scheme for blood flow simulation. Int. J. Numer. Methods Fluids. 2013, 72, 177–205. [Google Scholar] [CrossRef]
Tadmor, E. The numerical viscosity of entropy stable schemes for systems of conservation laws, I. Math. Comput. 1987, 49, 91–103. [Google Scholar] [CrossRef]
Fjordholm, U.S.; Mishra, S.; Tadmor, E. Arbitrary high-order essentially non-oscillatory entropy stable schemes for systems of conservation laws. SIAM J. Numer. Anal. 2012, 50, 544–573. [Google Scholar] [CrossRef]
Bürger, R.; Valbuena, S.; Vega, C. A well-balanced and entropy stable scheme for a reduced blood flow model. Numer. Meth. Part Differ. Equ. 2023, 39(3), 2491–2509. [Google Scholar] [CrossRef]
Fjordholm, U.S.; Mishra, S.; Tadmor, E. ENO reconstruction and ENO interpolation are stable. Foundations of Computational Mathemathics. 2013, 13, 139–159. [Google Scholar] [CrossRef]
Fjordholm, U.S.; Mishra, S.; Tadmor, E. Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography. J. Comput. Phys. 2011, 230, 5587–5609. [Google Scholar] [CrossRef]
Wu, K.; Shu, C.W. Entropy Symmetrization and High-Order Accurate Entropy Stable Numerical Schemes for Relativistic MHD Equations. SIAM J. Sci. Comput. 2020, 42(4), A2230–A2261. [Google Scholar] [CrossRef]
Guitti, B.; Berton, C.; Hoang, M.; Toro, E. A fully well-balanced scheme for the 1D blood flow equations with friction source term. J. Comput. Phys. 2020, 421, 1–33. [Google Scholar] [CrossRef]
Britton, J.; Xing, Y. Well-balanced discontinuous Galerkin methods for the one-dimensional blood flow through arteries model with man-at-eternal-rest and living-man equilibria. Computers & Fluids. 2020, 203, 1–32. [Google Scholar] [CrossRef]

Figure 1. Rarefaction waves: area and velocity at

t = 0.009 s

for the reference and approximate solutions.

Figure 1. Rarefaction waves: area and velocity at

t = 0.009 s

for the reference and approximate solutions.

Figure 2. Shock waves: area and velocity at

t = 0.012 s

for the references and approximate solutions.

Figure 2. Shock waves: area and velocity at

t = 0.012 s

for the references and approximate solutions.

Figure 3. Shock and rarefaction waves: area and velocity at

t = 0.012 s

for the references and approximate solutions.

Figure 3. Shock and rarefaction waves: area and velocity at

t = 0.012 s

for the references and approximate solutions.

Figure 4. Entropy vs time on a mesh with 200 uniform cells.

Figure 5. Numerical solutions of the zero pressure man-at-eternal-rest problem on a mesh with 200 cells at time

t = 5 s

. Area and area at rest (left) and comparison between solutions obtained with WB-ES and Lax-Friedrichs schemes (right).

Figure 5. Numerical solutions of the zero pressure man-at-eternal-rest problem on a mesh with 200 cells at time

t = 5 s

. Area and area at rest (left) and comparison between solutions obtained with WB-ES and Lax-Friedrichs schemes (right).

Figure 6. numerical solutions of the non-zero pressure man-at-eternal-rest problem on a mesh with 200 cells at simulated time

t = 1 s

. Area minus area at rest (left) and velocity (right).

Figure 6. numerical solutions of the non-zero pressure man-at-eternal-rest problem on a mesh with 200 cells at simulated time

t = 1 s

. Area minus area at rest (left) and velocity (right).

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

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

Abstract

1. Introduction

2. Numerical method

3. Theoretical results and numerical experiments

Discretización temporal

3.1. Numerical tests

3.1.1. Ejemplo 1.

3.1.2. Ejemplo 2.

3.1.3. Ejemplo 3

4. Conclusions

Author Contributions

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe