Semi-analytical method for blood flow model with source terms

Abstract

This article discusses a semi-analytical method for solving modelling equations for blood flow through compliant vessels. He's method, also called the variational iteration method (VIM), is used to find an approach to a solution in an explicit form. The chosen problem has been handled by an asymptotic reduction of the associated incompressible Navier–Stokes equations, to obtain a one-dimensional nonlocal blood flow model with source terms. Numerical examples were chosen based on experimentally valued parameters and were used to highlight and illustrate the benefits of VIM. The results obtained from the blood flow model were compared to a numerical reference method.

Keywords:

• VIM
• 1D blood model
• compliant vessels

2010 MSCs:

• 35C99
• 65D99
• 76Z05

1. Introduction

Differential equations play a crucial role in modelling natural phenomena in reality. Scientists deal with many properties and actions (collisions, chemical bonding, solitary wave, viscosity, temperature, pressure, etc.) making the equations that govern these problems increasingly difficult to solve it analytically [1–6]. For this purpose, researchers have resorted to finding solutions using many methods to approach it, as the direct algebraic methods [7–9], rational expansion method [10], reductive perturbation technique [11], extended mapping [12], Adomian decomposition method, homotopy perturbation method, differential transform method, and VIM [13–18]. VIM, for example, is a semi-analytical technique among these mentioned methods. It was invented by He at the end of the last century [19,20]. Many authors have developed VIM in the last decades, as Wazwaz [21,22], which has proposed to solve various types of nonlinear problems. The VIM solution is treated by an iterative process leading to an accurate solution. The advantage of this method is that it can provide an approximate solution in an explicit form, unlike classical numerical methods which can only provide it by interpolation. In this study, we treat the problem of blood flow behaviour through compliant vessels using the VIM. It is therefore necessary to reduce the number of variables involved in this problem. This approach is based on an asymptotic reduction process applied to the associated incompressible Navier–Stokes equations. It should be noted that recently this subject has been treated using effective numerical methods by certain researchers, one can quote for example [23,24]. However, in this work, we present a way to use the advantages of VIM to find a solution approach in an explicit form without using a linearization or discretization procedure. We suggest applying this method to study the blood flow in the artery because this study allows us to know the blood flow by measuring the blood pressure and determining its flow in the arteries. Blood flow is assumed to be Newtonian and is governed by the equation of continuity and the equilibrium equation of momentum (which are known as the Navier–Stokes equations). We apply VIM to the equations governing blood flow in arteries in the presence of a source term. We proposed a comparison between the VIM solution with a numerical method based on the MATLAB PDEPE function to evaluate this proposed method.

2. Global one-dimensional blood flow modelling

In recent years, a lot of works have proposed mathematical and numerical modelling of the human cardiovascular system. These models are based on Navier–Stokes equations [26,27]. The 3D model describes the flow motion in compliant vessels and its interaction with wall displacement [28,29]. To simplify these equations, many authors have studied a monodirectional flow system [25,29,31], where the 1D flow equations are generated from a single vessel system. This model was used to study systematic circulation since their activity is based on reasonable calculation costs. This model describes an incompressible viscous fluid flow in an elastic tube taking account the blood vessels and their interaction with wall displacement. Several approaches can be used to derive the equations by using the asymptotic analysis of incompressible Navier–Stokes equations in long and narrow channels [29,30]. Two approaches are distinguished. The first integrates the Navier–Stokes equations in a generic section. The second is based on adimensional underlying equations [30,31]. This approach seems more realist and gives global results.

2.1. Asymptotic reduction

In the asymptotic reduction, the equations will be treated with nondimensional variables, for more details see [30]. Basing on that the artery or aorta is an almost cylindrical vessel and the blood flows in the axial direction. Then, the system of equations used must be in cylindrical coordinates $(x,r,\theta )$, when the x-axis is aligned with the symmetry axis of the vessel. We denote by $V=(V_x,V_r,U_{\theta })$ the velocity. We will first assume that the angular velocity is zero to obtain the following equations of motion: (1) $\begin{array}{rl}& ({\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}+V_r{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }r}}+V_x{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}-\nu {\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{r}^2}}-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }r}}-{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}})V_x\\ & +{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }x}}=0,\end{array}$(1) (2) $\begin{array}{rl}& ({\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}+V_r{\displaystyle \frac{\mathrm{\partial }{V}_r}{\mathrm{\partial }r}}+V_x{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}-\nu {\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{r}^2}}-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }r}}-{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}})V_r\\ & +{\displaystyle \frac{V_r}{r^2}}+{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }r}}=0.\end{array}$(2) Translating the reduced nondimensional equations and using the incompressibility condition, we obtain (3) ${\displaystyle \frac{\mathrm{\partial }{V}_x}{\mathrm{\partial }x}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{r\mathrm{\partial }{V}_r}{\mathrm{\partial }r}}=0.$(3) We consider the characteristic quantities $U_0$ and $V_0$ that are the radial and axial characteristic velocity, L is the characteristic length, $R_0$ is the ship's inner radius, and the following appropriate dimensionless variables: $\begin{array}{rl}& r=R_0\tilde{r},x=L\tilde{x},t={\displaystyle \frac{L}{V_0}}\tilde{t},V_x=V_0\widetilde{V_x},\\ & V_r=V_0\widetilde{V_r},p=\rho V_0^2\tilde{p}.\end{array}$By a simple calculation [30], the first reduced motion equation will be as follows: (4) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{t}}}(\tilde{r}\widetilde{V_x})+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{r}}}(\tilde{r}\widetilde{V_r}\widetilde{V_x})+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}(\tilde{r}{\widetilde{V_x}}^2)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{x}}}(\tilde{r}\tilde{p})\\ & ={\displaystyle \frac{\mathrm{\nu \lambda }}{V_0{R}_0^2}}({\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{r}}}(\tilde{r}{\displaystyle \frac{\mathrm{\partial }\widetilde{V_x}}{\mathrm{\partial }\tilde{r}}})).\end{array}$(4) The second reduced motion equation will be as follows: (5) ${\displaystyle \frac{\mathrm{\partial }\tilde{p}}{\mathrm{\partial }\tilde{r}}}=0.$(5) This equation implies that the pressure is constant across the cross-section of the vessel. Finally, the reduced incompressibility condition is as follows: (6) ${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{r}}}(\tilde{r}\widetilde{V_r})+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{x}}}(\tilde{r}\widetilde{V_x})=0.$(6) In this step, we will express the equations in terms of the average quantity in the transverse zone. Let $\tilde{R}$ be the inner ship radius.

Therefore, (7) $\begin{array}{rl}& \tilde{U}={\displaystyle \frac{1}{{\tilde{R}}^2}}{\displaystyle {\int }_0^{\tilde{R}}2\widetilde{V_x}\tilde{r}\mathrm{d}\tilde{r},}\end{array}$(7) (8) $\begin{array}{rl}& \alpha ={\displaystyle \frac{1}{{\tilde{R}}^2{\tilde{U}}^2}}{\displaystyle {\int }_0^{\tilde{R}}2{\widetilde{V_x}}^2\tilde{r}\mathrm{d}\tilde{r},}\end{array}$(8) where $\tilde{U}$ is the axial mean velocity, and α is the correction term.

We integrate the equations on $\tilde{r}=0$ to $\tilde{R}$ which expresses the terms of the average quantity by specifying the boundary condition with $\tilde{r}=\tilde{R}$: $[{\tilde{V}}_r{]}_{\tilde{r}=\tilde{R}}={\displaystyle \frac{\mathrm{\partial }\tilde{R}}{\mathrm{\partial }\tilde{x}}}[{\tilde{V}}_x{]}_{\tilde{r}=\tilde{R}}+{\displaystyle \frac{\mathrm{\partial }\tilde{R}}{\mathrm{\partial }\tilde{t}}}.$For the incompressibility condition, we integrate Equations (6(6) ${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{r}}}(\tilde{r}\widetilde{V_r})+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{x}}}(\tilde{r}\widetilde{V_x})=0.$(6) ) and using the definition of $\tilde{U}$, then ${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}({\tilde{R}}^2)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{x}}}\left(\tilde{V}{\displaystyle \frac{{\tilde{R}}^2}{2}}\right)=0.$For the first equation of motion, we integrate the first reduced motion equation, we have (9) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{t}}}\left({\displaystyle \frac{{\tilde{R}}^2}{2}}\tilde{V}\right)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{x}}}\left({\displaystyle \frac{\alpha {\tilde{R}}^2{\tilde{V}}^2}{2}}\right)+{\displaystyle \frac{{\tilde{R}}^2}{2}}{\displaystyle \frac{\mathrm{\partial }\tilde{p}}{\mathrm{\partial }\tilde{x}}}\\ & ={\displaystyle \frac{\mathrm{\nu \lambda }}{U_0{R}_0^2}}\tilde{R}{\left({\displaystyle \frac{\mathrm{\partial }\widetilde{U_x}}{\mathrm{\partial }\tilde{r}}})\right)}_{\tilde{r}=\tilde{R}}.\end{array}$(9) The dimensionless reduced and average equations are as follows: (10) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}({\tilde{R}}^2)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{x}}}(\tilde{U}{\tilde{R}}^2)=0,\end{array}$(10) (11) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{t}}}({\tilde{R}}^2\tilde{U})+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\tilde{x}}}(\alpha {\tilde{R}}^2{\tilde{U}}^2)+{\tilde{R}}^2{\displaystyle \frac{\mathrm{\partial }\tilde{p}}{\mathrm{\partial }\tilde{x}}}\\ & =2{\displaystyle \frac{\mathrm{\nu \lambda }}{V_0{R}_0^2}}\tilde{R}{\left({\displaystyle \frac{\mathrm{\partial }\widetilde{V_x}}{\mathrm{\partial }\tilde{r}}})\right)}_{\tilde{r}=\tilde{R}}.\end{array}$(11) Finally, we define the average cross-sectional velocity U and the coefficient α, then we find their relation with dimensionless quantities (12) $\tilde{U}={\displaystyle \frac{2}{V_0{R}^2}}{\displaystyle {\int }_0^{R}r{V}_x\mathrm{d}r,}$(12) where R the internal ship radius in dimensional terms.

The dimensional axial velocity is (13) $U={\displaystyle \frac{1}{R^2}}{\displaystyle {\int }_0^{R}2r{V}_x\mathrm{d}r,}$(13) with $U=V_0\tilde{U},$ even for (14) $\alpha ={\displaystyle \frac{1}{R^2{U}^2}}{\displaystyle {\int }_0^{R}2V_x^{2}r\mathrm{d}r.}$(14) Using these definitions, the reduced equations are transformed into their dimensional form to obtain (15) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}(R^2)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}(U{R}^2)=0,\end{array}$(15) (16) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}(R^{2}U)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}(\alpha R^2{U}^2)+R^2{\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }x}}\\ & =2{\displaystyle \frac{\mathrm{\nu \lambda }}{V_0{R}_0^2}}R{\left({\displaystyle \frac{\mathrm{\partial }{V}_x}{\mathrm{\partial }r}})\right)}_{r=R}.\end{array}$(16) To obtain the equations written according to the average magnitude, one can take a typical approximation of the speed profile [30] (17) $V_x={\displaystyle \frac{\gamma +2}{\gamma }}U{(1-{\displaystyle \frac{r}{R}})}^{\gamma }.$(17) Equations (14(14) $\alpha ={\displaystyle \frac{1}{R^2{U}^2}}{\displaystyle {\int }_0^{R}2V_x^{2}r\mathrm{d}r.}$(14) ) and (17(17) $V_x={\displaystyle \frac{\gamma +2}{\gamma }}U{(1-{\displaystyle \frac{r}{R}})}^{\gamma }.$(17) ) involve the following relationship between the speed profile form determined by γ and the correction coefficient α: $\gamma ={\displaystyle \frac{2-\alpha }{\alpha -1}}.$By posing $X=R^2$ and $Y=R^{2}U$ denote the volumetric flow, we obtain (18) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }x}}=f_X,\\ {\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\left({\displaystyle \frac{Y^2}{X}}\right)+{\displaystyle \frac{X}{\rho }}{\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }x}}+K_R{\displaystyle \frac{Y}{X}}=f_Y,\end{array}$(18) where $f_X$ is the mass source per unit length, and $f_Y$ represents the acceleration of the flow due to several forces which are applied on the fluid (gravity, friction, etc.) [29], and $K_R$ is the viscous flow resistance per unit length of the vessel.

3. VIM solution of blood flow problem

The variational iteration method [21,22], which is a modified and generalized method of the Lagrange multipliers [19], demonstrates an effective, easy, and accurate resolution of a large class of nonlinear problems with approximations that converge quickly to a precise solution. To achieve our goal in the management of systems of nonlinear partial differential equations, we use VIM, and to illustrate this method we will consider the following system which writes in the form of an operator: (19) $\mathcal{L}X(t)+\mathcal{N}X(t)=h(t),\mathcal{L}Y(t)+\mathcal{N}Y(t)=g(t),$(19) where $\mathcal{L}$ and $\mathcal{N}$ are linear and nonlinear operators, respectively, h and g design the source terms. We can build the following functional correction: (20) $\begin{array}{rl}& X_{n+1}=X_n+{\displaystyle {\int }_0^t{\lambda }_1(\xi )(\mathcal{L}{X}_n(\xi )+\mathcal{N}{\tilde{X}}_n(\xi )-h(\xi ))\mathrm{d}\xi ,}\end{array}$(20) (21) $\begin{array}{rl}& Y_{n+1}=Y_n+{\displaystyle {\int }_0^t{\lambda }_2(\xi )(\mathcal{L}{Y}_n(\xi )+\mathcal{N}{\tilde{Y}}_n(\xi )-g(\xi ))\mathrm{d}\xi .}\end{array}$(21) The Lagrange multipliers ${\lambda }_1,{\lambda }_2$ can be recognized using the variational theory, $X_0,Y_0$ are appropriate initial functions, the index n indicates the nth-iteration , and $\tilde{X}(n),\tilde{Y}(n)$ are considered as restricted variations, i.e. $\delta {\tilde{X}}_n=0\text{and};\delta {\tilde{Y}}_n=0$ [14,22]. The successive approximations $X_n\text{and}{Y}_n$ solutions $X\text{and}Y$ can be easily obtained.

For the system Equation (18(18) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }x}}=f_X,\\ {\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\left({\displaystyle \frac{Y^2}{X}}\right)+{\displaystyle \frac{X}{\rho }}{\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }x}}+K_R{\displaystyle \frac{Y}{X}}=f_Y,\end{array}$(18) ), we have ${\lambda }_1={\lambda }_2=-1$. Then, the VIM method to solve the blood flow problem can be written as follows: (22) $\{\begin{array}{rl}{X}_0& =X(0,x);Y_0=Y(0,y),\\ X_{n+1}& =X_n-{\displaystyle {\int }_0^t({\displaystyle \frac{\mathrm{\partial }{X}_n}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }{Y}_n}{\mathrm{\partial }x}}-f_{X_n})\mathrm{d}\xi ,}\\ Y_{n+1}& =Y_n-{\displaystyle {\int }_0^t({\displaystyle \frac{\mathrm{\partial }{Y}_n}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\left({\displaystyle \frac{(Y_n{)}^2}{X_n}}\right)}\\ & +{\displaystyle \frac{X_n}{\rho }}\left({\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }x}}\right)+K_R{\displaystyle \frac{Y_n}{X_n}}-f_{Y_n})\mathrm{d}\xi .\end{array}$(22)

4. Numerical results

We will present in this section some examples to show the effectiveness of the described method.

4.1. Example 1

The first example is chosen with linear source terms to illustrate the VIM for a 1D-blood model. To verify the accuracy of the proposed solution, we relate the numerical results to a numerical solution using MATLAB's PDEPE function. Let the following parameters: $\begin{array}{rl}\rho & =1050{\mathrm{kg}/\mathrm{m}}^3,\\ \nu & =3.2\times {10}^{-6}{\mathrm{m}}^2/\mathrm{s}\mathrm{and}{K}_R=8\mathrm{\pi \nu }.\end{array}$We consider the initial condition (23) $X(x,0)=x+1;Y(x,0)=0.$(23) Let $f_X=X$, $f_Y=Y$ and $P(x)=x$, then the problem of the blood flow model (18(18) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }x}}=f_X,\\ {\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\left({\displaystyle \frac{Y^2}{X}}\right)+{\displaystyle \frac{X}{\rho }}{\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }x}}+K_R{\displaystyle \frac{Y}{X}}=f_Y,\end{array}$(18) ) is as follows: (24) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }x}}=X,\\ {\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\left({\displaystyle \frac{Y^2}{X}}\right)+{\displaystyle \frac{X}{\rho }}+K_R{\displaystyle \frac{Y}{X}}=Y.\end{array}$(24) Using Equations (22(22) $\{\begin{array}{rl}{X}_0& =X(0,x);Y_0=Y(0,y),\\ X_{n+1}& =X_n-{\displaystyle {\int }_0^t({\displaystyle \frac{\mathrm{\partial }{X}_n}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }{Y}_n}{\mathrm{\partial }x}}-f_{X_n})\mathrm{d}\xi ,}\\ Y_{n+1}& =Y_n-{\displaystyle {\int }_0^t({\displaystyle \frac{\mathrm{\partial }{Y}_n}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\left({\displaystyle \frac{(Y_n{)}^2}{X_n}}\right)}\\ & +{\displaystyle \frac{X_n}{\rho }}\left({\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }x}}\right)+K_R{\displaystyle \frac{Y_n}{X_n}}-f_{Y_n})\mathrm{d}\xi .\end{array}$(22) ), we get the following recursive process: (25) $\{\begin{array}{rl}{X}_0& =x+1;Y_0=0,\\ X_{n+1}& =X_n-{\displaystyle {\int }_0^t({\displaystyle \frac{\mathrm{\partial }{X}_n}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }{Y}_n}{\mathrm{\partial }x}}-X_n)\mathrm{d}\xi ,}\\ Y_{n+1}& =Y_n-{\displaystyle {\int }_0^t({\displaystyle \frac{\mathrm{\partial }{Y}_n}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\left({\displaystyle \frac{(Y_n{)}^2}{X_n}}\right)}\\ & +{\displaystyle \frac{X_n}{\rho }}+K_R{\displaystyle \frac{Y_n}{X_n}}-Y_n\phantom{\left({\displaystyle \frac{(Y_n{)}^2}{X_n}}\right)})\mathrm{d}\xi .\end{array}$(25) The first approach of VIM is written by $X_1=(t+1)(x+1);Y_1=-{\displaystyle \frac{t(x+1)}{\rho }}.$The second $\begin{array}{rl}{X}_2& =t^2{\displaystyle \frac{((x+1)\rho +1)}{(2\rho )}}+(t+1)(x+1),\\ Y_2& =-{\displaystyle \frac{1}{{\rho }^2}}(\log (t+1)-t\\ & +\rho (t^2(x+1)-\mathrm{KR}(t+\log (t+1))+{\displaystyle \frac{t^2}{2}})\phantom{(t^2(x+1)-\mathrm{KR}(t+\log (t+1))+{\displaystyle \frac{t^2}{2}})})\\ & -{\displaystyle \frac{(t(x+1))}{\rho }}.\end{array}$The approximation of averaged over the cross-section linear velocity is given by $U=\frac{Y_2}{X_2}$. It is a two-step VIM solution for the blood flow model.

Figure 1 shows the VIM solution (solid line) with the numerical solution (solid line with circle markers). Both solutions are almost identical. We note that the maximum relative error is on the order of ${10}^{-2}$ (Table 1).

Figure 1. The numerical solution obtained by PDEPE is indicated by the solid blue line marked with circle symbols. The red curve line presents a VIM solution.

Table 1. Relative error computed between the numerical solution and the VIM solution.

Time	0	$0.05$	0.1	0.15	0.2	0.25	0.3
Relative error	0	0.0007	0.0019	0.0032	0.0051	0.0076	0.0107

4.2. Example 2

In this example, we use a model that takes into account the elastic properties of the arteries described in [33]. This model is characterized by the following pressure expression: (26) $P(X)=\{\begin{array}{ll}\mathrm{\alpha exp}({\displaystyle \frac{X}{X_0}}-1)& \mathrm{if}X>X_0,\\ \mathrm{\alpha log}\left({\displaystyle \frac{X}{X_0}}\right)& \mathrm{if}X\le X_0,\end{array}$(26) where $\alpha =\rho c_0^2$, and $c_0$ is the velocity of small disturbances propagation along the vessel.

Let the initial condition (27) $\begin{array}{rl}& \{X(x,0)=1;Y(x,0)=x,\\ & f_X=\mathrm{\delta X},\delta >0;f_Y={\displaystyle \frac{-\beta }{X^3}}Y({\displaystyle \frac{X}{X_0}}+{\displaystyle \frac{X_0}{X}}),\beta >0.\end{array}$(27) Using Equation (27(27) $\begin{array}{rl}& \{X(x,0)=1;Y(x,0)=x,\\ & f_X=\mathrm{\delta X},\delta >0;f_Y={\displaystyle \frac{-\beta }{X^3}}Y({\displaystyle \frac{X}{X_0}}+{\displaystyle \frac{X_0}{X}}),\beta >0.\end{array}$(27) ), we get $X_0=1;Y_0=x$.

The first approximation is given $X_1=(\delta -1)t+1;Y_1=(1-(\mathrm{KR}+2\beta +2)t)x.$The second approximation $\begin{array}{rl}{X}_2(x,t)& =1+(\delta -1)t\\ & +(1+\beta +K_R/2+\delta (\delta -1)/2)t^2.\\ Y_2(x,t)& =x-a_1\mathrm{xt}+{\displaystyle \frac{2(2+K_R+2\beta {)}^2}{(1-\delta )}}x{\displaystyle \frac{t^2}{2}}\\ & -{\displaystyle \frac{\beta (1-(2+K_R+2\beta ))}{1-\delta }}{\displaystyle \frac{x}{1-(1-\delta )t}}\\ & -{\displaystyle \frac{\beta (2+K_R+2\beta )}{2(1-\delta {)}^2}}{\displaystyle \frac{x}{(1-(1-\delta )t{)}^2}}\\ & -a_2{\displaystyle \frac{x}{(1-(1-\delta )t{)}^3}}\\ & -a_3\mathrm{xlog}(\mid 1-(1-\delta )t\mid ).\end{array}$With $\begin{array}{rl}{a}_1& ={\displaystyle \frac{(2+K_R+2\beta )(K_R+4)}{1-\delta }}-{\displaystyle \frac{2(2+K_R+2\beta {)}^2}{(1-\delta {)}^2}},\\ a_2& ={\displaystyle \frac{\beta }{3(1-\delta )}}-{\displaystyle \frac{\beta (2+K_R+2\beta )}{(1-\delta {)}^2}},\\ a_3& ={\displaystyle \frac{\beta (2+K_R+2\beta )-(K_R+2)}{(1-\delta )}}\\ & -{\displaystyle \frac{(2+K_R+2\beta )(K_R+4)}{(1-\delta {)}^2}}+{\displaystyle \frac{2(2+K_R+2\beta {)}^2}{(1-\delta {)}^3}}.\end{array}$

Figure 2 illustrates the variation of the axial velocity for the time $t\in [0,1]$ and the position of the cross-section $x\in [0.3,0.9]$, for a the blood flow model taking into account the elastic properties of the arteries. Figure 3 shows the behaviour of the axial velocity as a function of time. Knowing that we took different fixed values of x for each solution to see the effect of the position of the artery cross-section on the axial velocity.

Figure 2. The VIM solution obtained by two steps of calculation for the blood flow model takes into account the elastic properties of the arteries.

Figure 3. The effect of the different position x of cross-section on the solution which was obtained by two calculation steps for the blood flow model taking into account the elastic properties of the arteries.

4.3. Example 3

We consider the mechanical model treated in [29,32,34] and describes the elastic properties of both arteries and veins given by the following pressure: (28) $P(X)=G_0[{\left({\displaystyle \frac{X}{X_0}}\right)}^m-{\left({\displaystyle \frac{X}{X_0}}\right)}^n]\mathrm{\forall }{X}_0\ne 0,$(28) where $G_0$ depends of Poisson's ratio, Young's modulus, and the wall thickness to radius ratio.

Considering the following source terms: $f_X=\mathrm{\delta X};f_Y={\displaystyle \frac{-\beta }{X^3}}Y({\displaystyle \frac{X}{X_0}}+{\displaystyle \frac{X_0}{X}}).$For the numerical study, we put $m=\frac{1}{2}$ and n=0.

Let the following initial condition: $X_0(x,t)=x;Y_0(x,t)=1.$We get the first step of VIM $X_1=x(1+\mathrm{\delta t});Y_1=1-({\displaystyle \frac{K_R}{x}}-{\displaystyle \frac{1}{x^2}}+{\displaystyle \frac{2\beta }{x^3}})t.$The second step is given by $\begin{array}{rl}{X}_2& =X_1-({\displaystyle \frac{K_R}{x^2}}-{\displaystyle \frac{2}{x^3}}+{\displaystyle \frac{6\beta }{x^4}}+x{\delta }^2){\displaystyle \frac{t^2}{2}},\\ Y_2& =1-c_1{\displaystyle \frac{t}{\delta }}-c_2{\displaystyle \frac{t^2}{2\delta }}+({\displaystyle \frac{1}{x^3}}+{\displaystyle \frac{1}{\delta }}){\displaystyle \frac{\beta }{1+\mathrm{\delta t}}}\\ & -{\displaystyle \frac{\beta }{2\delta }}({\displaystyle \frac{K_R}{x^4}}-{\displaystyle \frac{1}{x^5}}+{\displaystyle \frac{2\beta }{x^6}}){\displaystyle \frac{1}{(1+\mathrm{\delta t}{)}^2}}\\ & -{\displaystyle \frac{\beta }{3\delta x^3}}({\displaystyle \frac{K_R}{x}}+{\displaystyle \frac{1}{x^2}}+{\displaystyle \frac{2\beta }{x^3}}+\delta ){\displaystyle \frac{1}{(1+\mathrm{\delta t}{)}^3}}\\ & -c_3{\displaystyle \frac{1}{\delta }}\log (|1+\mathrm{\delta t}|).\end{array}$With $\begin{array}{rl}{c}_1& ={\displaystyle \frac{140\beta }{\delta x^8}}+{\displaystyle \frac{30(K_R-\beta )}{\delta x^6}}\\ & +{\displaystyle \frac{8[\delta (K_R+1)\beta -6\beta (3\beta +2K_R)-3]}{\delta x^5}}\\ & -{\displaystyle \frac{6K_R}{x^4}}+{\displaystyle \frac{2K_R^2(\delta -1)}{\delta x^3}},\\ c_2& =-{\displaystyle \frac{140\beta }{\delta x^8}}-{\displaystyle \frac{30K_R}{x^6}}+{\displaystyle \frac{24(6{\beta }^2+4\beta K_R+1)}{x^5}}+{\displaystyle \frac{2K_R^2}{x^3}},\\ c_3& =-{\displaystyle \frac{140\beta }{{\delta }^2{x}^8}}-{\displaystyle \frac{30(K_R-\mathrm{\delta \beta })}{{\delta }^2{x}^6}}\\ & +{\displaystyle \frac{2(48\beta K_R+72{\beta }^2+12-4\delta -\mathrm{\beta \delta }-4K_R\beta {\delta }^2)}{{\delta }^2{x}^5}}\\ & +{\displaystyle \frac{(3\delta K_R+3K_R+\delta )}{\delta x^4}}-{\displaystyle \frac{\beta K_R}{x^3}}\\ & +{\displaystyle \frac{(2K_R^2-\beta K_R{\delta }^2-2K_R^2{\delta }^2)}{{\delta }^2{x}^3}}-{\displaystyle \frac{K_R\delta }{x^2}}.\end{array}$

Figure 4 illustrates the change in axial velocity as a function of time t which varies in the interval $[0.2,0.7]$ and the position of the cross-section with $x\in [0.5,1]$.

Figure 4. The VIM solution for the mechanical model taking into account some elastic properties of both arteries and veins.

5. Discussion

In Example 4.1, the performance of the VIM solution is compared to a numerical method using the MATLAB PDEPE function. Note that this example is entirely illustrative and based on simple expressions of pressure and source terms. The outcomes show that after only two iterations of VIM, we obtain a reasonable result. Examples 4.2 and 4.3 were chosen on the basis of parameters evaluated experimentally. We can distinguish the flexibility and ability of VIM to solve the blood flow problem through compliant vessels, in the presence of the complexity of the added terms. Despite these performances, we distinguish some limitations of VIM. In these examples, going from one iteration to another can lead to more complicated integrals in the computation of the solution.

6. Conclusion

This article examined a model based on the Navier–Stokes system. This system describes the behaviour of blood flow in the arteries. Thanks to source terms and pressure expressions that take into account the properties of the elastic walls of the arteries, the model can cover certain interactions of the blood with the surfaces of the arteries. To simplify the study of the problem, a process of asymptotic reduction of the incompressible Navier–Stokes equations was applied while keeping the problem its nonlocal character. Numerically, we have proposed an approximate solution to the blood flow problem in an explicit form, using the variation iteration method. The simplicity and usability of VIM is demonstrated despite the complexity and diversity of this topic. We also used a numerical method based on Matlab's PDEPE function to discuss and evaluate the VIM. The main advantage of the VIM is that the results obtained are reasonable in two iteration steps. The VIM is perfectly capable of providing a solution despite the complexity of the pressure and/or the expression of the source terms.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Additional information

Funding

This project was supported by the University of Jeddah represented by the project UJ-21-DR-16, to which we extend our sincere thanks for its support.