- Email: [email protected]
Forced Korteweg–de Vries equation in an elastic tube filled with an inviscid fluid
International Journal of Engineering Science 44 (2006) 621–632 www.elsevier.com/locate/ijengsci
Forced Korteweg–de Vries equation in an elastic tube filled with an inviscid fluid Tay Kim Gaik
*
Panel of Mathematics, Sciences Studies Center, University College Technology of Tun Hussein Onn, 86400 Parit Raja, Batu Pahat, Johor, Malaysia Received 19 January 2006; accepted 13 April 2006 Available online 17 July 2006
Abstract In the present work, treating the arteries as a prestressed thin walled elastic tube with a stenosis and the blood as an inviscid fluid, we have studied the propagation of weakly nonlinear waves in such a composite medium, in the long wave approximation, by use of the reductive perturbation method [C.S. Gardner, G.K. Morikawa, Similarity in the asymptotic behavior of collision-free hydromagnetic waves and water waves, Courant Institute Math. Sci. Report, NYO-9082 (1960) 1–30, T. Taniuti, C.C. Wei, Reductive perturbation method in non-linear wave propagation I, J. Phys. Soc. Jpn., 24 (1968) 941–946]. We obtained the forced Korteweg–de Vries (FKdV) equation with variable coefficients as the evolution equation. By use of the coordinate transformation, it is shown that this type of evolution equation admits a progressive wave solution with variable wave speed. As might be expected from physical consideration, the wave speed reaches its maximum value at the center of stenosis and gets smaller and smaller as we go away from the center of the stenosis. The variations of radial displacement and the fluid pressure with the distance parameter are also examined numerically. The results seem to be consistent with Bernoulli’s law for inviscid fluid. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Solitary waves; Elastic tubes with stenosis
1. Introduction Due to its applications in arterial mechanics, the propagation of pressure pulses in fluid-filled distensible tubes has been studied by several researchers [1,2]. Most of the works on wave propagation in compliant tubes have considered small amplitude waves ignoring the nonlinear effects and focused on the dispersive character of waves (see [3–5]). However, when the nonlinear terms arising from the constitutive equations and kinematical relations are introduced, one has to consider either finite amplitude, or small-but-finite amplitude waves, depending on the order of nonlinearity.
*
Tel.: +07 4536721/4536701; fax: +07 4536051. E-mail address: [email protected]
0020-7225/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2006.04.008
622
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
The propagation of finite amplitude waves in fluid-filled elastic or viscoelastic tubes has been examined, for instance, by Rudinger [6], Ling and Atabek [7], Anliker et al. [8] and Tait and Moodie [9] by using the method of characteristics, in studying the shock formation. On the other hand, the propagation of smallbut-finite amplitude waves in distensible tubes has been investigated by Johnson [10], Hashizume [11,12], and Yomosa [13]. In all these works [10–13], the effect of initial deformation is neglected. Recently in a series of works of Demiray, Antar and Bakirtas (see [14–28]) in which they treated artery as incompressible prestressed thin isotropic elastic, thick viscoelastic or tapered elastic tube filled with inviscid, viscous or layered fluid as blood, using approximate method on fluid equations and reductive perturbation method in the long-wave approximation, they obtained various evolution equations of Korteweg–de Vries, Burgers and Korteweg–de Vries–Burgers type equations. In all previous works, they treated the arteries as circularly cylindrical long thin tubes with a constant cross-section. However due to decomposition of fat or cholesterol in artery over time, the artery become narrower and may have variable radius along the axis of the tube. Thus, in this work, treating the arteries as an incompressible prestressed thin walled elastic tube with a stenosis and the blood as an incompressible inviscid fluid, we have studied the propagation of weakly nonlinear waves in such a composite medium, in the long wave approximation, by use of the reductive perturbation method [29,30]. We obtained the forced Korteweg–de Vries (FKdV) equation with variable coefficients as the evolution equation. By use of the coordinate transformation, it is shown that this type of evolution equation admits a progressive wave solution with variable wave speed. As might be expected from physical consideration, the wave speed reaches its maximum value at the center of stenosis and gets smaller and smaller as one goes away from the center of the stenosis. The variations of radial displacement and the fluid pressure with the distance parameter are also examined numerically. The results seem to be consistent with Bernoulli’s law for inviscid fluid. 2. Basic equations and theoretical preliminaries In this section, we shall give the derivation of the field equations of an elastic tube, which is considered to be a model for an artery, and an inviscid fluid, which is considered to be a model for blood. 2.1. Equations of tube In this sub-section, we shall derive the governing equations of an elastic tube filled with an inviscid fluid. Such a combination of a solid and a fluid is considered to be a model for blood flow in arteries. For a healthy human being, the systolic pressure is about 120 mm Hg, and the diastolic pressure is around 80 mm Hg. This means that the arteries are subjected to a mean pressure P0 = 100 mm Hg, and in the course of blood flow, a dynamical pressure increment DP = ±20 mm Hg is added on this initial field. Moreover, experimental studies [2] revealed that the arteries are also subjected to an initial axial stretch kz, which is about kz = 1.6. These observations show that the arteries are initially subjected to static deformation both in the radial and the axial directions, and a dynamical pressure (or a radial displacement u*) is superimposed on this initial deformation. Due to the external tethering in the axial direction, the effect of axial displacement is neglected. Now, we consider a thin and long tube of circular cross-section with radius R*(Z*) in the cylindrical polar coordinates (R*, H, Z*). Then, the position vector of a point on the tube may be described by R ¼ R ðZ Þer þ Z ez ;
ð1Þ
where er, eh and ez are the unit base vectors in the cylindrical polar coordinates and Z* is the axial coordinates of a material point in the natural state. The arclengths along the meridional and circumferential curves are given by " 2 #1=2 dR dS Z ¼ 1 þ dZ ; dS H ¼ R ðZ Þ dH: ð2Þ dZ
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
623
Motivated with the experimental observations [2], we shall assume that the elastic tube is subjected to an axial stretch ratio kz, and the static pressure P 0 ðZ Þ. Then, the deformation may be described by r0 ¼ r ðz Þer þ z ez ;
z ¼ kz Z ;
ð3Þ
where z* is the axial coordinate after static deformation, and r*(z*) characterizes the variable radius after this static deformation. Then the arclengths after static deformation along the meridional and circumferential directions are given by 0
2 1=2
ds0z ¼ ½1 þ ðr Þ
dz ;
dsh ¼ r ðz Þ dh;
where a prime denotes the differentiation of the corresponding field variable with respect to Then, the stretch ratios along the meridional and circumferential curves are given by 0
k01 ¼
ds0z ½1 þ ðr Þ2 1=2 ¼ kz ; 0 2 1=2 dS Z ½1 þ k2z ðR Þ
k02 ¼
ð4Þ z *.
ds0h r ðz Þ ¼ : dS H R ðz Þ
ð5Þ
Upon this initial static deformation, we shall superimpose a finite dynamical radial displacement u*(z*, t*), where t* is the time parameter, but, in view of the external tethering in the axial direction, the axial displacement is assumed to be negligible. Then, the position vector r of a generic point on the tube may be described by r ¼ ½r ðz Þ þ u ðz ; t Þer þ z ez : The arc lengths along the deformed meridional and circumferential curves are respectively given by " #1=2 2 ou 0 dsz ¼ 1 þ r þ dz ; dsh ¼ ½r ðz Þ þ u ðz ; t Þ dh: oz
ð6Þ
ð7Þ
Then, the stretch ratios along the meridional and circumferential curves in the final configuration read, respectively, by h i1=2 0 2 1 þ r þ ou oz r ðz Þ þ u ðz ; t Þ k1 ¼ k z : ð8Þ ; k ¼ 2 0 2 1=2 R ðz Þ ½1 þ k2z ðR Þ The unit tangent vector t along the deformed meridional curve and the unit exterior normal vector n to the deformed membrane are given by 0 ou 0 r þ oz er þ ez er r þ ou ez oz t¼ ; n¼ ; ð9Þ K K where the function K is defined by " 2 #1=2 ou 0 K¼ 1þ r þ dz : oz The material that we shall consider is assumed to be incompressible. This condition imposes the following restriction on the thickness H, and h, before and after deformation respectively h¼
H : k1 k2
ð10Þ
Let T1 and T2 be the membrane forces along the meridional and circumferential curves, respectively. Then, the equation of the radial motion of a small tube element placed between the planes z* = const, z* + dz* = const, h = const and h + dh = const may be given by 2 o T1 ou H 2 0 0 2 1=2 o u ðr þ u Þ r þ K þ ðr þ u ÞP ¼ q R ðz Þ½1 þ k ðR Þ ; T 2 0 z oz K kz oz ot2
ð11Þ
624
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
where q0 is the mass density of the membrane material, and P* is the fluid pressure on the inner surface of the tube. Let lR be the strain energy density function of the tube material, where l is the shear modulus. Then, the membrane forces T1 and T2 may be expressed in terms of the stretch ratios as T1 ¼
lH oR ; k2 ok1
T2 ¼
lH oR : k1 ok2
ð12Þ
Introducing Eq. (12) into Eq. (11), the equation of motion of the tube in the radial direction takes the following form o lHR 0 ou oR lH 0 2 1=2 oR r þ ½1 þ k2z ðR Þ oz kz ok2 K oz ok1 þ ðr þ u ÞP q0
2 H 0 2 1=2 o u R ðz Þ½1 þ k2z ðR Þ ¼ 0: kz ot2
ð13Þ
2.2. Equations of fluid In general, blood is known to be an incompressible non-Newtonian fluid. However, in the course of flow in large arteries, the red blood cells in the vicinity of arterial wall move to the central region of the artery so that hematocrit ratio becomes quite low near the arterial wall, which results in lower viscosity in this region. Moreover, due to high shear rate near the arterial wall the viscosity of blood is further reduced. Therefore, for flow problems in large blood vessels, as a first approximation, blood may be treated as an inviscid fluid. The axially symmetrical equation of motion of an inviscid fluid in the cylindrical polar coordinates (r, h, z*) may be given by oV r oV r oV 1 oP ¼ 0; þ V z r þ þ V r qf or ot or oz
ð14Þ
oV z oV z oV 1 oP þ V z z þ þ V r ¼ 0; qf oz ot or oz
ð15Þ
oV r V r oV z þ þ ¼0 or r oz
ð16Þ
ðincompressibilityÞ;
with the boundary conditions V
r jr¼rf
ou ou 0 ¼ þ r þ V z jr¼rf ; ot oz
P jr¼rf ¼ P ;
ð17Þ
where V r and V z are the components of the fluid velocity vector in the radial and axial direction, respectively, qf is the mass density of the fluid, P is the pressure function and rf = r* + u* is the final inner radius of the tube. At this stage it is convenient to introduce the following non-dimensional quantities R0 t ¼ t; z ¼ R0 z; u ¼ R0 u; V r ¼ c0 V r ; V z ¼ c0 V z ; r ¼ R0 x; P ¼ qf c20 p c0 ð18Þ q0 H lH 2 2 P ¼ qf c0 ; c0 ¼ ; R ðz Þ ¼ R0 ½1 F ðzÞ; r ðz Þ ¼ R0 ½kh f ðzÞ; p; m ¼ qf R 0 qf R0 where kh = r0/R0 is the initial stretch ratio, r0 is the radius of the origin after finite static deformation, R0 is the initial reference radius and c0 is the Moens–Korteweg wave speed.
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
625
Introducing (18) into Eqs. (13)–(17), the following non-dimensional equations are obtained p¼
m½1 F ðzÞ½1 þ k2z ðF 0 Þ2 1=2 o2 u ½1 þ k2z ðF 0 Þ2 1=2 oR þ ot2 kz ðkh f ðzÞ þ uÞ ok2 kz ðkh f ðzÞ þ uÞ ( ) 1 o ½1 F ðzÞðf 0 þ ou=ozÞ oR ; ðkh f ðzÞ þ uÞ oz ½1 þ ðf 0 þ ou=ozÞ2 1=2 ok1
ð19Þ
oV r oV r oV r o p ¼ 0; þVr þVz þ ox ot ox oz oV z oV z oV z o p ¼ 0; þVr þVz þ oz ot ox oz oV r V r oV z þ þ ¼ 0; ox x oz with the boundary conditions ou ou þ f 0 þ V r jx¼kh f ðzÞþu ¼ V z jx¼kh f ðzÞþu ; ot oz
ð20Þ ð21Þ ð22Þ
pjx¼kh f ðzÞþu ¼ p:
ð23Þ
The Eqs. (19)–(23) give sufficient relations to determine the field quantities u, Vr, Vz and p completely. 3. Long wave approximation In this section, we shall examine the propagation of small-but-finite amplitude waves in a fluid-filled thin elastic tube whose non-dimensional governing equations are given in Eqs. (19)–(23). For this, we adopt the long wave approximation and employ the reductive perturbation method [29,30]. For this type of problems, it is convenient to introduce the following type of stretched coordinates n ¼ 1=2 ðz ctÞ;
s ¼ 3=2 z;
ð24Þ
where is a small parameter measuring the weakness of nonlinearity and dispersion, whereas c is the scale parameter to be determined from the solution. Solving z in terms of s we get z ¼ 3=2 s:
ð25Þ
Introducing (25) into the expression of F(z) and f(z), we obtain ^ sÞ; F ð3=2 sÞ ¼ Gðn;
f ð3=2 sÞ ¼ g^ðn; sÞ:
In order to take the effect of stenosis into account, F(z) and f(z) must be of order of ^ sÞ and g^ðn; sÞ have the following form work, we shall assume that Gðn; ^ sÞ ¼ GðsÞ; Gðn;
ð26Þ 5/2
. For the present
g^ðn; sÞ ¼ gðsÞ:
Introducing the following differential relations o o o o 1=2 o 1=2 ¼ c ; ¼ þ ot on oz on os
ð27Þ
ð28Þ
into the Eqs. (19)–(23), we obtain 2 1=2 2 1=2 mc2 ½1 GðsÞ½1 þ 5 k2z ðdG Þ o2 u ½1 þ 5 k2z ðdG Þ oR o o ds ds þ þ os kz ðkh g þ uÞ kz ðkh g þ uÞ ok2 ðkh g þ uÞ on on2 ( ) ½1 GðsÞ dg ou ou oR þ þ 2 ; 1=2 ou þ 3=2 ou Þ2 1=2 ds on os ok1 ½1 þ ð5=2 dg þ ds on os oV r oV r oV r op 1=2 oV r 1=2 ¼ 0; þVr þ Vz þ c þ ox on ox on os p¼
ð29Þ ð30Þ
626
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
oV z oV z oV z oV z p op 1=2 1=2 o þ ¼ 0; þVr þ Vz þ þ c on os on ox on os oV r V r oV z oV z þ þ 1=2 þ ¼ 0; ox x on os 1=2
ð31Þ ð32Þ
with the boundary conditions V r jx¼kh gþu ¼ 1=2 c
ou dg ou ou þ 1=2 2 þ þ V z jx¼kh gþu ; on ds on os
pjx¼kh gþu ¼ p:
ð33Þ
For the long wave limit, it is assumed that the field quantities may be expanded into asymptotic series of as V r ¼ 1=2 ðV rð1Þ þ 2 V rð2Þ þ Þ;
u ¼ u1 þ 2 u2 þ ;
p ¼ p0 þ p1 þ 2 p2 þ ;
V z ¼ V zð1Þ þ 2 V zð2Þ þ ;
p¼ p0 þ p1 þ 2 p2 þ
ð34Þ
Introducing the expansions (34) into the Eqs. (30)–(33), the following sets of differential equations are obtained O() equations o p1 ¼ 0; ox
c
oV ð1Þ o p1 z þ ¼ 0; on on
oV ð1Þ V ð1Þ oV ð1Þ r þ r þ z ¼ 0; ox x on
ð35Þ
with the boundary conditions V ð1Þ r jx¼kh ¼ c
ou1 ; on
p1 jx¼kh ¼ p1 :
ð36Þ
O(2) equations c
oV ð1Þ o p2 oV ð2Þ oV ð1Þ oV zð1Þ op2 op1 oV rð2Þ V rð2Þ oV zð2Þ oV zð2Þ r þ ¼ 0; c z þ V rð1Þ z þ V ð1Þ þ þ ¼ 0; þ þ þ ¼ 0; z on ox on ox on on os ox x on os ð37Þ
with the boundary conditions
oV ð1Þ ou2 ou1 ð1Þ ð2Þ r
V r þ ½u1 gðsÞ þ V j ; ¼ c ox x¼kh on on z x¼kh
o p1
p2 þ ½u1 gðsÞ ¼ p2 : ox
ð38Þ
x¼kh
In order to complete the equations, one must know the expressions of p1 and p2 in terms of the radial displacement u. For that purpose, we need the series expansion of the stretch ratios k1 and k2, which read k1 ’ kz ;
k2 ¼ kh þ ðu1 g þ kh GÞ þ 2 ðkh G2 gG þ Gu1 þ u2 Þ:
ð39Þ
Using the expansion (39) in the expression of p, given in (29), we have p1 ¼ b1 ðu1 gÞ þ c1 kh G; 2 2 mc o u1 p2 ¼ b0 þ b1 u2 þ b2 u21 þ b3 ðsÞu1 þ pðsÞ; kh kz on2
ð40Þ ð41Þ
where the function b3(s) and p(s) are defined by b3 ðsÞ ¼ 2½kh c2 G b2 g;
pðsÞ ¼ ðk2h c2 þ kh c1 ÞG2 þ b2 g2 2c2 kh Gg:
ð42Þ
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
627
Here the coefficients of c0, c1, c2, b0, b1 and b2 are defined by c0 ¼
1 oR ; kh kz okh
1 oR ; b0 ¼ kh okz
1 o2 R 1 o3 R ; c ¼ ; 2 kh kz ok2h 2kh kz ok3h c b b1 ¼ c1 0 ; b2 ¼ c2 1 : kh kh c1 ¼
ð43Þ
3.1. Solution of the field equations From the solution of the Eqs. (35) and (40) under the boundary conditions (36), we have u1 ¼ U ðn; sÞ;
p1 ¼ b1 ðU gÞ þ c1 kh G;
V ð1Þ z ¼
b1 U; c
V rð1Þ ¼
b1 oU x; 2c on
ð44Þ
provided that the following condition holds true b1 ¼ 2c2 =kh :
ð45Þ
Here U(n, s) is an unknown function whose governing equation will be obtained later and c is the phase velocity in the long wave approximation. To obtain the solution to O(2) equations, we introduce (44) into the Eqs. (37) and (41) and the boundary conditions (38) which results in b1 o2 U o p2 ¼ 0; xþ 2 on2 ox 2 oV ð2Þ b1 oU o p2 oU þ b1 g0 ðsÞ þ c1 kh G0 ðsÞ ¼ 0; c z þ þ b1 U on os on c on oV rð2Þ V ð2Þ oV ð2Þ b oU ¼ 0; þ r þ z þ 1 ox x on c os 2 2 mc oU p2 ¼ b0 þ b1 u2 þ b2 U 2 þ b3 ðsÞU þ pðsÞ; kh kz on2
ð46Þ ð47Þ ð48Þ ð49Þ
and the boundary conditions V rð2Þ jx¼kh ¼
3b1 oU b1 oU ou2 gðsÞ c U ; on 2c on 2c on
p2 jx¼kh ¼ p2 :
From the integration of (46) and the use of the boundary condition (50)2 we have 2 b o2 U 2 mc b1 k2h o2 U x þ b þ þ b1 u2 þ b2 U 2 þ b3 ðsÞu1 þ pðsÞ: p2 ¼ 1 0 4 on2 kh kz 4 on2 Introducing Eq. (48) and solution of (51) into Eq. (47), we obtain 2 oV rð2Þ V rð2Þ b1 o3 U 2 2b1 oU b1 2b2 oU mc b0 b1 k2h o3 U þ 3þ þ þ ¼ x þ U on kh kz ox x 4c on3 c os c c c 4c on3 b ou2 b3 ðsÞ oU b1 g0 ðsÞ c1 kh G0 ðsÞ þ þ þ 1 : c on c c c on From the integration of Eq. (52), one can get 2 b1 o3 U 3 x 2b1 oU b1 2b2 oU mc b0 b1 k2h o3 U ð2Þ þ 3þ þ x þ Vr ¼ U 2 c os on k h kz 16c on3 c c c 4c on3 0 0 b ou2 b3 ðsÞ oU b1 g ðsÞ c1 kh G ðsÞ þ þ þ 1 : c on c c c on
ð50Þ
ð51Þ
ð52Þ
ð53Þ
628
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
Applying boundary condition (50)1 in (53), we get 2 k3h b1 o3 U kh 2b1 oU b1 2b2 oU mc b0 b1 k2h o3 U þ 3þ þ þ U on kh kz 16c on3 2 c os c c c 4c on3 0 b ou2 b3 ðsÞ oU b1 g0 ðsÞ c1 kh G ðsÞ 3b oU b1 oU ou2 þ gðsÞ c þ 1 þ : ¼ 1U c on c c on 2c on c on 2c on
ð54Þ
Eliminating u2 in Eq. (54), we obtain the following forced Korteweg–de Vries (FKdV) equation with variable coefficients oU oU o3 U oU þ l1 U þ l2 3 þ l3 ðsÞ ¼ lðsÞ; os on on on
ð55Þ
where the coefficients l1, l2, l3(s) and l(s) are defined by 5 b m k2 b þ 2 ; l2 ¼ þ h 0 ; 2kh b1 4kz 16 2b1 kh c 2 b2 1 1 c kh l3 ðsÞ ¼ GðsÞ þ gðsÞ; lðsÞ ¼ g0 ðsÞ 1 G0 ðsÞ: 2 b1 b1 2kh 2b1 l1 ¼
ð56Þ
3.2. Progressive wave solution to FKdV equation In this sub-section, we shall present a progressive wave solution to the evolution equation given in (55). For that purpose, we introduce the following new dependent variable V as Z s 1 kh c 1 U ðn; sÞ ¼ V ðn; sÞ þ lðsÞ ds ¼ V ðn; sÞ þ gðsÞ GðsÞ : ð57Þ 2 b1 0 Introducing this expression of U into the Eq. (55), we have oV oV o3 V l1 kh c 1 oV þ l1 V þ l2 3 þ ¼ 0: gðsÞ GðsÞ þ l3 ðsÞ os on on 2 b on 1 Now, we introduce the following coordinate transformation Z s l1 kh c 1 0 0 s ¼ s; n ¼ n gðsÞ GðsÞ þ l3 ðsÞ ds: 2 b1 0
ð58Þ
ð59Þ
In the new coordinate system, the evolution equation reduces to the conventional Korteweg–de Vries equation oV oV o3 V þ l V þ l ¼ 0: 1 2 os0 on0 on03 As is well known, this evolution equation has the following solitary wave solution 1=2 l1 a la 2 V ¼ asech f; f ¼ n0 1 s 0 ; 12l2 3
ð60Þ
ð61Þ
where a is the amplitude of the solitary wave. The phase function f can be expressed as in terms of the coordinates (n,s) 1=2 Z s l1 a la l1 kh c 1 f¼ gðsÞ n 1 s GðsÞ þ l3 ðsÞ ds : ð62Þ 12l2 3 2 b1 0 Inserting Eqs. (61), (62) and (56)3 into (57), the solution of FKdV equation with variable coefficients (55) is given by
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
1 kh c 1 U ¼ asech f þ gðsÞ GðsÞ ; 2 b1 2
629
ð63Þ
where f¼
l1 a 12l2
1=2 n
l1 a s 3
Z s 0
l1 kh c 1 kh c 2 b 1 GðsÞ þ GðsÞ 2 þ gðsÞ ds : gðsÞ 2 b1 b1 b1 2kh
ð64Þ
From Eq. (44)2 and (63), the fluid pressure p1 is given by b1 kh c 1 gðsÞ þ GðsÞ: ð65Þ 2 2 As is seen from the expression of the phase function f, (64), the trajectory of the wave is not a straight line anymore, it is rather a curve in the (n, s) plane. This is the result of the stenosis in the tube. As a matter of fact, the existence of a stenosis causes the variable wave speed. Noting that s is the space variable and n is a temporal variable, the wave speed is defined by p1 ¼ ab1 sech2 f
vp ¼
ds ¼ dn l1 a þ 3
3 4kh
b2 2b1
1
: gðsÞ þ bkh1 c2 l12c1 GðsÞ
ð66Þ
As is seen from Eq. (66), due to the existence of a stenosis, the wave speed is variable along the tube axis. 4. Numerical results and discussions For the numerical calculation, we need the values of the coefficients c0, c1, c2, b0, b1, b2, l1, l2, l3(s) and l(s). In order to do that, one has to know the constitutive relation of the tube material. In this work, we shall utilize the constitutive relation proposed by Demiray [31] for soft biological tissues. Following Demiray [31], the strain energy density function may be expressed as R¼
1 fexp½aðI 1 3Þ 1g; 2a
ð67Þ
where a is a material constant and I1 is the first invariant of Finger deformation tensor defined by I1 = kz2 + kh2 + 1/(kz2kh2). Introducing (67) into Eq. (43), the coefficients c0, c1, c2, b0, b1 and b2 are obtained as 2 ! ! !2 3 1 1 1 4 3 1 c0 ¼ kh 3 2 F ðkh ; kz Þ; c1 ¼ 1 þ 4 2 þ 2a kh 3 2 5F ðkh ; kz Þ; kh kz kh kz kh kz kh kz kh kz 2 3 ! ! !3 1 4 12 1 3 1 2 ð68Þ 5 2 þ 6a kh 3 2 c2 ¼ 1 þ 4 2 þ 4a kh 3 2 5F ðkh ; kz Þ; 2kh kz kh kz kh kz kh kz kh kz ! 1 1 c b b0 ¼ kz 2 3 F ðkh ; kz Þ; b1 ¼ c1 0 ; b2 ¼ c2 1 ; kh kh kh kh kz where the function F(kh, kz) is defined by " !# 1 2 2 F ðkh ; kz Þ ¼ exp a kh þ kz þ 2 2 3 : kh k z
ð69Þ
Right now, we need the value of material constant a. For the static case, the present model was compared by Demiray [32] with the experimental measurements by Simon et al. [33] on canine abdominal artery with the characteristics Ri = 0.31 cm, R0 = 0.38 cm and kz = 1.53, and the value of the material constant a was found to be a = 1.948. Using this numerical value of the coefficient a, and for the initial deformation kh = kz = 1.6, we obtain c0 = 49.183, c1 = 326.844, c2 = 1176.561, b0 = 78.692, b1 = 296.105, b2 = 991.496, l1 = 4.911, and
630
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
l2 = 0.043 provided m = 0.1. In order to get the coefficient of l3(s) and l(s), we have to specify the function G(s) and g(s) which characterize the shape of the stenosis in the undeformed and deformed state. In this work, we consider G(s) = 0, g(s) = sech (ds), and a = 1. Utilizing these numerical values, we obtain the solution of FKdV equation with variable coefficients (55) as 1 U ¼ sech2 f þ sechðdsÞ; 2 where
ð70Þ
1:2055 1 tan ½sinhðdsÞ : f ¼ 3:085 n 1:637s þ d
ð71Þ
The fluid pressure function is as p1 ¼ 296:105sech2 f 148:053sechðdsÞ;
ð72Þ
and the wave speed is
a
1.6 Radial displacement, U( ξ,τ)
150
ξ=0 ξ=3 ξ=5
Fluid pressure function, p( ξ,τ)
1.8
1.4 1.2 1 0.8 0.6 0.4 –5
0
5
10
15
ξ=0 ξ=3 ξ=5
b
100 50 0 –50 –100 –150 –5
0
5
Space, τ
10
15
Space, τ
Fig. 1. Graph of (a) radial displacement U, (b) pressure function p1 versus space s for d = 0.01 at time n = 0, n = 3 and n = 5.
250
a
δ=0.1 δ=1
Fluid pressure function, p(ξ ,τ)
Radial displacement, U(ξ ,τ)
1.5
1
0.5
0 –5
0 Space, τ
5
b
δ=0.1 δ=1
200 150 100 50 0 –50 –100 –150 –6
–4
–2
0 Space, τ
2
4
6
Fig. 2. Graph of (a) radial displacement at time n = 0, (b) pressure function at time n = 1 versus space s for d = 0.1 and d = 1.
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
631
2.6 2.4
Wave speed, Vp
2.2 2 δ=0.01 δ=0.02 δ=0.05
1.8 1.6 1.4 1.2 1 0.8 –30
–20
–10
0
10
20
30
Space, τ Fig. 3. Wave speed, Vp versus space s for d = 0.01, d = 0.02, d = 0.05.
vp ¼
1 : 1:637 1:205sechðdsÞ
ð73Þ
Since the stenosis takes its maximum value (the radius is minimum) at the center, as is seen from the Eq. (73), the speed of wave reaches to its maximum value at the center and it becomes smaller and smaller as we go away from the center. Graph of radial displacement U and fluid pressure function p1 versus space s for d = 0.01 at time n = 0, n = 3 and n = 5 are shown in Fig. 1. It shows radial displacement and fluid pressure function admit a solitary wave solution and propagate with a constant amplitude to the right as time n increases. Fig. 2(a) shows the radial displacement with respect to space s at time n = 0. It is seen that when d is reduced from 1 to 0.1, the solitary wave solution becomes broader and shifted up. While in Fig. 2(b) shows the pressure function versus space s at time n = 1, we observe that the solitary wave becomes narrower as d is increased. Besides, it shows the fluid pressure function is minimum at the center of the stenosis, whereas the fluid velocity is maximum at the same point. Such a result is consistent with Bernoulli’s law. Fig. 3 shows wave speed Vp versus space variable s. The outmost curve is at d = 0.01, the middle is at d = 0.02 and the inmost is at d = 0.05. It shows the wave speed takes its maximum value at the center of the stenosis and it gets smaller and smaller as we go away from the center of the stenosis approaches a constant value of 0.637. Such a result seems to be reasonable from physical considerations. Note that if the shape of the stenosis is sharp, the wave speed decays rapidly. References [1] T.J. Pedley, Fluid Mechanics of Large Blood Vessels, Cambridge University Press, Cambridge, 1980. [2] Y.C. Fung, Biodynamics: Circulation, Springer-Verlag, New York, 1981. [3] H.B. Atabek, H.S. Lew, Wave propagation through a viscous incompressible fluid contained in an initially stressed elastic tube, Biophys. J. 7 (1966) 486–503. [4] A.J. Rachev, Effects of transmural pressure and muscular activity on pulse waves in arteries, J. Biomech. Eng., ASME 102 (1980) 119– 123. [5] H. Demiray, Wave propagation through a viscous fluid contained in a prestressed thin elastic tube, Int. J. Eng. Sci. 30 (1992) 1607– 1620. [6] G. Rudinger, Shock waves in a mathematical model of aorta, J. Appl. Mech. 37 (1970) 34–37. [7] S.C. Ling, H.B. Atabek, A nonlinear analysis of pulsatile blood flow in arteries, J. Fluid Mech. 55 (1972) 492–511. [8] M. Anliker, R.L. Rockwell, E. Ogden, Nonlinear analysis of flow pulses and shock waves in arteries, Z. Angew. Math. Phys. 22 (1968) 217–246. [9] R.J. Tait, T.B. Moodie, Waves in nonlinear fluid filled tubes, Wave Motion 6 (1984) 197–203.
632 [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632 R.S. Johnson, A nonlinear equation incorporating damping and dispersion, J. Fluid Mech. 42 (1970) 49–60. Y. Hashizume, Nonlinear pressure waves in a fluid-filled elastic tube, J. Phys. Soc. Jpn. 54 (1985) 3305–3312. Y. Hashizume, Nonlinear pressure waves arteries, J. Phys. Soc. Jpn. 57 (1988) 4160–4168. S. Yomosa, Solitary waves in large blood vessels, J. Phys. Soc. Jpn. 56 (1987) 506–520. H. Demiray, Solitary waves in a prestressed elastic tube, Bull. Math. Biol. 58 (1996) 939–955. H. Demiray, Solitary waves in initially stressed thin elastic tubes, Int. J. Non-linear Mech. 32 (1997) 1165–1176. H. Demiray, N. Antar, Nonlinear waves in an inviscid fluid contained in a prestressed viscoelastic thin tube, Z. Angew. Math. Phys. 48 (1997) 325–340. H. Demiray, Propagation of weakly nonlinear waves in fluid-filled thick viscoelastic tubes, Appl. Math. Modell. 23 (1999). N. Antar, H. Demiray, Weakly nonlinear waves in a prestressed thin elastic tube containing a viscous fluid, Int. J. Eng. Sci. 37 (1999) 1859–1876. N. Antar, H. Demiray, The boundary layer approximation and nonlinear waves in elastic tubes, Int. J. Eng. Sci. 38 (2000) 1441–1457. N. Antar, The Korteweg–de Vries-Burgers hierarchy in fluid-filled elastic tubes, Int. J. Eng. Sci. 40 (2002) 1179–1198. H. Demiray, Propagation of weakly nonlinear waves in fluid-filled thin elastic tubes, Appl. Math. Comp. 133 (2002). H. Demiray, Nonlinear waves in a prestressed elastic tube filled with a layered fluid, Int. J. Eng. Sci. 40 (2002) 713–726. H. Demiray, On the derivation of some non-linear evolution equations and their progressive wave solutions, Int. J. Non-linear Mech. 38 (2003) 63–70. I. Bakirtas, N. Antar, Evolution equaitons for nonlinear waves in a tapered elastic tube filled with a viscous fluid, Int. J. Eng. Sci. 41 (2003) 1163–1176. H. Demiray, On the existence of some evolution equations in fluid-filled elastic tubes and their progressive wave solutions, Int. J. Eng. Sci. 42 (2004) 1693–1706. H. Demiray, Weakly nonlinear waves in a lineary tapered elastic tube filled with a fluid, Math. Comp. Modell. 39 (2004) 151–162. I. Bakirtas, H. Demiray, Weakly non-linear waves in a tapered elastic tube filled with an inviscid fluid, Int. J. Non-linear Mech. 40 (2005) 785–793. H. Demiray, On some nonlinear waves in fluid-filled viscoelastic tubes : weakly dispersive case, Commun. Non-linear Sci. Num. Sim. 10 (2005) 425–440. C.S. Gardner, G.K. Morikawa, Similarity in the asymptotic behavior of collision-free hydromagnetic waves and water waves, Courant Institute Math. Sci. Report, NYO-9082, (1960) 1–30. T. Taniuti, C.C. Wei, Reductive perturbation method in non-linear wave propagation I, J. Phys. Soc. Jpn. 24 (1968) 941–946. H. Demiray, On the elasticity of soft biological tissues, J. Biomechan. 5 (1972) 309–311. H. Demiray, Large deformation analysis of some basic problems in biophysics, Bull. Math. Biol. 38 (1976) 701–711. B.R. Simon, A.S. Kobayashi, D.E. Stradness, C.A. Wiederhielm, Re-evaluation of arterial constitutive laws, Circ. Res. 30 (1972) 491–500.
Forced Korteweg–de Vries equation in an elastic tube filled with an inviscid fluid Tay Kim Gaik
*
Panel of Mathematics, Sciences Studies Center, University College Technology of Tun Hussein Onn, 86400 Parit Raja, Batu Pahat, Johor, Malaysia Received 19 January 2006; accepted 13 April 2006 Available online 17 July 2006
Abstract In the present work, treating the arteries as a prestressed thin walled elastic tube with a stenosis and the blood as an inviscid fluid, we have studied the propagation of weakly nonlinear waves in such a composite medium, in the long wave approximation, by use of the reductive perturbation method [C.S. Gardner, G.K. Morikawa, Similarity in the asymptotic behavior of collision-free hydromagnetic waves and water waves, Courant Institute Math. Sci. Report, NYO-9082 (1960) 1–30, T. Taniuti, C.C. Wei, Reductive perturbation method in non-linear wave propagation I, J. Phys. Soc. Jpn., 24 (1968) 941–946]. We obtained the forced Korteweg–de Vries (FKdV) equation with variable coefficients as the evolution equation. By use of the coordinate transformation, it is shown that this type of evolution equation admits a progressive wave solution with variable wave speed. As might be expected from physical consideration, the wave speed reaches its maximum value at the center of stenosis and gets smaller and smaller as we go away from the center of the stenosis. The variations of radial displacement and the fluid pressure with the distance parameter are also examined numerically. The results seem to be consistent with Bernoulli’s law for inviscid fluid. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Solitary waves; Elastic tubes with stenosis
1. Introduction Due to its applications in arterial mechanics, the propagation of pressure pulses in fluid-filled distensible tubes has been studied by several researchers [1,2]. Most of the works on wave propagation in compliant tubes have considered small amplitude waves ignoring the nonlinear effects and focused on the dispersive character of waves (see [3–5]). However, when the nonlinear terms arising from the constitutive equations and kinematical relations are introduced, one has to consider either finite amplitude, or small-but-finite amplitude waves, depending on the order of nonlinearity.
*
Tel.: +07 4536721/4536701; fax: +07 4536051. E-mail address: [email protected]
0020-7225/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2006.04.008
622
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
The propagation of finite amplitude waves in fluid-filled elastic or viscoelastic tubes has been examined, for instance, by Rudinger [6], Ling and Atabek [7], Anliker et al. [8] and Tait and Moodie [9] by using the method of characteristics, in studying the shock formation. On the other hand, the propagation of smallbut-finite amplitude waves in distensible tubes has been investigated by Johnson [10], Hashizume [11,12], and Yomosa [13]. In all these works [10–13], the effect of initial deformation is neglected. Recently in a series of works of Demiray, Antar and Bakirtas (see [14–28]) in which they treated artery as incompressible prestressed thin isotropic elastic, thick viscoelastic or tapered elastic tube filled with inviscid, viscous or layered fluid as blood, using approximate method on fluid equations and reductive perturbation method in the long-wave approximation, they obtained various evolution equations of Korteweg–de Vries, Burgers and Korteweg–de Vries–Burgers type equations. In all previous works, they treated the arteries as circularly cylindrical long thin tubes with a constant cross-section. However due to decomposition of fat or cholesterol in artery over time, the artery become narrower and may have variable radius along the axis of the tube. Thus, in this work, treating the arteries as an incompressible prestressed thin walled elastic tube with a stenosis and the blood as an incompressible inviscid fluid, we have studied the propagation of weakly nonlinear waves in such a composite medium, in the long wave approximation, by use of the reductive perturbation method [29,30]. We obtained the forced Korteweg–de Vries (FKdV) equation with variable coefficients as the evolution equation. By use of the coordinate transformation, it is shown that this type of evolution equation admits a progressive wave solution with variable wave speed. As might be expected from physical consideration, the wave speed reaches its maximum value at the center of stenosis and gets smaller and smaller as one goes away from the center of the stenosis. The variations of radial displacement and the fluid pressure with the distance parameter are also examined numerically. The results seem to be consistent with Bernoulli’s law for inviscid fluid. 2. Basic equations and theoretical preliminaries In this section, we shall give the derivation of the field equations of an elastic tube, which is considered to be a model for an artery, and an inviscid fluid, which is considered to be a model for blood. 2.1. Equations of tube In this sub-section, we shall derive the governing equations of an elastic tube filled with an inviscid fluid. Such a combination of a solid and a fluid is considered to be a model for blood flow in arteries. For a healthy human being, the systolic pressure is about 120 mm Hg, and the diastolic pressure is around 80 mm Hg. This means that the arteries are subjected to a mean pressure P0 = 100 mm Hg, and in the course of blood flow, a dynamical pressure increment DP = ±20 mm Hg is added on this initial field. Moreover, experimental studies [2] revealed that the arteries are also subjected to an initial axial stretch kz, which is about kz = 1.6. These observations show that the arteries are initially subjected to static deformation both in the radial and the axial directions, and a dynamical pressure (or a radial displacement u*) is superimposed on this initial deformation. Due to the external tethering in the axial direction, the effect of axial displacement is neglected. Now, we consider a thin and long tube of circular cross-section with radius R*(Z*) in the cylindrical polar coordinates (R*, H, Z*). Then, the position vector of a point on the tube may be described by R ¼ R ðZ Þer þ Z ez ;
ð1Þ
where er, eh and ez are the unit base vectors in the cylindrical polar coordinates and Z* is the axial coordinates of a material point in the natural state. The arclengths along the meridional and circumferential curves are given by " 2 #1=2 dR dS Z ¼ 1 þ dZ ; dS H ¼ R ðZ Þ dH: ð2Þ dZ
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
623
Motivated with the experimental observations [2], we shall assume that the elastic tube is subjected to an axial stretch ratio kz, and the static pressure P 0 ðZ Þ. Then, the deformation may be described by r0 ¼ r ðz Þer þ z ez ;
z ¼ kz Z ;
ð3Þ
where z* is the axial coordinate after static deformation, and r*(z*) characterizes the variable radius after this static deformation. Then the arclengths after static deformation along the meridional and circumferential directions are given by 0
2 1=2
ds0z ¼ ½1 þ ðr Þ
dz ;
dsh ¼ r ðz Þ dh;
where a prime denotes the differentiation of the corresponding field variable with respect to Then, the stretch ratios along the meridional and circumferential curves are given by 0
k01 ¼
ds0z ½1 þ ðr Þ2 1=2 ¼ kz ; 0 2 1=2 dS Z ½1 þ k2z ðR Þ
k02 ¼
ð4Þ z *.
ds0h r ðz Þ ¼ : dS H R ðz Þ
ð5Þ
Upon this initial static deformation, we shall superimpose a finite dynamical radial displacement u*(z*, t*), where t* is the time parameter, but, in view of the external tethering in the axial direction, the axial displacement is assumed to be negligible. Then, the position vector r of a generic point on the tube may be described by r ¼ ½r ðz Þ þ u ðz ; t Þer þ z ez : The arc lengths along the deformed meridional and circumferential curves are respectively given by " #1=2 2 ou 0 dsz ¼ 1 þ r þ dz ; dsh ¼ ½r ðz Þ þ u ðz ; t Þ dh: oz
ð6Þ
ð7Þ
Then, the stretch ratios along the meridional and circumferential curves in the final configuration read, respectively, by h i1=2 0 2 1 þ r þ ou oz r ðz Þ þ u ðz ; t Þ k1 ¼ k z : ð8Þ ; k ¼ 2 0 2 1=2 R ðz Þ ½1 þ k2z ðR Þ The unit tangent vector t along the deformed meridional curve and the unit exterior normal vector n to the deformed membrane are given by 0 ou 0 r þ oz er þ ez er r þ ou ez oz t¼ ; n¼ ; ð9Þ K K where the function K is defined by " 2 #1=2 ou 0 K¼ 1þ r þ dz : oz The material that we shall consider is assumed to be incompressible. This condition imposes the following restriction on the thickness H, and h, before and after deformation respectively h¼
H : k1 k2
ð10Þ
Let T1 and T2 be the membrane forces along the meridional and circumferential curves, respectively. Then, the equation of the radial motion of a small tube element placed between the planes z* = const, z* + dz* = const, h = const and h + dh = const may be given by 2 o T1 ou H 2 0 0 2 1=2 o u ðr þ u Þ r þ K þ ðr þ u ÞP ¼ q R ðz Þ½1 þ k ðR Þ ; T 2 0 z oz K kz oz ot2
ð11Þ
624
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
where q0 is the mass density of the membrane material, and P* is the fluid pressure on the inner surface of the tube. Let lR be the strain energy density function of the tube material, where l is the shear modulus. Then, the membrane forces T1 and T2 may be expressed in terms of the stretch ratios as T1 ¼
lH oR ; k2 ok1
T2 ¼
lH oR : k1 ok2
ð12Þ
Introducing Eq. (12) into Eq. (11), the equation of motion of the tube in the radial direction takes the following form o lHR 0 ou oR lH 0 2 1=2 oR r þ ½1 þ k2z ðR Þ oz kz ok2 K oz ok1 þ ðr þ u ÞP q0
2 H 0 2 1=2 o u R ðz Þ½1 þ k2z ðR Þ ¼ 0: kz ot2
ð13Þ
2.2. Equations of fluid In general, blood is known to be an incompressible non-Newtonian fluid. However, in the course of flow in large arteries, the red blood cells in the vicinity of arterial wall move to the central region of the artery so that hematocrit ratio becomes quite low near the arterial wall, which results in lower viscosity in this region. Moreover, due to high shear rate near the arterial wall the viscosity of blood is further reduced. Therefore, for flow problems in large blood vessels, as a first approximation, blood may be treated as an inviscid fluid. The axially symmetrical equation of motion of an inviscid fluid in the cylindrical polar coordinates (r, h, z*) may be given by oV r oV r oV 1 oP ¼ 0; þ V z r þ þ V r qf or ot or oz
ð14Þ
oV z oV z oV 1 oP þ V z z þ þ V r ¼ 0; qf oz ot or oz
ð15Þ
oV r V r oV z þ þ ¼0 or r oz
ð16Þ
ðincompressibilityÞ;
with the boundary conditions V
r jr¼rf
ou ou 0 ¼ þ r þ V z jr¼rf ; ot oz
P jr¼rf ¼ P ;
ð17Þ
where V r and V z are the components of the fluid velocity vector in the radial and axial direction, respectively, qf is the mass density of the fluid, P is the pressure function and rf = r* + u* is the final inner radius of the tube. At this stage it is convenient to introduce the following non-dimensional quantities R0 t ¼ t; z ¼ R0 z; u ¼ R0 u; V r ¼ c0 V r ; V z ¼ c0 V z ; r ¼ R0 x; P ¼ qf c20 p c0 ð18Þ q0 H lH 2 2 P ¼ qf c0 ; c0 ¼ ; R ðz Þ ¼ R0 ½1 F ðzÞ; r ðz Þ ¼ R0 ½kh f ðzÞ; p; m ¼ qf R 0 qf R0 where kh = r0/R0 is the initial stretch ratio, r0 is the radius of the origin after finite static deformation, R0 is the initial reference radius and c0 is the Moens–Korteweg wave speed.
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
625
Introducing (18) into Eqs. (13)–(17), the following non-dimensional equations are obtained p¼
m½1 F ðzÞ½1 þ k2z ðF 0 Þ2 1=2 o2 u ½1 þ k2z ðF 0 Þ2 1=2 oR þ ot2 kz ðkh f ðzÞ þ uÞ ok2 kz ðkh f ðzÞ þ uÞ ( ) 1 o ½1 F ðzÞðf 0 þ ou=ozÞ oR ; ðkh f ðzÞ þ uÞ oz ½1 þ ðf 0 þ ou=ozÞ2 1=2 ok1
ð19Þ
oV r oV r oV r o p ¼ 0; þVr þVz þ ox ot ox oz oV z oV z oV z o p ¼ 0; þVr þVz þ oz ot ox oz oV r V r oV z þ þ ¼ 0; ox x oz with the boundary conditions ou ou þ f 0 þ V r jx¼kh f ðzÞþu ¼ V z jx¼kh f ðzÞþu ; ot oz
ð20Þ ð21Þ ð22Þ
pjx¼kh f ðzÞþu ¼ p:
ð23Þ
The Eqs. (19)–(23) give sufficient relations to determine the field quantities u, Vr, Vz and p completely. 3. Long wave approximation In this section, we shall examine the propagation of small-but-finite amplitude waves in a fluid-filled thin elastic tube whose non-dimensional governing equations are given in Eqs. (19)–(23). For this, we adopt the long wave approximation and employ the reductive perturbation method [29,30]. For this type of problems, it is convenient to introduce the following type of stretched coordinates n ¼ 1=2 ðz ctÞ;
s ¼ 3=2 z;
ð24Þ
where is a small parameter measuring the weakness of nonlinearity and dispersion, whereas c is the scale parameter to be determined from the solution. Solving z in terms of s we get z ¼ 3=2 s:
ð25Þ
Introducing (25) into the expression of F(z) and f(z), we obtain ^ sÞ; F ð3=2 sÞ ¼ Gðn;
f ð3=2 sÞ ¼ g^ðn; sÞ:
In order to take the effect of stenosis into account, F(z) and f(z) must be of order of ^ sÞ and g^ðn; sÞ have the following form work, we shall assume that Gðn; ^ sÞ ¼ GðsÞ; Gðn;
ð26Þ 5/2
. For the present
g^ðn; sÞ ¼ gðsÞ:
Introducing the following differential relations o o o o 1=2 o 1=2 ¼ c ; ¼ þ ot on oz on os
ð27Þ
ð28Þ
into the Eqs. (19)–(23), we obtain 2 1=2 2 1=2 mc2 ½1 GðsÞ½1 þ 5 k2z ðdG Þ o2 u ½1 þ 5 k2z ðdG Þ oR o o ds ds þ þ os kz ðkh g þ uÞ kz ðkh g þ uÞ ok2 ðkh g þ uÞ on on2 ( ) ½1 GðsÞ dg ou ou oR þ þ 2 ; 1=2 ou þ 3=2 ou Þ2 1=2 ds on os ok1 ½1 þ ð5=2 dg þ ds on os oV r oV r oV r op 1=2 oV r 1=2 ¼ 0; þVr þ Vz þ c þ ox on ox on os p¼
ð29Þ ð30Þ
626
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
oV z oV z oV z oV z p op 1=2 1=2 o þ ¼ 0; þVr þ Vz þ þ c on os on ox on os oV r V r oV z oV z þ þ 1=2 þ ¼ 0; ox x on os 1=2
ð31Þ ð32Þ
with the boundary conditions V r jx¼kh gþu ¼ 1=2 c
ou dg ou ou þ 1=2 2 þ þ V z jx¼kh gþu ; on ds on os
pjx¼kh gþu ¼ p:
ð33Þ
For the long wave limit, it is assumed that the field quantities may be expanded into asymptotic series of as V r ¼ 1=2 ðV rð1Þ þ 2 V rð2Þ þ Þ;
u ¼ u1 þ 2 u2 þ ;
p ¼ p0 þ p1 þ 2 p2 þ ;
V z ¼ V zð1Þ þ 2 V zð2Þ þ ;
p¼ p0 þ p1 þ 2 p2 þ
ð34Þ
Introducing the expansions (34) into the Eqs. (30)–(33), the following sets of differential equations are obtained O() equations o p1 ¼ 0; ox
c
oV ð1Þ o p1 z þ ¼ 0; on on
oV ð1Þ V ð1Þ oV ð1Þ r þ r þ z ¼ 0; ox x on
ð35Þ
with the boundary conditions V ð1Þ r jx¼kh ¼ c
ou1 ; on
p1 jx¼kh ¼ p1 :
ð36Þ
O(2) equations c
oV ð1Þ o p2 oV ð2Þ oV ð1Þ oV zð1Þ op2 op1 oV rð2Þ V rð2Þ oV zð2Þ oV zð2Þ r þ ¼ 0; c z þ V rð1Þ z þ V ð1Þ þ þ ¼ 0; þ þ þ ¼ 0; z on ox on ox on on os ox x on os ð37Þ
with the boundary conditions
oV ð1Þ ou2 ou1 ð1Þ ð2Þ r
V r þ ½u1 gðsÞ þ V j ; ¼ c ox x¼kh on on z x¼kh
o p1
p2 þ ½u1 gðsÞ ¼ p2 : ox
ð38Þ
x¼kh
In order to complete the equations, one must know the expressions of p1 and p2 in terms of the radial displacement u. For that purpose, we need the series expansion of the stretch ratios k1 and k2, which read k1 ’ kz ;
k2 ¼ kh þ ðu1 g þ kh GÞ þ 2 ðkh G2 gG þ Gu1 þ u2 Þ:
ð39Þ
Using the expansion (39) in the expression of p, given in (29), we have p1 ¼ b1 ðu1 gÞ þ c1 kh G; 2 2 mc o u1 p2 ¼ b0 þ b1 u2 þ b2 u21 þ b3 ðsÞu1 þ pðsÞ; kh kz on2
ð40Þ ð41Þ
where the function b3(s) and p(s) are defined by b3 ðsÞ ¼ 2½kh c2 G b2 g;
pðsÞ ¼ ðk2h c2 þ kh c1 ÞG2 þ b2 g2 2c2 kh Gg:
ð42Þ
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
627
Here the coefficients of c0, c1, c2, b0, b1 and b2 are defined by c0 ¼
1 oR ; kh kz okh
1 oR ; b0 ¼ kh okz
1 o2 R 1 o3 R ; c ¼ ; 2 kh kz ok2h 2kh kz ok3h c b b1 ¼ c1 0 ; b2 ¼ c2 1 : kh kh c1 ¼
ð43Þ
3.1. Solution of the field equations From the solution of the Eqs. (35) and (40) under the boundary conditions (36), we have u1 ¼ U ðn; sÞ;
p1 ¼ b1 ðU gÞ þ c1 kh G;
V ð1Þ z ¼
b1 U; c
V rð1Þ ¼
b1 oU x; 2c on
ð44Þ
provided that the following condition holds true b1 ¼ 2c2 =kh :
ð45Þ
Here U(n, s) is an unknown function whose governing equation will be obtained later and c is the phase velocity in the long wave approximation. To obtain the solution to O(2) equations, we introduce (44) into the Eqs. (37) and (41) and the boundary conditions (38) which results in b1 o2 U o p2 ¼ 0; xþ 2 on2 ox 2 oV ð2Þ b1 oU o p2 oU þ b1 g0 ðsÞ þ c1 kh G0 ðsÞ ¼ 0; c z þ þ b1 U on os on c on oV rð2Þ V ð2Þ oV ð2Þ b oU ¼ 0; þ r þ z þ 1 ox x on c os 2 2 mc oU p2 ¼ b0 þ b1 u2 þ b2 U 2 þ b3 ðsÞU þ pðsÞ; kh kz on2
ð46Þ ð47Þ ð48Þ ð49Þ
and the boundary conditions V rð2Þ jx¼kh ¼
3b1 oU b1 oU ou2 gðsÞ c U ; on 2c on 2c on
p2 jx¼kh ¼ p2 :
From the integration of (46) and the use of the boundary condition (50)2 we have 2 b o2 U 2 mc b1 k2h o2 U x þ b þ þ b1 u2 þ b2 U 2 þ b3 ðsÞu1 þ pðsÞ: p2 ¼ 1 0 4 on2 kh kz 4 on2 Introducing Eq. (48) and solution of (51) into Eq. (47), we obtain 2 oV rð2Þ V rð2Þ b1 o3 U 2 2b1 oU b1 2b2 oU mc b0 b1 k2h o3 U þ 3þ þ þ ¼ x þ U on kh kz ox x 4c on3 c os c c c 4c on3 b ou2 b3 ðsÞ oU b1 g0 ðsÞ c1 kh G0 ðsÞ þ þ þ 1 : c on c c c on From the integration of Eq. (52), one can get 2 b1 o3 U 3 x 2b1 oU b1 2b2 oU mc b0 b1 k2h o3 U ð2Þ þ 3þ þ x þ Vr ¼ U 2 c os on k h kz 16c on3 c c c 4c on3 0 0 b ou2 b3 ðsÞ oU b1 g ðsÞ c1 kh G ðsÞ þ þ þ 1 : c on c c c on
ð50Þ
ð51Þ
ð52Þ
ð53Þ
628
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
Applying boundary condition (50)1 in (53), we get 2 k3h b1 o3 U kh 2b1 oU b1 2b2 oU mc b0 b1 k2h o3 U þ 3þ þ þ U on kh kz 16c on3 2 c os c c c 4c on3 0 b ou2 b3 ðsÞ oU b1 g0 ðsÞ c1 kh G ðsÞ 3b oU b1 oU ou2 þ gðsÞ c þ 1 þ : ¼ 1U c on c c on 2c on c on 2c on
ð54Þ
Eliminating u2 in Eq. (54), we obtain the following forced Korteweg–de Vries (FKdV) equation with variable coefficients oU oU o3 U oU þ l1 U þ l2 3 þ l3 ðsÞ ¼ lðsÞ; os on on on
ð55Þ
where the coefficients l1, l2, l3(s) and l(s) are defined by 5 b m k2 b þ 2 ; l2 ¼ þ h 0 ; 2kh b1 4kz 16 2b1 kh c 2 b2 1 1 c kh l3 ðsÞ ¼ GðsÞ þ gðsÞ; lðsÞ ¼ g0 ðsÞ 1 G0 ðsÞ: 2 b1 b1 2kh 2b1 l1 ¼
ð56Þ
3.2. Progressive wave solution to FKdV equation In this sub-section, we shall present a progressive wave solution to the evolution equation given in (55). For that purpose, we introduce the following new dependent variable V as Z s 1 kh c 1 U ðn; sÞ ¼ V ðn; sÞ þ lðsÞ ds ¼ V ðn; sÞ þ gðsÞ GðsÞ : ð57Þ 2 b1 0 Introducing this expression of U into the Eq. (55), we have oV oV o3 V l1 kh c 1 oV þ l1 V þ l2 3 þ ¼ 0: gðsÞ GðsÞ þ l3 ðsÞ os on on 2 b on 1 Now, we introduce the following coordinate transformation Z s l1 kh c 1 0 0 s ¼ s; n ¼ n gðsÞ GðsÞ þ l3 ðsÞ ds: 2 b1 0
ð58Þ
ð59Þ
In the new coordinate system, the evolution equation reduces to the conventional Korteweg–de Vries equation oV oV o3 V þ l V þ l ¼ 0: 1 2 os0 on0 on03 As is well known, this evolution equation has the following solitary wave solution 1=2 l1 a la 2 V ¼ asech f; f ¼ n0 1 s 0 ; 12l2 3
ð60Þ
ð61Þ
where a is the amplitude of the solitary wave. The phase function f can be expressed as in terms of the coordinates (n,s) 1=2 Z s l1 a la l1 kh c 1 f¼ gðsÞ n 1 s GðsÞ þ l3 ðsÞ ds : ð62Þ 12l2 3 2 b1 0 Inserting Eqs. (61), (62) and (56)3 into (57), the solution of FKdV equation with variable coefficients (55) is given by
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
1 kh c 1 U ¼ asech f þ gðsÞ GðsÞ ; 2 b1 2
629
ð63Þ
where f¼
l1 a 12l2
1=2 n
l1 a s 3
Z s 0
l1 kh c 1 kh c 2 b 1 GðsÞ þ GðsÞ 2 þ gðsÞ ds : gðsÞ 2 b1 b1 b1 2kh
ð64Þ
From Eq. (44)2 and (63), the fluid pressure p1 is given by b1 kh c 1 gðsÞ þ GðsÞ: ð65Þ 2 2 As is seen from the expression of the phase function f, (64), the trajectory of the wave is not a straight line anymore, it is rather a curve in the (n, s) plane. This is the result of the stenosis in the tube. As a matter of fact, the existence of a stenosis causes the variable wave speed. Noting that s is the space variable and n is a temporal variable, the wave speed is defined by p1 ¼ ab1 sech2 f
vp ¼
ds ¼ dn l1 a þ 3
3 4kh
b2 2b1
1
: gðsÞ þ bkh1 c2 l12c1 GðsÞ
ð66Þ
As is seen from Eq. (66), due to the existence of a stenosis, the wave speed is variable along the tube axis. 4. Numerical results and discussions For the numerical calculation, we need the values of the coefficients c0, c1, c2, b0, b1, b2, l1, l2, l3(s) and l(s). In order to do that, one has to know the constitutive relation of the tube material. In this work, we shall utilize the constitutive relation proposed by Demiray [31] for soft biological tissues. Following Demiray [31], the strain energy density function may be expressed as R¼
1 fexp½aðI 1 3Þ 1g; 2a
ð67Þ
where a is a material constant and I1 is the first invariant of Finger deformation tensor defined by I1 = kz2 + kh2 + 1/(kz2kh2). Introducing (67) into Eq. (43), the coefficients c0, c1, c2, b0, b1 and b2 are obtained as 2 ! ! !2 3 1 1 1 4 3 1 c0 ¼ kh 3 2 F ðkh ; kz Þ; c1 ¼ 1 þ 4 2 þ 2a kh 3 2 5F ðkh ; kz Þ; kh kz kh kz kh kz kh kz kh kz 2 3 ! ! !3 1 4 12 1 3 1 2 ð68Þ 5 2 þ 6a kh 3 2 c2 ¼ 1 þ 4 2 þ 4a kh 3 2 5F ðkh ; kz Þ; 2kh kz kh kz kh kz kh kz kh kz ! 1 1 c b b0 ¼ kz 2 3 F ðkh ; kz Þ; b1 ¼ c1 0 ; b2 ¼ c2 1 ; kh kh kh kh kz where the function F(kh, kz) is defined by " !# 1 2 2 F ðkh ; kz Þ ¼ exp a kh þ kz þ 2 2 3 : kh k z
ð69Þ
Right now, we need the value of material constant a. For the static case, the present model was compared by Demiray [32] with the experimental measurements by Simon et al. [33] on canine abdominal artery with the characteristics Ri = 0.31 cm, R0 = 0.38 cm and kz = 1.53, and the value of the material constant a was found to be a = 1.948. Using this numerical value of the coefficient a, and for the initial deformation kh = kz = 1.6, we obtain c0 = 49.183, c1 = 326.844, c2 = 1176.561, b0 = 78.692, b1 = 296.105, b2 = 991.496, l1 = 4.911, and
630
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
l2 = 0.043 provided m = 0.1. In order to get the coefficient of l3(s) and l(s), we have to specify the function G(s) and g(s) which characterize the shape of the stenosis in the undeformed and deformed state. In this work, we consider G(s) = 0, g(s) = sech (ds), and a = 1. Utilizing these numerical values, we obtain the solution of FKdV equation with variable coefficients (55) as 1 U ¼ sech2 f þ sechðdsÞ; 2 where
ð70Þ
1:2055 1 tan ½sinhðdsÞ : f ¼ 3:085 n 1:637s þ d
ð71Þ
The fluid pressure function is as p1 ¼ 296:105sech2 f 148:053sechðdsÞ;
ð72Þ
and the wave speed is
a
1.6 Radial displacement, U( ξ,τ)
150
ξ=0 ξ=3 ξ=5
Fluid pressure function, p( ξ,τ)
1.8
1.4 1.2 1 0.8 0.6 0.4 –5
0
5
10
15
ξ=0 ξ=3 ξ=5
b
100 50 0 –50 –100 –150 –5
0
5
Space, τ
10
15
Space, τ
Fig. 1. Graph of (a) radial displacement U, (b) pressure function p1 versus space s for d = 0.01 at time n = 0, n = 3 and n = 5.
250
a
δ=0.1 δ=1
Fluid pressure function, p(ξ ,τ)
Radial displacement, U(ξ ,τ)
1.5
1
0.5
0 –5
0 Space, τ
5
b
δ=0.1 δ=1
200 150 100 50 0 –50 –100 –150 –6
–4
–2
0 Space, τ
2
4
6
Fig. 2. Graph of (a) radial displacement at time n = 0, (b) pressure function at time n = 1 versus space s for d = 0.1 and d = 1.
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632
631
2.6 2.4
Wave speed, Vp
2.2 2 δ=0.01 δ=0.02 δ=0.05
1.8 1.6 1.4 1.2 1 0.8 –30
–20
–10
0
10
20
30
Space, τ Fig. 3. Wave speed, Vp versus space s for d = 0.01, d = 0.02, d = 0.05.
vp ¼
1 : 1:637 1:205sechðdsÞ
ð73Þ
Since the stenosis takes its maximum value (the radius is minimum) at the center, as is seen from the Eq. (73), the speed of wave reaches to its maximum value at the center and it becomes smaller and smaller as we go away from the center. Graph of radial displacement U and fluid pressure function p1 versus space s for d = 0.01 at time n = 0, n = 3 and n = 5 are shown in Fig. 1. It shows radial displacement and fluid pressure function admit a solitary wave solution and propagate with a constant amplitude to the right as time n increases. Fig. 2(a) shows the radial displacement with respect to space s at time n = 0. It is seen that when d is reduced from 1 to 0.1, the solitary wave solution becomes broader and shifted up. While in Fig. 2(b) shows the pressure function versus space s at time n = 1, we observe that the solitary wave becomes narrower as d is increased. Besides, it shows the fluid pressure function is minimum at the center of the stenosis, whereas the fluid velocity is maximum at the same point. Such a result is consistent with Bernoulli’s law. Fig. 3 shows wave speed Vp versus space variable s. The outmost curve is at d = 0.01, the middle is at d = 0.02 and the inmost is at d = 0.05. It shows the wave speed takes its maximum value at the center of the stenosis and it gets smaller and smaller as we go away from the center of the stenosis approaches a constant value of 0.637. Such a result seems to be reasonable from physical considerations. Note that if the shape of the stenosis is sharp, the wave speed decays rapidly. References [1] T.J. Pedley, Fluid Mechanics of Large Blood Vessels, Cambridge University Press, Cambridge, 1980. [2] Y.C. Fung, Biodynamics: Circulation, Springer-Verlag, New York, 1981. [3] H.B. Atabek, H.S. Lew, Wave propagation through a viscous incompressible fluid contained in an initially stressed elastic tube, Biophys. J. 7 (1966) 486–503. [4] A.J. Rachev, Effects of transmural pressure and muscular activity on pulse waves in arteries, J. Biomech. Eng., ASME 102 (1980) 119– 123. [5] H. Demiray, Wave propagation through a viscous fluid contained in a prestressed thin elastic tube, Int. J. Eng. Sci. 30 (1992) 1607– 1620. [6] G. Rudinger, Shock waves in a mathematical model of aorta, J. Appl. Mech. 37 (1970) 34–37. [7] S.C. Ling, H.B. Atabek, A nonlinear analysis of pulsatile blood flow in arteries, J. Fluid Mech. 55 (1972) 492–511. [8] M. Anliker, R.L. Rockwell, E. Ogden, Nonlinear analysis of flow pulses and shock waves in arteries, Z. Angew. Math. Phys. 22 (1968) 217–246. [9] R.J. Tait, T.B. Moodie, Waves in nonlinear fluid filled tubes, Wave Motion 6 (1984) 197–203.
632 [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]
K.G. Tay / International Journal of Engineering Science 44 (2006) 621–632 R.S. Johnson, A nonlinear equation incorporating damping and dispersion, J. Fluid Mech. 42 (1970) 49–60. Y. Hashizume, Nonlinear pressure waves in a fluid-filled elastic tube, J. Phys. Soc. Jpn. 54 (1985) 3305–3312. Y. Hashizume, Nonlinear pressure waves arteries, J. Phys. Soc. Jpn. 57 (1988) 4160–4168. S. Yomosa, Solitary waves in large blood vessels, J. Phys. Soc. Jpn. 56 (1987) 506–520. H. Demiray, Solitary waves in a prestressed elastic tube, Bull. Math. Biol. 58 (1996) 939–955. H. Demiray, Solitary waves in initially stressed thin elastic tubes, Int. J. Non-linear Mech. 32 (1997) 1165–1176. H. Demiray, N. Antar, Nonlinear waves in an inviscid fluid contained in a prestressed viscoelastic thin tube, Z. Angew. Math. Phys. 48 (1997) 325–340. H. Demiray, Propagation of weakly nonlinear waves in fluid-filled thick viscoelastic tubes, Appl. Math. Modell. 23 (1999). N. Antar, H. Demiray, Weakly nonlinear waves in a prestressed thin elastic tube containing a viscous fluid, Int. J. Eng. Sci. 37 (1999) 1859–1876. N. Antar, H. Demiray, The boundary layer approximation and nonlinear waves in elastic tubes, Int. J. Eng. Sci. 38 (2000) 1441–1457. N. Antar, The Korteweg–de Vries-Burgers hierarchy in fluid-filled elastic tubes, Int. J. Eng. Sci. 40 (2002) 1179–1198. H. Demiray, Propagation of weakly nonlinear waves in fluid-filled thin elastic tubes, Appl. Math. Comp. 133 (2002). H. Demiray, Nonlinear waves in a prestressed elastic tube filled with a layered fluid, Int. J. Eng. Sci. 40 (2002) 713–726. H. Demiray, On the derivation of some non-linear evolution equations and their progressive wave solutions, Int. J. Non-linear Mech. 38 (2003) 63–70. I. Bakirtas, N. Antar, Evolution equaitons for nonlinear waves in a tapered elastic tube filled with a viscous fluid, Int. J. Eng. Sci. 41 (2003) 1163–1176. H. Demiray, On the existence of some evolution equations in fluid-filled elastic tubes and their progressive wave solutions, Int. J. Eng. Sci. 42 (2004) 1693–1706. H. Demiray, Weakly nonlinear waves in a lineary tapered elastic tube filled with a fluid, Math. Comp. Modell. 39 (2004) 151–162. I. Bakirtas, H. Demiray, Weakly non-linear waves in a tapered elastic tube filled with an inviscid fluid, Int. J. Non-linear Mech. 40 (2005) 785–793. H. Demiray, On some nonlinear waves in fluid-filled viscoelastic tubes : weakly dispersive case, Commun. Non-linear Sci. Num. Sim. 10 (2005) 425–440. C.S. Gardner, G.K. Morikawa, Similarity in the asymptotic behavior of collision-free hydromagnetic waves and water waves, Courant Institute Math. Sci. Report, NYO-9082, (1960) 1–30. T. Taniuti, C.C. Wei, Reductive perturbation method in non-linear wave propagation I, J. Phys. Soc. Jpn. 24 (1968) 941–946. H. Demiray, On the elasticity of soft biological tissues, J. Biomechan. 5 (1972) 309–311. H. Demiray, Large deformation analysis of some basic problems in biophysics, Bull. Math. Biol. 38 (1976) 701–711. B.R. Simon, A.S. Kobayashi, D.E. Stradness, C.A. Wiederhielm, Re-evaluation of arterial constitutive laws, Circ. Res. 30 (1972) 491–500.