Numerical modeling of the pulse wave propagation in large blood vessels based on liquid and wall interaction

The purpose of this article is to develop a non-linear, one-dimensional model of pulse wave propagation in the arterial cardiovascular system. The model includes partial differential equations resulting from the balance of mass and momentum for the fluid-filled area and the balance equation for the area of the wall and vessels. The considered mathematical model of pulse wave propagation in the thoracic aorta section takes into account the viscous dissipation of fluid energy, realistic values of parameters describing the physicochemical properties of blood and vessel wall. Boundary and initial conditions contain the appropriate information obtained from in vivo measurements. As a result of the numerical solution of the mass and momentum balance equations for the blood and the equilibrium equation for the arterial wall area, time- dependent deformation, respective velocity profiles and blood pressure were determined.


Introduction
The blood flow in the human body is mostly caused by the deformation of the vessels and other flow conduits. The interaction between the forces acting on fluid and vessel walls often leads to a non-linear pressure drop, wave-flow phenomena and the complex instabilities in the circulation. An existing need to possess a knowledge of these phenomena is some kind of a challenge for experimental research, analytic and numerical calculations [1,2,3,4]. Description of a number of numerical simulation works on aortic blood flow, assuming both rigid and elastically deformable walls as well are shown in [3]. The results of numerical calculation on the self-excited oscillation in the velocity and pressure of the selected sections of the artery, taking into account its local obstruction, is shown in article [4]. In [2] the flow of blood as a Newtonian fluid in the thoracic aorta is analysed, assuming that the wall material is characterized by isotropic linear elasticity. The known dependence of the volume flow rate in time was used in the numerical calculations as the initial condition. The purpose of this article is to develop a non-linear, one-dimensional model of pulse wave propagation in the arterial cardiovascular system. The main idea of this work involves an analysis of the average blood flow perturbation in the thoracic aorta in human, caused by a cyclically varying pressure in the inlet cross-section. Arterial segments (human thoracic aorta) are modeled using pipe sections undergoing linear elastic deformations. The model allows to perform an analysis of the performance of a healthy or partially obstructed cardiovascular system.

Problem formulation
We will consider an unsteady, laminar blood flow through the thoracic aorta segment in human. Figure  1 shows the considered flow diagram. Our considerations take into account the following assumptions:  blood vessel has the shape of a circular cylinder with changing cross-sectional area S=S(s,t), constant wall thickness b and fixed length L,  the vessel wall undergoes substantial elastic deformations in the radial direction only,  in a vessel filled with blood pressure p(s,t) and velocity υ=υ(s,t) wave propagation occur,  the blood filling the vessel is a Newtonian fluid demonstrating a constant viscosity and a slight compressibility. Given the above assumptions, the considered blood flow in the thoracic aorta is described by the following equations : Equations (1) and (2) result from the analysis of the overall mass and momentum balance equations for blood and from the equilibrium equations for an elastic blood vessel wall material. Symbols in equations (1) and (2) mean: For the thick blood vessel, constrained at both ends, the parameter C takes the form of [5]: where μ is the Poisson's ratio for the material of the arterial wall. The system of equations (1) and (2) is known in the literature for modeling water hammer [5]. To solve this system of equations the method of characteristics is most commonly used [5]. The main idea behind the method of characteristics is an effective reduction of a system of partial differential equations to an equivalent system of ordinary differential equations. It should be noted that equation (6) is valid only for a C + characteristic, described by a relationship: while equation (7) is valid along the Ccharacteristic: In equations (6) and (7) the value H(s,t) is the piezometric height of the liquid column expressing the pressure changes: From equations (8) and (9) it is apparent that the slope of the lines of characteristics depends on the fluid velocity υ(s,t).

Solving the problem
We will write differential equations (6) and (7) in the form of finite differences. Assuming small values of spatial Δs and time Δt steps, new values of H and υ at the point P of characteristics C + and Cintersection can be obtained, as shown in Figure 2.
It is worth noting that the characteristics C + and Cin Figure 2 are straight line segments due to the small value of the time step Δt. Given the above, equations (6) and (7) will be written in the form of finite differences: Equations (11) and (12) are valid along the characteristics C + and Crespectively. From the differential equations (11) and (12) we calculate υp and Hp for each time-step: (13) In equations (13) and (14) sinθ=dz/ds is positive when a blood vessel rises up the z axis ( Figure 1). In order to ensure the stability of solutions, the relationship between the time and spatial step takes the form of inequality: Initial and boundary conditions are based on measurements made in vivo with a piezoelectric tonometer (SphygmoCor, AtCor Medical). In particular, the dependence of the blood pressure p(t) at the opening of the aortic valve was obtained. Selected pressure waveforms for persons aged 27 and 78 have been presented in Figure 3.  For the purposes of the calculation, the pressure changes were converted according to equation (10): Taking into account the relation (16) in equation (12), we obtain the initial condition for velocity: In order to perform numerical integration of the equations (1) and (2) with conditions (16)-(19), in accordance with the idea of the method of characteristics, the numerical code in FORTRAN was developed. The measured changes in the pressure p(t) (H(t)) were interpolated using the spline function in the calculation process.

Results and discussion
The numerical calculation includes parameter values and physical properties of the thoracic aortic wall at T = 37ºC adopted from the literature [6]:       The time-dependent local changes in blood pressure calculated in the tenth heart cycle for the 78-year old man are shown in Figure 4. Blood pressure plots are assigned to the respective segments along the thoracic aorta. Square-shaped symbols in Figure 4 indicate a change in blood pressure in the inlet crosssection (x=s/L=0), triangle-shaped points represent the changes in blood pressure in the second segment (x=16.8/160=0.105). It can be easily observed that the predetermined profile of numerical changes in blood pressure in the second aortic segment is considered very similar to the corresponding shape preset in the inlet cross-section, characterized by a suitable, lower values of pressure. Respective local changes in blood pressure profiles in the case of 27-year old man, also calculated during the tenth cardiac cycle, are shown in Figure 5.
The time-dependent changes in blood velocity calculated in the tenth heart cycle for the 78-year old man are shown in Figure 6. Similar waveforms of blood velocity changes in equally distant segments for a 27-year old man are presented in Figure 7. From a comparison of the blood velocity waveforms in figures 6 and 7, a greater degree of severity of the flow perturbation is present in the stiffer vessel section in the 78-year old man (Eθ(78) = 18.7·10 5 Pa) compared to the corresponding flow through more elastic part of the aorta in younger man (Eθ(27) = 5.5·10 5 Pa). Figures 8 and 9 represent the spatial changes in blood pressure through the thoracic aorta segment in 78-and 27-year old man respectively, at selected moments of the tenth cardiac cycle. Due to lower heart rate in the younger man and longer cardiac cycle, the spatial pressure waveforms are selected at different time points. In figures 10 and 11 changes in blood velocity along the thoracic aorta are presented for 78-and 27-year old man respectively. The diagrams refer to the same moments of the tenth cardiac cycle which are analyzed in the spatial pressure variations in Figures 8 and 9. It can be seen that the blood velocity at the end of the considered aortic segment is achieved about the same value for the three randomly selected moments of the tenth cardiac cycle.

Final remarks and conclusions
In this paper we present numerical results of the integration of the one-dimensional differential equations of an unsteady blood flow in the human thoracic aorta. The calculations account for the pressure waveforms in the inlet cross section of the considered aortic segment, measured in vivo in two people of different age. The realistic, taken from the literature, values of parameters describing the blood and arterial wall properties are also included. Elastic deformation of the wall caused by pulsatile blood flow was taken into account by varying the value of the modulus of elasticity E, dependent on the human age. The analysis obtained from the numerical results have been limited to the tenth of the cardiac cycle. As a result from the analysis it can be concluded that in cardiac cycle a strong stabilization of the blood flow parameters is found in the considered aorta segment, both in the elderly and the young man. Elastic deformation of the considered aortic wall also stabilizes. The maximum relative changes in the aortic diameter are not greater than ± 5%.