Propagation of Nonlinear Pressure Waves in Blood

Abstract

The propagation of weakly nonlinear pressure waves in a fluid-filled elastic tube has been investigated. The reductive perturbation method has been employed to derive the Korteweg-de Vries equation for small but finite amplitude. The effect of the final inner radius of the tube 𝑟_𝑓 on the basic properties of the soliton wave was discussed. Moreover, the conditions of stability and the soliton existence via the potential and the corresponding phase portrait were computed. The applicability of the present investigation to flow problems in arteries is discussed.

1. Introduction

The intermittent ejection of blood from the left ventricle produces pressure waves that flow in the arterial tree. Experimental data found that the flow velocity in blood vessels largely depends on the elastic properties of the vessel wall and they propagate towards the periphery with a characteristic diagram [1]. In arterial mechanics, the propagation of pressure waves in fluid-filled distensible tubes has been theoretically studied by several researchers [2–4]. Experimental observation for the simultaneous changes in amplitude of the pressure waves at five sites from the ascending aorta to the saphenous artery in dogs showed that the pulsatile character of the blood is soliton waves [5]. Yomosa [6] investigated the nonlinear propagation of solitary waves in large blood vessels. He found that the pulse waves of pressure and flow propagating through the arteries can be described as solitary waves excited by cardiac ejections of blood, and the features of the pulse wave such as “peaking” and “steepening” are interpreted in the viewpoint of soliton. Later, R. M. Shoucri and M. M. Shoucri studied the application of the method of characteristics of shock waves in models of blood flow in the Aorta [7]. Recently, Gaik and Demiray [8] treated the arteries as an incompressible prestressed thin walled elastic tube with a stenosis and the blood as a Newtonian fluid with variable viscosity, which vanishes on the arterial wall, and it takes the maximum value at the center of the artery. They studied the propagation of weakly nonlinear waves in the long wave approximation by the use of the perturbation methods [9]. Many authors examined the stability and the soliton existence condition via the potential and the corresponding phase portrait [10–12]. Elwakil et al. [10] studied the propagation of solitary electron acoustic waves in unmagnetized collisionless plasma. They found that, there are saddle and two-center equilibrium state in the phase plane and there is two-finite separatrix going from a saddle returning to it in the phase portrait. The evolution of small but finite-amplitude propagating solitary structures was studied the by means of Korteweg-de Vries equation (KdV). Many authors have been using the KdV equation to study the properties of solitary waves [13, 14]. The major topic of this work is to study the propagation of pressure waves in weakly nonlinear waves in a fluid-filled elastic tube. This paper is organized as follows: in Section 2, we present the basic set of fluid equations governing our model. In Section 3, long wave approximation is used to drive KdV equation, and the solution for KdV equation is obtained. In Section 4, some discussions and conclusions are given.

2. Basic Equations

To drive the equation of motion of the tube, let us consider a circular cylindrical long tube of radius𝑅₀with a uniform inner pressure 𝑃₀, let the axial stretch ratio𝜆_𝑧and𝑟₀be the radius of the cylindrical tube after such a finite static deformation. The position vector of a generic point of the tube may be described by the following: 𝑟=𝑟₀+𝑢^∗𝑒_𝑟+𝑧^∗𝑒_𝑧,𝑧^∗=𝜆_𝑧𝑍, (1) where𝑢^∗(𝑧^∗,𝑡^∗),𝑒_𝑟,𝑒_𝜃,𝑒_𝑧,𝑧^∗ and𝑍are a finite time dependent radial displacement, unit base vectors in the cylindrical polar coordinates, the spatial coordinate in the intermediate configuration and the axial coordinate of a point in the undeformed configuration, respectively. The axial displacement is neglected in view of the external tethering. The unit tangent vector𝑡to the meridional curve and the unit exterior normal vector 𝑛 to the deformed membrane are given by 1𝑡=Λ𝜕𝑢^∗𝜕𝑧^∗𝑒_𝑟+𝑒_𝑧1,𝑛=Λ𝑒_𝑟−𝜕𝑢^∗𝜕𝑧^∗𝑒_𝑧, (2) whereΛis defined by Λ=1+𝜕𝑢^∗𝜕𝑧^∗21/2. (3)

The stretches in the axial and circumferential directions may be given as follows: 𝜆₁=𝜆_𝑧Λ,𝜆₂=𝜆_𝜃+𝑢^∗𝑅₀, (4) where𝜆_𝜃=𝑟₀/𝑅₀is the stretch ratio in the circumferential direction after finite static deformation.

Let 𝐹₁ and 𝐹₂ be the membrane forces acting along each unit length of the meridional and circumferential curves of the tube, respectively. The equation of motion of the tube in the radial direction is given by 𝜕𝜕𝑧^∗𝜇1+(𝜕𝑢^∗/𝜕𝑧^∗)^21/2𝜕Σ𝜕𝜆₁𝜕𝑢^∗𝜕𝑧^∗−𝜇𝜆_𝑧𝑅₀𝜕Σ𝜕𝜆₂+𝜆_𝜃+𝑢^∗𝑅₀𝑃^∗𝐻=𝜌₀𝜆_𝑧𝜕²𝑢^∗𝜕𝑡^∗2, (5) where𝜇is the shear modulus of the tube material,𝜇Σis the strain energy density function,𝐻is the initial tube thickness,𝑃^∗is the fluid pressure, and𝜌₀is the mass density of the tube material. In order to complete the field equations one must know the value of the fluid pressure𝑃^∗. Therefore, (5) is to be complemented with the equations governing the blood fluid. Blood is known to be an incompressible non-Newtonian fluid. The main factor for blood to behave like a non-Newtonian fluid is the deformability of red blood cells and the level of cell concentration (hematocrit ratio). When blood flows in arteries the red cells move to the central region of the artery and, thus, the hematocrit ratio is reduced near the arterial wall, where the shear rate is quite high, as can be seen from Poiseuille flow. In other words, experimental observations indicate that when the shear rate is high, blood behaves like a Newtonian fluid. The ratio of the viscous terms to the nonlinear term is(𝜇_𝑣/𝜌_𝑓)(𝜕²𝜐^∗/𝜕𝑧²)/𝜐^∗(𝜕𝜐^∗/𝜕𝑧)≈𝜇_𝑣/𝜌_𝑓𝑉^′𝑇^′𝜐^∗≈𝑂(10⁻⁵), considering𝜇_𝑣=0.04pand𝜌_𝑓=1.05g/cm³. Therefore, the viscous effect in comparison to the nonlinear effect can be neglected. Based on these observations, we assume that blood is an incompressible inviscid fluid whose equations of axially symmetrical motion in the cylindrical polar coordinates are given by 𝜕𝑉^∗_𝑟+𝑉𝜕𝑟𝑟^∗𝑟+𝜕𝑉^∗_𝑧∗𝜕𝑧^∗=0,𝜕𝑉^∗_𝑟𝜕𝑡^∗+𝑉^∗_𝑟𝜕𝑉^∗_𝑟𝜕𝑟+𝑉^∗_𝑧∗𝜕𝑉^∗_𝑟𝜕𝑧^∗+1𝜌_𝑓𝜕𝑃𝜕𝑟=0,𝜕𝑉^∗_𝑧∗𝜕𝑡^∗+𝑉^∗_𝑟𝜕𝑉^∗_𝑧∗𝜕𝑟+𝑉^∗_𝑧∗𝜕𝑉^∗_𝑧∗𝜕𝑧^∗+1𝜌_𝑓𝜕𝑃𝜕𝑧^∗=0, (6) where𝑉^∗_𝑟and𝑉^∗_𝑧∗are the fluid velocity components in the radial and axial directions, respectively,𝜌_𝑓is the mass density of the fluid, and𝑃is the fluid pressure function. These field equations must satisfy the following boundary conditions: 𝑉^∗_{𝑟𝑟=𝑟𝑓}=𝜕𝑢^∗𝜕𝑡^∗+𝜕𝑢^∗𝜕𝑧^∗𝑉^∗_{𝑧∗𝑟=𝑟𝑓},𝑃_{𝑟=𝑟𝑓}=𝑃^∗. (7)

Here𝑃is the fluid pressure function,𝜌_𝑓is the fluid mass density, and𝑟_𝑓is the final inner radius of the tube. where𝜆_𝜃=(𝑟₀/𝑅₀)is the stretch ratio in the circumferential direction after the finite static deformation, where from (5) it is obtained that 𝑃^∗=𝜌₀𝑅₀𝐻𝜆_𝑧𝑅₀𝜆_𝜃+𝑢^∗𝜕²𝑢^∗𝜕𝑡^∗2+𝜇𝐻𝜆_𝑧𝑅₀𝜆_𝜃+𝑢^∗𝜕Σ𝜕𝜆₂−𝑅₀𝐻𝑅₀𝜆_𝜃+𝑢^∗𝜕𝜕𝑧^∗×𝜇1+(𝜕𝑢^∗/𝜕𝑧^∗)^21/2𝜕Σ𝜕𝜆₁𝜕𝑢^∗𝜕𝑧^∗. (8)

In general, the strain energy densityΣis a function of𝜆₁and𝜆₂. For our purposes, we shall assume thatΣis analytic in𝜆₁and𝜆₂and can be expanded into power series of the following form: 𝑃^∗=𝛽₁𝜇𝐻𝑅²₀𝑢^∗+𝜌₀𝐻𝜆_𝜃𝜆_𝑧𝜕²𝑢^∗𝜕𝑡^∗2−𝛼₀𝜕𝜇𝐻²𝑢^∗𝜕𝑧^∗2+𝛽₂𝜇𝐻𝑅³₀𝑢^∗2−𝜌₀𝐻𝜆²_𝜃𝜆_𝑧𝑅₀𝑢^∗𝜕²𝑢^∗𝜕𝑡^∗2−𝛼₁𝜇𝐻𝑅₀𝜕𝑢^∗𝜕𝑧^∗2+𝛼₀𝜇𝐻𝜆_𝜃𝑅₀−2𝛼₁𝜇𝐻𝑅₀𝑢^∗𝜕²𝑢^∗𝜕𝑧^∗2+𝜇𝐻𝑅₀𝑃₀, (9) where the coefficients𝛼₀,…,𝛽₂are defined by 𝛼₀=1𝜆_𝜃𝜆_𝑧𝜕Σ𝜕Λ_𝑢∗=0,𝛽₁=𝑅₀𝜆_𝜃𝜆_𝑧𝑅₀𝜕²Σ𝜕𝑢^∗2−1𝜆_𝜃𝜕Σ𝜕𝑢^∗_𝑢∗=0,𝛽₂=𝑅³₀2𝜆_𝜃𝜆_𝑧𝜕³Σ𝜕𝑢^∗3_𝑢∗=0−𝛽₁𝜆_𝜃. (10)

Equations (6)–(9) give sufficient relations to determine the unknowns𝑢^∗,𝑉^∗_𝑟,𝑉^∗_𝑧∗, and 𝑃.

3. Long Wave Approximation

The reductive perturbation method is used to study the propagation of small but finite amplitude. Let us introduce the following types of stretched coordinates [15]: 𝜉=𝜖^1/2𝑧^∗−𝑔𝑡^∗,𝜏=𝜖^3/2𝑧^∗, (11) where𝜖is a small parameter measuring the smallness of nonlinearity, dissipation, and dispersion,𝑔is the phase velocity in the longwave approximation to be determined later. All physical quantities appearing in (6) and (7) are expanded as power series in 𝜖 about their equilibrium values as: 𝑉^∗_𝑟=^∞_𝑛=1𝜖^1/2+𝑛𝑉_𝑟^∗(𝑛)(),𝑉𝜉,𝜏,𝑟_∗𝑧∗=^∞_𝑛=1𝜖^𝑛𝑉_𝑧^∗(𝑛)_∗(𝜉,𝜏,𝑟),𝑃=^∞_𝑛=1𝜖^𝑛𝑃_(𝑛)(),𝑢𝜉,𝜏,𝑟_∗=^∞_𝑛=1𝜖^𝑛𝑢^∗_𝑛(𝜉,𝜏,𝑟). (12)

We impose the boundary conditions as follows: 𝑃_{𝑟=𝑟𝑓}𝜌=𝜖⁰𝐻𝑔²𝜆_𝜃𝜆_𝑧−𝛼₀𝜕𝜇𝐻²𝑢^∗𝜕𝜉²+𝛽₁𝜇𝐻𝑅²₀𝑢^∗−𝜖³𝛼₀𝜕𝜇𝐻²𝑢^∗𝜕𝜏²−2𝜖²𝛼₀𝜕𝜇𝐻²𝑢^∗+𝛽𝜕𝜉𝜕𝜏²𝜇𝐻𝑅³₀𝑢^∗2𝜌−𝜖⁰𝐻𝑔²𝜆²_𝜃𝜆_𝑧𝑅₀𝑢^∗𝜕²𝑢^∗𝜕𝜉²−𝛼₁𝜇𝐻𝑅₀𝜖𝜕𝑢^∗𝜕𝜉²+2𝜖²𝜕𝑢^∗𝜕𝜉𝜕𝑢^∗𝜕𝜏+𝜖³𝜕𝑢^∗𝜕𝜏²+𝛼₀𝜇𝐻𝜆_𝜃𝑅₀−2𝛼₁𝜇𝐻𝑅₀×𝑢^∗𝜖𝜕²𝑢^∗𝜕𝜉²+𝜖³𝜕²𝑢^∗𝜕𝜏²+𝜖²𝜕²𝑢^∗+𝜕𝜉𝜕𝜏𝜇𝐻𝑅₀𝑃₀,𝑉^∗_{𝑟𝑟=𝑟𝑓}=−𝜖^1/2𝑔𝜕𝑢^∗𝜕𝜉+𝜖^1/2𝜕𝑢^∗𝜕𝜉+𝜖𝜕𝑢^∗𝑉𝜕𝜏^∗𝑧∗_{𝑟=𝑟𝑓}. (13)

Substituting (11) and (12) into (6) and (9) and equating coefficients of like powers of𝜖. Then, from the lowest-order equations in 𝜖 the following results are obtained: 1𝜌_𝑓𝜕𝑃₁𝜕𝑟=0,−𝑔𝜕𝑉_𝑧^∗(1)_∗+1𝜕𝜉𝜌_𝑓𝜕𝑃₁𝜕𝜉=0,𝜕𝑉_𝑟^∗(1)+𝑉𝜕𝑟𝑟^𝑟∗(1)+𝜕𝑉_𝑧^∗(1)_∗𝜕𝜉=0. (14)

And the boundary conditions can be written as follows: 𝑃_{1𝑟=𝑟𝑓}=𝛽₁𝜇𝐻𝑅²₀𝑢^∗₁,𝑉_𝑟^∗(1)_{𝑟=𝑟𝑓}=−𝑔𝜕𝑢^∗₁.𝜕𝜉 (15)

From the solution of the sets (14) and (15) we obtain that 𝑢^∗₁()()=𝑅=𝑈𝜉,𝜏,𝑈𝜉,𝜏²⁰𝛽₁𝜇𝐻𝑃₁(),𝑉𝜉,𝜏_{𝑧∗(1)∗}=1𝑔𝜌_𝑓𝑃₁(𝜉,𝜏),𝑉_𝑟^∗(1)1=−2𝑔𝜌_𝑓𝜕𝑃₁(𝜉,𝜏)𝑔𝜕𝜉𝑟,₂=𝛽₁𝜇𝐻𝑟_𝑓2𝜌_𝑓𝑅²₀, (16) where 𝑈(𝜉,𝜏) is an unknown function whose governing equation will be obtained later and𝑔is the phase velocity in the longwave approximation.

Considering now the coefficients of𝑂(𝜖²), we derive with the aid of (16) the following set of equations 12𝜌_𝑓𝜕²𝑃₁(𝜉,𝜏)𝜕𝜉²1𝑟+𝜌_𝑓𝜕𝑃₂𝜕𝑟=0,−𝑔𝜕𝑉_𝑧^∗(2)_∗+1𝜕𝜉𝑔²𝜌²_𝑓𝑃₁𝜕𝑃₁+1𝜕𝜉𝜌_𝑓𝜕𝑃₂+1𝜕𝜉𝜌_𝑓𝜕𝑃₁𝜕𝜏=0,𝜕𝑉_𝑟^∗(2)+𝑉𝜕𝑟𝑟^𝑟∗(2)+𝜕𝑉_𝑧^∗(2)_∗+1𝜕𝜉𝑔𝜌_𝑓𝜕𝑃₁𝜕𝜏=0. (17)

And the boundary conditions can be written as follows 𝑃_{(2)𝑟=𝑟𝑓}=𝜌₀𝐻𝑔²𝜆_𝜃𝜆_𝑧−𝛼₀𝑅𝜇𝐻²⁰𝛽₁𝜕𝜇𝐻²𝑃₁𝜕𝜉²+𝛽₁𝜇𝐻𝑅²₀𝑢^∗₂+𝑅₀𝛽₂𝜇𝐻𝛽²₁𝑃₁²,𝑉_𝑟^∗(2)_{𝑟=𝑟𝑓}=−𝑔𝜕𝑢^∗₂+𝑅𝜕𝜉𝑔𝜌_𝑓𝛽²⁰₁𝜇𝐻𝑃₁𝜕𝑃₁𝜕𝜉_{𝑟=𝑟𝑓}. (18)

From (17-a) we have 𝑃₂1=−4𝜕²𝑃₁𝜕𝜉²𝑟². (19)

Eliminating the second order perturbed quantities𝑃₂and𝑉_𝑧^∗(2)_∗, the desired KdV equation is obtained as follows: 𝜕𝑃₁𝜕𝜏+𝐴𝑃₁𝜕𝑃₁𝜕𝜕𝜉+𝐵³𝑃₁𝜕𝜉³=0, (20) where the coefficients𝐴and𝐵are defined by 𝑅𝐴=²⁰𝛽₁𝛽𝜇𝐻¹𝜇𝐻2𝑔²𝜌_𝑓𝑅²₀+2𝑔²𝜌_𝑓𝛽₂𝑅₀𝛽²₁𝜇𝐻𝑟_𝑓+1𝑟_𝑓,𝑔𝐵=²𝜌_𝑓𝑅²₀𝑟_𝑓4𝛽₁+𝜌𝜇𝐻⁰𝜌_𝑓𝑔⁴𝑅⁴₀𝜆_𝜃𝜆_𝑧𝛽²₁𝜇²𝐻𝑟_𝑓−𝑟²_𝑓−𝛼16⁰𝑔²𝑅⁴₀𝜌_𝑓𝛽²₁𝜇𝐻𝑟_𝑓. (21) Our system can support two kinds of potential structure depending on the sign of the coefficient of the nonlinear term (𝐴). A stationary solitary wave solution of the KdV equation can be obtained by transforming the space variable to 𝜂=(𝜁−𝑣𝜏), (22) where𝑣is velocity of the wave. This has been done by imposing the boundary conditions for localized perturbations, namely, 𝑃₁=0, 𝑑𝑃₁/𝑑𝜂=0, and𝑑²𝑃₁/𝑑𝜂²=0 for 𝜂→±∞. Thus, the steady state solution of (22) can be expressed as 𝑃₁=𝑃₀sech²𝜂Δ, (23) where the soliton amplitude𝑃₀and the soliton widthΔare given by 𝑃₀=3𝑣𝐴,Δ=2𝐵𝑣. (24)

4. Numerical Results and Discussion

The weakly nonlinear pressure waves in a fluid-filled elastic tube have been investigated. To make our result physically relevant, numerical studies have been made using parameters close to those values corresponding to actual biologically relevant parameters for experimental data in dogs [6, 16]. The effect of the final inner radius of the tube𝑟_𝑓on the basic properties of the amplitude and width of the soliton is shown in Figure 1. It is obvious that the magnitude of the soliton amplitude and width increases with𝑟_𝑓. From (20), the nonlinear equation of motion can be obtained as 12𝑑𝑃₁𝑑𝜂²+𝑉𝑃₁=0. (25)

Figure 1 

The soliton profile for expression (23) for different values of𝑟𝑓for𝛽1=296.105,𝛽2=991.496,𝜌𝑓=1.05gm/cm3,𝑅0=0.38 cm,𝜌0=1.03 gm/cm3,𝜆𝑧=𝜆𝜃=1.6,𝛼0=78.692,𝐻=2×10−2cm,𝜇=0.4, and𝑣=8cm/s.

This equation can be regarded as an “energy integral” of an oscillating particle of unit mass, with a velocity𝑑𝑃₁/𝑑𝜂and position𝑃₁in a potential 𝑉𝑃₁=𝑃₁−𝑣6𝐵2𝐵𝑃²₁. (26)

The kinetic energy of oscillating particle is shown in Figure 2. A necessary condition for the existence of the solitary waves is that𝑑²𝑉/𝑑𝑃²₁<0for𝑃₁=0. A value of 𝑑²𝑉/𝑑𝑃²₁>0predicts the formation of a shock wave. The potential and the corresponding phase portrait are shown in Figure 3. Stable solitons will exist when𝑣/2𝐵>0; otherwise stable solitons will not be present. The application of our model might be particularly interesting in the new observations for the biological experimental data.

Figure 2 

The kinetic energy for𝑟𝑓=0.75,𝛽1=296.105,𝛽2=991.496,𝜌𝑓=1.05gm/cm3,𝑅0=0.38 cm,𝜌0=1.03 gm/cm3,𝜆𝑧=𝜆𝜃=1.6,𝛼0=78.692,𝐻=2×10−2cm,𝜇=0.4, and𝑣=8 cm/s.

Figure 3 

The potential and the corresponding phase portraitfor𝑟𝑓=0.75,𝛽1=296.105,𝛽2=991.496,𝜌𝑓=1.05gm/cm3,𝑅0=0.38 cm,𝜌0=1.03 gm/cm3,𝜆𝑧=𝜆𝜃=1.6,𝛼0=78.692,𝐻=2×10−2cm,𝜇=0.4, and𝑣=8 cm/s.