Using a Separable Mathematical Entropy to Construct Entropy-Stable Schemes for a Reduced Blood Flow Model

: The aim of this paper is to derive a separable entropy for a one-dimensional reduced blood ﬂow model, which will be used to treat the symmetrizability of the model in full generality and for constructing entropy conservative ﬂuxes, which are one of the essential building blocks of designing entropy-stable schemes. Time discretization is conducted by implicit–explicit (IMEX) Runge–Kutta schemes, but solutions for nonlinear systems will not be required due to the particular form of the source term. To validate the numerical schemes obtained, some numerical tests are presented.


Introduction
This paper is concerned with a well-known reduced blood flow model described by a nonlinear hyperbolic system of conservation laws in one space dimension [1], which is used to model the flow of blood in axisymmetric vessels with compliant walls. The governing equations in terms of the vessel cross-sectional area A(x, t) and the mean blood velocity U(x, t) in the axial direction x are given by where ρ is the blood density, assumed to be constant for blood, which is essentially incompressible, C f is the skin friction coefficient and P = P(A) is the internal pressure, which is taken here as (see, for instance, [2]) Here A 0 is the vessel cross-sectional at rest and P 0 is the pressure when A = A 0 . Hereafter, it is assumed that β and A 0 are constants, but in reality, they may depend on x in the case of some pathologies.
The system of Equation (1) is known as the (A, U)-system. The velocity U is not a conservative quantity, in contrast to Q = AU within the so-called (A, Q)-system (see [3] for a complete discussion). For continuous solutions, the two formulations are equivalent. However, we will restrict our attention to the (A, U)-system, which can be written in the vector form as and the corresponding quasilinear form is where The Jacobian matrix H has two real eigenvalues, namely λ 1 = U − c and λ 2 = U + c with corresponding right eigenvectors Let us define the Shapiro number S h as [4] S h = U c .
The quantity S h is the analog of the Froude number for the shallow water equations. A state u is said to be subcritical if S h < 1, critical if S h = 1, and supercritical if S h > 1. The system is strictly hyperbolic in subcritical and supercritical regimes. In physiological conditions, blood flow is almost always subcritical. Nevertheless, very specific pathologies may lead to supercritical flows [5].
One-dimensional models are notably recognized to be computationally inexpensive in comparison with 3D models. In addition, 1D models are not suitable for describing blood flow in complicated morphological regions, and they can be coupled with 3D models to obtain a considerable reduction in the computational complexity [6].
The system (1) without friction, that is, with C f = 0, admits the following steady-state solution, known as the (non-zero pressure) man-at-eternal-rest steady state or dead-man equilibrium [3] (by analogy to the lake at rest in the shallow water equations): In particular, the (zero pressure) man-at-eternal-rest steady state is given by Let us also recall that a convex scalar function η = η(u) is the entropy for the system of conservation laws with associated entropy flux where is the vector of entropy variables. (η, G) is called an entropy pair for the conservation law (11). When η is strictly convex, the entropy variables v symmetrize the system (11) by making the change of variables u = u(v) [7], which puts the system into its equivalent symmetric form Note that the Jacobian of g(v) is the Hessian of the function The function ψ is called entropy potential and plays an important role in the construction of entropy conservative fluxes.
In Reference [6], an entropy pair for the inviscid (A, U) system (1) was derived, namely This entropy pair was used in [8] for constructing a well-balanced and entropy-stable scheme for the inviscid (A, U)-system but with A 0 depending on x.
The remainder of the paper is organized as follows: in Section 2, a separable entropy pair for the inviscid model (1) is derived and then employed to prove the symmetrizability of the last-mentioned model. In Section 3, another application of the last-mentioned entropy pair is obtained, namely, the construction of entropy conservative fluxes, which in turn, are used to obtain entropy-stable schemes by adding numerical diffusion [9]. After dealing with spatial discretization, we end this section with the treatment of the friction source term by using IMEX schemes. In Section 4, the obtained numerical schemes are validated with some benchmark tests taken from the literature. At last, some conclusions are drawn in Section 5.

Entropy Pair and Simmetrizability
It is well-known (see [10]) that symmetrizability is equivalent to the existence of a convex entropy function. Using the entropy given by (15), a result on the symmetrizability of the inviscid form of the system (1) was conducted in [11] to the subcritical case, that is, under the assumption U(x, t) < c(A(x, t)) for (x, t) ∈ [0, L] × [0, T]. We next follow the idea used in [12] to construct a separable entropy function for the inviscid (A, U) system, which allows us to obtain symmetrizability without the aforementioned assumption.
If η(u) is an additively separable function, that is, η(u) = e 1 (A) + e 2 (U), then the Hessian matrix of η, denoted by η uu (u), is a diagonal matrix. Now, it is fairly easy to see that η uu (u)H(u) is symmetric if Thus, e 1 (A) = −2β √ A and e 2 (U) = ρU 2 /2. Accordingly, is an entropy function for the inviscid (A, U) system, and the associated entropy flux is From the entropy pair (16) and (17), we obtain the entropy potential Next, we use the entropy function given above to demonstrate the symmetrizability of the inviscid (A, U) system (compare with the hypothesis in [11] (Lemma 2)).
Proof. From (16) it follows that the entropy variables are . Then the Jacobian is clearly symmetric positive definite and is symmetric, which proves the lemma.
In [11,13,14], the authors point out the fundamental importance of the numerical analysis with symmetrizability, in particular, to study the error estimates of the Runge-Kutta discontinuous Galerkin method.
Using the entropy function (16) and the eigenvector rescaling theorem [15] (Theorem 4), we obtain the lemma below, which provides a scaling of the eigenvectors R = [r 1 Then we have where u v is the symmetric positive definite change-of-variables matrix given by (19).

Proof. The result can be obtained directly by insertion.
This lemma will be used in the next section to construct a numerical diffusion operator.

Entropy Conservative and Entropy-Stable Numerical Schemes
For the homogeneous system (11) (until Section 3.2, we take C f = 0), a semi-discrete conservative scheme on a uniform spatial mesh x j = j∆x, j ∈ Z writes as where u j (t) denotes the numerical approximation of u(x j , t) and the numerical flux F j+1/2 is an approximation of the flux function at the cell interface j + 1/2.
The scheme (21) is called entropy stable with respect to the entropy pair (η, G) if it satisfies a discrete entropy inequality for some numerical entropy fluxG j+1/2 consistent with the entropy flux G. If equality holds in (22), then the scheme (21) is called entropy conservative. We focus on entropy-stable numerical fluxes of the form [9] j being the cell interface values of a reconstructed function v j (x) and D j+1/2 is a suitable numerical diffusion matrix, which will be specified later.
then the scheme is second-order accurate and entropy conservative.
Observe that the existence of an explicitly given entropy pair is an important ingredient in designing entropy conservative schemes. For the scalar case, the solution of (24) is unique. However, for systems of conservation laws, this is no longer true.
Tadmor also proposed (see [16] for more details) the following solution of (24): where v j+1/2 (ξ) denotes the straight line connecting v j and v j+1 , i.e., The flux (25) is sometimes called Averaged Energy Conservative (AEC for short) flux [17]. A straightforward computation of the integral in (25) for the (A, U) system yields the following components of the numerical fluxF j+1/2 : It is easy to check that the flux (26) is consistent. To do this, it is sufficient to recall that Another way of solving (24) is based on different paths in the phase space of the entropy variables [18]. This is described as follows: Let {r i } n i=1 be an arbitrary set of n linearly independent vectors, and let {l i } n i=1 be the corresponding orthogonal set. At an interface x j+1/2 , we define the paths Then the entropy conservative flux is given bỹ This flux is termed the Pathwise Energy Conservative (PEC) flux. In [17], it was reported that the computation of (27) may be numerically unstable.
A third strategy to construct the entropy conservative flux at the interface x j+1/2 was proposed in [17]. In this reference, an explicit solution of (24) for the shallow water equations was obtained by using the identity The strategy is called Explicit Energy Conservative (EEC) flux, and we employ this approach to the (A, u) blood flow model. Using (28), the jump of the entropy potential (18) Under the assumptions that A 0 and β are constants, the jump of entropy variables can be written as Writing down the desired flux componentwise asF j+1/2 = [F (1) j+1/2 ,F j+1/2 ] T , inserting the above two quantities into (24), equating jumps in U and P(A) and then solving the resulting system, we obtaiñ This flux is clearly consistent, very simple to code and computationally inexpensive.
The two-point entropy conservative fluxes obtained from (24) are only second-order accurate. However, high-order entropy conservative fluxes can be constructed by linear combinations of two-point entropy conservative fluxesF [19]. In this work, we use the fourth-order entropy conservative flux given bỹ To deal with the numerical diffusion part in (23), the diffusion matrix D j+1/2 is taken as [9] D j+1/2 = R j+1/2 Λ j+1/2 R T j+1/2 , where Λ j+1/2 = diag(|λ 1 |, |λ 2 |) is a Roe-type diagonal matrix, and R j+1/2 is the matrix of scaled right eigenvectors of the flux Jacobian H(u j+1/2 ) that is evaluated at the average state u j+1/2 := (u j + u j+1 )/2. To complete the description of (23), it only remains to perform a suitable reconstruction of the entropy variables v. Let v j (x) be a p−th reconstruction function of the entropy variables v. Denoting and defining the scaled entropy variables Instead of reconstructing the entropy variables, we reconstruct the scaled entropy variables such that the so-called sign property sign z l j j+1/2 = sign z l j j+1/2 will be satisfied. Here, z l j andz l j denote the l-th component of z j andz j , respectively. The advantage of using reconstruction procedures satisfying the sign property lies in the fact that they are entropy-stable (see Lemma 3.2 in [9]). On the other hand, the use of high-order nonoscillatory reconstruction is needed in order to avoid large oscillations around shocks; in particular, we use the fourth-order ENO. The crucial fact is that the ENO method satisfies the sign property [20], which in turn, guarantees that the reconstruction does not destroy entropy stability. The combination of entropy conservative fluxes and ENO reconstruction is termed TeCNO schemes [9].

Friction Source Term Discretization
Up to now, we have restricted our attention to the system (1)-(3) without friction, that is, the homogeneous case. For a treatment of the non-homogeneous case, the spatial semi-discretization of system (3) can be written as where is the spatial discretization of the convective term, and S(u) corresponds to the source term s(u). To solve (35), an implicit-explicit Runge-Kutta (IMEX-RK) method (see [21] and the references therein) will be used. Let us first recall that an IMEX-RK scheme consists of applying an implicit discretization to the source term and an explicit one to the convective part. An m-stage IMEX-RK scheme applied to system (35) takes the form a il S(u (l) ), i = 1, ..., m, This scheme is characterized by m × m matricesÃ = (ã il ) (withã il = 0 for l ≥ i ) and A = (a il ) that correspond to the explicit (ERK) and (diagonally) implicit (DIRK) parts of the method, respectively, whileb = (b 1 , . . . ,b m ) T and b = (b 1 , . . . , b m ) T are m-dimensional vectors of real coefficients.
In this work, we employ a second-order IMEX-RK scheme based on the Heun method coupled with the L-stable DIRK method, namely the scheme H-LDIRK3(2,2,2) [21,22], which is defined bỹ The particular value of γ guarantees that the implicit part is a third-order DIRK scheme with the best dampening properties (see [22,23] and references given there).
Applying the H-LDIRK3(2,2,2) scheme to (35) and taking into account the particular form of the source term s(u) = (0, −C f U/A) T , we obtain that the required computations to advance u n from time t n to t n+1 = t n + ∆t are given by where u = u n + ∆tã 21K1 + ∆ta 21 K 1 , It is worth pointing out that solutions for nonlinear systems are not required, reducing computational costs and making this method much simpler to implement. Let us also mention that (36) corresponds to the formula obtained in [3] (with a 11 = 1) by using the so-called semi-implicit treatment (SI).
Notice that when C f = 0, the method described above reduces to the Heun method, which is an explicit strong stability preserving Runge-Kutta method (SSPRK) [24]: (2) .
The time step ∆t is computed adaptatively in order for the CFL condition to be satisfied. We use the value ∆t = κ * ∆x/a, where a is an estimate of the maximal characteristic velocity. Here, the CFL number κ is taken as 0.5.

Numerical Results
In this section, we present some numerical examples by testing the numerical flux (23) with F given by (26) and (30) for the inviscid (A, U) system as well as for the viscous case. In the latter case, the scheme H-LDIRK3(2,2,2) is employed and compared with the SI scheme. For all the tests, the blood density is taken as ρ = 1060 kg/m 3 .

Example 1: The Ideal Tourniquet
This example is proposed by Delestre and Lagrée [3], and it resembles the dam break problem in shallow water equations. A tourniquet is applied, and we remove it instantaneously. When the tourniquet is removed, the blood flows from upstream to downstream in the vessel. The initial conditions are

Example 2: Wave Equation
We consider the system without friction and constant cross-section at rest to validate the capability of the proposed scheme to approximate the perturbed steady-state solutions. The cross section at rest is given by A 0 (x) = πR 2 0 with R 0 = 4 × 10 −3 m, and the initial conditions are with L = 0.16 m, x 1 = 0.2L, x 2 = 0.4L, x 3 = 0.6L and ε = 5 × 10 −3 . The computational domain is [0, L] and β = π −1 × 10 8 Pa/m. The exact solution (see [3] for more details) can be expressed as

Example 3: Wave Damping
In this last test [3], the viscous damping term is investigated in the linearized momentum equation. This is the analog of the Womersley problem [25], and a periodic signal at the inflow is considered with a constant section at rest. System (1), in terms of the variables (R, U) with the friction term, takes the form In the above model, the skin friction coefficient C f is 8πν, with ν being the viscosity of the blood. We consider this example on the computational domain [0, 3] subject to the given initial conditions U(x, 0) = 0 m/s, along with the following parameters: β = π −1 × 10 8 Pa/m, R 0 = 4 × 10 −3 m. The incoming discharge is Q b (0, t) = Q amp sin(ωt) m 3 /s, where Q amp = 3.45 × 10 −7 m 3 /s is the amplitude of the inflow discharge and ω = 2π/T pulse being T pulse = 0.5 s the time length of a pulse. A damping wave is obtained in the domain (see [3] for more details) with  Exact and numerical solutions of the damping of a discharge wave are depicted in Figure 4. Enlarged views are included to compare the performance of the H-LDIRK3(2,2,2) method and the SI method used for the discretization of the friction source term.

Conclusions
In this paper, we have presented a separable entropy pair for a 1D reduced blood flow model. We have used the proposed entropy pair to show the symmetrizability of the model and for constructing explicit, computationally inexpensive and easy-to-implement entropy conservative fluxes. After adding numerical diffusion to the entropy conservative fluxes (as recommended in [9]) along with a suitable reconstruction of the entropy variables, we obtained entropy-stable schemes. No significant differences were found between the numerical solutions obtained with the EEC flux and the AEC flux. However, small differences observed in the second numerical example were favorable for the EEC flux. Finally, we have solved the friction source term by IMEX-RK methods, with the advantage that it was unnecessary to solve nonlinear systems. In fact, the particular form of the friction term allowed us to derive a fully explicit scheme.

Conflicts of Interest:
The authors declare no conflict of interest.