On the Derivation of a New One-Dimensional Model for Blood Flows and Its Numerical Approximation

Yolhan Mannes¹, Mehmet Ersoy¹, Ömer Faruk Eker², Aimed Ajroud³

¹Institut de Mathématiques de Toulon (IMATH), Université de Toulon, La Garde, France.
²CHU de Bordeaux, Université Claude Bernard Lyon 1, Lyon, France.
³Ecole d’Ingénieurs Seatech, Université de Toulon, Toulon, France.

DOI: 10.4236/jamp.2025.1310198                

Abstract

We propose a new section-averaged one-dimensional model for blood flows in deformable arteries. The model is derived from the three-dimensional Navier-Stokes equations, written in cylindrical coordinates, under the “thin-artery” assumption (similar to the “shallow-water” assumption for free surface models). The blood flow/artery interaction is taken into account through suitable boundary conditions. The obtained equations enter the scope of the non-linear convection-diffusion problems. We show that the resulting model is energetically consistent. The proposed model extends most extant models by adding more scope, depending on an additional viscous term. We compare both models computationally based on an Incomplete Interior Penalty Galerkin (IIPG) method for the parabolic part, and on a Runge Kutta Discontinuous Galerkin (RKDG) method for the hyperbolic part. The time discretization explicit/implicit is based on the well-known Additive Runge-Kutta (ARK) method. Moreover, through a suitable change of variables, by construction, we show that the numerical scheme is well-balanced, i.e., it preserves exactly still-steady state solutions. To end, we numerically investigate its efficiency through several test cases with a confrontation to an exact solution.

Keywords

Blood Flow, Convection-Diffusion Problems, Model Reduction, Asymptotic Analysis, Well-Balanced Scheme, DG, IIPG, RKDG, ARK

1. Introduction

Modeling the cardiovascular system in arteries holds a central place in medical science, particularly in connection with cardiovascular diseases such as coronary heart disease, stroke, peripheral artery disease, aneurysms, and among others. This is especially important today in understanding and forecasting the impact of developed countries’ way of life on people’s healthcare (around 30% of cardiovascular disease deaths are developed from countries). Therefore, it is of major interest to develop an accurate mathematical model. In this part, we derive a new one-dimensional model for blood flow in the spirit of the work by [1]-[3].

The dynamic of such flow is mainly influenced by the fluid-structure interaction with the artery wall. The forecast to predict the motion of blood through the artery is a difficult task to which substantial effort has been devoted [4]-[11].

One of the most widely used models to describe the motion of blood through the artery is the one-dimensional (1D) Blood Flow equation derived, for instance, in [4] [8] [9] [12]-[14]. This classical model takes the form of a hyperbolic system of Partial Differential Equations (PDE) describing the conservation laws linking the wall elasticity to the fluid dynamics. This model is derived from an ansatz for the velocity profile (3) and reads,

$\{\begin{array}{l}{\partial }_{t}A+{\partial }_{x}Q=0,\\ {\partial }_{t}Q+{\partial }_x\left(\alpha \frac{Q^2}{A}+\frac{1}{\rho }AP\left(A,x\right)\right)=\frac{1}{\rho }P\left(A,x\right){\partial }_{x}A-K\frac{Q}{A},\end{array}$ (1)

where the unknowns $A\left(t,x\right)$ stand for artery’s section area (assumed to be cylindrical), $Q\left(t,x\right)=A\left(t,x\right)\overline{u_x}\left(t,x\right)$ is the flow rate and $\overline{u_x}$ is the mean speed over the artery’s section (see [4] [9] for further details). The function $P\left(A,x\right)$ denotes the pressure of blood at the wall and reads,

$P\left(A,x\right)=b\left(x\right)\frac{\sqrt{A}-\sqrt{A_0}}{A_0},$ (2)

where $b$ encompasses the elastic behavior of the artery, i.e., $b\left(x\right)=\frac{E\left(x\right)}{1-\xi }h\sqrt{\pi }$ , where $E$ is the Young’s modulus, $\xi$ is the Poisson ratio, and $h$ is the wall thickness. The velocity profile is given by

$u_x\left(t,r,x\right)=\frac{Q\left(t,x\right)}{A\left(t,x\right)}\frac{\gamma +2}{\gamma }\left[1-{\left(\frac{r}{R\left(t,x\right)}\right)}^{\gamma }\right],$ (3)

where $\gamma$ is an integer (often set to 9, see for instance [9] [15]). The friction term $K$ is defined as a function of $\gamma$ by $K=2\pi \nu \left(\gamma +2\right)$ , where $\nu$ is the kinematic viscosity of the blood. In Equation (1), the Bousinessq coefficient is given by $\alpha =\frac{\gamma +2}{\gamma +1}$ (see [9] [10] [15]).

Our main goal in this section is to derive, starting from the incompressible Navier-Stokes equations with suitable Navier boundary conditions¹ to account for the dynamic of the artery wall and the friction generated, a model akin to (1) via section-averaging (see [1]) under a thin-artery assumption. The averaged model that we obtain differs from the model (1) in that it has an additional diffusion term, and the velocity profile is a direct consequence of the performed asymptotic approximation (as in [3]). More precisely, one has

$u_x\left(t,r,x\right)=\frac{Q}{A}\left[1-\frac{kR\left(t,x\right)}{4\nu }\left(1-2{\left(\frac{r}{R\left(t,x\right)}\right)}^2\right)\right]$ (4)

where $k<0$ is a friction term.

The new model we propose is

$\{\begin{array}{l}{\partial }_{t}A+{\partial }_{x}Q=0,\\ {\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{A}+\frac{1}{\rho }AP\left(A,x\right)\right)-{\partial }_x\left(3\nu A{\partial }_x\left(\frac{Q}{A}\right)\right)=\frac{1}{\rho }P\left(A,x\right){\partial }_{x}A+\frac{2\pi Rk}{1-\frac{Rk}{4\nu }}\frac{Q}{A},\end{array}$ (5)

where $P$ is defined by (2). We show that system (5) possesses a mathematical entropy given by,

$E\left(t,x\right)=\widehat{E}\left(A,\overline{u_x},x\right)=\frac{A{\overline{u_x}}^2}{2}+\frac{1}{\rho }AP\left(A,x\right)-\frac{\beta \left(x\right)}{3\rho A_0\left(x\right)}{A}^{\frac{3}{2}}$

where $\overline{u_x}$ stands here for the mean speed over the artery’s section of the quantity (4).

We show that this entropy satisfies the following entropy relation for smooth solutions:

${\partial }_{t}E+{\partial }_x\left(\left(E+\frac{\beta \left(x\right)}{3\rho A_0\left(x\right)}{A}^{\frac{3}{2}}\right)\overline{u_x}\right)={\partial }_x\left(3\nu A{\partial }_x\left(\frac{Q}{A}\right)\right)\overline{u_x}+\frac{2\pi Rk}{1-\frac{Rk}{4\nu }}{\overline{u_x}}^2.$

Under the assumptions, $Q\left(t,0\right)=Q\left(t,L\right)=0$ at the inlet and outlet respectively, we show that the global energy decreases:

${\partial }_t\left({\displaystyle {\int }_0^{L}E}\right)=-3\nu {\displaystyle {\int }_0^L\left(A{\left({\partial }_x\overline{u_x}\right)}^2\right)}+\frac{2\pi Rk}{1-\frac{Rk}{4\nu }}{\displaystyle {\int }_0^L{\overline{u_x}}^2}<0.$

We outline the rest of the article as follows: in Section 2, as our starting point, we present the Navier-Stokes equations and the boundary conditions, including friction and the wall law. We derive the section-averaged one-dimensional equations. In Section 3, we use a Discontinuous Galerkin (DG) method from [16] called the Incomplete Interior Penalty Galerkin (IIPG) method together with the Runge-Kutta Discontinuous Galerkin method (RKDG) from [17]. In Section 4, we use additive implicit/explicit Runge-Kutta for time integration for the parabolic and hyperbolic parts, respectively. We provide extensive numerical testing in Section 5 of the resulting code².

2. Derivation of a New One-Dimensional Model for Blood Flow in Arteries

In this section, we present the full derivation of the new one-dimensional model for blood flow.

2.1. Navier-Stokes Equations in Cylindrical Coordinates

Our aim is to construct a mathematical model for blood flow in an artery consistent with the phenomena that can affect its motion. To this purpose, we propose a model reduction of the three-dimensional Navier-Stokes equations leading to a new one-dimensional model following the technique in [2] [3] [18]. We study the case of an axisymmetric artery (with a circular cross-section) and consider suitable boundary conditions to account for artery wall radial deformation and friction. More precisely, the displacements are only in the radial direction, for all angles $\theta$ .

Figure 1. Artery shape (DALL-E generated picture).

We start in Section 2.1.1 by reviewing the Navier-Stokes equations in cylindrical coordinates, describing the physics with the artery wall boundary. We then introduce the boundary condition at the wall in Section 2.1.2.

2.1.1. Geometric Set-Up and the Three-Dimensional Navier-Stokes Equations in Cylindrical Coordinates

With reference to Figure 1, we consider an incompressible fluid moving in the time-space domain

$\text{Ω}=\left\{\left(t,x,y,z\right)\in \left[0,T\right]\times {ℝ}^3|\sqrt{y^2+z^2}\le R\left(t,x\right),x\in \left[0,L\right]\right\},$ (6)

where $R$ is the radius of a cross-section of the artery (supposed circular), $L$ is its length, and $T>0$ is an arbitrary final time.

We assume that the viscous flow is axisymmetric (following [6] [9] [15]) and its velocity $u$ satisfies, on the domain Ω, the three-dimensional incompressible Navier-Stokes equations

$\{\begin{array}{l}\text{div}u=0,\\ {\partial }_{t}u+\text{div}\left(u\otimes u\right)+\nabla p-\text{div}\sigma =0,\end{array}$ (7)

where $u=u_{x}i+u_{y}j+u_{z}k$ with $\left(i,j,k\right)$ the cartesian basis, $p=\frac{P}{\rho }$ where $P$ is the pressure and $\rho$ , the density of the fluid. Finally, the stress tensor is

$\sigma =\nu \left(\nabla u+{\left(\nabla u\right)}^t\right),$

where $\nu$ is the kinematic viscosity.

Consider the following change of variable,

$e_r\left(\theta \right)=\cos \theta i+\sin \theta j,e_{\theta }\left(\theta \right)=-\sin \theta i+\cos \theta j,$ (8)

then the domain Ω (6) can be expressed as

$\Omega =\left\{\left(t,xi+r{e}_r\left(\theta \right)\right)\in \left[0,T\right]\times {ℝ}^3;r\in \left[0,R\left(t,x\right)\right],\theta \in \left[0,2\pi \right],x\in \left[0,L\right]\right\}.$ (9)

The velocity $u$ , thanks to the axisymmetry assumption, is written $u=u_r{e}_r+u_{x}i$ where $u_r$ is the radial speed and $u_x$ is called, later on, the axial speed.

Remark. As done in [9], thanks to the axisymmetry assumption, the velocity $u$ has no angular component. For the same reason, radial and axial components of $u$ don’t depend on $\theta$ .

The coordinate transformation reads,

$M\left(r,\theta ,x\right)=xi+r{e}_r\left(\theta \right),$

leading to the following Jacobian matrix,

$A^{-1}={\left(\begin{array}{ccc}1& 0& 0\\ 0& r& 0\\ 0& 0& 1\end{array}\right)}_{e_r,e_{\theta },i}\text{and its inverse}A=\frac{1}{J}{\left(\begin{array}{ccc}r& 0& 0\\ 0& 1& 0\\ 0& 0& r\end{array}\right)}_{e_r,e_{\theta },i},$

where $J=r=\text{det}\left(A^{-1}\right)$ . Coordinate change (8) together with axisymmetry assumption led to

$\nabla p={\left(\begin{array}{c}{\partial }_{r}p\\ \frac{{\partial }_{\theta }p}{r}\\ {\partial }_{x}p\end{array}\right)}_{e_r,e_{\theta },i},\text{div}u=\frac{1}{r}{\partial }_r\left(r{u}_r\right)+{\partial }_x{u}_x,\nabla u={\left(\begin{array}{ccc}{\partial }_r{u}_r& 0& {\partial }_r{u}_x\\ 0& \frac{u_r}{r}& 0\\ {\partial }_x{u}_r& 0& {\partial }_x{u}_x\end{array}\right)}_{e_r,e_{\theta },i}$

$\text{div}\sigma ={\left(\begin{array}{c}\frac{1}{r}{\partial }_r\left(r{\sigma }_{rr}\right)+{\partial }_x{\sigma }_{xr}-\frac{1}{r}{\sigma }_{\theta \theta }\\ \frac{1}{r}{\partial }_r\left(r{\sigma }_{r\theta }\right)+{\partial }_x{\sigma }_{x\theta }+\frac{1}{r}{\sigma }_{\theta r}\\ \frac{1}{r}{\partial }_r\left(r{\sigma }_{rx}\right)+{\partial }_x{\sigma }_{xx}\end{array}\right)}_{e_r,e_{\theta },i},$

where,

$\sigma ={\left(\begin{array}{ccc}{\sigma }_{rr}& {\sigma }_{r\theta }& {\sigma }_{rx}\\ {\sigma }_{\theta r}& {\sigma }_{\theta \theta }& {\sigma }_{\theta x}\\ {\sigma }_{xr}& {\sigma }_{x\theta }& {\sigma }_{xx}\end{array}\right)}_{e_r,e_{\theta },i}=\nu {\left(\begin{array}{ccc}2{\partial }_r{u}_r& 0& {\partial }_r{u}_x+{\partial }_x{u}_r\\ //& 2\frac{u_r}{r}& 0\\ //& //& 2{\partial }_x{u}_x\end{array}\right)}_{e_r,e_{\theta },i}.$ (10)

Finally, the Navier Stokes Equations (7) in cylindrical coordinates (8) are

$\{\begin{array}{l}\frac{1}{r}{\partial }_r\left(r{u}_r\right)+{\partial }_x{u}_x=0,\\ {\partial }_t{u}_r+\frac{1}{r}{\partial }_r\left(r{u}_r^2\right)+{\partial }_x\left(u_x{u}_r\right)+{\partial }_{r}p=\nu \left[\frac{2}{r}{\partial }_r\left(r{\partial }_r{u}_r\right)+{\partial }_x\left({\partial }_r{u}_x+{\partial }_x{u}_r\right)-\frac{2u_r}{r^2}\right],\\ \frac{1}{r}{\partial }_{\theta }p=0,\\ {\partial }_t{u}_x+\frac{1}{r}{\partial }_r\left(r{u}_r{u}_x\right)+{\partial }_x\left(u_x^2\right)+{\partial }_{x}p=\nu \left[\frac{1}{r}{\partial }_r\left(r\left({\partial }_r{u}_x+{\partial }_x{u}_r\right)\right)+2{\partial }_x^2{u}_x\right].\end{array}$ (11)

2.1.2. The Artery Wall Boundary

Crucial to our model derivation is the particular situation at the wall boundary, where the effect of wall deformation plays a central role. The wall boundary is the set of points

$\Gamma =\left\{xi+R\left(t,x\right)e_r\left(\theta \right);\theta \in \left[0,2\pi \right[,x\in \left[0,L\right],\left[t\in 0,T\right]\right\}.$ (12)

We define Γ’s tangential and outward normal vectors by

$t_x^w=\frac{1}{G}\left(i+{\partial }_{x}R\left(t,x\right)e_r\right),n^w=\frac{1}{G}\left(e_r-{\partial }_{x}R\left(t,x\right)i\right),G=\sqrt{1+{\left({\partial }_{x}R\right)}^2}.$

where $G$ is the axial arclength.

Since the wall is assumed to be rough, it produces friction, and due to its elastic behavior, it may get deformed. We take friction into account by considering the following Navier boundary condition on the wall Γ,

$\left(\sigma n^w\right)\cdot t_x^w=ku\cdot t_x^w,$ (13)

where $k\le 0$ is a friction term. Fluid-structure interaction is modeled with the condition

$u\cdot n^w={\partial }_{t}R{e}_r\cdot n^w.$ (14)

Thanks to the following hypothesis:

H1. Small thickness and plain stresses: The vessel wall thickness $h$ is assumed to be, a constant, small enough to allow a shell-type representation of the artery geometry. The vessel structure is subjected to plain stresses.

H2. Cylindrical reference geometry and radial displacement: As described before in (9) and (12), the artery is described by a circular cylindrical surface with straight axes and its displacements are only in the radial direction.

H3. Small deformation gradients and linear elastic behavior: We suppose that the artery wall behaves like a linear elastic solid where ${\partial }_{x}R$ is assumed to be bounded in time.

H4. Incompressibility: The wall tissue is incompressible [15].

H5. Dominance of circumferential stresses: Stresses acting along the axial direction can be neglected compared to circumferential ones.

One can derive the wall dynamic law (see [9] for further details):

${\partial }_t^2\eta +\frac{\rho }{h{\rho }_w}b\frac{\eta }{R_0^2}=\frac{R}{R_0}\frac{\rho }{h{\rho }^w}\left[p-p_{\text{ext}}-G\sigma n^w\cdot e_r\right],$ (15)

where $\eta =R-R_0$ is the displacement of the wall, $R_0=R\left(t=0,x\right)$ is the initial radius of the artery, ${\rho }^w$ is the wall tissues density, $h$ is the wall thickness, $p_{ext}$ is the external pressure, $b\left(x\right)=\frac{E\left(x\right)h}{\rho \left(1-{\xi }^2\right)}$ is a function of the Young modulus $E$ , $\sigma$ (10) is the fluid stress tensor at $r=R$ , and $p$ is the fluid pressure at $r=R$ .

Gathering (13), (14) and (15), boundary conditions can be written:

$\{\begin{array}{l}{u}_r-{\partial }_{x}R{u}_x={\partial }_{t}R,\\ \nu \left[2{\partial }_{x}R\left({\partial }_r{u}_r-{\partial }_x{u}_x\right)+\left(1-{({\partial }_{x}R)}^2\right)\left({\partial }_r{u}_x+{\partial }_x{u}_r\right)\right]=Gk\left(u_x+{\partial }_{x}R{u}_r\right),\\ \frac{h{\rho }_w}{\rho }\frac{R_0}{R}{\partial }_t^{2}R+b\frac{R-R_0}{R{R}_0}+\nu \left(2{\partial }_{r}u-{\partial }_{x}R\left({\partial }_r{u}_x+{\partial }_x{u}_r\right)\right)=p-p_{ext}.\end{array}$ (16)

2.2. Section-Averaging

We now proceed in Section 2.2.1 to write the Navier-Stokes equations with boundary conditions in non-dimensional form. Next, under a thin-artery assumption, we consider the radius of the artery to be small compared to the length, introducing a small parameter $\epsilon$ . We formally make an asymptotic expansion of the Navier-Stokes system in first (in Section 2.2.2) and second (in Section 2.2.4) order with respect to $\epsilon$ . Finally, we derive a section-averaged first- and second-order one-dimensional model for blood flow. Mathematical properties of both models are presented in Section 2.2.3 and in Section 2.2.5. Our approach is similar to those used in [2] [3] [9] [18] [19].

2.2.1. Dimensionless Navier-Stokes Equations

To derive the blood flow model, we assume that the artery’s radius is small compared to its length and that radial variations in velocity are small compared to axial ones. For large systemic arteries, typical values of the radius-to-length ratio are $R/L~{10}^{-2}$ to ${10}^{-1}$ , which justifies the thin-artery assumption. For smaller vessels (e.g., arterioles), this ratio can increase up to $O\left({10}^{-1}\right)$ , conserving the validity of the asymptotic expansion. This is achieved by postulating a small parameter ratio

$\epsilon :=\frac{\overline{R}}{L}=\frac{U_r}{U_x}\ll 1,$

where, $\overline{R}$ , $L$ , $U_r$ , and $U_x$ are the scales of, respectively, radius, length, radial velocity, and axial velocity. As a consequence, the time scale $T$ is such that

$T=\frac{L}{U_x}=\frac{\overline{R}}{U_r}.$

We also choose the pressure scale to be

$\overline{p}={\overline{U_x}}^2.$ (17)

It is convenient to define $L,U_x$ and $T$ , as finite constants with respect to $\epsilon$ , while $\overline{R}=\epsilon L$ and $U_r=\epsilon U_x$ . This allows us to introduce the dimensionless quantities of time $\tilde{t}$ , space ($\tilde{x}$ , $\tilde{r}$ ), pressure $\tilde{p}$ and velocity field (${\tilde{u}}_x$ , ${\tilde{u}}_r$ ) via the following scaling relations

$\begin{array}{ll}\tilde{t}:=\frac{t}{T},\hfill & \tilde{p}\left(\tilde{t},\tilde{r},\tilde{x}\right)=\frac{p\left(t,r,x\right)}{\overline{p}},\hfill \\ \tilde{x}:=\frac{x}{L},\hfill & {\tilde{u}}_x\left(\tilde{t},\tilde{r},\tilde{x}\right)=\frac{u_x\left(t,r,x\right)}{U_x},\hfill \\ \tilde{r}:=\frac{r}{\overline{R}}=\frac{r}{\epsilon L},\hfill & {\tilde{u}}_r\left(\tilde{t},\tilde{r},\tilde{x}\right)=\frac{u_r\left(t,r,x\right)}{U_r}.\hfill \end{array}$ (18)

We also rescale the following coefficients

$\begin{array}{ll}\tilde{k}=\frac{k}{U_r}=\frac{k}{\epsilon U_x},\hfill & \tilde{R}\left(\tilde{t},\tilde{x}\right)=\frac{R\left(t,x\right)}{\overline{R}},\hfill \\ \tilde{h}=\frac{h}{\epsilon \overline{R}},\hfill & \tilde{E}\left(\tilde{x}\right)=\epsilon \frac{E\left(x\right)}{\rho U_x^2},\hfill \\ {\tilde{R}}_0\left(\tilde{x}\right)=\frac{R_0\left(x\right)}{\overline{R}}.\hfill & \hfill \end{array}$ (19)

Finally, we define the non-dimensional Reynolds number,

$R_e=\frac{L{U}_x}{\nu },$

and

${\nu }_0={\left(\epsilon R_e\right)}^{-1}$

yielding to the asymptotic regime

$R_e^{-1}={\nu }_0\epsilon .$ (20)

Using these dimensionless variables (17), (18), (19) and (20) in the Navier-Stokes Equations (11) and the boundary conditions (16), we get

$\{\begin{array}{l}\frac{1}{\tilde{r}}{\partial }_{\tilde{r}}\left(\tilde{r}{\tilde{u}}_r\right)+{\partial }_{\tilde{x}}{\tilde{u}}_x=0,\\ {\partial }_{\tilde{r}}\tilde{p}=\epsilon {\nu }_0\left[\frac{2}{\tilde{r}}{\partial }_{\tilde{r}}\left(\tilde{r}{\partial }_{\tilde{r}}{\tilde{u}}_r\right)+{\partial }_{\tilde{x}}\left({\partial }_{\tilde{r}}{\tilde{u}}_x\right)-2\frac{{\tilde{u}}_r}{{\tilde{r}}^2}\right]+{\epsilon }^2{\delta }_{r,\epsilon ,\tilde{u}},\\ {\partial }_{\tilde{t}}{\tilde{u}}_x+\frac{1}{\tilde{r}}{\partial }_{\tilde{r}}\left(\tilde{r}{\tilde{u}}_r{\tilde{u}}_x\right)+{\partial }_{\tilde{x}}\left({\tilde{u}}_x^2\right)+{\partial }_{\tilde{x}}\tilde{p}=\frac{{\nu }_0}{\epsilon }\frac{1}{\tilde{r}}{\partial }_{\tilde{r}}\left(\tilde{r}{\partial }_{\tilde{r}}{\tilde{u}}_x\right)+\epsilon {\nu }_0\left[\frac{1}{\tilde{r}}{\partial }_{\tilde{r}}\left(\tilde{r}{\partial }_{\tilde{x}}{\tilde{u}}_r\right)+2{\partial }_{\tilde{x}}^2{\tilde{u}}_x\right],\\ {\tilde{u}}_r-{\partial }_{\tilde{x}}\tilde{R}{\tilde{u}}_x={\partial }_{\tilde{t}}\tilde{R},\\ \frac{1}{\epsilon }{\nu }_0{\partial }_{\tilde{r}}{\tilde{u}}_x-G\tilde{k}{\tilde{u}}_x=-\epsilon {\nu }_0\left[2{\partial }_{\tilde{x}}\tilde{R}\left({\partial }_{\tilde{r}}{\tilde{u}}_r-{\partial }_{\tilde{x}}{\tilde{u}}_x\right)+{\partial }_{\tilde{x}}{\tilde{u}}_r-{\left({\partial }_{\tilde{x}}\tilde{R}\right)}^2{\partial }_{\tilde{r}}{\tilde{u}}_x\right]+{\epsilon }^2{\delta }_{R,\epsilon ,\tilde{u}},\\ \tilde{b}\frac{\tilde{R}-{\tilde{R}}_0}{\tilde{R}{\tilde{R}}_0}+\epsilon {\nu }_0\left(2{\partial }_{\tilde{r}}{\tilde{u}}_r-{\partial }_{\tilde{x}}\tilde{R}{\partial }_{\tilde{r}}{\tilde{u}}_x\right)=\tilde{p}+{\epsilon }^2{\delta }_{p,\epsilon ,\tilde{u}}.\end{array}$ (21)

where,

$\begin{array}{l}{\delta }_{r,\epsilon ,\tilde{u}}=-\left({\partial }_{\tilde{t}}{\tilde{u}}_r+\frac{1}{\tilde{r}}{\partial }_{\tilde{r}}\left(\tilde{r}{\tilde{u}}_r{}^2\right)+{\partial }_{\tilde{x}}\left({\tilde{u}}_x{\tilde{u}}_r\right)\right)+\epsilon {\nu }_0{\partial }_{\tilde{x}}^2{\tilde{u}}_r,\\ {\delta }_{R,\epsilon ,\tilde{u}}=G\tilde{k}{\partial }_{\tilde{x}}\tilde{R}{\tilde{u}}_r+{\nu }_0\epsilon {\left({\partial }_{\tilde{x}}\tilde{R}\right)}^2{\partial }_{\tilde{x}}{\tilde{u}}_r,\\ {\delta }_{p,\epsilon ,\tilde{u}}=\epsilon {\nu }_0{\partial }_{\tilde{x}}\tilde{R}{\partial }_{\tilde{x}}{\tilde{u}}_r-\epsilon \frac{\tilde{h}{\rho }_w}{\rho }\frac{{\tilde{R}}_0}{\tilde{R}}{\partial }_{\tilde{t}}^2\tilde{R},\\ \tilde{b}\left(\tilde{x}\right)=\frac{\tilde{E}\left(\tilde{x}\right)\tilde{h}}{1-{\xi }^2}.\end{array}$

2.2.2. First-Order Approximation of the Dimensionless Navier-Stokes Equations

Omitting $\tilde{\cdot }$ , remarking that $G=\sqrt{1+{\epsilon }^2{\left({\partial }_{\tilde{x}}\tilde{R}\right)}^2}=1+O\left({\epsilon }^2\right)$ , gathering all first-order terms in $O\left(\epsilon \right)$ , Equations (21) become

$\{\begin{array}{l}\frac{1}{r}{\partial }_r\left(r{u}_r\right)+{\partial }_x{u}_x=0,\\ {\partial }_{r}p=O\left(\epsilon \right),\\ {\partial }_t{u}_x+\frac{1}{r}{\partial }_r\left(r{u}_r{u}_x\right)+{\partial }_x\left(u_x^2\right)+{\partial }_{x}p=\frac{{\nu }_0}{\epsilon }\frac{1}{r}{\partial }_r\left(r{\partial }_r{u}_x\right)+O\left(\epsilon \right),\\ u_r-{\partial }_{x}R{u}_x={\partial }_{t}R,\\ \frac{{\nu }_0}{\epsilon }{\partial }_r{u}_x=k{u}_x+O\left(\epsilon \right),\\ b\frac{R-R_0}{R{R}_0}=p+O\left(\epsilon \right).\end{array}$ (22)

Integrating from $r$ and $R$ the radial momentum equation (the second equation in (22)), we obtain the following pressure

$p\left(t,r,x\right)=p\left(t,R,x\right)+O\left(\epsilon \right)=b\left(x\right)\frac{R\left(t,x\right)-R_0\left(t,x\right)}{R\left(t,x\right)R_0\left(t,x\right)}+O\left(\epsilon \right).$ (23)

As done in [9], we linearize the Equation (23) with respect to $R$ to obtain

$p\left(t,r,x\right)=b\left(x\right)\frac{R\left(t,x\right)-R_0\left(t,x\right)}{R_0^2\left(t,x\right)}+O\left(\epsilon \right).$ (24)

Moreover, by identifying terms at order $\frac{1}{\epsilon }$ in the axial momentum equation (the third equation in (22)), thanks to the boundary conditions (the fifth equation in (22)), we obtain the “motion by slice” decomposition

${\nu }_0\frac{1}{r}{\partial }_r\left(r{\partial }_r{u}_x\right)=O\left(\epsilon \right)\Rightarrow r{\partial }_r{u}_x=O\left(\epsilon \right)\Rightarrow u_x\left(t,r,x\right)=u_{x,0}\left(t,x\right)+O\left(\epsilon \right),$

for some function $u_{x,0}$ . Noting

$\overline{u_x}\left(t,x\right)=\frac{1}{\left|S_{t,x}\right|}{\displaystyle {\int }_{S_{t,x}}{u}_x\left(t,r,x\right)\text{d}{S}_{t,x}}=\frac{1}{A\left(t,x\right)}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R\left(t,x\right)}r{u}_x\left(t,r,x\right)\text{d}r\text{d}\theta }}$

as the mean speed of the fluid over a cross-section of the artery $S_{t,x}$ where

$S_{t,x}=\left\{xi+r{e}_r\left(\theta \right);r\in \left[0,R\left(t,x\right)\right[,\left[\theta \in 0,2\pi \right[\right\}$ (25)

and

$\left|S_{t,x}\right|=A\left(t,x\right)=\pi R{\left(t,x\right)}^2,$ (26)

we have the properties

$u_x\left(t,r,x\right)=\overline{u_x}\left(t,x\right)+O\left(\epsilon \right)\text{and}\overline{u_x^2\left(t,r,x\right)}={\overline{u_x}}^2\left(t,x\right)+O\left(\epsilon \right).$(27)

Integrating the divergence equation (the first equation in (22)) over the section $S_{t,x}$ (25), we obtain

$\begin{array}{c}0={\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^R\left[{\partial }_r\left(r{u}_r\right)+{\partial }_x\left(r{u}_x\right)\right]\text{d}r\text{d}\theta }}\\ =2\pi R{u}_r\left(r=R\right)+{\partial }_x\left({\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R}r{u}_x}}\right)-2\pi {\partial }_{x}RR{u}_x\left(r=R\right).\end{array}$

In view of the normal boundary condition (the fourth equation in (22)), we get the mass conservation equation

${\partial }_{t}A+{\partial }_{x}Q=0,$ (28)

where

$Q\left(t,x\right)=A\left(t,x\right)\overline{u_x}\left(t,x\right)={\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R}r{u}_x\text{d}r\text{d}\theta }}$(29)

is the blood flow rate through the section $S_{t,x}$ of the artery.

Then, integrating the axial momentum equation (the third equation in (22)), we get

$\begin{array}{l}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R}r\left[{\partial }_t{u}_x+\frac{1}{r}{\partial }_r\left(r{u}_r{u}_x\right)+{\partial }_x\left(u_x^2\right)+{\partial }_{x}p\right]\text{d}r\text{d}\theta }}\\ ={\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R}r\left[\frac{{\nu }_0}{\epsilon }\frac{1}{r}{\partial }_r\left(r{\partial }_r{u}_x\right)\right]\text{d}r\text{d}\theta }}+O\left(\epsilon \right)\end{array}$

yielding to

$\begin{array}{l}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^R{\partial }_t\left(r{u}_x\right)\text{d}r\text{d}\theta }}+2\pi R{u}_r\left(r=R\right)u_x\left(r=R\right)+{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^R{\partial }_x\left(r{u}_x^2\right)\text{d}r\text{d}\theta }}\\ +{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^R{\partial }_x\left(rp\right)\text{d}r\text{d}\theta }}=\frac{{\nu }_0}{\epsilon }2\pi R\left({\partial }_r{u}_x\right)\left(r=R\right)+O\left(\epsilon \right).\end{array}$

Using the definition of $A$ (26) and $Q$ (29), thanks to the normal boundary (the fourth, fifth, and sixth equation in (22)), we have

${\partial }_{t}Q+{\partial }_x\left(A\overline{u_x^2}\right)+{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^R{\partial }_x\left(rp\right)\text{d}r\text{d}\theta }}=2\pi Rk{u}_x\left(r=R\right)+O\left(\epsilon \right)$

and thanks to the radial momentum equation (the second equation in (22)) and (27), we deduce the equation of momentum

${\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{A}+Ap\right)=p{\partial }_{x}A+2\pi Rk\frac{Q}{A}+O\left(\epsilon \right)$ (30)

Finally, from Equations (24), (28) and (30), dropping all terms of the first order, we obtain the following section-averaged one-dimensional model for blood flow

$\{\begin{array}{l}{\partial }_{t}A+{\partial }_{x}Q=0,\\ {\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{A}+Ap\right)=p{\partial }_{x}A+2\pi Rk\frac{Q}{A},\\ p=b\frac{R-R_0}{R_0^2}.\end{array}$ (31)

2.2.3. Mathematical Properties for the First-Order Model

We have the following results.

Theorem 1. Let $\left(A,\overline{u_x}\right)$ and $Q=A\overline{u_x}$ satisfy the one-dimensional blood flow system (31), we have:

1) System (31) is strictly hyperbolic on the set $\left\{A>0\right\}$ .

2) For smooth $\left(A,\overline{u_x}\right)$ , in the region where $A>0$ , the mean velocity $\overline{u_x}$ satisfies the head equation

${\partial }_t\overline{u_x}+{\partial }_x\psi \left(\overline{u_x},p\right)=2\pi Rk\frac{\overline{u_x}}{A},$ (32)

where

$\psi \left(\overline{u_x},p\right)=\frac{{\overline{u_x}}^2}{2}+p$ (33)

is the total head.

3) For smooth $\left(A,\overline{u_x}\right)$ , the steady state reads

$\overline{u_x}=0,p=p_0\text{for}\text{some}\text{constant}{p}_0.$

In particular, for $A=A_0$ , we have $p_0=0$ .

4) The pair of function $\left(E,E+\tilde{p}\right)$ with $E:=A\psi -\tilde{p}$ forms a mathematical entropy pair for system (31), in that they satisfy the following entropy relation for smooth $\left(A,\overline{u_x}\right)$ :

${\partial }_{t}E+{\partial }_x\left(\left(E+\tilde{p}\right)\overline{u_x}\right)=2\pi Rk{\overline{u_x}}^2\le 0$ (34)

where $\tilde{p}=\frac{b\sqrt{\pi }}{3A_0}\left(A^{\frac{3}{2}}-A_0^{\frac{3}{2}}\right)$ .

5) For $\overline{u_x}$ such that $\overline{u_x}\left(x=0\right)=\overline{u_x}\left(x=L\right)=0$ , the total energy in Ω, $E_{tot}\left(t\right)={\displaystyle {\int }_0^{L}E\left(t,x\right)\text{d}x}$ decreases and satisfies the inequality

$\frac{\text{d}}{\text{d}t}{E}_{tot}=\frac{2\pi Rk}{1-\epsilon \frac{Rk}{4{\nu }_0}}{\Vert \overline{u_x}\Vert }_{L^2\left(\left[0,L\right]\right)}^2\le 0.$

Proof.

1) To prove the hyperbolicity of system (31), we write

${\partial }_{t}U+H{\partial }_{x}U=S$

where $U=\left(\begin{array}{l}A\hfill \\ Q\hfill \end{array}\right)$ , $H=\left(\begin{array}{cc}0& 1\\ -\frac{Q^2}{A^2}+A{\partial }_{A}p& 2\frac{Q}{A}\end{array}\right)=\left(\begin{array}{cc}0& 1\\ -\frac{Q^2}{A^2}+\frac{b\sqrt{\pi }}{2A_0}\sqrt{A}& 2\frac{Q}{A}\end{array}\right)$ , $S=\left(\begin{array}{c}0\\ 2\pi Rk\frac{Q}{A}\end{array}\right)$ . The eigenvalues of $H$ read ${\lambda }_1\left(A,Q\right)=\frac{Q}{A}-\sqrt{\frac{b\sqrt{\pi }}{2A_0}}{A}^{\frac{1}{4}}$ , ${\lambda }_2\left(A,Q\right)=\frac{Q}{A}+\sqrt{\frac{b\sqrt{\pi }}{2A_0}}{A}^{\frac{1}{4}}$ and are real and distinct for $A>0$ . Therefore, system (31) is strictly hyperbolic on the set $\left\{A>0\right\}$ .

2) We rewrite the momentum (second) equation in system (31) in terms of $\left(A,\overline{u_x}\right)$ , with $\overline{u_x}=\frac{Q}{A}$ , as

${\partial }_t\left(A\overline{u_x}\right)+{\partial }_x\left(A{\overline{u_x}}^2+Ap\right)=p{\partial }_{x}A+2\pi Rk\overline{u_x}.$

Applying the product rule to the first term of (30) and substituting in the conservation of mass equation, we get

$A{\partial }_t\overline{u_x}-\overline{u_x}{\partial }_x\left(A\overline{u_x}\right)+{\partial }_x\left(A{\overline{u_x}}^2\right)+A{\partial }_{x}p=2\pi Rk\overline{u_x}.$

Again, using the product rule, it follows

$A{\partial }_t\overline{u_x}+A\overline{u_x}{\partial }_x\left(\overline{u_x}\right)+A{\partial }_{x}p=2\pi Rk\overline{u_x}.$

Finally, dividing by $A>0$ , we get the head equation

${\partial }_t\overline{u_x}+{\partial }_x\psi \left(\overline{u_x},p\right)=2\pi Rk\frac{\overline{u_x}}{A}.$

3) Looking for still steady states, in the Equation (32) for instance (or (31)), we have, ${\partial }_x\psi \left(\overline{u_x},p\right)=0$ meaning that ${\partial }_{x}p=0$ , i.e., $p=p_0$ for some constant $p_0$ . Recalling that $p=b\sqrt{\pi }\frac{\sqrt{A}-\sqrt{A_0}}{A_0}$ , if $A=A_0$ , then we immediately deduce that $p_0=0$ . This steady state can be easily preserved from a numerical point of view.

4) To derive the entropy relation for system (31), we multiply the conservation of mass by $\psi$ which leads to

$\psi {\partial }_{t}A+\psi {\partial }_x\left(A\overline{u_x}\right)=0,$

then,

${\partial }_t\left(A\psi \right)+{\partial }_x\left(A\psi \overline{u_x}\right)-A\left[{\partial }_t\psi +\overline{u_x}{\partial }_x\psi \right]=0,$

now, using the definition of $\psi$ in (33), we get,

${\partial }_t\left(A\psi \right)+{\partial }_x\left(A\psi \overline{u_x}\right)-A\left[\overline{u_x}\left({\partial }_t\overline{u_x}+{\partial }_x\psi \right)+{\partial }_{t}p\right]=0.$

We use the momentum Equation (32) to get,

${\partial }_t\left(A\psi -\tilde{p}\right)+{\partial }_x\left(A\psi \overline{u_x}\right)=2\pi Rk{\overline{u_x}}^2,$

where $\tilde{p}={\displaystyle \int A{\partial }_{A}p}=\frac{b\sqrt{\pi }}{3A_0}\left(A^{\frac{3}{2}}-A_0^{\frac{3}{2}}\right)$ . We define the mathematical entropy $E=A\psi -\tilde{p}$ and we obtain

${\partial }_{t}E+{\partial }_x\left(\left(E+\tilde{p}\right)\overline{u_x}\right)=2\pi Rk{\overline{u_x}}^2.$

Finally, using that $k<0$ , we get consistent energy decay, meaning,

${\partial }_{t}E+{\partial }_x\left(\left(E+\tilde{p}\right)\overline{u_x}\right)\le 0.$

5) Integrating (34) between 0 and $L$ , we get,

${\partial }_t{\displaystyle {\int }_0^{L}E}+{\left[\left(E+\tilde{p}\right)\overline{u_x}\right]}_0^L=2\pi Rk{\displaystyle {\int }_0^L{\overline{u_x}}^2},$

which, when we suppose $\overline{u_x}\left(x=0\right)=\overline{u_x}\left(x=L\right)=0$ leads to,

$\frac{\text{d}}{\text{d}t}{E}_{tot}\left(t\right)=\frac{2\pi Rk}{1-\epsilon \frac{Rk}{4{\nu }_0}}{\Vert \overline{u_x}\Vert }_{L^2\left(\left[0,L\right]\right)}^2\le 0,$

where $E_{tot}\left(t\right)={\displaystyle {\int }_0^{L}E\text{d}x}$ is the total energy in Ω.

□

2.2.4. Second-Order Approximation of the Dimensionless Navier-Stokes Equations

We can improve the order of accuracy of the section-averaged one-dimensional blood flow system (31) by determining the first-order correction depending on $r$ in the expansion of $u_x\left(t,r,s\right)$ (see (27)). To do so, we come back to Equation (21) and gather all second-order terms in $O\left({\epsilon }^2\right)$ , leading to,

$\{\begin{array}{l}\frac{1}{r}{\partial }_r\left(r{u}_r\right)+{\partial }_x{u}_x=0,\\ {\partial }_{r}p=\epsilon {\nu }_0\left[\frac{2}{r}{\partial }_r\left(r{\partial }_r{u}_r\right)+{\partial }_x\left({\partial }_r{u}_x\right)-2\frac{u_r}{r^2}\right]+O\left({\epsilon }^2\right),\\ {\partial }_t{u}_x+\frac{1}{r}{\partial }_r\left(r{u}_r{u}_x\right)+{\partial }_x\left(u_x^2\right)+{\partial }_{x}p=\frac{{\nu }_0}{\epsilon }\frac{1}{r}{\partial }_r\left(r{\partial }_r{u}_x\right)+\epsilon {\nu }_0\left[\frac{1}{r}{\partial }_r\left(r{\partial }_x{u}_r\right)+2{\partial }_x^2{u}_x\right]+O\left({\epsilon }^2\right),\\ u_r-{\partial }_{x}R{u}_x={\partial }_{t}R,\\ \frac{1}{\epsilon }{\nu }_0{\partial }_r{u}_x-Gk{u}_x+O\left({\epsilon }^2\right)=-\epsilon {\nu }_0\left[2{\partial }_{x}R\left({\partial }_r{u}_r-{\partial }_x{u}_x\right)+{\partial }_x{u}_r-{\left({\partial }_{x}R\right)}^2{\partial }_r{u}_x\right],\\ b\frac{R-R_0}{R{R}_0}+\epsilon {\nu }_0\left(2{\partial }_r{u}_r-{\partial }_{x}R{\partial }_r{u}_x\right)=p+O\left({\epsilon }^2\right).\end{array}$ (35)

Then, we remark on the axial momentum equation

$\begin{array}{c}\frac{{\nu }_0}{\epsilon }\frac{1}{r}{\partial }_r\left(r{\partial }_r{u}_x\right)={\partial }_t{u}_x+\frac{1}{r}{\partial }_r\left(r{u}_r{u}_x\right)+{\partial }_x\left(u_x^2\right)+{\partial }_{x}p+O\left(\epsilon \right)\\ ={\partial }_t{u}_x+r{\partial }_r\frac{1}{r}{u}_r{u}_x+u_r{\partial }_r{u}_x+2u_x{\partial }_x\left(u_x\right)+{\partial }_{x}p+O\left(\epsilon \right)\\ ={\partial }_t{u}_x+u_r{\partial }_r{u}_x+u_x{\partial }_x\left(u_x\right)+{\partial }_{x}p+O\left(\epsilon \right)\\ ={\partial }_t\overline{u_x}+{\partial }_x\left(\frac{{\overline{u_x}}^2}{2}+p\right)+O\left(\epsilon \right)\\ =2\pi Rk\frac{u_x\left(r=0\right)}{A}+O\left(\epsilon \right).\end{array}$

Multiplying by $r$ and integrating from 0 to $r$ , we get

$\frac{{\nu }_0}{\epsilon }r{\partial }_r{u}_x=2\pi Rk\frac{r^2}{2}\frac{u_x\left(r=0\right)}{A}+O\left(\epsilon \right),$

and we obtain the following formula which gives a more detailed view of the axial velocity through a parabolic correction

$u_x=u_x\left(r=0\right)\left[1+\epsilon \frac{Rk}{2{\nu }_0}{\left(\frac{r}{R}\right)}^2\right]+O\left({\epsilon }^2\right).$

Then, integrating over a section $S_{t,x}$ as defined in (25), we obtain

$\overline{u_x}=\left[1+\epsilon \frac{Rk}{4{\nu }_0}\right]u_x\left(r=0\right)+O\left({\epsilon }^2\right),$

meaning,

$\begin{array}{c}{u}_x=\overline{u_x}\left[1-\epsilon \frac{Rk}{4{\nu }_0}\right]\left[1+\epsilon \frac{Rk}{2{\nu }_0}{\left(\frac{r}{R}\right)}^2\right]+O\left({\epsilon }^2\right)\\ =\overline{u_x}\left[1-\epsilon \frac{Rk}{4{\nu }_0}\left(1-2{\left(\frac{r}{R}\right)}^2\right)\right]+O\left({\epsilon }^2\right).\end{array}$ (36)

Moreover,

$\begin{array}{c}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R}r{u}_x^2\text{d}r\text{d}\theta }}={\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R}r{\overline{u_x}}^2}}{\left[1-\epsilon \frac{Rk}{4{\nu }_0}\left(1-2{\left(\frac{r}{R}\right)}^2\right)\right]}^2+O\left({\epsilon }^2\right)\\ =A{\overline{u_x}}^2+O\left({\epsilon }^2\right)\end{array}$(37)

We write the pressure as follows

$\begin{array}{c}{\partial }_{r}p=\epsilon \nu \left[\frac{2}{r}{\partial }_r\left(r{\partial }_r{u}_r\right)+{\partial }_x\left({\partial }_r{u}_x\right)-2\frac{u_r}{r^2}\right]+O\left({\epsilon }^2\right)\\ =\epsilon \nu \left[\frac{2}{r}{\partial }_r{u}_r+2{\partial }_r^2{u}_r+{\partial }_x\left({\partial }_r{u}_x\right)-2\frac{u_r}{r^2}\right]+O\left({\epsilon }^2\right)\\ =\epsilon \nu \left[2{\partial }_r\frac{u_r}{r}+2{\partial }_r^2{u}_r+{\partial }_x\left({\partial }_r{u}_x\right)\right]+O\left({\epsilon }^2\right)\\ ={\partial }_r\left(\epsilon {\nu }_0\left(2\frac{u_r}{r}+2{\partial }_r{u}_r+{\partial }_x{u}_x\right)\right)+O\left({\epsilon }^2\right).\end{array}$

Integrating from $R$ to $r$ , we obtain

$p=p\left(R\right)+\epsilon {\nu }_0\left[2\frac{u_r}{r}+2{\partial }_r{u}_r+{\partial }_x{u}_x-2\frac{u_r\left(R\right)}{R}-2{\partial }_r{u}_r\left(R\right)-{\partial }_x{u}_x\left(R\right)\right]+O\left({\epsilon }^2\right),$

We now use the dynamic wall equation (the sixth equation in (35)) replacing $p\left(R\right)$ by its expression,

$p=b\frac{R-R_0}{R{R}_0}+\epsilon {\nu }_0\left[-{\partial }_{x}R{\partial }_r{u}_x\left(R\right)+\frac{2}{r}{\partial }_r\left(r{u}_r\right)+{\partial }_x{u}_x-2\frac{u_r\left(R\right)}{R}-{\partial }_x{u}_x\left(R\right)\right]+O\left({\epsilon }^2\right),$

and using the divergence equation (the first equation in (35)),

$p=b\frac{R-R_0}{R{R}_0}+\epsilon {\nu }_0\left[-{\partial }_{x}R{\partial }_r{u}_x\left(R\right)-{\partial }_x{u}_x-2\frac{u_r\left(R\right)}{R}-{\partial }_x{u}_x\left(R\right)\right]+O\left({\epsilon }^2\right).$

We now use the tangential Navier boundary condition (the fifth equation in (35)), mainly that ${\partial }_r{u}_x\left(r=R\right)=O\left(\epsilon \right)$ and the normal boundary condition at the wall (the fourth equation in (35)) yielding to

$p=b\frac{R-R_0}{R{R}_0}+\epsilon {\nu }_0\left[-{\partial }_x{u}_x-\frac{{\partial }_{t}A+{\partial }_{x}A{u}_x\left(R\right)}{A}-{\partial }_x{u}_x\left(R\right)\right]+O\left({\epsilon }^2\right),$

finally giving,

$p=b\frac{R-R_0}{R{R}_0}-\epsilon {\nu }_0{\partial }_x{u}_x\left(R\right)+O\left({\epsilon }^2\right).$

which can be approximated by

$p=b\frac{R-R_0}{R_0^2}-\epsilon {\nu }_0{\partial }_x{u}_x\left(R\right)+O\left({\epsilon }^2\right)$

as done in Section 2.2.2 (see also Equation (24)).

The left-hand side of the axial momentum equation (the third equation in (35)) can be integrated, keeping in mind (36) and (37), as follows

$\begin{array}{l}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R}r}}\left[{\partial }_t{u}_x+\frac{1}{r}{\partial }_r\left(r{u}_r{u}_x\right)+{\partial }_x\left(u_x^2\right)+{\partial }_{x}p\right]\text{d}r\text{d}\theta \\ ={\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{A}\right)+2\pi {\partial }_x{\displaystyle {\int }_0^{R}rp\text{d}r}-2\pi {\partial }_{x}RRp\left(R\right)+O\left({\epsilon }^2\right)\\ ={\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{A}+Ab\frac{R-R_0}{R_0^2}\right)-\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\left(\overline{u_x}\right)\right)-2\pi {\partial }_{x}RRp\left(R\right)+O\left({\epsilon }^2\right)\end{array}$

The right-hand side of the axial momentum equation (the third equation in (35)) together with the pressure term provides

$\begin{array}{l}2\pi {\partial }_{x}RRp\left(R\right)+{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{R}r}}\left[\frac{{\nu }_0}{\epsilon }\frac{1}{r}{\partial }_r\left(r{\partial }_r{u}_x\right)+\epsilon {\nu }_0\left[\frac{1}{r}{\partial }_r\left(r{\partial }_x{u}_r\right)+2{\partial }_x^2{u}_x\right]\right]\text{d}r\text{d}\theta +O\left({\epsilon }^2\right)\\ =2\pi {\partial }_{x}RR\left[b\frac{R-R_0}{R{R}_0}+\epsilon {\nu }_0\left(2{\partial }_r{u}_r\left(R\right)-{\partial }_{x}R{\partial }_r{u}_x\left(R\right)\right)\right]\\ +\left[2\pi \frac{{\nu }_0}{\epsilon }R{\partial }_r{u}_x\left(R\right)+2\pi \epsilon {\nu }_0\left[R{\partial }_x{u}_r\left(R\right)+R^2{\partial }_x^2{u}_x\right]\right]+O\left({\epsilon }^2\right)\end{array}$

$\begin{array}{l}={\partial }_{x}Ab\frac{R-R_0}{R{R}_0}+2\pi Rk{u}_x\left(R\right)+2\pi R\epsilon {\nu }_0\left[2{\partial }_{x}R{\partial }_r{u}_r\left(R\right)-{\left({\partial }_{x}R\right)}^2{\partial }_r{u}_x\left(R\right)\right]\\ -2\pi R\epsilon {\nu }_0\left[2{\partial }_{x}R\left({\partial }_r{u}_r\left(R\right)-{\partial }_x{u}_x\left(R\right)\right)+{\partial }_x{u}_r\left(R\right)-{\left({\partial }_{x}R\right)}^2{\partial }_r{u}_x\left(R\right)\right]\\ +2\pi R\epsilon {\nu }_0\left[{\partial }_x{u}_r\left(R\right)+R{\partial }_x^2{u}_x\right]+O\left({\epsilon }^2\right)\\ ={\partial }_{x}Ab\frac{R-R_0}{R{R}_0}+2\pi Rk{u}_x\left(R\right)+2\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\overline{u_x}\right)+O\left({\epsilon }^2\right)\end{array}$

Finally, the section-averaged one-dimensional viscous model for blood flow reads,

$\{\begin{array}{l}{\partial }_{t}A+{\partial }_{x}Q=0,\\ {\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{A}+Ap\right)-3\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\left(\frac{Q}{A}\right)\right)={\partial }_{x}Ap+\frac{2\pi Rk}{1-\epsilon \frac{Rk}{4{\nu }_0}}\frac{Q}{A},\end{array}$ (38)

resulting from an approximation in $O\left({\epsilon }^2\right)$ of the Navier Stokes equations. From now on, in these equations, $p$ stands for the pressure $p=b\frac{R-R_0}{R_0^2}$ . We emphasize that, at this order, the Coriolis-Boussinesq coefficient $\alpha =\frac{\overline{u_x^2}}{{\overline{u_x}}^2}$ is equal to 1. Our model differs from the existing ones from this parameter, see for instance [9].

2.2.5. Mathematical Properties for the Second-Order Model

We have the following results.

Theorem 2. Let $\left(A,\overline{u_x}\right)$ and $Q=A\overline{u_x}$ satisfy the one-dimensional blood flow system (38).

1) For smooth $\left(A,\overline{u_x}\right)$ , in the region where $A>0$ , the mean velocity $\overline{u_x}$ satisfies the head equation

${\partial }_t\overline{u_x}+{\partial }_x\psi \left(A,\overline{u_x},x\right)=3\epsilon {\nu }_0\frac{1}{A}{\partial }_x\left(A{\partial }_x\overline{u_x}\right)+2\pi Rk\frac{\overline{u_x}}{A},$

where $\psi \left(\overline{u_x},p\right)$ is the total head defined in (33).

2) For smooth $\left(A,\overline{u_x}\right)$ , the steady state reads

$\overline{u_x}=0,p=p_0\text{for}\text{some}\text{constant}{p}_0.$

In particular, for $A=A_0$ , we have $p_0=0$ .

3) The pair of function $\left(E,E+\tilde{p}\right)$ with $E:=A\psi -\tilde{p}$ forms a mathematical entropy pair for system (38), in that they satisfy the following entropy relation for smooth $\left(A,\overline{u_x}\right)$ :

${\partial }_{t}E+{\partial }_x\left(\left(E+\tilde{p}\right)\overline{u_x}\right)=3\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\overline{u_x}\right)\overline{u_x}+2\pi Rk{\overline{u_x}}^2.$

where $\tilde{p}=\frac{b\sqrt{\pi }}{3A_0}\left(A^{\frac{3}{2}}-A_0^{\frac{3}{2}}\right)$ .

4) For $\overline{u_x}$ such that $\overline{u_x}\left(x=0\right)=\overline{u_x}\left(x=L\right)=0$ , the total energy in our domain $E_{tot}\left(t\right)={\displaystyle {\int }_0^{L}E\text{d}x}$ decreases and satisfies

$\frac{\text{d}}{\text{d}t}{E}_{tot}\left(t\right)=-3\epsilon {\nu }_0{\displaystyle {\int }_0^{L}A{\left({\partial }_x\overline{u_x}\right)}^2}+\frac{2\pi Rk}{1-\epsilon \frac{Rk}{4{\nu }_0}}{\Vert \overline{u_x}\Vert }_{L^2\left(\left[0,L\right]\right)}^2\le 0.$

Proof.

1) We rewrite the momentum (second) equation in system (31) in terms of $\left(A,\overline{u_x}\right)$ , with $\overline{u_x}=\frac{Q}{A}$ , as

${\partial }_t\left(A\overline{u_x}\right)+{\partial }_x\left(A{\overline{u_x}}^2+Ap\right)=3\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\overline{u_x}\right)+p{\partial }_{x}A+2\pi Rk\overline{u_x}.$

Applying the product rule to the first term of (30) and substituting in the conservation of mass equation, we get

$A{\partial }_t\overline{u_x}-\overline{u_x}{\partial }_x\left(A\overline{u_x}\right)+{\partial }_x\left(A{\overline{u_x}}^2\right)+A{\partial }_{x}p=3\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\overline{u_x}\right)+2\pi Rk\overline{u_x}.$

Again, using the product rule, it follows

$A{\partial }_t\overline{u_x}+A\overline{u_x}{\partial }_x\left(\overline{u_x}\right)+A{\partial }_{x}p=3\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\overline{u_x}\right)+2\pi Rk\overline{u_x}.$

Finally, dividing by $A>0$ , we get the head equation associated with system (31),

${\partial }_t\overline{u_x}+{\partial }_x\psi \left(\overline{u_x},p\right)=3\epsilon {\nu }_0\frac{1}{A}{\partial }_x\left(A{\partial }_x\overline{u_x}\right)+2\pi Rk\frac{\overline{u_x}}{A}.$

2) We refer to the proof of Theorem 1.

3) To derive the entropy relation for system (31), we multiply the conservation of mass by $\psi$ which leads to

$\psi {\partial }_{t}A+\psi {\partial }_x\left(A\overline{u_x}\right)=0,$

then,

${\partial }_t\left(A\psi \right)+{\partial }_x\left(A\psi \overline{u_x}\right)-A\left[{\partial }_t\psi +\overline{u_x}{\partial }_x\psi \right],$

now, using the definition of $\psi$ in (33), we get,

${\partial }_t\left(A\psi \right)+{\partial }_x\left(A\psi \overline{u_x}\right)-A\left[\overline{u_x}\left({\partial }_t\overline{u_x}+{\partial }_x\psi \right)+{\partial }_{t}p\right].$

We use the momentum Equation (32) to get,

${\partial }_t\left(A\psi -\tilde{p}\right)+{\partial }_x\left(A\psi \overline{u_x}\right)=3\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\overline{u_x}\right)\overline{u_x}+2\pi Rk{\overline{u_x}}^2,$

where,

$\tilde{p}={\displaystyle \int A{\partial }_{A}p}=\frac{b\sqrt{\pi }}{3A_0}\left(A^{\frac{3}{2}}-A_0^{\frac{3}{2}}\right).$

We define the mathematical entropy $E=A\psi -\tilde{p}$ and remark

${\partial }_{t}E+{\partial }_x\left(\left(E+\tilde{p}\right)\overline{u_x}\right)=3\epsilon {\nu }_0{\partial }_x\left(A{\partial }_x\overline{u_x}\right)\overline{u_x}+2\pi Rk{\overline{u_x}}^2.$

4) Integrating (34) between 0 and $L$ , we get,

${\partial }_t{\displaystyle {\int }_0^{L}E\text{d}x}+{\left[\left(E+\tilde{p}\right)\overline{u_x}\right]}_0^L=3\epsilon {\nu }_0{\displaystyle {\int }_0^L{\partial }_x\left(A{\partial }_x\overline{u_x}\right)\overline{u_x}}+2\pi Rk{\displaystyle {\int }_0^L{\overline{u_x}}^2},$

which, when we suppose $\overline{u_x}\left(x=0\right)=\overline{u_x}\left(x=L\right)=0$ and use integration by parts leads to,

$\frac{\text{d}}{\text{d}t}{E}_{tot}\left(t\right)=-3\epsilon {\nu }_0{\displaystyle {\int }_0^{L}A{\left({\partial }_x\overline{u_x}\right)}^2}+\frac{2\pi Rk}{1-\epsilon \frac{Rk}{4{\nu }_0}}{\Vert \overline{u_x}\Vert }_{L^2\left(\left[0,L\right]\right)}^2\le 0.$

□

Remark. Although the velocity profile remains parabolic, the second-order asymptotics enforce a Boussinesq coefficient $\alpha =1$ . This may seem restrictive, but it is a direct consequence of the modeling and is required to ensure energetic consistency with the full Navier–Stokes equations. Moreover, in contrast with classical 1D models where wall friction is prescribed through an empirical coefficient, here the velocity profile depends explicitly on friction, which enhances the physical consistency of the model. Finally, we stress that one-dimensional models are intrinsically limited when applied to highly complex arterial geometries, such as the cerebral circulation. In these cases, multidimensional models (2D or 3D) are better suited to capture flow separation, recirculation, and strong curvature effects.

3. A Discontinuous Galerkin Method for Convection-Diffusion Equations

In this section, we present a discontinuous Galerkin method to solve a general one-dimensional convection-diffusion system of equations using the IIPG [16] [20]-[22] and the RKDG methods [17] [23]-[25].

3.1. Model Problem

Our aim is to construct a high-order numerical method for the blood flow problem (38) derived in Section 2. To this purpose, we propose a Discontinuous Galerkin approach for one-dimensional convection-diffusion problem following [16] [17] [21]. Let us consider the following non-linear one-dimensional convection-diffusion problem:

$\{\begin{array}{l}{\partial }_{t}u+{\partial }_{x}f-{\partial }_x\left(g{\partial }_{x}u\right)=s,\forall \left(t,x\right)\in \left]0,T\right]\times \left]a,b\right[,\\ u\left(t=0,x\right)=u_0\left(x\right),\forall x\in \left[a,b\right]\\ u\left(t,x=a\right)=u_a\left(t\right),\forall t\in \left[0,T\right]\\ u\left(t,x=b\right)=u_b\left(t\right),\forall t\in \left[0,T\right]\end{array}$ (39)

where $\left(u_a\left(t\right),u_b\left(t\right)\right)\in {ℝ}^d\times {ℝ}^d$ are Dirichlet boundary conditions and $u_0\left(x\right)\in {ℝ}^d$ is the initial condition. Here, $f\left(t,x,u\right)\in {ℝ}^d$ is the convection component, $g\left(t,x,u\right)\in {ℝ}^{d\times d}$ is the diffusive component and $s\left(t,x,u,{\partial }_{x}u\right)\in {ℝ}^d$ is the source term, where $d=1,2$ or 3.

Remark. For Neumann boundary conditions, one can refer to [16] and [17] for more details.

3.2. Space Discretization

Let $a=x_0<x_1<\cdots <x_N=b$ be a partition of the interval $\left[a,b\right]$ and $I_n=\left(x_n,x_{n+1}\right)$ where we define $h_n=x_{n+1}-x_n$ . For the sake of simplicity, in what follows, we consider $h=h_n$ , $\forall n\in 〚0,N-1〛$ . The space of piecewise discontinuous polynomials of degree $p$ is

$V_p\left(\left[a,b\right]\right)=\left\{v:{v|}_{I_n}\in {\mathbb{P}}_p\left(I_n\right)\forall n\in 0,\cdots ,N-1\right\},$

where ${\mathbb{P}}_p\left(I_n\right)$ is the space of polynomials of degree $p$ on the interval $I_n$ . Defining $v\left(x_n^{+}\right)={\lim }_{\begin{array}{l}\epsilon \to 0\\ \epsilon >0\end{array}}v\left(x_n+\epsilon \right)$ and $v\left(x_n^{-}\right)={\lim }_{\begin{array}{l}\epsilon \to 0\\ \epsilon >0\end{array}}v\left(x_n-\epsilon \right)$ , the jump and the average of $v$ at the endpoints of $I_n$ , $\forall n\in 〚0,N-1〛$ , are defined by:

$\left[v\left(x_n\right)\right]=v\left(x_n^{-}\right)-v\left(x_n^{+}\right),\left\{v\left(x_n\right)\right\}=\frac{1}{2}\left(v\left(x_n^{-}\right)+v\left(x_n^{+}\right)\right).$

By convention, we also extend the definition of jump and average at the endpoints of the interval $\left[a,b\right]$ by

$\begin{array}{lll}\left[v\left(x_a\right)\right]=-v\left(x_a^{+}\right),\hfill & \hfill & \left\{v\left(x_a\right)\right\}=v\left(x_a^{+}\right),\hfill \\ \left[v\left(x_b\right)\right]=v\left(x_b^{-}\right),\hfill & \hfill & \left\{v\left(x_b\right)\right\}=v\left(x_b^{-}\right).\hfill \end{array}$ (40)

For convenience, we define the numerical flux functions

$\begin{array}{l}\widehat{f}\left(t,x,v\right)=\left\{f\left(t,x,v\right)\right\}+\frac{c}{2}\left[v\right],\\ \widehat{g}\left(t,x,v,{\partial }_{x}v\right)=\left\{g\left(t,x,v\right){\partial }_{x}v\right\}-\frac{\sigma }{h}\left[g\left(t,x,v\right)v\right],\end{array}$ (41)

where $c={\max }_{i=1,\cdots ,d}\left|{\lambda }_i\left(t,x_n,u\right)\right|$ , with ${\lambda }_i$ , the eigenvalues of ${\nabla }_{u}f$ , and $\sigma >0$ a penalty parameter. The numerical flux $\widehat{f}\left(t,x,v\right)$ used for the advective component is the classical Rusanov flux [17] while $\widehat{g}\left(t,x,v,{\partial }_{x}v\right)$ yields to the well-known IIPG (Interior Incomplete Penalty Galerkin) [16].

3.3. Weak Problem

To obtain, the weak formulation of Equations (39), we multiply by $v\in V_p\left(\left[a,b\right]\right)$ and we integrate by parts on each interval $I_n$ , to get

$\begin{array}{l}{\displaystyle {\int }_{x_n}^{x_{n+1}}{\partial }_{t}uv}-{\displaystyle {\int }_{x_n}^{x_{n+1}}f}\left(t,x,u\right){\partial }_{x}v+f\left(t,x_{n+1},u\right)v\left(x_{n+1}^{-}\right)-f\left(t,x_n,u\right)v\left(x_n^{+}\right)\\ +{\displaystyle {\int }_{x_n}^{x_{n+1}}g}\left(t,x,u\right){\partial }_{x}u{\partial }_{x}v-g\left(t,x_{n+1},u\right){\partial }_{x}uv\left(x_{n+1}^{-}\right)+g\left(t,x_n,u\right){\partial }_{x}uv\left(x_n^{+}\right)\\ ={\displaystyle {\int }_{x_n}^{x_{n+1}}s}\left(t,x,u,{\partial }_{x}u\right)v.\end{array}$

By adding all $N$ equations above, we obtain

$\begin{array}{l}{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}{\partial }_{t}uv}}-{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}f}}\left(t,x,u\right){\partial }_{x}v+{\displaystyle {\sum }_{n=0}^N\left[f\left(t,x_n,u\right)v\left(x_n\right)\right]}\\ +{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}g\left(t,x,u\right){\partial }_{x}u{\partial }_{x}v}}-{\displaystyle {\sum }_{n=0}^N\left[g\left(t,x_n,u\right){\partial }_{x}uv\left(x_n\right)\right]}\\ ={\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}s}}\left(t,x,u,{\partial }_{x}u\right)v,\end{array}$ (42)

where we used convention (40). The discrete problem is finding $u$ solution of (42) for any $v$ with proper boundary conditions from (39). For the sake of simplicity, in what follows, we consider $d=1$ . The generalization to $d>1$ can be easily done and is left to the reader.

ODE System

Consider the following monomial basis functions of ${\mathbb{P}}_p\left(I_n\right)$ reading $\forall i\in 0,p$ ,

${\varphi }_i^n\left(x\right)={\left[\frac{2}{h}\left(x-x_{n+\frac{1}{2}}\right)\right]}^i,$ (43)

where $x_{n+\frac{1}{2}}=\frac{1}{2}\left(x_n+x_{n+1}\right)$ . Their derivative is $\frac{\text{d}}{\text{d}x}{\varphi }_i^n\left(x\right)=\frac{2}{h}i{\left[\frac{2}{h}\left(x-x_{n+\frac{1}{2}}\right)\right]}^{i-1}$ . We look for an approximate solution $u_{DG}$ of (42) in $V_p\left(\left[a,b\right]\right)$ that we write, in the basis (43),

$u_{DG}\left(t,x\right)={\displaystyle \sum }_{m=0}^{N-1}{\displaystyle \sum }_{j=0}^p{\varphi }_j^m\left(x\right)U_j^m\left(t\right)1{\text{l}}_{I_n}\left(x\right),$ (44)

where the coefficients $U_j^m$ are unknown time-dependent functions and $1{\text{l}}_{I_n}\left(x\right)$ is the characteristic function of the set $I_n$ .

We use the following notation for boundary conditions:

$\left[u_{DG}\left(x_a\right)\right]=u_a\left(t\right)-u_{DG}\left(x_a^{+}\right),\left[u_{DG}\left(x_b\right)\right]=u_{DG}\left(x_b^{-}\right)-u_b\left(t\right),$

where $u_a$ and $u_b$ are the boundary conditions defined in (39), while, for $v$ , the conventions described in (40) are used.

Then, considering the approximations

$\left[f\left(t,x_n,u_{DG}\right)v\left(x_n\right)\right]\approx \widehat{f}\left(t,x_n,u_{DG}\right)\left[v\left(x_n\right)\right],$

and,

$\left[g\left(t,x_n,u_{DG}\right){\partial }_x{u}_{DG}v\left(x_n\right)\right]\approx \widehat{g}\left(t,x_n,u_{DG},{\partial }_x{u}_{DG}\right)\left[v\left(x_n\right)\right],$

in (42) where $\widehat{f}$ and $\widehat{g}$ are the numerical flux defined in (41), we look for $u_{DG}$ solution of

$\begin{array}{l}{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}{\partial }_t{u}_{DG}v}}-{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}f}}\left(t,x,u_{DG}\right){\partial }_{x}v+{\displaystyle {\sum }_{n=0}^N\widehat{f}}\left(t,x_n,u_{DG}\right)\left[v\left(x_n\right)\right]\\ +{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}g}}\left(t,x,u_{DG}\right){\partial }_x{u}_{DG}{\partial }_{x}v-{\displaystyle {\sum }_{n=0}^N\widehat{g}}\left(t,x_n,u_{DG},{\partial }_x{u}_{DG}\right)\left[v\left(x_n\right)\right]\\ ={\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}s}}\left(t,x,u_{DG},{\partial }_x{u}_{DG}\right)v.\end{array}$(45)

including the boundary conditions of (39) in $\widehat{g}\left(t,x_k,u_{DG},{\partial }_x{u}_{DG}\right)$ and $\widehat{f}\left(t,x_k,u_{DG}\right)$ for $k=0$ and $k=N$ .

Keeping in mind the definition of $u_{DG}$ (44), Equation (45) becomes, for $i\in 0,p$ ,

$\begin{array}{l}{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\sum }_{j=0}^p\left({\displaystyle {\int }_{x_n}^{x_{n+1}}{\varphi }_j^n\left(x\right){\varphi }_i^n\left(x\right)}\right){\partial }_t{U}_j^n\left(t\right)}}-{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}f\left(t,x,u_{DG}\right){\partial }_x{\varphi }_i^n}}\\ +{\displaystyle {\sum }_{n=0}^N\widehat{f}\left(t,x_n,u_{DG}\right)\left[v_i\left(x_n\right)\right]}+{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}g\left(t,x,u_{DG}\right){\partial }_x{u}_{DG}{\partial }_x{\varphi }_i^n\left(x\right)}}\\ -{\displaystyle {\sum }_{n=0}^N\widehat{g}\left(t,x_n,u_{DG},{\partial }_x{u}_{DG}\right)\left[v_i\left(x_n\right)\right]}\\ ={\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}s}}\left(t,x,u_{DG},{\partial }_x{u}_{DG}\right){\varphi }_i^n\left(x\right)\end{array}$(46)

where for $1\le n\le N$ ,

$\left[v_i\left(x_n\right)\right]=v_i\left(x_n^{-}\right)-v_i\left(x_n^{+}\right)={\varphi }_i^{n-1}\left(x_n\right)-{\varphi }_i^n\left(x_n\right),$

$\left[v_i\left(x_0=x_a\right)\right]=-v_i\left(x_a^{+}\right)=-{\varphi }_i^0\left(x_a\right)\text{and}\left[v_i\left(x_N=x_b\right)\right]=v_i\left(x_b^{-}\right)={\varphi }_i^{N-1}\left(x_b\right).$

Using those in (46), we obtain

$\begin{array}{l}{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\sum }_{j=0}^p\left({\displaystyle {\int }_{x_n}^{x_{n+1}}{\varphi }_j^n\left(x\right){\varphi }_i^n\left(x\right)}\right){\partial }_t{U}_j^n\left(t\right)}}-{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}f}}\left(t,x,u_{DG}\right){\partial }_x{\varphi }_i^n\\ +{\displaystyle {\sum }_{n=1}^N\widehat{f}}\left(t,x_n,u_{DG}\right){\varphi }_i^{n-1}\left(x_n\right)-{\displaystyle {\sum }_{n=1}^N\widehat{f}}\left(t,x_n,u_{DG}\right){\varphi }_i^n\left(x_n\right)\\ +{\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}g}}\left(t,x,u_{DG}\right){\partial }_x{u}_{DG}{\partial }_x{\varphi }_i^n\left(x\right)-{\displaystyle {\sum }_{n=1}^N\widehat{g}}\left(t,x_n,u_{DG},{\partial }_x{u}_{DG}\right){\varphi }_i^{n-1}\left(x_n\right)\\ +{\displaystyle {\sum }_{n=0}^{N-1}\widehat{g}}\left(t,x_n,u_{DG},{\partial }_x{u}_{DG}\right){\varphi }_i^n\left(x_n\right)\\ ={\displaystyle {\sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}s}}\left(t,x,u_{DG},{\partial }_x{u}_{DG}\right){\varphi }_i^n\left(x\right),\end{array}$

which we rewrite as,

${\displaystyle {\sum }_{n=0}^{N-1}{M}_{ij}^n{\partial }_t{U}_j^n\left(t\right)}-{\displaystyle {\sum }_{n=0}^{N-1}{F}_i^n\left(t,u_{DG}\right)}+{\displaystyle {\sum }_{n=0}^{N-1}{G}_i^n\left(t,u_{DG}\right)}={\displaystyle {\sum }_{n=0}^{N-1}{S}_i^n}\left(t,u_{DG}\right),$

where $M_{ij}^n={\displaystyle {\int }_{x_n}^{x_{n+1}}{\varphi }_j^n\left(x\right){\varphi }_i^n\left(x\right)}$ ,

$\begin{array}{c}{F}_i^n\left(t,u_{DG}\right)={\displaystyle {\int }_{x_n}^{x_{n+1}}f}\left(t,x,u_{DG}\right){\partial }_x{\varphi }_i^n-\widehat{f}\left(t,x_{n+1},u_{DG}\right){\varphi }_i^n\left(x_{n+1}\right)\\ +\widehat{f}\left(t,x_n,u_{DG}\right){\varphi }_i^n\left(x_n\right),\end{array}$(47)

$\begin{array}{c}{G}_i^n\left(t,u_{DG}\right)={\displaystyle {\int }_{x_n}^{x_{n+1}}g}\left(t,x,u_{DG}\right){\partial }_x{u}_{DG}{\partial }_x{\varphi }_i^n\left(x\right)-\widehat{g}\left(t,x_{n+1},u_{DG},{\partial }_x{u}_{DG}\right){\varphi }_i^n\left(x_{n+1}\right)\\ +\widehat{g}\left(t,x_n,u_{DG},{\partial }_x{u}_{DG}\right){\varphi }_i^n\left(x_n\right),\end{array}$ (48)

and,

$S_i^n\left(t,u_{DG},{\partial }_x{u}_{DG}\right)={\displaystyle {\int }_{x_n}^{x_{n+1}}s}\left(t,x,u_{DG},{\partial }_x{u}_{DG}\right){\varphi }_i^n\left(x\right).$ (49)

Finally, we define,

$\begin{array}{l}F\left(t,u_{DG}\right)={\left(F_i^n\left(t,u_{DG}\right)\right)}_{i\in 〚0,p〛,n\in 〚0,N-1〛},\hfill \\ S\left(t,u_{DG},{\partial }_x{u}_{DG}\right)={\left(S_i^n\left(t,u_{DG},{\partial }_x{u}_{DG}\right)\right)}_{i\in 〚0,p〛,n\in 〚0,N-1〛},\hfill \\ G\left(t,u_{DG}\right)={\left(G_i^n\left(t,u_{DG}\right)\right)}_{i\in 〚0,p〛,n\in 〚0,N-1〛},\hfill \\ M={\left(M_{ij}^n\right)}_{i\in 〚0,p〛,j\in 〚0,p〛,n\in 〚0,N-1〛,m\in 〚0,N-1〛},\hfill \\ U={\left(U_j^n\right)}_{i\in 〚0,p〛,n\in 〚0,N-1〛}.\hfill \end{array}$ (50)

This leads to the following system of ODE’s,

$M{\partial }_{t}U=F\left(t,U\right)+S\left(t,U\right)-G\left(t,U\right).$ (51)

Remark.

1) All integrals are computed using $2p+1$ Gauss-Legendre quadrature points.

2) The initial solution $u_0$ of problem (39) is projected to get $U_0$ , i.e.,

${\displaystyle \sum }_{n=0}^{N-1}{\displaystyle \sum }_{j=0}^p{\displaystyle {\int }_{x_n}^{x_{n+1}}{\left(U_0\right)}_j^n{\varphi }_j\left(x\right){\varphi }_i\left(x\right)}={\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{x_n}^{x_{n+1}}{u}_0}\left(x\right){\varphi }_i\left(x\right),$

for all $i$ in $〚0,p〛$ , this leads to a linear problem, $M{U}_0={\left({\displaystyle {\int }_{x_n}^{x_{n+1}}{u}_0\left(x\right){\varphi }_i\left(x\right)}\right)}_{i\in 〚0,p〛,n\in 〚0,N-1〛}.$

4. Time Discretization

In this section, we propose a time discretization based on the Additive Runge-Kutta (ARK) methods applied to the ODE system (51) as done in [26]-[28] for instance.

ARK Method

We recall the additive ARK methods derived in [26]. ARK methods are used to solve equations of the form

$U_t={\displaystyle \sum }_{m=1}^D{F}^{\left[m\right]}\left(t,U\right),$ (52)

where $F^{\left[m\right]}$ are some functions. Let $t_k\in \left[0,T\right]$ and $\text{Δ}t$ a well-chosen step size, s-stage additive Runge-Kutta methods reads, for $i\in 1,s$ ,

$\{\begin{array}{l}{K}^{\left(i\right)}=U^{\left(k\right)}+\text{Δ}t{\displaystyle {\sum }_{m=1}^D{\displaystyle {\sum }_{j=1}^s{a}_{ij}^{\left[m\right]}{F}^{\left[m\right]}}}\left(t_k+c_j^{\left[m\right]}\text{Δ}t,K^{\left(j\right)}\right),\\ U^{\left(k+1\right)}=U^{\left(k\right)}+\text{Δ}t{\displaystyle {\sum }_{m=1}^D{\displaystyle {\sum }_{i=1}^s{b}_i^{\left[m\right]}{F}^{\left[m\right]}}}\left(t_k+c_i^{\left[m\right]}\text{Δ}t,K^{\left(i\right)}\right),\\ {\widehat{U}}^{\left(k+1\right)}=U^{\left(k\right)}+\text{Δ}t{\displaystyle {\sum }_{m=1}^D{\displaystyle {\sum }_{i=1}^s{\widehat{b}}_i^{\left[m\right]}{F}^{\left[m\right]}}}\left(t_k+c_i^{\left[m\right]}\text{Δ}t,K^{\left(i\right)}\right),\end{array}$

where $U^{\left(k\right)}\approx U\left(t_k\right)$ is the approximate solution of (52) at time $t_k$ , $U^{\left(k+1\right)}\approx U\left(t_k+\text{Δ}t\right)$ the approximate solution at time $t_{k+1}=t_k+\text{Δ}t$ at a certain order of accuracy and ${\widehat{U}}^{\left(k+1\right)}\approx U\left(t_k+\text{Δ}t\right)$ another approximate solution at time $t_{k+1}=t_k+\text{Δ}t$ at a lower order of accuracy. Each of the respective Butcher coefficients ${\left(a_{ij}^{\left[m\right]}\right)}_{i\in 〚1,s〛,j\in 〚1,s〛,m\in 〚1,D〛}$ , ${\left(b_i^{\left[m\right]}\right)}_{i\in 〚1,s〛,m\in 〚1,D〛}$ , ${\left(c_i^{\left[m\right]}\right)}_{i\in 〚1,s〛,m\in 〚1,D〛}$ and ${\left({\widehat{b}}_i^{\left[m\right]}\right)}_{i\in 〚1,s〛,m\in 〚1,D〛}$ are constrained, at a minimum, by certain order of accuracy and stability considerations discussed in [26]. $s$ is the number of steps associated with the ARK method.

Application to the ODE System

We rewrite the ODE (51) in the form of (52) leading to,

${\partial }_t\left(U\right)=M^{-1}{F}^{\left[1\right]}\left(t,U\right)+M^{-1}{F}^{\left[2\right]}\left(t,U\right),$

where

$F^{\left[1\right]}\left(t,U\right)=F\left(t,U\right)+S\left(t,U\right)\text{and}{F}^{\left[2\right]}\left(t,U\right)=-G\left(t,U\right).$ (53)

This leads to the following scheme,

$\{\begin{array}{l}M{K}^{\left(i\right)}=M{U}^{\left(k\right)}+\Delta t{\displaystyle {\sum }_{j=1}^s\left[a_{ij}^{\left[1\right]}{F}^{\left[1\right]}\left(t_k+c_j^{\left[1\right]}\Delta t,K^{\left(j\right)}\right)+a_{ij}^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_j^{\left[2\right]}\Delta t,K^{\left(j\right)}\right)\right]},\\ M{U}^{\left(k+1\right)}=M{U}^{\left(k\right)}+\Delta t{\displaystyle {\sum }_{i=1}^s\left[b_i^{\left[1\right]}{F}^{\left[1\right]}\left(t_k+c_i^{\left[1\right]}\Delta t,K^{\left(i\right)}\right)+b_i^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_i^{\left[2\right]}\Delta t,K^{\left(i\right)}\right)\right]},\\ M{\widehat{U}}^{\left(k+1\right)}=M{U}^{\left(k\right)}+\Delta t{\displaystyle {\sum }_{i=1}^s\left[{\widehat{b}}_i^{\left[1\right]}{F}^{\left[1\right]}\left(t_k+c_i^{\left[1\right]}\Delta t,K^{\left(i\right)}\right)+{\widehat{b}}_i^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_i^{\left[2\right]}\Delta t,K^{\left(i\right)}\right)\right]}.\end{array}$

We use the PID-controller from [26] to adapt $\text{Δ}t$ with an upper-bound taken from [17] as in the hyperbolic case. The maximal time step is given by

$\text{Δ}t=\frac{\text{cfl}}{2p+1}\frac{h}{c^k},$

where $\text{cfl}\in \left]0,1\right]$ , $p$ is the polynomial degree, $h$ , the spacial step size and $c$ the characteristic speed define as $c^k={\max }_{n\in 〚0,N-1〛}\left|\lambda \left(t_k,U^k,x_n\right)\right|$ , with $\lambda \left(t,u,x\right)$ the eigenvalues of ${\nabla }_{u}f\left(t,u,x\right)$ . An explicit Runge-Kutta scheme is used for $f^{\left[1\right]}$ and a diagonally implicit one for $f^{\left[2\right]}$ .

In this case, the scheme can be written as,

$\{\begin{array}{l}M{K}^{\left(i\right)}-\text{Δ}t{a}_{ii}^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_i^{\left[2\right]}\text{Δ}t,K^{\left(i\right)}\right)\\ =M{U}^{\left(k\right)}+\text{Δ}t{\displaystyle {\sum }_{j=1}^{i-1}\left[a_{ij}^{\left[1\right]}{F}^{\left[1\right]}\left(t_k+c_j^{\left[1\right]}\text{Δ}t,K^{\left(j\right)}\right)+a_{ij}^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_j^{\left[2\right]}\text{Δ}t,K^{\left(j\right)}\right)\right]},\\ M{U}^{\left(k+1\right)}=M{U}^{\left(k\right)}+\text{Δ}t{\displaystyle {\sum }_{i=1}^s\left[b_i^{\left[1\right]}{F}^{\left[1\right]}\left(t_k+c_i^{\left[1\right]}\text{Δ}t,K^{\left(i\right)}\right)+b_i^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_i^{\left[2\right]}\text{Δ}t,K^{\left(i\right)}\right)\right]},\\ M{\widehat{U}}^{\left(k+1\right)}=M{U}^{\left(k\right)}+\text{Δ}t{\displaystyle {\sum }_{i=1}^s\left[{\widehat{b}}_i^{\left[1\right]}{F}^{\left[1\right]}\left(t_k+c_i^{\left[1\right]}\text{Δ}t,K^{\left(i\right)}\right)+{\widehat{b}}_i^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_i^{\left[2\right]}\text{Δ}t,K^{\left(i\right)}\right)\right]}.\end{array}$ (54)

This means we have $s$ non-linear equations to solve. We use a fix-point method here to solve those, meaning, we define the following recursion on $r$ ,

$\begin{array}{l}M{K}_{r+1}^{\left(i\right)}-\text{Δ}t{a}_{ii}^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_i^{\left[2\right]}\text{Δ}t,K_r^{\left(i\right)}\right)\\ =M{U}^{\left(k\right)}+\text{Δ}t{\displaystyle {\sum }_{j=1}^{i-1}\left[a_{ij}^{\left[1\right]}{F}^{\left[1\right]}\left(t_k+c_j^{\left[1\right]}\text{Δ}t,K^{\left(j\right)}\right)+a_{ij}^{\left[2\right]}{F}^{\left[2\right]}\left(t_k+c_j^{\left[2\right]}\text{Δ}t,K^{\left(j\right)}\right)\right]},\end{array}$ (55)

with $K_1^{\left(i\right)}=U^{\left(k\right)}$ while $\Vert K_{r+1}^{\left(i\right)}-K_r^{\left(i\right)}\Vert >\epsilon$ for a given tolerance, $\epsilon$ (set to 10⁻¹² in practice).

We use the following time scheme depending on $p$ (see Table 1) and butcher tables which can be found in [26] for instance.

Table 1. Time schemes for different values of $p$ .

$p$	Time scheme
1	ARK3
2	ARK3
3	ARK4
4	ARK5

Well-Balanced Scheme for Still-Steady-States Solutions

We show in this section that one can easily obtain a well-balanced scheme without modifying the numerical method contrary to the well-known existing methods [29]-[31]. It can be achieved through a simple change of variables and the presented method can also be applied to the Saint-Venant equations with classical approximate Riemann Solver.

Let us introduce

$a=A-A_0.$

One can transform the system (38) into the form

${\partial }_{t}a+{\partial }_{x}Q=0,$ (56)

and,

$\begin{array}{l}{\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{a+A_0}+\left(a+A_0\right)P\right)-{\partial }_x\left(3\nu \left(a+A_0\right){\partial }_x\left(\frac{Q}{a+A_0}\right)\right)\\ =P{\partial }_x\left(a+A_0\right)+\gamma \frac{Q}{a+A_0}.\end{array}$ (57)

Then, we have the following property:

Proposition 1. The numerical scheme (54) preserves exactly the still-steady states solutions (see Theorem 1 and Theorem 2). Let us note $U^{\left(k\right)}$ be an approximation of $U\left(t_k\right)$ . If $U^{\left(0\right)}=0$ then $U^{\left(k\right)}=0$ for all $k>0$ . More precisely, the approximation of $a$ and $Q$ at time $t_k$ are 0.

Proof.

We want to prove that if ${\left(U^{\left(k\right)}\right)}_k$ is a sequence verifying Equation (54) such that $U^{\left(0\right)}=0$ , then, $U^{\left(k\right)}=0$ , $\forall k>0$ . In this proof, we justify only the case $U^{\left(0\right)}=0\Rightarrow U^{\left(1\right)}=0$ since the proof at rank $k$ , $U^{\left(k\right)}=0\Rightarrow U^{\left(k+1\right)}=0$ , is a straightforward consequence.

Let us also remark that if $u=0$ and ${\partial }_{x}u=0$ , then $F^{\left[1\right]}\left(t,0\right)=F^{\left[2\right]}\left(t,0\right)=0$ for all $t>0$ which is the case when we consider the model (56) and (57) with $f\left(t,x,u\right)=\frac{Q^2}{a+A_0}+\left(a+A_0\right)P$ , $g\left(t,x,u\right)=3\nu \left(a+A_0\right)$ , $s\left(t,x,u,{\partial }_{x}u\right)=P{\partial }_x\left(a+A_0\right)+\gamma \frac{Q}{a+A_0}$ with $u={\left(a,Q\right)}^t$ through the definition of (41), (47), (48), (49), (50) and (53).

Assume that $U^{\left(0\right)}=0$ , then in (54), for $i=1$ , we have

$M{K}^{\left(1\right)}-\text{Δ}t{a}_{11}^{\left[2\right]}{F}^{\left[2\right]}\left(t_0+c_1^{\left[2\right]}\text{Δ}t,K^{\left(1\right)}\right)=0.$

To obtain the value of $K^{\left(1\right)}$ , we perform a fixed-point iteration (55). For $r=0$ , we have $M{K}_1^{\left(1\right)}=0$ yielding to $K^{\left(1\right)}=0$ . For $i=2$ , $K^{\left(2\right)}$ satisfies

$M{K}^{\left(2\right)}-\text{Δ}t{a}_{11}^{\left[2\right]}{F}^{\left[2\right]}\left(t_1+c_1^{\left[2\right]}\text{Δ}t,K^{\left(1\right)}\right)=0$

yielding to, as $K^{\left(1\right)}$ , to $K^{\left(2\right)}=0$ . Then, recursively, we get $K^{\left(i\right)}=0$ for $i=1,\cdots ,s$ . As a consequence, the second and third equation in (54), provide $U^{\left(1\right)}=0$ and ${\widehat{U}}^{\left(1\right)}=0$ .

□

5. Numerical Test Cases

In this section, we present some test cases for the Discontinuous Galerkin (DG) method on convection-diffusion problems as presented in Section 3. We deal with a scalar case in Section 5.1 and the blood flow models (31) and (38) in Section 5.2. The main objectives of these test cases are to:

1) Examine the behavior of the numerical scheme in the presence of degenerating parabolic terms. This will involve testing the scheme’s performance when the diffusion component becomes negligible or zero.

2) Assess the numerical convergence of the DG method. We aim to evaluate the accuracy and stability of the method under various scenarios, including different grid resolutions and complexities.

3) Apply the DG method to the 1D blood flow models (31) and (38). This step will demonstrate the method’s practical applicability in simulating real-world scenarios in blood flow dynamics.

These test cases will provide insights into the efficacy and robustness of the DG method.

5.1. Convergence Order for Scalar Convection-Diffusion Equations

To test our numerical scheme, we consider two test cases. The first one, presented in Section 5.1.1 will highlight the robustness of the scheme in scenarios where the exact solution gets close to a discontinuous solution for diffusion coefficient smaller and smaller, with respect to the penalty parameters.

The second test case, presented in Section 5.1.2, however is made to make the diffusion tend to 0 in time to highlight how the scheme acts when the model has degenerate parabolicity.

5.1.1. Quasi-Discontinuous Solution

For this first test case, the following one-dimensional convection-diffusion problem is considered:

$\{\begin{array}{ll}{\partial }_{t}u+{\partial }_{x}u-\nu {\partial }_{xx}u=\frac{2}{\nu }u\left(1-u^2\right),\hfill & \hfill \\ u\left(x,0\right)=\text{tanh}\left(\frac{x}{\nu }\right)\hfill & x\in \left[0,1\right],\hfill \\ u\left(0,t\right)=\text{tanh}\left(-\frac{t}{\nu }\right)\hfill & t\in \left[0,0.5\right],\hfill \\ u\left(1,t\right)=\text{tanh}\left(\frac{1-t}{\nu }\right)\hfill & t\in \left[0,0.5\right].\hfill \end{array}$ (58)

The analytical solution is $u\left(x,t\right)=\text{tanh}\left(\frac{x-t}{\nu }\right)$ (shown in Figure 2). Following the notations from Section 3, we consider:

• $f\left(u\right)=u$ : The convection term;

• $g\left(u\right)=\nu$ : The constant diffusion coefficient;

• $s\left(u\right)=\frac{2}{\nu }u\left(1-u^2\right)$ : The source term;

where $\nu$ is a strictly positive parameter.

Numerical Solutions

In Figure 2, for $\sigma =100$ , we observe different behavior of the numerical solution when the diffusion coefficient decreases. Indeed, as the diffusion coefficient $\nu$ decreases, the solution approaches a discontinuity leading to high error in the numerical solution and the accuracy of the solution can be improved using slope limiter [17].

Figure 2. Solution of (58) with degree $p=1$ and at time $t=\frac{1}{2}$ .

In the case of $\nu =0.001$ , we observe, as shown in Figure 3 (to compare with Figure 2), the smaller $\sigma$ is, the better the results are accurate meaning that the automatization of the choice of the parameter $\sigma$ is relevant (as in [16]). Moreover, in the case of large viscosity, we will see that the parameter $\sigma$ is crucial to get accurate solution (see Figure 5).

Figure 3. Solution of (58) with degree $p=1$ and at time $t=\frac{1}{2}$ with $\nu =0.001$ and $\sigma =0$ .

Numerical Order

We proceed with the computation of the error $ϵ$ to get the numerical order of convergence where $ϵ={\Vert u_{\text{exact}}-u_{\text{DG}}\Vert }_{l^2\left(\left[0,1\right]\right)}$ with $u_{\text{exact}}$ the exact solution to the problem (58) and $u_{\text{DG}}$ the solution obtained from the numerical scheme presented in Section 3 at time $t=\frac{1}{2}$ .

In our convergence study, the numerical solution is computed for $N=16,32$ , and 64. and the results are given in Figure 4 for $\sigma =100$ . We also study the order with respect to the parameter $\sigma$ in Figure 5.

In Figure 4, we recover the usual numerical order $no$ of the IIPG method (as in [16] for instance) whenever the viscosity $\nu$ is large enough, i.e.,

$no=p+1$ if $p$ is odd

$no=p$ if $p$ is even

Whenever the viscosity is small, as shown in Figure 4, the numerical order of convergence is almost $0.5$ , i.e., we recover the classical order for convection problem in the presence of discontinuity or very sharp numerical solution.

Figure 4. Numerical convergence order for different values of $p$ and $\nu$ with $\sigma =100$ .

Figure 5. Numerical convergence order for different values of $\sigma$ and $\nu$ with $p=1$ .

Finally, to show the importance of the choice of the penalty parameter $\sigma$ , Figure 5 show the numerical order of convergence for different values of $\sigma$ and $\nu$ with $p=1$ . For a high value of $\nu$ , a high penalty parameter is required while, for the low value of $\nu$ , the penalty parameter should be small. Automatization of the choice of the parameter $\sigma$ is then relevant to get accurate numerical simulation (as in [16], done in the case of a linear problem).

5.1.2. Quasi-Degenerate Diffusion

For this second test case, the following one-dimensional convection-diffusion problem is considered:

$\{\begin{array}{ll}{\partial }_{t}u+{\partial }_x\left(\frac{u^2}{2}\right)-\nu {\text{e}}^{-t}{\partial }_{xx}u=20\left(1-u^2\right)\left(x-\frac{1}{2}+\left(1+40t\nu {\text{e}}^{-t}\right)tu\right),\hfill & \hfill \\ u\left(x,0\right)=0\hfill & \text{for}x\in \left[0,1\right],\hfill \\ u\left(0,t\right)=\tanh \left(-10t\right)\hfill & \text{for}t\in \left[0,10\right],\hfill \\ u\left(1,t\right)=\tanh \left(10t\right)\hfill & \text{for}t\in \left[0,10\right].\hfill \end{array}$ (59)

The analytical solution is $u\left(x,t\right)=\text{tanh}\left(20t\left(x-\frac{1}{2}\right)\right)$ (as shown in Figure 6). Following the notations from Section 3, the functions are defined as:

• $f\left(u\right)=\frac{u^2}{2}$ : The convection term.

• $g\left(u\right)=\nu {\text{e}}^{-t}$ : The diffusion term.

• $s\left(u\right)=20\left(1-u^2\right)\left(x-\frac{1}{2}+\left(1+40t\nu {\text{e}}^{-t}\right)tu\right)$ : The source term.

For this test case, we set $\sigma =100$ and $\nu =0.01$ . Remark that, in (59), the parabolic part vanishes exponentially, which leads to the hyperbolic part being dominant. Moreover, due to the definition of the viscous numerical flux (see (41)), the penalty term vanishes, meaning that we expect to obtain an accurate solution given the observation in Section 5.1.1.

Figure 6. Solution of (59) with degree $p=1$ .

As already mentioned in the previous test case Section 5.1.1, these results could be improved using slope limiters (see [17]). Finally, in Figure 7, we compute the numerical order of convergence and we obtain the order of the RKDG scheme (see [17]) (i.e., if we approximate the solution of problem (59) by a piecewise polynomial of order $p$ , we obtain $p+1$ order of convergence).

Figure 7. Errors for the test case (59) at $t=0.5$ .

5.2. Test Case for Blood Flow Systems

To solve the second-order (viscous) one-dimensional blood flow model (see Section 2, Equations (38)), we use the DG method presented in Section 3. Our aim is to compute the numerical order of convergence for a given exact solution and to show that, under a suitable change of variable, the numerical scheme captures exactly the still-steady states solutions. Moreover, we will also compare the non-viscous model to the viscous one.

5.2.1. Convergence Order

We consider the blood flow model

$\{\begin{array}{l}{\partial }_{t}A+{\partial }_{x}Q=0,\\ {\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{A}+AP\right)-P{\partial }_{x}A-3\nu {\partial }_x\left(A{\partial }_x\left(\frac{Q}{A}\right)\right)=K_r\frac{Q}{A}+S_Q,\end{array}$ (60)

for which the exact solution is

$A\left(x,t\right)=2+\cos \left(x\right)\sin \left(t\right),Q\left(x,t\right)=-\sin \left(x\right)\cos \left(t\right).$ (61)

In these equations, we have $K_r=2\pi \frac{Rk}{1-\frac{Rk}{4\nu }}=4\frac{22}{15}\pi \nu$ for the friction term and the source term is given by

$\begin{array}{l}{S}_Q=-2\cos \left(t\right)\cos \left(x\right)\frac{Q}{A}-\frac{\beta \sin \left(x\right)\sin \left(t\right)}{2A_0}\sqrt{A}+\frac{KrQ}{A}+\sin \left(x\right)\sin \left(t\right)\\ +\sin \left(x\right)\sin \left(t\right)\frac{Q^2}{A^2}-3\nu \left(-Q+2\frac{\cos \left(x\right)\sin \left(t\right)Q}{A}+\frac{Q^3}{A^2}{\tan }^2\left(t\right)\right),\end{array}$

where $P=\beta \frac{\sqrt{A}-\sqrt{A_0}}{A_0}$ is the pressure.

For this test case, we consider the following initial and boundary conditions:

$\begin{array}{ll}A\left(x,0\right)=A_0,\hfill & Q\left(x,0\right)=-\sin \left(x\right),\hfill \\ A\left(0,t\right)=2+\sin \left(t\right),\hfill & Q\left(0,t\right)=0,\hfill \\ A\left(2\pi ,t\right)=2+\sin \left(t\right),\hfill & Q\left(2\pi ,t\right)=0,\hfill \end{array}$

and the following physical and numerical parameters in Table 2.

Table 2. Physical and numerical parameters.

Parameter	Value
$\beta$	1
$A_0$	2
$\nu$	0 or 1
$T$	0.5
$L$	2π
$\sigma$	100
$\text{cfl}$	0.5

We compute the error $ϵ={\Vert u_{\text{exact}}-u_{\text{DG}}\Vert }_{l^2\left(\left[0,1\right]\right)}$ with $u_{\text{exact}}$ the exact solution (61) and $u_{\text{DG}}$ the solution obtained from the numerical scheme (presented in Section 3) at time $t=T$ .

Non-viscous case $\nu =0$ (see Equation (60)): In Figure 8, we show the pressure and speed for different values of $N$ (4, 8, 16, 32) in the non-viscous case. In Figure 9, we compute the numerical order of convergence for the pressure and speed for different polynomial degree $p$ . We obtain numerical order between $p+\frac{1}{2}$ and $p+1$ for both the pressure and the speed. As expected, those results are similar to those obtained in the RKDG case [17].

Figure 8. Convergence study for the exact solution 61 in the non-viscous case at time $t=T$ .

Figure 9. Convergence order in the non-viscous case.

Viscous case $\nu >0$ (see Equation (60)): As done before, in Figure 10, we show the pressure and speed for different values of $N$ (4, 8, 16, 32) for the viscous model (60). In Figure 11, we obtain numerical order between $p+\frac{1}{2}$ and $p+1$ when $p$ is odd and from $p-\frac{1}{2}$ to $p+\frac{1}{2}$ when $p$ is even for both the pressure and the speed. Those results are similar to those obtained in the IIPG case [16].

On average, the non-linear solver required between 2 and 4 Newton iterations per stage. The fixed-point iterations were similarly efficient. The tolerance 10⁻¹² ensured robust convergence, while the CFL restriction limited excessive time-step growth. Overall, these settings balanced stability with efficiency, leading to competitive runtimes.

Figure 10. Convergence study for the exact solution (61) in the viscous case at time $t=T$ .

Figure 11. Convergence order in the viscous case.

5.2.2. Still-Steady States Solutions

As done in subsection of Section 4, one can transform the previous system (60) into the equivalent form

${\partial }_{t}a+{\partial }_{x}Q=0,$

and,

$\begin{array}{l}{\partial }_{t}Q+{\partial }_x\left(\frac{Q^2}{a+A_0}+\left(a+A_0\right)P\right)-{\partial }_x\left(3\nu \left(a+A_0\right){\partial }_x\left(\frac{Q}{a+A_0}\right)\right)\\ =P{\partial }_x\left(a+A_0\right)+\gamma \frac{Q}{a+A_0},\end{array}$

where

$a=A-A_0.$

The main advantage of this formulation is that, as shown in subsection of Section 4, we capture exactly the still-steady states solutions. To show the ability of the scheme to capture those states, we consider the following initial and boundary conditions:

$\begin{array}{ll}A\left(x,0\right)=A_0\left(x\right),\hfill & Q\left(x,0\right)=0,\hfill \\ A\left(0,t\right)=A_0\left(0\right),\hfill & Q\left(0,t\right)=0,\hfill \\ A\left(L,t\right)=A_0\left(L\right),\hfill & Q\left(L,t\right)=0\hfill \end{array}$

where $A_0\left(x\right)=\pi {\left(R_0+\left(R_1-R_0\right){\text{e}}^{-2{\left(x-\frac{L}{2}\right)}^2}\right)}^2$ and

$\beta \left(x\right)=\left[E_0+\frac{E_1-E_0}{2}\left(\tanh \left(10\left(x-5\right)\right)-\tanh \left(10\left(x-10\right)\right)\right)\right]\frac{h\sqrt{\pi }}{\rho \left(1-{\xi }^2\right)}.$

From a physical viewpoint, $A_0$ represents a stenosis (as represented in Figure 12) and $\beta$ represents the elasticity distribution within the artery. We have used the physical and numerical parameters in Table 3.

Table 3. Physical and numerical parameters for the 1D steady-state solutions.

Parameter	Value	Unit	Description
$E_0$	$3\times {10}^6$	kg∙cm−1∙s−2	Minimum young modulus
$E_1$	$3\times {10}^8$	kg∙cm−1∙s−2	Maximum young modulus
$h$	0.05	cm	Thickness
$\rho$	1	kg∙cm−3	Density
$\xi$	0.0		Poisson coefficient
$R_0$	0.5	cm	Maximum radius
$R_1$	0.3	cm	Minimum radius
$\nu$	0.03	cm2∙s−1	Kinematic viscosity
$T$	0.25	s	Final time
$L$	15	cm	Artery’s length
$\sigma$	1		Penalty parameter
$\text{cfl}$	0.5		Cfl
$N$	32		Number of sub-intervals

As expected and as shown in Figure 13, the still-steady states are exactly preserved.

Figure 12. Artery Geometry.

Figure 13. Steady state solutions at $t=T=0.25$ .

6. Concluding Remarks and Perspectives

In this work, we have presented the derivation of a new one-dimensional model for blood flows, including accounting for additional viscous terms in comparison with existing models in the literature. The model is derived from the axisymmetric assumption, for which the mean axis is assumed straight and the effects of the curvature and torsion are neglected. In the case of the brain arteries, those effects play an important role in the fluid dynamics of the blood and cannot be omitted and more accurate models are required. This is the main topic of a forthcoming work. Furthermore, we have presented a numerical method based on a Discontinuous Galerkin method (ARK-IIPG-RKDG) to show its robustness and accuracy. However, as done in [16], the automatization of the penalty parameter is crucial and the construction of such a method is still in progress in the case of convection-diffusion problems. It is the topic of a forthcoming paper. Finally, we stress that one-dimensional models are intrinsically limited when applied to highly complex arterial geometries, such as the cerebral circulation. In these cases, multidimensional models (2D or 3D) are better suited to capture flow separation, recirculation, and strong curvature effects.

Acknowledgements

This work has been supported by the ADEN-MED project (Adaptability to Extreme Events and Natural Risks-Application to the Mediterranean and Djibouti), funded by the Région Sud Provence-Alpes-Côte d’Azur under the AAP MEDCLIMAT “Natural Risks and Food Sovereignty”.

NOTES

¹Contrary to , where they use Dirichlet boundary conditions.

²A Julia implementation of this code, written by Y. Mannes and M. Ersoy, is freely available on request, see also https://julialang.org/.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ersoy, M. (2015) Dimension Reduction for Incompressible Pipe and Open Channel Flow Including Friction. Conference Applications of Mathematics 2015, Prague, November 2015, 17-33.
[2]	Ersoy, M., Lakkis, O. and Townsend, P. (2020) A Saint-Venant Model for Overland Flows with Precipitation and Recharge. Mathematical and Computational Applications, 26, Article 1.[CrossRef]
[3]	Gerbeau, J. and Perthame, B. (2001) Derivation of Viscous Saint-Venant System for Laminar Shallow Water; Numerical Validation. Discrete & Continuous Dynamical Systems—B, 1, 89-102.[CrossRef]
[4]	Formaggia, L., Nobile, F. and Quarteroni, A. (2002) A One Dimensional Model for Blood Flow: Application to Vascular Prosthesis. In: Babuška, I., Ciarlet, P.G. and Miyoshi, T., Eds., Mathematical Modeling and Numerical Simulation in Continuum Mechanics, Springer, 137-153.[CrossRef]
[5]	Formaggia, L., Gerbeau, J.F., Nobile, F. and Quarteroni, A. (2001) On the Coupling of 3D and 1D Navier-Stokes Equations for Flow Problems in Compliant Vessels. Computer Methods in Applied Mechanics and Engineering, 191, 561-582.[CrossRef]
[6]	Čanić, S. and Kim, E.H. (2003) Mathematical Analysis of the Quasilinear Effects in a Hyperbolic Model Blood Flow through Compliant Axi-Symmetric Vessels. Mathematical Methods in the Applied Sciences, 26, 1161-1186.[CrossRef]
[7]	Čanić, S. (2002) Blood Flow through Compliant Vessels after Endovascular Repair: Wall Deformations Induced by the Discontinuous Wall Properties. Computing and Visualization in Science, 4, 147-155.[CrossRef]
[8]	Berntsson, F., Ghosh, A., Kozlov, V.A. and Nazarov, S.A. (2018) A One Dimensional Model of Blood Flow through a Curvilinear Artery. Applied Mathematical Modelling, 63, 633-643.[CrossRef]
[9]	Quarteroni, A. and Formaggia, L. (2004) Mathematical Modelling and Numerical Simulation of the Cardiovascular System. Handbook of Numerical Analysis, 12, 3-127.[CrossRef]
[10]	Ventre, J. (2020) Reduced-Order Models for Blood Flow in the Large Arteries: Applications to Cardiovascular Pathologies. Ph.D. Thesis, Sorbonne Université.
[11]	Krivovichev, G.V. (2022) Computational Analysis of One-Dimensional Models for Simulation of Blood Flow in Vascular Networks. Journal of Computational Science, 62, Article ID: 101705.[CrossRef]
[12]	Azer, K. and Peskin, C.S. (2007) A One-Dimensional Model of Blood Flow in Arteries with Friction and Convection Based on the Womersley Velocity Profile. Cardiovascular Engineering, 7, 51-73.[CrossRef] [PubMed]
[13]	Nardinocchi, P., Pontrelli, G. and Teresi, L. (2005) A One-Dimensional Model for Blood Flow in Prestressed Vessels. European Journal of Mechanics—A/Solids, 24, 23-33.[CrossRef]
[14]	Huo, Y. and Kassab, G.S. (2007) A Hybrid One-Dimensional/Womersley Model of Pulsatile Blood Flow in the Entire Coronary Arterial Tree. American Journal of Physiology-Heart and Circulatory Physiology, 292, H2623-H2633.[CrossRef] [PubMed]
[15]	Fung, Y.C. (2013) Biomechanics: Mechanical Properties of Living Tissues. Springer.
[16]	Rivière, B. (2008) Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations. Society for Industrial and Applied Mathematics.[CrossRef]
[17]	Cockburn, B. (2003) Discontinuous Galerkin Methods. ZAMM—Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 83, 731-754.[CrossRef]
[18]	Ersoy, M. (2016) Dimension Reduction for Compressible Pipe Flows Including Friction. Asymptotic Analysis, 98, 237-255.[CrossRef]
[19]	Debyaoui, M.A. and Ersoy, M. (2021) A Generalised Serre-Green-Naghdi Equations for Variable Rectangular Open Channel Hydraulics and Its Finite Volume Approximation. In: Muñoz-Ruiz, M.L., Parés, C. and Russo, G., Eds., Recent Advances in Numerical Methods for Hyperbolic PDE Systems, Springer, 251-268.[CrossRef]
[20]	Dolejší, V. (2008) Incomplete Interior Penalty Galerkin Method for a Nonlinear Convection-Diffusion Equation. In: Kunisch, K., Of, G. and Steinbach, O., Eds., Numerical Mathematics and Advanced Applications, Springer, 457-464.[CrossRef]
[21]	Clément, J.B. (2021) Numerical Simulation of Flows in Unsaturated Porous Media by an Adaptive Discontinuous Galerkin Method: Application to Sandy Beaches. Ph.D. Thesis, Université de Toulon.
[22]	Poussel, C., Ersoy, M., Golay, F. and Mannes, Y. (2023) Runge-Kutta Discontinuous Galerkin Method Applied to Shallow Water Equations. Topical Problems of Fluid Mechanics 2023, Prague, 22-24 February 2023, 152-159.[CrossRef]
[23]	Cockburn, B. and Shu, C. (2001) Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems. Journal of Scientific Computing, 16, 173-261.[CrossRef]
[24]	Qiu, J., Khoo, B.C. and Shu, C. (2006) A Numerical Study for the Performance of the Runge-Kutta Discontinuous Galerkin Method Based on Different Numerical Fluxes. Journal of Computational Physics, 212, 540-565.[CrossRef]
[25]	Qiu., J. and Shu, C. (2005) Runge-Kutta Discontinuous Galerkin Method Using WENO Limiters. SIAM Journal on Scientific Computing, 26, 907-929.[CrossRef]
[26]	Kennedy, C.A. and Carpenter, M.H. (2003) Additive Runge-Kutta Schemes for Convection-Diffusion-Reaction Equations. Applied Numerical Mathematics, 44, 139-181.[CrossRef]
[27]	Xia, Y., Xu, Y. and Shu, C. (2007) Efficient Time Discretization for Local Discontinuous Galerkin Methods. Discrete & Continuous Dynamical Systems—B, 8, 677-693.[CrossRef]
[28]	Li, L. and Wu, S. (2020) A Hybrid Time Integration Scheme for the Discontinuous Galerkin Discretizations of Convection-Dominated Problems. Energies, 13, Article 1870.[CrossRef]
[29]	Perthame, B. and Simeoni, C. (2001) A Kinetic Scheme for the Saint-Venant System with a Source Term. Calcolo, 38, 201-231.[CrossRef]
[30]	Bourdarias, C., Ersoy, M. and Gerbi, S. (2014) Unsteady Mixed Flows in Non Uniform Closed Water Pipes: A Full Kinetic Approach. Numerische Mathematik, 128, 217-263.[CrossRef]
[31]	Greenberg, J.M. and Leroux, A.Y. (1996) A Well-Balanced Scheme for the Numerical Processing of Source Terms in Hyperbolic Equations. SIAM Journal on Numerical Analysis, 33, 1-16.[CrossRef]