On the Effects of Boundary Conditions in One-Dimensional Models of Hemodynamics

Abstract

The paper is devoted to the theoretical analysis of the effects of boundary conditions on the solutions of the system of one-dimensional (1D) hemodynamics. The integral inequalities, which realize the energy inequalities for the solutions of initial-boundary-value problems, are obtained. It is demonstrated that the unphysical unbounded solutions can take place for the case of bounded functions from boundary conditions. For the periodic boundary conditions, the integral estimation illustrates the correct behavior of the solution. For this case of boundary conditions, the effective Fourier method for the analytical solution is proposed. The analytical solutions, obtained by this approach, can be used for the comparison of different 1D blood-flow models. The results obtained in the paper allow for an the alternatively view of the stated boundary conditions and can explain some problems, which can arise in numerical simulations. They expand the possibilities of the application of analytical methods in the field of blood-flow simulation. The results can be useful for the specialists on blood-flow modeling.

1. Introduction

Nowadays, the mathematical models of blood flow are widely used in applied medical investigations as a powerful tool for the prediction of the results of surgeries and different vascular defects (such as stenoses and aneurysms) [1,2]. For the simulation of processes in large vascular systems, 1D models are widely used in practice [3,4]. These models are obtained by the averaging of the Navier—Stokes equations on the vessel cross-section [5,6,7]. As a result, the model based on the hyperbolic nonlinear system of 1D equations on the cross-sectional area and flow rate (or mean velocity) is constructed. These models can be coupled with 0D and 3D models [8], and the obtained results are close to the averaged solutions of 3D models [9].

From the physical viewpoint, blood is considered as an incompressible viscous fluid [10]. In most works, the viscosity of blood is described by the Newtonian model with the proper representation of the velocity profile [5,7,8,11]. However, in many works [12,13,14,15,16,17,18,19] the viscosity is ignored, and blood is considered as an inviscid fluid. This approximation can be used for the case of large arteries and veins, where the viscosity is not as important as it is in small vessels.

In recent decades, 1D models have been used for the simulation of many processes in physiology and medicine. For example, they have been applied to the simulation of system partial hepatectomy [3]; brain hemodynamics predictions [6,20]; the simulation of flow in large arteries [9] and near bifurcations [8]; the modelling of flows in systems with stenosed vessels [21]; the analysis of blood-flow effects on cerebrospinal fluid dynamics [22]; the analysis of flows in tapered vessels [23]; and the simulation of flows in arterial systems [24,25,26]. It must be noted that the results of the simulations are close to the results of the experiments (e.g., see [24,25]). Unfortunately, according to the nonlinear nature of the main equations, the different problems for such models in the general case can be solved only numerically. However, analytical methods are also actual in qualitative investigations of 1D hemodynamics because they can be used for the analysis of the qualitative properties of the solution (such as stability, dispersion, dissipation, etc.); for the comparison of different rheological models of blood; and for providing an effective tool for testing programs, which implement the algorithms of numerical methods [27].

The presented paper is devoted to the analytical investigation of possible qualitative properties of the solutions of 1D models of blood flow. It must be noted that investigations of this kind are presented in a relatively small number of works. For example, in the papers of Ashmetkov et al. [28,29], the linearized system for pressure and velocity is considered, and the method for the obtaining of analytical solutions, based on the reflection of waves, is proposed. Paquerot and Remoissenet [30] present a nonlinear model for the description of pressure-wave propagation in arteries. The model leads to the solution of the Boussinesq-type scalar equation. The one-soliton and two-soliton solutions, obtained by the Hirota method, are analyzed. Ilyin in [31] construct the model, based on the inviscid Burgers equation, and apply the method of characteristics for the solution of an initial-boundary-value problem on the semi-infinite interval in order to obtain the nonlinear pressure–velocity relations. Canic et al. [7] present the proof of the theorem, which provides the condition of the existence of a smooth global solution of a nonlinear conservative law system, and provide the conditions for the shock formation. The recent work [32] is devoted to the analysis of the well-posedness of the system of 1D equations in a network in the case of the vessels with viscoelastic properties. The conditions of the existence and uniqueness of the maximal strong solution are obtained, such as the stability estimations.

In the presented paper, we try to expand the field of application of analytical methods in the analysis of 1D blood-flow models and to demonstrate that the problems for linearized 1D system of inviscid hemodynamics can have physically irrelevant (unbounded) solutions for the physically relevant (bounded) boundary conditions. It is demonstrated that for the periodic conditions, this problem is not observed. For this case of boundary conditions, the simple and effective Fourier method is proposed.

The paper has the following structure. In Section 2, the mathematical model and problem statement are presented. In Section 3, the integral estimations of the solution in cases of different boundary conditions are obtained. Some concluding remarks are made in Section 4.

2. Model

The 1D model of the inviscid blood flow in a distensible vessel with impermeable walls is based on the averaged incompressibility condition and momentum equation [6,31]:

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial z}=0,\frac{\partial Q}{\partial t}+\frac{\partial }{\partial z}\left(\frac{Q^2}{A}\right)+\frac{A}{\rho }\frac{\partial P}{\partial z}=0,$$

(1)

where t is a time, z is a cylindrical coordinate, $A(t,z)$ is the cross-sectional area, $Q(t,z)$ is the flow rate ($Q=AU$, where U is the mean axial velocity), $P(t,z)$ is the pressure, and $\rho$ is a constant density.

System (1) is closed by the equation-of-state $P=P\left(A\right)$. In the case of the arteries, the following equation is used [6,9]:

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

(2)

where $P_{ext}$ is the external pressure; $P_d$ and $A_d$ are the diastolic pressure and cross-sectional area; $\beta =\frac{4}{3}\sqrt{\pi }Eh$, where E is the Young’s modulus; and h is the vessel wall thickness. In the paper, the case of constant values of $P_d$, $A_d$, E, and h is considered.

After the substitution of (2) into (1), the following system is obtained:

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial z}=0,\frac{\partial Q}{\partial t}+\frac{\partial }{\partial z}\left(\frac{Q^2}{A}\right)+\gamma \sqrt{A}\frac{\partial A}{\partial z}=0,$$

(3)

where $\gamma =\frac{\beta }{2\rho A_d}$.

Let the following dimensionless variables be introduced:

$$\tilde{z}=\frac{z}{L_c},\tilde{t}=\frac{t}{T_c},\tilde{A}=\frac{A}{A_c},\tilde{P}=\frac{P}{\rho U_c^2},\tilde{Q}=\frac{Q}{A_c{U}_c},$$

(4)

where $L_c$, $T_c$, $A_c$, $U_c$ are the characteristic length, time, cross-sectional area, and velocity. The value of $L_c$ is considered as a vessel length, $U_c$ can be estimated by the maximal wave speed, $T_c$ is computed from $L_c$ and $U_c$, and $A_c$ can be estimated, for example, from the diastolic cross-sectional radius. The tilde sign in (4) will be ignored in the text below.

After the substitution of (4) into (3), the following equations are obtained:

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial z}=0,\frac{\partial Q}{\partial t}+\frac{\partial }{\partial z}\left(\frac{Q^2}{A}\right)+\xi \sqrt{A}\frac{\partial A}{\partial z}=0,$$

(5)

where $\xi =\gamma \frac{\sqrt{A_c}}{U_c^2}$.

The interesting results on the solutions of the initial-boundary-value problems for (5) can be demonstrated for the linearized version of this system. For the linearization, we represent the solution of (5) as $A=\overline{A}+a(t,z)$, $Q=\overline{Q}+q(t,z)$, where $(\overline{A},\overline{Q})$ is the constant steady-state solution of (5), and $a(t,z)$, $q(t,z)$ are the small perturbations, as their first derivatives.

The linearized version of (5) is written as

$$\frac{\partial a}{\partial t}+\frac{\partial q}{\partial z}=0,\frac{\partial q}{\partial t}+d_1\frac{\partial a}{\partial z}+d_2\frac{\partial q}{\partial z}=0,$$

(6)

where

$$d_1=\xi \sqrt{\overline{A}}-\frac{{\overline{Q}}^2}{{\overline{A}}^2},d_2=\frac{2\overline{Q}}{\overline{A}}.$$

It must be noted that the speed of sound is presented as [9]: $c=\sqrt{\frac{\beta }{2\rho A_d}}{A}^{\frac{1}{4}}$. So, $c^2=\gamma \sqrt{A}$, and in the dimensionless variables it is written as $c^2=\xi \sqrt{A}$. At the same time, $\overline{U}=\overline{Q}/\overline{A}$. The blood flow in physiological conditions is realized when $\left|U\right|<c$ [33], so the coefficient $d_1$ in (6) is strictly positive.

Let the single blood vessel, corresponding to the spatial interval $z\in [0,l]$ be considered, and the following initial and boundary conditions are stated for (6):

$$a(0,z)=\phi \left(z\right),q(0,z)=\psi \left(z\right),$$

(7)

$${\delta }_{1}a(t,0)+{\delta }_{2}q(t,0)={\chi }_1\left(t\right),{\beta }_{1}a(t,l)+{\beta }_{2}q(t,l)={\chi }_2\left(t\right),$$

(8)

where ${\delta }_1$, ${\delta }_2$, ${\beta }_1$, ${\beta }_2$ are the known coefficients, and ${\delta }_1^2+{\delta }_2^2\ne 0$, ${\beta }_1^2+{\beta }_2^2\ne 0$, and $\phi \left(z\right)$, $\psi \left(z\right)$, ${\chi }_1\left(t\right)$, ${\chi }_2\left(t\right)$ are the known functions, which are smooth enough. To ensure the smoothness of solution in $C^1([0,+\infty )\times [0,l])$, the following compatibility conditions are imposed:

$${\delta }_1\phi \left(0\right)+{\delta }_2\psi \left(0\right)={\chi }_1\left(0\right),{\beta }_1\phi \left(l\right)+{\beta }_2\psi \left(l\right)={\chi }_2\left(0\right),$$

$${\delta }_2{d}_1{\phi }^{\prime }\left(0\right)+({\delta }_1+{\delta }_2{d}_2){\psi }^{\prime }\left(0\right)=-{\dot{\chi }}_1\left(0\right),{\beta }_2{d}_1{\phi }^{\prime }\left(l\right)+({\beta }_1+{\beta }_2{d}_2){\psi }^{\prime }\left(l\right)=-{\dot{\chi }}_2\left(0\right).$$

The physiologically adequate solutions of (6)–(8) are expected to be bounded at $\left\{\right(t,z)$: $t\in [0,+\infty ),z\in [0,l]\}$:

$$\underset{z\in [0,l]}{max}\left(\left|a\right(t,z\left)\right|,\left|q\right(t,z\left)\right|\right)\le C=const,t\to +\infty .$$

(9)

As it will be demonstrated in Section 3, for the bounded and physically relevant functions ${\chi }_1\left(t\right)$ and ${\chi }_2\left(t\right)$, the situation, when the realistic condition (9) is not realized, can be observed.

3. Results

In this section, the integral inequalities, which provide the energy estimates of the solutions of (6)–(8), are obtained in order to demonstrate that the condition (9) can be broken.

Let we rewrite the system (6) in the matrix-vector form:

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{A}\frac{\partial \mathbf{u}}{\partial z}=\mathbf{0},$$

where $\mathbf{u}={(a,q)}^T$. Matrix $\mathbf{A}$ has two real eigenvalues:

$${\lambda }_{1,2}=\frac{d_2\pm \sqrt{d_2^2+4d_1}}{2}.$$

As can be seen, ${\lambda }_1>0$ and ${\lambda }_2<0$. Without the loss of generality, let the case of $\overline{U}>0$ be considered, so $d_2>0$, and we obtain that ${\lambda }_1>\left|{\lambda }_2\right|$.

The matrix $\mathbf{A}$ can be presented as $\mathbf{A}=\mathbf{R}\Lambda \mathbf{L}$, where $\Lambda$ is the diagonal matrix of eigenvalues and $\mathbf{L}$, $\mathbf{R}$ are the matrices of left and right eigenvectors of $\mathbf{A}$, respectively.

Theorem 1.

Let the linear function $p(z)$ be presented as

$$p\left(z\right)=p\left(0\right)\frac{l-z}{l}+p\left(l\right)\frac{z}{l},$$

(10)

where $p\left(0\right)>0$, $p\left(l\right)>0$ and the following inequalities are realized:

$$|p\left(0\right)e_1|\le \frac{1}{2},\left|\frac{e_2}{p\left(l\right)}\right|\le \frac{1}{2},$$

(11)

where

$$e_1=\frac{{\delta }_1{R}_{12}+{\delta }_2{R}_{22}}{{\delta }_1{R}_{11}+{\delta }_2{R}_{21}},e_2=\frac{{\beta }_1{R}_{11}+{\beta }_2{R}_{21}}{{\beta }_1{R}_{12}+{\beta }_2{R}_{22}},$$

where $R_{ij}$ are the components of matrix $\mathbf{R}$.

Then, the following inequality takes place:

$$\underset{0}{\overset{l}{\int }}\left(p^2\left(z\right){\left(L_{11}a(t,z)+L_{12}q(t,z)\right)}^2+{\left(L_{21}a(t,z)+L_{22}q(t,z)\right)}^2\right)dz\le$$

$$exp(\Lambda t)\underset{0}{\overset{l}{\int }}\left(p^2\left(z\right){\left(L_{11}\phi \left(z\right)+L_{12}\psi \left(z\right)\right)}^2+{\left(L_{21}\phi \left(z\right)+L_{22}\psi \left(z\right)\right)}^2\right)dz+$$

$$2{\lambda }_{1}exp(\Lambda t)\underset{0}{\overset{t}{\int }}\left(p^2\left(0\right)\frac{{\chi }_1^2\left(\tau \right)}{{({\delta }_1{R}_{11}+{\delta }_2{R}_{21})}^2}+\frac{{\chi }_2^2\left(\tau \right)}{{({\beta }_1{R}_{12}+{\beta }_2{R}_{22})}^2}\right)d\tau ,$$

(12)

where

$$\Lambda =\frac{2{\lambda }_1|p\left(l\right)-p\left(0\right)|}{lmin\left(p\right(0),p(l\left)\right)},$$

and $L_{ij}$ are the components of matrix $\mathbf{L}$.

Proof 

(Proof). We introduce the new variables $w_1$, $w_2$ as: $\mathbf{u}=\mathbf{R}\mathbf{w}$, and according to the fact that $\mathbf{L}={\mathbf{R}}^{-1}$, system (6) is rewritten as:

$$\frac{\partial w_1}{\partial t}+{\lambda }_1\frac{\partial w_1}{\partial z}=0,\frac{\partial w_2}{\partial t}+{\lambda }_2\frac{\partial w_2}{\partial z}=0,$$

(13)

where $w_1=L_{11}a+L_{12}q$, $w_2=L_{21}a+L_{22}q$.

The initial and boundary conditions are rewritten as

$$w_1(0,z)=\tilde{\phi }\left(z\right)=L_{11}\phi \left(z\right)+L_{12}\psi \left(z\right),w_2(0,z)=\tilde{\psi }\left(z\right)=L_{21}\phi \left(z\right)+L_{22}\psi \left(z\right),$$

$$w_1(t,0)+e_1{w}_2(t,0)={\tilde{\chi }}_1\left(t\right),w_2(t,l)+e_2{w}_1(t,l)={\tilde{\chi }}_2\left(t\right),$$

where

$${\tilde{\chi }}_1\left(t\right)=\frac{{\chi }_1\left(t\right)}{{\delta }_1{R}_{11}+{\delta }_2{R}_{21}},{\tilde{\chi }}_2\left(t\right)=\frac{{\chi }_2\left(t\right)}{{\beta }_1{R}_{12}+{\beta }_2{R}_{22}}.$$

The following new variables are introduced:

$${\overline{w}}_1(t,z)=p\left(z\right)w_1(t,z),{\overline{w}}_2(t,z)=w_2(t,z).$$

It must be noted that function $p\left(z\right)$ with properties (11) can be selected for any finite values of $e_1$ and $e_2$.

So, the Equation (13) are rewritten as

$$\frac{\partial {\overline{w}}_1}{\partial t}+{\lambda }_1\frac{\partial {\overline{w}}_1}{\partial z}-g\left(z\right){\overline{w}}_1=0,\frac{\partial {\overline{w}}_2}{\partial t}+{\lambda }_2\frac{\partial {\overline{w}}_2}{\partial z}=0,$$

(14)

where $g\left(z\right)=\frac{{\lambda }_1{p}^{\prime }\left(z\right)}{p\left(z\right)}$. The initial and boundary conditions are rewritten as

$${\overline{w}}_1(0,z)=\overline{\phi }\left(z\right)=p\left(z\right)\tilde{\phi }\left(z\right),{\overline{w}}_2(0,z)=\tilde{\psi }\left(z\right),$$

(15)

$${\overline{w}}_1(t,0)+{\overline{e}}_1{\overline{w}}_2(0,z)={\overline{\chi }}_1\left(t\right),{\overline{w}}_2(t,l)+{\overline{e}}_2{\overline{w}}_1(t,l)={\overline{\chi }}_2\left(t\right),$$

(16)

where

$${\overline{e}}_1=e_{1}p\left(0\right),{\overline{e}}_2=\frac{e_2}{p\left(l\right)},{\overline{\chi }}_1\left(t\right)=p\left(0\right){\tilde{\chi }}_1\left(t\right),{\overline{\chi }}_2\left(t\right)={\tilde{\chi }}_2\left(t\right).$$

We consider some functions $f\left(z\right)$, $k\left(z\right)$, $h\left(t\right)$ and let the following relation take place for $z=a$: $f\left(a\right)=-sk\left(a\right)+h\left(t\right)$, where $s\in \mathrm{R}$. It is easy to demonstrate that $f^2\left(a\right)\le 2h^2\left(t\right)+2s^2{k}^2\left(a\right)$, so $f^2\left(a\right)-k^2\left(a\right)\le 2h^2\left(t\right)-(1-2s^2)k^2\left(a\right)$. With the use of this inequality we obtain that

$${\overline{w}}_1^2(t,0)-{\overline{w}}_2^2(t,0)\le 2{\overline{\chi }}_1^2\left(t\right)-(1-2{\overline{e}}_1^2){\overline{w}}_2^2(t,0),$$

and according to the condition $|{\overline{e}}_1|\le \frac{1}{2}$ (see (11)), the following inequality can be written:

$${\overline{w}}_1^2(t,0)-{\overline{w}}_2^2(t,0)\le 2{\overline{\chi }}_1^2\left(t\right),$$

(17)

and by the same way, we obtain that

$${\overline{w}}_2^2(t,l)-{\overline{w}}_1^2(t,l)\le 2{\overline{\chi }}_2^2\left(t\right).$$

(18)

It is easy to demonstrate that

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{\overline{w}}_1^{2}dz=-{\lambda }_1{\overline{w}}_1^2{|}_0^l+2\underset{0}{\overset{l}{\int }}g{\overline{w}}_1^{2}dz,\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{\overline{w}}_2^{2}dz=-{\lambda }_1{\overline{w}}_2^2{|}_0^l.$$

So, the following inequality is written:

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}({\overline{w}}_1^2+{\overline{w}}_2^2)dz=-{\lambda }_1{\overline{w}}_1^2{|}_0^l-{\lambda }_2{\overline{w}}_2^2{|}_0^l+2\underset{0}{\overset{l}{\int }}g{\overline{w}}_1^{2}dz=-{\lambda }_1{\overline{w}}_1^2{|}_0^l+$$

$$|{\lambda }_2|{\overline{w}}_2^2{|}_0^l+2\underset{0}{\overset{l}{\int }}g{\overline{w}}_1^{2}dz\le {\lambda }_1({\overline{w}}_2^2-{\overline{w}}_1^2){|}_0^l+2\underset{0}{\overset{l}{\int }}g{\overline{w}}_1^{2}dz.$$

From (17), (18) we obtain that

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}({\overline{w}}_1^2+{\overline{w}}_2^2)dz\le 2{\lambda }_1({\overline{\chi }}_1^2\left(t\right)+{\overline{\chi }}_2^2\left(t\right))+2\underset{0}{\overset{l}{\int }}g{\overline{w}}_1^{2}dz.$$

With the use of the expression for $p\left(z\right)$ and the integral mean value theorem, the following inequality can be obtained:

$$2\underset{0}{\overset{l}{\int }}g{\overline{w}}_1^{2}dz\le \Lambda \underset{0}{\overset{l}{\int }}{\overline{w}}_1^{2}dz.$$

So, the following estimation takes place:

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}({\overline{w}}_1^2+{\overline{w}}_2^2)dz\le 2{\lambda }_1({\overline{\chi }}_1^2\left(t\right)+{\overline{\chi }}_2^2\left(t\right))+\Lambda \underset{0}{\overset{l}{\int }}({\overline{w}}_1^2+{\overline{w}}_2^2)dz.$$

(19)

After the integration of (19) in time, we obtain the following inequality:

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}({\overline{w}}_1^2+{\overline{w}}_2^2)dz\le exp(\Lambda t)\underset{0}{\overset{l}{\int }}({\overline{\phi }}^2\left(z\right)+{\overline{\psi }}^2\left(z\right))dz+$$

$$2{\lambda }_{1}exp(\Lambda t)\underset{0}{\overset{t}{\int }}({\overline{\chi }}_1^2\left(\tau \right)+{\overline{\chi }}_2^2\left(\tau \right))d\tau ,$$

and from this inequality, after returning to the variables a and q, we can obtain (12). □

So, as we can see, according to (12), the problem (6)–(8) can have solutions, which are unbounded at $t\to +\infty$, even though the bounded functions ${\chi }_1\left(t\right)$ and ${\chi }_2\left(t\right)$ are considered, and the physically relevant condition (9) may not be realized.

3.1. Case of Boundary Conditions for q

It must be noted that for the special classes of boundary conditions, the integral extimations, which are different from (12), can be obtained. LWe consider system (6) with initial conditions (7) and the following boundary conditions:

$$q(t,0)=q_1\left(t\right),q(t,l)=q_2\left(t\right).$$

(20)

Theorem 2.

Let $q_1\left(t\right)$, $q_2\left(t\right)\in {\mathrm{C}}^1[0,+\infty )$. Then, the following integral inequality takes place:

$$\underset{0}{\overset{l}{\int }}\left(a^2(t,z)+\frac{1}{d_1}{\left(q(t,z)-\frac{1}{l}(q_2\left(t\right)z-q_1\left(t\right)(z-l))\right)}^2\right)dz\le$$

$$2\underset{0}{\overset{l}{\int }}\left({\phi }^2\left(z\right)+\frac{1}{d_1}{\left(\psi \left(z\right)-\frac{1}{l}(q_2\left(0\right)z-q_1\left(0\right)(z-l))\right)}^2\right)dz+\frac{2t}{l}\underset{0}{\overset{t}{\int }}{(q_2\left(\tau \right)-q_1\left(\tau \right))}^{2}d\tau +$$

$$\frac{2t}{d_1{l}^2}\underset{0}{\overset{t}{\int }}\underset{0}{\overset{l}{\int }}{(d_2(q_2\left(\tau \right)-q_1\left(\tau \right))+{\dot{q}}_2\left(\tau \right)z-{\dot{q}}_1\left(\tau \right)(z-l))}^{2}dzd\tau .$$

(21)

Proof. 

The new function, which satisfies the boundary condition (20), is introduced:

$$\mathcal{Q}(t,z)=q(t,z)-\frac{1}{l}(q_2\left(t\right)z-q_1\left(t\right)(z-l)).$$

So, the system (6) is rewritten as

$$\frac{\partial a}{\partial t}+\frac{\partial \mathcal{Q}}{\partial z}=F\left(t\right),\frac{\partial \mathcal{Q}}{\partial t}+d_1\frac{\partial a}{\partial z}+d_2\frac{\partial \mathcal{Q}}{\partial z}=G(t,z),$$

(22)

where

$$F\left(t\right)=\frac{q_1\left(t\right)-q_2\left(t\right)}{l},G(t,z)=\frac{d_2(q_1\left(t\right)-q_2\left(t\right))-({\dot{q}}_2\left(t\right)z-{\dot{q}}_1\left(t\right)(z-l))}{l}.$$

System (22) is considered with the following initial and boundary conditions:

$$a(0,z)=\phi \left(z\right),\mathcal{Q}(0,z)=\zeta \left(z\right),\mathcal{Q}(t,0)=\mathcal{Q}(t,l)=0,$$

where $\zeta \left(z\right)=\psi \left(z\right)-\frac{q_2\left(0\right)z-q_1\left(0\right)(z-l)}{l}$. So, the inhomogeneous system (22) with the homogeneous boundary conditions for $\mathcal{Q}(t,z)$ is considered. According to the Duhamel’s principle [34], the solution of this problem can be presented as

$$a(t,z)=\tilde{a}(t,z)+\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau ,\mathcal{Q}(t,z)=\tilde{\mathcal{Q}}(t,z)+\underset{0}{\overset{t}{\int }}{\mathcal{Q}}_{\tau }(t,z;\tau )d\tau ,$$

(23)

where $(\tilde{a},\tilde{\mathcal{Q}})$ is the solution of the following problem:

$$\frac{\partial \tilde{a}}{\partial t}+\frac{\partial \tilde{\mathcal{Q}}}{\partial z}=0,\frac{\partial \tilde{\mathcal{Q}}}{\partial t}+d_1\frac{\partial \tilde{a}}{\partial z}+d_2\frac{\partial \tilde{\mathcal{Q}}}{\partial z}=0,t\in (0,+\infty ),z\in (0,l),$$

(24)

$$\tilde{a}(0,z)=\phi \left(z\right),\tilde{\mathcal{Q}}(t,z)=\zeta \left(z\right),$$

$$\tilde{\mathcal{Q}}(t,0)=\tilde{\mathcal{Q}}(t,l)=0,$$

(25)

and $(a_{\tau },{\mathcal{Q}}_{\tau })$ is the solution of the following problem:

$$\frac{\partial a_{\tau }}{\partial t}+\frac{\partial {\mathcal{Q}}_{\tau }}{\partial z}=0,\frac{\partial {\mathcal{Q}}_{\tau }}{\partial t}+d_1\frac{\partial a_{\tau }}{\partial z}+d_2\frac{\partial {\mathcal{Q}}_{\tau }}{\partial z}=0,t\in (\tau ,+\infty ),z\in (0,l),$$

$$a_{\tau }(\tau ,z;\tau )=F\left(\tau \right),{\mathcal{Q}}_{\tau }(\tau ,z;\tau )=G(\tau ,z),{\mathcal{Q}}_{\tau }(t,0;\tau )={\mathcal{Q}}_{\tau }(t,l;\tau )=0.$$

We obtain the expression for the time derivative of $\underset{0}{\overset{l}{\int }}\left({\tilde{a}}^2(t,z)+{\tilde{\mathcal{Q}}}^2(t,z)\right)dz$:

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{\tilde{a}}^2(t,z)dz=-2\underset{0}{\overset{l}{\int }}\tilde{a}(t,z)\frac{\partial \tilde{\mathcal{Q}}}{\partial z}(t,z)dz=-2(\tilde{a}(t,z)\tilde{\mathcal{Q}}(t,z){|}_0^l-$$

$$-\underset{0}{\overset{l}{\int }}\tilde{\mathcal{Q}}(t,z)\frac{\partial \tilde{a}}{\partial z}(t,z)dz)=-\frac{2}{d_1}\underset{0}{\overset{l}{\int }}\mathcal{Q}(t,z)\frac{\partial \mathcal{Q}}{\partial t}(t,z)dz-\frac{2d_2}{d_1}\underset{0}{\overset{l}{\int }}\tilde{\mathcal{Q}}(t,z)\frac{\partial \tilde{\mathcal{Q}}}{\partial z}(t,z)dz=$$

$$=-\frac{1}{d_1}\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{\tilde{\mathcal{Q}}}^2(t,z)dz,$$

where the boundary conditions (25) are taken into account.

So, we can obtain the following:

$$\underset{0}{\overset{l}{\int }}\left({\tilde{a}}^2(t,z)+{\tilde{\mathcal{Q}}}^2(t,z)\right)dz=\left(1-\frac{1}{d_1}\right)\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{\tilde{\mathcal{Q}}}^2(t,z)dz.$$

After the integration, the following equality is obtained:

$$\underset{0}{\overset{l}{\int }}\left({\tilde{a}}^2(t,z)+\frac{1}{d_1}{\tilde{\mathcal{Q}}}^2(t,z)\right)dz=\underset{0}{\overset{l}{\int }}\left({\phi }^2\left(z\right)+\frac{1}{d_1}{\zeta }^2\left(z\right)\right)dz.$$

(26)

By the same operations we obtain

$$\underset{0}{\overset{l}{\int }}\left(a_{\tau }^2(t,z;\tau )+\frac{1}{d_1}{\mathcal{Q}}_{\tau }^2(t,z;\tau )\right)dz=F^2\left(\tau \right)l+\frac{1}{d_1}\underset{0}{\overset{l}{\int }}{G}^2(\tau ,z)dz.$$

(27)

According to Equation (23):

$$\underset{0}{\overset{l}{\int }}{a}^2(t,z)dz=\underset{0}{\overset{l}{\int }}\left({\tilde{a}}^2(t,z)+2\tilde{a}(t,z)\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau +{\left(\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau \right)}^2\right)dz.$$

From $\underset{0}{\overset{l}{\int }}{\left(\tilde{a}(t,z)-\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau \right)}^{2}dz\ge 0$ we can obtain

$$2\underset{0}{\overset{l}{\int }}\tilde{a}(t,z)\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau dz\le \underset{0}{\overset{l}{\int }}\left({\tilde{a}}^2(t,z)+{\left(\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau \right)}^2\right)dz,$$

and the following inequality can be written:

$$\underset{0}{\overset{l}{\int }}{a}^2(t,z)dz\le 2\underset{0}{\overset{l}{\int }}\left({\tilde{a}}^2(t,z)+{\left(\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau \right)}^2\right)dz.$$

(28)

The same inequality takes place for $\mathcal{Q}(t,z)$.

So, the following estimation is obtained:

$$\underset{0}{\overset{l}{\int }}\left(a^2(t,z)+\frac{1}{d_1}{\mathcal{Q}}^2(t,z)\right)dz\le 2\underset{0}{\overset{l}{\int }}\left({\tilde{a}}^2(t,z)+\frac{1}{d_1}{\tilde{\mathcal{Q}}}^2(t,z)\right)dz+$$

$$+2\underset{0}{\overset{l}{\int }}\left({\left(\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau \right)}^2+\frac{1}{d_1}{\left(\underset{0}{\overset{t}{\int }}{\tilde{\mathcal{Q}}}_{\tau }(t,z;\tau )d\tau \right)}^2\right)dz.$$

(29)

According to the Cauchy–Schwarz inequality, we obtain

$${\left(\underset{0}{\overset{t}{\int }}{a}_{\tau }(t,z;\tau )d\tau \right)}^2\le t\underset{0}{\overset{t}{\int }}{a}_{\tau }^2(t,z;\tau )d\tau ,{\left(\underset{0}{\overset{t}{\int }}{\tilde{\mathcal{Q}}}_{\tau }(t,z;\tau )d\tau \right)}^2\le t\underset{0}{\overset{t}{\int }}{\tilde{\mathcal{Q}}}_{\tau }^2(t,z;\tau )d\tau .$$

So, the following integral inequality can be obtained from (26)–(29):

$$\underset{0}{\overset{l}{\int }}\left(a^2(t,z)+\frac{1}{d_1}{\mathcal{Q}}^2(t,z)\right)dz\le 2\underset{0}{\overset{l}{\int }}\left({\phi }^2\left(z\right)+\frac{1}{d_1}{\zeta }^2\left(z\right)\right)+2tl\underset{0}{\overset{t}{\int }}{F}^2\left(\tau \right)d\tau +$$

$$+\frac{2t}{d_1}\underset{0}{\overset{l}{\int }}\underset{0}{\overset{t}{\int }}{G}^2(\tau ,z)d\tau dz,$$

and after the return to the function $q(t,z)$ we can see, that inequality (21) is correct. □

As we can see, the last terms in (21) provide the possibility of the time-growing solutions for the bounded physically relevant functions $q_1\left(t\right)$ and $q_2\left(t\right)$. We construct the simple example, which illustrates this behavior. Let $q_1$ and $q_2$ be the constants $q_1=C_1$, $q_2=C_2$, $C_1>C_2$. Let $\psi \left(z\right)=Dz+C_1$, where D = $\frac{C_2-C_1}{l}<0$. Let $\phi \left(z\right)=-\frac{d_2}{d_1}Dz$. This problem has the following solution:

$$q(t,z)=Dz+C_1,a(t,z)=-Dt+\phi \left(z\right).$$

So we can see that constant boundary conditions and bounded initial conditions lead to the unbounded unphysical ($a(t,z)\to +\infty$) solution.

In order to demonstrate the similar behavior in the fully nonlinear model (1), which is widely used in practice, the physiological example on the simulation of blood flow in the part of the human aorta can be considered. For the simulation, the segment Ao. V from aortic model, presented in [9], is chosen. For this case, the variable mechanical properties take place: $A_d=A_d\left(z\right)$, $\beta =\beta \left(z\right)$. For this vessel, $l=$ 15.2 cm, the radius of the inlet cross-section $R_{in}$ is equal to 1.25 cm, and the radius of the outlet cross-section $R_{out}$ is equal to 0.99 cm. The $A_d\left(z\right)$ is chosen as

$$A_d\left(z\right)=\pi {\left(R_{in}+\frac{R_{out}-R_{in}}{l}z\right)}^2.$$

The values of $\beta$ for inlet and outlet cross-sections are computed from the values of the speed of sound, presented in Table IV from [9], and $\beta \left(z\right)$ is considered as a linear function: $\beta \left(z\right)={\beta }_{in}+\frac{{\beta }_{out}-{\beta }_{in}}{l}z$.

The following boundary conditions are stated: $Q(t,0)=q_1=500$${\mathrm{cm}}^3/\mathrm{s}$, $Q(t,l)=q_2=300$${\mathrm{cm}}^3/\mathrm{s}$, which are close to the typical values of Q from [9]. The initial conditions are stated as

$$A(0,z)=A_d\left(z\right),Q(0,z)=\frac{q_2-q_1}{l}z+q_1.$$

The numerical simulations are realized on the time interval from 0 to 1 s. The second-order Lax–Wendroff scheme is used for the computations on the spatial grid with 100 nodes and the time grid with 10,000 nodes. The values of A in the boundary points are computed from the compatibility conditions (for details, see [26]).

In Figure 1, the results of the numerical simulations are presented. As can be seen from the plots of A in the left, middle, and right vessel points, for this type of bounded boundary and the initial conditions, the obtained solutions do not have any physical or physiological meaning because they lead to the unbounded values of A. So, for this practical example we have the same behavior as for the linear theory.

3.2. Case of Periodic Boundary Conditions

We consider the following periodic boundary conditions for the system (6):

$$a(t,l)=a(t,0),q(t,l)=q(t,0),$$

(30)

with the same initial conditions, but it is proposed that $\phi \left(l\right)=\phi \left(0\right)$ and $\psi \left(l\right)=\psi \left(0\right)$.

3.2.1. Integral Estimation

Theorem 3.

The following equation takes place for the solution of problem (6), (7), (30):

$$\underset{0}{\overset{l}{\int }}\left(a^2(t,z)+\frac{1}{d_1}{q}^2(t,z)\right)dz=\underset{0}{\overset{l}{\int }}\left({\phi }^2\left(z\right)+\frac{1}{d_1}{\psi }^2\left(z\right)\right)dz,$$

(31)

Proof. 

We obtain the expressions for $\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}\left(a^2(t,z)+q^2(t,z)\right)dz$:

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{a}^2(t,z)dz=-2\underset{0}{\overset{l}{\int }}\tilde{a}(t,z)\frac{\partial \tilde{q}}{\partial z}(t,z)dz=-2(a(t,z)q(t,z){|}_0^l-$$

$$-\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial a}{\partial z}(t,z)dz)=2\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial a}{\partial z}(t,z)dz.$$

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{q}^2(t,z)dz=2\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial q}{\partial z}(t,z)dz=-2\underset{0}{\overset{l}{\int }}q(t,z)\left(d_1\frac{\partial a}{\partial z}(t,z)+d_2\frac{\partial q}{\partial z}(t,z)\right)dz.$$

$$\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial q}{\partial z}(t,z)dz=q^2(t,z){|}_0^l-\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial q}{\partial z}(t,z)dz,$$

so $\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial q}{\partial z}(t,z)dz=0$,

$$d_1\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial a}{\partial z}(t,z)dz=-\frac{1}{2}\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{q}^2(t,z)dz.$$

So, we obtain that

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}\left(a^2(t,z)+q^2(t,z)\right)dz=2\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial a}{\partial z}(t,z)dz-2d_1\underset{0}{\overset{l}{\int }}q(t,z)\frac{\partial a}{\partial z}(t,z)dz=$$

$$=-\frac{1}{d_1}\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{q}^2(t,z)dz+\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}{q}^2(t,z)dz.$$

So, the following equation is obtained:

$$\frac{\partial }{\partial t}\underset{0}{\overset{l}{\int }}\left(a^2(t,z)+\frac{1}{d_1}{q}^2(t,z)\right)dz=0,$$

and after its integration on t we obtain (31). □

So, for the periodic conditions (30), the condition (9) is realized for any appropriate functions $\phi \left(z\right)$ and $\psi \left(z\right)$.

3.2.2. Fourier Method

As is mentioned in the Introduction, the approach to obtain the solution of the initial-boundary-value problem for system (6) in the general form is proposed in [28,29], and it is based on the reflection of waves. In the case of periodic conditions (30), the Fourier method can be easily and effectively used.

Let the solution of system (6) be presented as

$$a(t,z)=\alpha \left(t\right)Z\left(z\right),q(t,z)=\kappa \left(t\right)Z\left(z\right).$$

(32)

After the substitution of (32) in the first equation of (6), the following relation is obtained:

$$\frac{\dot{\alpha }\left(t\right)}{\kappa \left(t\right)}=-\frac{Z^{\prime }\left(z\right)}{Z\left(z\right)}=-\lambda =const.$$

So, we obtain the following eigenvalue problem:

$$Z^{\prime }\left(z\right)-\lambda Z\left(z\right)=0,Z\left(0\right)=Z\left(l\right).$$

(33)

Problem (33) has the following eigenfunctions and corresponding eigenvalues:

$$Z_m\left(z\right)=e^{\frac{2\pi m\mathrm{i}}{l}z},{\lambda }_m=\frac{2\pi m\mathrm{i}}{l},m\in \mathrm{Z}.$$

The solution to the problem with periodic conditions (30) is presented as Fourier expansions on these eigenfunctions:

$$a(t,z)=\sum _{m=-\infty }^{+\infty }{\alpha }_m\left(t\right)e^{\frac{2\pi m\mathrm{i}}{l}z},q(t,z)=\sum _{m=-\infty }^{+\infty }{\kappa }_m\left(t\right)e^{\frac{2\pi m\mathrm{i}}{l}z}.$$

(34)

After the substitution of (34) into (6), the following equations are obtained for ${\alpha }_m\left(t\right)$ and ${\kappa }_m\left(t\right)$:

$${\ddot{\alpha }}_m+d_2{\lambda }_m\dot{\alpha }-d_1{\lambda }_m^2{\alpha }_m=0,{\kappa }_m=-\frac{{\dot{\alpha }}_m}{{\lambda }_m},$$

(35)

for $m\ne 0$.

The solutions of (35) are presented as

$${\alpha }_m\left(t\right)=C_{1m}{e}^{\mathrm{i}\frac{\pi m{r}_1}{l}t}+C_{2m}{e}^{\mathrm{i}\frac{\pi m{r}_2}{l}t},$$

where $r_{1,2}=2{\lambda }_{1,2}$, and $C_{1m}$ and $C_{2m}$ are the constants. For $m=0$ we obtain that ${\alpha }_0\left(t\right)={\sigma }_0$, ${\kappa }_0\left(t\right)={\sigma }_1$, where ${\sigma }_0$ and ${\sigma }_1$ are the constants.

The constants in the expressions for ${\alpha }_m\left(t\right)$ and ${\kappa }_m\left(t\right)$ are defined after the substitution of (34) into initial conditions (7), expanded in the Fourier series on eigenfunctions of (33):

$$\phi \left(z\right)=\sum _{m=-\infty }^{+\infty }{\phi }_m{e}^{\frac{2\pi m\mathrm{i}}{l}z},\psi \left(z\right)=\sum _{m=-\infty }^{+\infty }{\psi }_m{e}^{\frac{2\pi m\mathrm{i}}{l}z},$$

where

$${\phi }_m=\frac{1}{l}\underset{0}{\overset{l}{\int }}\phi \left(z\right)e^{-\frac{2\pi m\mathrm{i}}{l}z}dz,{\psi }_m=\frac{1}{l}\underset{0}{\overset{l}{\int }}\psi \left(z\right)e^{-\frac{2\pi m\mathrm{i}}{l}z}dz.$$

After this operation, the following expressions for constants are obtained:

$$C_{1m}=-\frac{2{\psi }_m+{\phi }_m{r}_2}{r_1-r_2},C_{2m}=\frac{r_1{\phi }_m+2{\psi }_m}{r_1-r_2},{\sigma }_0={\phi }_0,{\sigma }_1={\psi }_0.$$

We consider the following example: let $l=1$, and the initial functions are presented as $\phi \left(z\right)=A_0+\frac{sin\left(2\pi z\right)}{2}$, $\psi \left(z\right)=0$, where $A_0\ge 1$. So, we obtain the following values of constants:

$${\sigma }_0=A_0,{\sigma }_1=0,C_{1m}=C_{2m}=0,m\ne \pm 1,$$

$$C_{11}=\frac{\mathrm{i}{r}_2}{2(r_1-r_2)},C_{21}=-\frac{\mathrm{i}{r}_1}{2(r_1-r_2)},C_{1-1}=-C_{11},C_{2-1}=-C_{21}.$$

For this example, $a(t,z)$ and $q(t,z)$ are the real functions, written as

$$a(t,z)=A_0+\frac{1}{2(r_2-r_1)}(r_{2}sin\left(\pi (2z+r_{1}t)\right)-r_{1}sin\left(\pi (2z+r_{2}t)\right)),$$

$$q(t,z)=\frac{r_1{r}_2}{2(r_2-r_1)}(sin\left(\pi (2z+r_{1}t)\right)-sin\left(\pi (2z+r_{2}t)\right)).$$

4. Discussion and Conclusions

In the presented paper, the author tries to expand the field of application of analytical methods in the simulation of blood flow and to analyze the effects of boundary conditions on the solutions of the linearized system of equations of 1D hemodynamics. The main results can be formulated as follows:

• The integral estimation (12) for the solution in the case of the general form of the boundary conditions is obtained. As can be seen, the unphysical unbounded solutions can take place for the case of bounded functions ${\chi }_1\left(t\right)$, ${\chi }_2\left(t\right)$.
• For the case of the boundary conditions on the flow rate, the estimation (21) is obtained, which is more accurate than (12), and it illustrates the similar possible behavior of the solutions.
• The general theory of the integral inequalities for hyperbolic equations is presented in [35,36]. However, in the presented paper, the more exact estimations are presented for the specific example of the hyperbolic system, where the specific constants in exponentials are obtained. Moreover, for the boundary conditions only for q, the estimation (21), which is more accurate than the exponential function (the linear function of t), is obtained.
  As it is mentioned in [36], the energetic norm is presented by the square of the $L_2$-norm of the solution. In the left parts of (12), (21), and (32), the integrals, which can be considered as squares of norms of functions, linearly related to the solutions, are presented. So, our estimations can be considered as some kind of energy inequalities for this specific type of hyperbolic system with the appropriate boundary conditions.
• For the periodic boundary conditions, the exact integral estimation (31) is obtained, which illustrates the correct behavior of the solution—it is bounded at $t\to +\infty$. For this case of boundary conditions, the Fourier method for the analytical solution can be applied. Such analytical solutions can be used for the comparison of different 1D blood-flow models [27].
• The results obtained in the paper can be useful for the specialists on blood-flow modeling because they allow for an alternative view of the stated boundary conditions and can explain some of the problems that can arise in numerical simulations.
• In the numerical experiment, where the fully nonlinear model, used in many works, is considered, it is demonstrated that the situation described by the author’s theoretical results can be observed in practice, when the bounded initial and boundary conditions lead to the incorrect results from the physical point of view. So, it is important to correctly impose the boundary conditions for the practical predictive simulations. From the medical point of view, this means that the users of the software must choose such conditions carefully because in the opposite case, it can lead to the incorrect results.