- Email: [email protected]
Waves in nonlinear fluid filled tubes
197
WAVE MOTION 6 (1984) 197-203 NORTH-HOLLAND
WAVES IN NONLINEAR
FLUID FILLED TUBES
R.J. TAIT and T. Bryant MOODIE Department of Mathematics, University of Alberta, Edmonton, T6G 2G1, Canada Received 9 February 1983 Wave propagation and shock formation in nonlinear elastic and viscoelastic fluid filled tubes is discussed. For a Mooney-Rivlin material a simple exact solution exhibiting distortionless propagation is found.
equation given by
1. Introduction In this paper we investigate the propagation of a pressure pulse along an initially uniform thin walled semi-infinite nonlinear elastic or viscoelastic tube filled with an incompressible, inviscid fluid. It is possible, in the model developed below, to incorporate a compressible fluid in the analysis, but if the fluid is a liquid with say the density and compressibility of water this effect will be negligible. Again we might include the effect of fluid friction, but in a number of practical applications, say for parameters appropriate to biological problems involving the larger blood vessels, this effect is also negligible [l]. We shall derive a one dimensional model dealing with variations of the quantities involved in the axial direction only. The radial variations are not completely neglected, however, as shown in [2], and the equations used result from a perturbation of the two dimensional axially symmetric equations. The governing equations are the conservation of mass equation given by
aA Mu) _.
at+
ax
’
(1.1)
where A denotes the internal area of the tube, t is time, x denotes axial distance, and tr is the fluid velocity in the x direction, and the momentum 0165-2125/84/%3.00
g+
$+l
ap=
ax p ax
0,
(1.2)
where p is the fluid density and P the pressure difference between the inside and outside of the tube. In addition we require a relationship between A, P, x, t appropriate to the non-linear elastic or viscoelastic situation and one such relation is derived below. For the perfectly elastic case dealt with in [3] this relationship is taken in the form @(A, P, x) = 0,
(I.31
where @ is a function determined by the elastic properties and the initial geometry of the tube. It is asserted in [3] that such a relation will be valid provided the radius of the tube varies slowly with distance along the tube. We make a similar assumption in the derivation of our model which includes that discussed in [3]. In the elastic case it is possible, for a MooneyRivlin material, to exhibit closed-form solutions to our problem, which show that no shocks occur. If more complicated strain energy functions are adopted it can be shown, using the method of characteristics, that shocks can occur; that is the solution will break down. In the viscoelastic case the situation is more difficult. However, information can be found on the propagation of the
@ 1984, Elsevier Science Publishers B.V. (North-Holland)
R.J. Tait, TB. Moodie / Waves in fluid filled tubes
198
wavefront, on the behavior of initial derivative discontinuities there, and on the possibility of shock formation.
2. Pressure/area relation In order to derive an appropriate pressure/area relation we discuss briefly the theory of non-linear viscoelasticity. Our notation is standard. Assuming an underlying Cartesian coordinate system we denote by x, X, the spatial and material coordinates of a particle X F denotes the deformation gradient tensor, which we may also write as the product RU where R is an orthogonal tensor, and U a positive definite symmetric tensor. Then B = FFT and C = FTF are the right and left Cauchy tensors, respectively. The constitutive equation for an incompressible, isotropic, simple material may be written as T=--+I+R(t)S;(C)R=(t),
(2.1)
T=-jil+F(t)~~(C)F’(t),
(2.2)
or
where T is the Cauchy or present stress tensor, I is the identity and ST, %z are functionals of the history of C, that is 90”(C)=
Consider now the purely radial motion given by r’=R’+@(t),
(2.5)
where (r, 8, z), (R, 0,Z) are polar coordinates at r, X, respectively. Relative to the underlying Cartesian coordinate system the orthonormal triads associated with these polar systems coincide for a given particle X for all time and we use physical components for all tensors referred to this triad. Then F is symmetric, R = I, and the components of B and C coincide, and T=-#I+.%;(C),
(2.6)
where C ={(Rlr)‘,
(r/R)‘,
1).
(2.7)
Using the methods of Fosdick [4], and Carroll [5], we set Q in (2.4) to be one of the matrices
in turn to show that T is also diagonal. Taking # independent of 0, z we have the equation of motion in the form aT(ll)+
T(ll)-
dr
7-(22)= ki’ r
,
(2.8)
k4+R(t)
T(ll)= (2.3)
-
and # is a scalar pressure. Isotropy implies that
Q%;Q= = %J={QCQ=}.
.Z=Z
where k is the density of the material, which may be integrated to give
; {C(t-S)}, s=o
%2(C) = G {Cft-S)}, s=o
6 = 0,
(2.4)
We will use the form given by equation (2.1), although the argument below applied equally to (2.2). It should be noted however that when the linear approximation to the functionals S, % are introduced they will, in general, lead to entirely different constitutive equations.
{S;(ll>-
I
.%$‘(22)}dr/r
(2.9)
where 0 (t) is arbitrary and we have assumed i’= a+/ar. If now we apply eq. (2.9) to a thin tube of thickness H, initial internal radius R. in the reference configuration and assume that we may neglect radial acceleration terms, we find
P=-5(5f> ~~o”(ll>(Co)-~~(22)(Co)~, 0
0
(2.10)
199
R.J. Tail, T.B. Moadie / Waves in fluid filled tubes
where r0 denotes the current inner radius of the tube, and C, = diag{Ri/r&
rz/Ri, 1).
(2.11)
The motion given by (2.5) is suitable for an axially constrained tube, which is appropriate in a number of practical applications. Major blood vessels, in situ, operate under considerable longitudinal constraint, called tethering [6], and a number of experiments carried out on fluid filled elastic and viscoelastic tubes appear to have the same property. If we consider an incompressible hyperelastic material, strain energy function W(I,, I*), where II, I2 are the principal invariants of C, and so of B, the constitutive eq. (2.6) reduces to T=-#I+rjB-@B-l
(2.12)
so that on substituting for 9 in (2.10) we have
P=#(g++)(l-A;/A’),
(2.13)
0
terms are f = _m{+r(t-X) I
G(C)
tr C(s)1
+~//~(t--s)C(~)}ds.
(2.15)
The first term may be included with the term -41, and it has been suggested in [8] that by choosing coefficients of higher order terms in the expansion as suitable combinations of step functions that the term C(S) may be replaced by (d/ds)R[C(s)], where 0 is a polynomial in C For a Mooney-Rivlin elastic material the strain energy function is taken in the form w=;{6(I,-3)+(1-6)(1,-3)},
(2.16)
where b is a constant and p is taken as the modulus of rigidity for infinitesimal deformation from the undeformed state. Then (2.13) reduces to P=+(l-A;/A2). 0
in agreement with [3]. We note that in the present case I,=I,=l+A,/A+A/A,,
G2(t) = ~{m + (1 - m) e-“‘},
where we have written A0 = mR& A = 1~rgand 4=2$,
$=2$Y 1
This suggests that by analogy with the linear situation we take
2
In [3] it was proposed that if A varies slowly with x that equation (2.13) will afford a good approximation to the pressure-area relation when P, A are considered as functions of x, t. We extend this approximation to the expression (2.10) and rewrite this expression in the form P=-~~I~o(ll)(c)-~o(22)(C)}, 0
(2.14)
(2.17)
where 0 < m d 1, and r is a relaxation constant. In this way, using (2.14), (2.15) and (2.17) we have a simple possible generalisation to nonlinear viscoelasticity of the Mooney-Rivlin material, with a model exhibiting impact and long time non-linear elastic response. With these assumptions the governing equations of our model are
aA off
a
=o
ap=
E+ g!+l ax
(2.18)
o
p
ax
(2.19)
where C = diag{Ao/A, A/Ao, 1) and
A = A(x, t).
In general 9 is given as a multiple integral expansion in terms of C, see Lockett [7]. The linear
AP=py
’ {m+(l-m)e-“-S)“} I -al
+I’(/i) ds,
(2.20)
200
R.J. Tair. T B. Moodie / Waves in fluid filled tubes
where we have written 2 = AlAo, y = H/R,,, and @ is determined by R. In the linear case (2.15) 0 will be given by (A - a-r>. More generally since C is diagonal and R a polynomial in C, @ will have a common factor (A- A-‘) so that @,(1) = 0.
?=7
the equations by setting
i=t
P=P/p, ’ Jl P
Ro, D=v
.f =x/R,,
X=/L
A = AlAo,
(3.1)
= 0,
(3.2)
v,+vv,+P,=0, AP=y
(3.3)
’ {m+(l-m) e-“-““‘} J-02
x$@(A)
ds,
(3.4)
in general, and in the elastic case equation (3.4) is replaced by P= yd&(A)A-‘,
(3.5)
where &(A)
=
?(A-A-r),
(3.6)
and 4, + are as before. Differentiating (3.4) with respect to t gives (;+t)(AP)=y(;+;)@(A),
(3.7)
so that the system of equations (3.2), (3.3), and
(3.8)
in which case the characteristic directions are given by v*Q.
(3.9)
In the elastic case, 0 is omitted and P, @ replaced by (3.5), (3.6). Consider now a medium initially at rest with initial velocity of the fluid being zero, initial pressure PO, and initial cross sectional area of the tube A,, and suppose a pressure is imposed at x = 0, that is, P(0, t)=f(t),
Inserting these quantities and dropping the bars gives A, + (Au),
Q*= @‘(A)-P>O,
dx z=O,
3. Properties of solutions We non-dimensionalize
(3.7) is hyperbolic provided
t>0,
(3.10)
with f(0) = PO. It follows that any discontinuities in the derivatives of P, assuming no shocks are formed, will propagate into the medium with velocity Qo, where Q; = y@‘(Ao) -PO.
(3.11)
In the viscoelastic case we assume that the initial state has persisted for a long time so that PO= ymA;r@(A,J,
(3.12)
and is determined by the relaxed modulus. In the case governed by eq. (2.15) this gives Qi =(l-m)+(l+m)Ao*, which in the elastic case m = 1, reduces to the wave speed for a Mooney-Rivlin material. This predicts that the speed of propagation of the wave is dependent on the initial ambient pressure. Equation (3.12) may also imply a restriction on the possible magnitude of the pressure. Thus in the MooneyRivlin case, m = 1, @ = (A - A-‘), it is clear that PO< y, and in the generalised viscoelastic case that P,
(3.13)
201
R.J. Tait, T.B. Moodie / Waves in fluid filled tubes
where u is a column vector of order 2 or 3 as is g. D is a square matrix of corresponding order. We assume that initially we have a uniform state uO,with g( u,,) = 0. We assume also that one of the entries of u say P, or in the elastic case A, is known at x=0 for ta0. We take D to have distinct eigenvalues A, and denote the largest positive eigenvalue by A+corresponding to the characteristic curve C+, and assume that D has a full set of left and right eigenvectors I, r respectively with I,, r+ corresponding to A+. Along the wavefront we denote the appropriate quantities by C”,, A:, r:, 1% We are concerned with the propagation of discontinuities in derivatives along C”,, passing through the origin of coordinates. Discontinuities in first derivatives are simpler to deal with, though second order derivatives would be more physically realistic. In the present case the equations to be derived below are essentially the same whether we deal with C’ discontinuities, or, in the absence of C’ discontinuities, with C2 discontinuities. We shall then confine our attention to the simpler situation. We use the work of Jeffrey [9] in the following. Introducing new variables 5=5(x,t),
5(x,0)=x;
with 5 constant along C+ characteristics, ==
[ql,
(3.14)
v=t
we set
x = [x,17
(3.15)
where [q] denotes the jump in q across C”,, that is [ql= q(O-) -q(O+). The analysis in [9] shows that 7r is a multiple of the right eigenvector of D on Ct and that II=ut+,
0
(3.16) (3.17)
Consider, first of all, the elastic situation. Here
u=($ D=&yA)t),
b=O, (3.19)
with P given by (3.5), (3.6). The wavefront propagates along C”, : x =JA,P’(A,)t,
(3.20)
and rr is constant there, given by 1 *=uo
( JP’(A,)/A,
>’
(3.21)
where a0 is the initial jump in A,. Since A+==++,
(3.22)
eq. (3.18) reads
dx uo -=dq 2
J+P”(Ao)+3P’(AoYAo~. 0
(3.23) It is clear that the Mooney-Rivlin material is an exceptional case since (VA+)‘n = 0 for any Ao. It is in fact clear that in this case we have an explicit solution P(x, t) =f(t--X/Co), A(x, t) =[l-P/Y]-“~,
u(x, t) = C,[l -AoA-‘I,
(3.24)
where Ci = 2rA02 provided P/y < 1. This indicates that waves of this type travel down the tube, without distortion, with velocity depending on the initial pressure. For biological applications it is suggested in [3] that a more appropriate model is provided by a strain energy function of the form w=;{b(&-3)+S(&-3)2+(I-b)(12-3)}. (3.25)
$=
(V.A+)‘lr,
along C”, : dxfdt = AT, and V, denotes gradient operator with respect to u.
(3.18) the
In this case (3.23) becomes
dx co -=dv 2
J
$j{
I- A04)6SyA;‘, 0
(3.26)
202
R.J. Tait, T.B. Moodie / Waves in f7uid filled tubes
where P=r(l-A-2){1+2S(A+A-1-2)},
(3.27)
so that P’(A,)>O if A,,2 1. If w denotes the positive jump in #l&(0, 0) then on computing a0 we see that xc = 0 at time fB given by te = Ao(l -A,4)-1[P’(Ao)]2/306y,
(3.28)
and breakdown occurs along the wave front. In general the initial breakdown occurs behind the front. Thus if A0 = 1, no breakdown occurs along the front but does so behind the front. For further details we refer to [3] where typographical errors should be noted in eqs. (23) and (28). Next consider the viscoelastic case. We have
and if w denotes the positive jump in aP/at (0,O) we may replace a0 by -o/Q& On integrating (3.33), we find that if ‘J’=(l-dQ;[Qf+Po] -
,T{Q; + ‘yAo@“(Ao)~,
(3.34)
is positive or zero, no breakdown can occur along the front. On the other hand if ?P < 0, breakdown will occur at time ti given by t”, = -7 In{-V/WT[Qi +
~Ao@“(Ao)I~
x(l+PolQ;).
(3.35)
Consider again the simplest case when @(A) = A-A-‘. Then F=y’(l-m){4-(4+2m+G)i + (1 + m)P’},
with Q as in (3.11) and
(3.30) g= LLlJ
*
(3.37)
then Cy2 0. IAs pointed out earlier 0 4 P < 1 with
(3.31)
A0 + co as P + 1. Thus for low initial pressures large jumps can be tolerated in the derivative aP/at at time t =0 without shocks developing on the
wavefront. As PO increases the allowable range of c3 is decreased.
where o=o,exp{
-(y)(l+P,io$)t/T}, (3.32)
along Here initial front.
where we have taken G = UT/-~, p = PO/my. It follows that if GS4/b+(l+m)F-2(2+m),
It now follows that
(3.36)
C”, : x T+jmot,where Q. is given by (3.10). a0 is tl. _Iinitial jump. It follows that any discontinuity decays exponentially along the In addition we have
dx= a{Q: + ‘yAo@‘(Ao)V2Qo, dT
(3.33)
4. Conclusion In the preceding paragraphs we have developed a possible model for the propagation of pressure pulses in nonlinear elastic and viscoelastic fluid filled tubes. In the elastic case, for a MooneyRivlin material, we have displayed a simple exact solution exhibiting waves propagating without distortion. If the quantities involved are infinitesimal
R.J. Tait, T.B. Moodie / Waves in @id filled tubes
this reduces to a simple progressing wave of fixed velocity, whereas the general case shows that the velocity is a function of the initial pressure. For other types of elastic materials one encounters shock formation. One would not expect the model to be appropriate in the neighbourhood of a shock. In the viscoelastic case the model predicts that any initial discontinuities in derivatives will decay in magnitude along the wave front. The condition for shock formation then depends on several parameters, w, 7, m, PO. In the simple generalisation of a Mooney-Rivlin material which we have chosen, if the parameters fall in a reasonable range no shock will form on the wavefront and we conjecture that none will form behind the wavefront. To proceed further will apparently require numerical calculation of the solution. Whether or not for more complicated viscoelastic materials shock formation will occur will apparently have to be answered by computational methods.
203
References 111G. Rudinger,
“Shock waves in mathematical models of the aorta”, J. Appl. Me& 37, 34 (1970). 121J.W. Lambert, “Fluid flow in a nonrigid tube”, Doctoral Dissertafion Series 19,418, University Microfilms Inc., Ann Arbor, Mich. (1956). [31 T.B. Moodie and J.B. Haddow, “Waves in thin walled elastic tubes containing an incompressible inviscid fluid”, Inr. J. Non-Linear Mechanics 12, 223 (1977). possible motions of incomr41 R.L. Fosdick, “Dynamically pressible, isotropic, simple materials”, Arch. Rat. Mech. Anal. 29, 272 (1968). deformations of incompress[51 M.M. Carroll, “Controllable ible simple materials”, In?. J. Eng. SC. 5, 515 (1967). [61 W.G. Frasher, “What is known about the physiology of larger blood vessels”, Biomechanics Symposium, ed. Y.C. Fung, New York, ASME, App. Mech. Div., 1 (1966). [71 F.J. Lockett, Nonlinear Viscoelastic Solids, Academic Press (1972). PI D.A. Spence, Nonlinear Wave Propagation in Viscoelastic Materials in Nonlinear Elasticify, ed. R.W. Dickey, Academic Press (1973). hyperbolic systems and waves”, [91 A. Jeffrey, “Quasilinear Research Notes in Math., Pitman (1976).
WAVE MOTION 6 (1984) 197-203 NORTH-HOLLAND
WAVES IN NONLINEAR
FLUID FILLED TUBES
R.J. TAIT and T. Bryant MOODIE Department of Mathematics, University of Alberta, Edmonton, T6G 2G1, Canada Received 9 February 1983 Wave propagation and shock formation in nonlinear elastic and viscoelastic fluid filled tubes is discussed. For a Mooney-Rivlin material a simple exact solution exhibiting distortionless propagation is found.
equation given by
1. Introduction In this paper we investigate the propagation of a pressure pulse along an initially uniform thin walled semi-infinite nonlinear elastic or viscoelastic tube filled with an incompressible, inviscid fluid. It is possible, in the model developed below, to incorporate a compressible fluid in the analysis, but if the fluid is a liquid with say the density and compressibility of water this effect will be negligible. Again we might include the effect of fluid friction, but in a number of practical applications, say for parameters appropriate to biological problems involving the larger blood vessels, this effect is also negligible [l]. We shall derive a one dimensional model dealing with variations of the quantities involved in the axial direction only. The radial variations are not completely neglected, however, as shown in [2], and the equations used result from a perturbation of the two dimensional axially symmetric equations. The governing equations are the conservation of mass equation given by
aA Mu) _.
at+
ax
’
(1.1)
where A denotes the internal area of the tube, t is time, x denotes axial distance, and tr is the fluid velocity in the x direction, and the momentum 0165-2125/84/%3.00
g+
$+l
ap=
ax p ax
0,
(1.2)
where p is the fluid density and P the pressure difference between the inside and outside of the tube. In addition we require a relationship between A, P, x, t appropriate to the non-linear elastic or viscoelastic situation and one such relation is derived below. For the perfectly elastic case dealt with in [3] this relationship is taken in the form @(A, P, x) = 0,
(I.31
where @ is a function determined by the elastic properties and the initial geometry of the tube. It is asserted in [3] that such a relation will be valid provided the radius of the tube varies slowly with distance along the tube. We make a similar assumption in the derivation of our model which includes that discussed in [3]. In the elastic case it is possible, for a MooneyRivlin material, to exhibit closed-form solutions to our problem, which show that no shocks occur. If more complicated strain energy functions are adopted it can be shown, using the method of characteristics, that shocks can occur; that is the solution will break down. In the viscoelastic case the situation is more difficult. However, information can be found on the propagation of the
@ 1984, Elsevier Science Publishers B.V. (North-Holland)
R.J. Tait, TB. Moodie / Waves in fluid filled tubes
198
wavefront, on the behavior of initial derivative discontinuities there, and on the possibility of shock formation.
2. Pressure/area relation In order to derive an appropriate pressure/area relation we discuss briefly the theory of non-linear viscoelasticity. Our notation is standard. Assuming an underlying Cartesian coordinate system we denote by x, X, the spatial and material coordinates of a particle X F denotes the deformation gradient tensor, which we may also write as the product RU where R is an orthogonal tensor, and U a positive definite symmetric tensor. Then B = FFT and C = FTF are the right and left Cauchy tensors, respectively. The constitutive equation for an incompressible, isotropic, simple material may be written as T=--+I+R(t)S;(C)R=(t),
(2.1)
T=-jil+F(t)~~(C)F’(t),
(2.2)
or
where T is the Cauchy or present stress tensor, I is the identity and ST, %z are functionals of the history of C, that is 90”(C)=
Consider now the purely radial motion given by r’=R’+@(t),
(2.5)
where (r, 8, z), (R, 0,Z) are polar coordinates at r, X, respectively. Relative to the underlying Cartesian coordinate system the orthonormal triads associated with these polar systems coincide for a given particle X for all time and we use physical components for all tensors referred to this triad. Then F is symmetric, R = I, and the components of B and C coincide, and T=-#I+.%;(C),
(2.6)
where C ={(Rlr)‘,
(r/R)‘,
1).
(2.7)
Using the methods of Fosdick [4], and Carroll [5], we set Q in (2.4) to be one of the matrices
in turn to show that T is also diagonal. Taking # independent of 0, z we have the equation of motion in the form aT(ll)+
T(ll)-
dr
7-(22)= ki’ r
,
(2.8)
k4+R(t)
T(ll)= (2.3)
-
and # is a scalar pressure. Isotropy implies that
Q%;Q= = %J={QCQ=}.
.Z=Z
where k is the density of the material, which may be integrated to give
; {C(t-S)}, s=o
%2(C) = G {Cft-S)}, s=o
6 = 0,
(2.4)
We will use the form given by equation (2.1), although the argument below applied equally to (2.2). It should be noted however that when the linear approximation to the functionals S, % are introduced they will, in general, lead to entirely different constitutive equations.
{S;(ll>-
I
.%$‘(22)}dr/r
(2.9)
where 0 (t) is arbitrary and we have assumed i’= a+/ar. If now we apply eq. (2.9) to a thin tube of thickness H, initial internal radius R. in the reference configuration and assume that we may neglect radial acceleration terms, we find
P=-5(5f> ~~o”(ll>(Co)-~~(22)(Co)~, 0
0
(2.10)
199
R.J. Tail, T.B. Moadie / Waves in fluid filled tubes
where r0 denotes the current inner radius of the tube, and C, = diag{Ri/r&
rz/Ri, 1).
(2.11)
The motion given by (2.5) is suitable for an axially constrained tube, which is appropriate in a number of practical applications. Major blood vessels, in situ, operate under considerable longitudinal constraint, called tethering [6], and a number of experiments carried out on fluid filled elastic and viscoelastic tubes appear to have the same property. If we consider an incompressible hyperelastic material, strain energy function W(I,, I*), where II, I2 are the principal invariants of C, and so of B, the constitutive eq. (2.6) reduces to T=-#I+rjB-@B-l
(2.12)
so that on substituting for 9 in (2.10) we have
P=#(g++)(l-A;/A’),
(2.13)
0
terms are f = _m{+r(t-X) I
G(C)
tr C(s)1
+~//~(t--s)C(~)}ds.
(2.15)
The first term may be included with the term -41, and it has been suggested in [8] that by choosing coefficients of higher order terms in the expansion as suitable combinations of step functions that the term C(S) may be replaced by (d/ds)R[C(s)], where 0 is a polynomial in C For a Mooney-Rivlin elastic material the strain energy function is taken in the form w=;{6(I,-3)+(1-6)(1,-3)},
(2.16)
where b is a constant and p is taken as the modulus of rigidity for infinitesimal deformation from the undeformed state. Then (2.13) reduces to P=+(l-A;/A2). 0
in agreement with [3]. We note that in the present case I,=I,=l+A,/A+A/A,,
G2(t) = ~{m + (1 - m) e-“‘},
where we have written A0 = mR& A = 1~rgand 4=2$,
$=2$Y 1
This suggests that by analogy with the linear situation we take
2
In [3] it was proposed that if A varies slowly with x that equation (2.13) will afford a good approximation to the pressure-area relation when P, A are considered as functions of x, t. We extend this approximation to the expression (2.10) and rewrite this expression in the form P=-~~I~o(ll)(c)-~o(22)(C)}, 0
(2.14)
(2.17)
where 0 < m d 1, and r is a relaxation constant. In this way, using (2.14), (2.15) and (2.17) we have a simple possible generalisation to nonlinear viscoelasticity of the Mooney-Rivlin material, with a model exhibiting impact and long time non-linear elastic response. With these assumptions the governing equations of our model are
aA off
a
=o
ap=
E+ g!+l ax
(2.18)
o
p
ax
(2.19)
where C = diag{Ao/A, A/Ao, 1) and
A = A(x, t).
In general 9 is given as a multiple integral expansion in terms of C, see Lockett [7]. The linear
AP=py
’ {m+(l-m)e-“-S)“} I -al
+I’(/i) ds,
(2.20)
200
R.J. Tair. T B. Moodie / Waves in fluid filled tubes
where we have written 2 = AlAo, y = H/R,,, and @ is determined by R. In the linear case (2.15) 0 will be given by (A - a-r>. More generally since C is diagonal and R a polynomial in C, @ will have a common factor (A- A-‘) so that @,(1) = 0.
?=7
the equations by setting
i=t
P=P/p, ’ Jl P
Ro, D=v
.f =x/R,,
X=/L
A = AlAo,
(3.1)
= 0,
(3.2)
v,+vv,+P,=0, AP=y
(3.3)
’ {m+(l-m) e-“-““‘} J-02
x$@(A)
ds,
(3.4)
in general, and in the elastic case equation (3.4) is replaced by P= yd&(A)A-‘,
(3.5)
where &(A)
=
?(A-A-r),
(3.6)
and 4, + are as before. Differentiating (3.4) with respect to t gives (;+t)(AP)=y(;+;)@(A),
(3.7)
so that the system of equations (3.2), (3.3), and
(3.8)
in which case the characteristic directions are given by v*Q.
(3.9)
In the elastic case, 0 is omitted and P, @ replaced by (3.5), (3.6). Consider now a medium initially at rest with initial velocity of the fluid being zero, initial pressure PO, and initial cross sectional area of the tube A,, and suppose a pressure is imposed at x = 0, that is, P(0, t)=f(t),
Inserting these quantities and dropping the bars gives A, + (Au),
Q*= @‘(A)-P>O,
dx z=O,
3. Properties of solutions We non-dimensionalize
(3.7) is hyperbolic provided
t>0,
(3.10)
with f(0) = PO. It follows that any discontinuities in the derivatives of P, assuming no shocks are formed, will propagate into the medium with velocity Qo, where Q; = y@‘(Ao) -PO.
(3.11)
In the viscoelastic case we assume that the initial state has persisted for a long time so that PO= ymA;r@(A,J,
(3.12)
and is determined by the relaxed modulus. In the case governed by eq. (2.15) this gives Qi =(l-m)+(l+m)Ao*, which in the elastic case m = 1, reduces to the wave speed for a Mooney-Rivlin material. This predicts that the speed of propagation of the wave is dependent on the initial ambient pressure. Equation (3.12) may also imply a restriction on the possible magnitude of the pressure. Thus in the MooneyRivlin case, m = 1, @ = (A - A-‘), it is clear that PO< y, and in the generalised viscoelastic case that P,
(3.13)
201
R.J. Tait, T.B. Moodie / Waves in fluid filled tubes
where u is a column vector of order 2 or 3 as is g. D is a square matrix of corresponding order. We assume that initially we have a uniform state uO,with g( u,,) = 0. We assume also that one of the entries of u say P, or in the elastic case A, is known at x=0 for ta0. We take D to have distinct eigenvalues A, and denote the largest positive eigenvalue by A+corresponding to the characteristic curve C+, and assume that D has a full set of left and right eigenvectors I, r respectively with I,, r+ corresponding to A+. Along the wavefront we denote the appropriate quantities by C”,, A:, r:, 1% We are concerned with the propagation of discontinuities in derivatives along C”,, passing through the origin of coordinates. Discontinuities in first derivatives are simpler to deal with, though second order derivatives would be more physically realistic. In the present case the equations to be derived below are essentially the same whether we deal with C’ discontinuities, or, in the absence of C’ discontinuities, with C2 discontinuities. We shall then confine our attention to the simpler situation. We use the work of Jeffrey [9] in the following. Introducing new variables 5=5(x,t),
5(x,0)=x;
with 5 constant along C+ characteristics, ==
[ql,
(3.14)
v=t
we set
x = [x,17
(3.15)
where [q] denotes the jump in q across C”,, that is [ql= q(O-) -q(O+). The analysis in [9] shows that 7r is a multiple of the right eigenvector of D on Ct and that II=ut+,
0
(3.16) (3.17)
Consider, first of all, the elastic situation. Here
u=($ D=&yA)t),
b=O, (3.19)
with P given by (3.5), (3.6). The wavefront propagates along C”, : x =JA,P’(A,)t,
(3.20)
and rr is constant there, given by 1 *=uo
( JP’(A,)/A,
>’
(3.21)
where a0 is the initial jump in A,. Since A+==++,
(3.22)
eq. (3.18) reads
dx uo -=dq 2
J+P”(Ao)+3P’(AoYAo~. 0
(3.23) It is clear that the Mooney-Rivlin material is an exceptional case since (VA+)‘n = 0 for any Ao. It is in fact clear that in this case we have an explicit solution P(x, t) =f(t--X/Co), A(x, t) =[l-P/Y]-“~,
u(x, t) = C,[l -AoA-‘I,
(3.24)
where Ci = 2rA02 provided P/y < 1. This indicates that waves of this type travel down the tube, without distortion, with velocity depending on the initial pressure. For biological applications it is suggested in [3] that a more appropriate model is provided by a strain energy function of the form w=;{b(&-3)+S(&-3)2+(I-b)(12-3)}. (3.25)
$=
(V.A+)‘lr,
along C”, : dxfdt = AT, and V, denotes gradient operator with respect to u.
(3.18) the
In this case (3.23) becomes
dx co -=dv 2
J
$j{
I- A04)6SyA;‘, 0
(3.26)
202
R.J. Tait, T.B. Moodie / Waves in f7uid filled tubes
where P=r(l-A-2){1+2S(A+A-1-2)},
(3.27)
so that P’(A,)>O if A,,2 1. If w denotes the positive jump in #l&(0, 0) then on computing a0 we see that xc = 0 at time fB given by te = Ao(l -A,4)-1[P’(Ao)]2/306y,
(3.28)
and breakdown occurs along the wave front. In general the initial breakdown occurs behind the front. Thus if A0 = 1, no breakdown occurs along the front but does so behind the front. For further details we refer to [3] where typographical errors should be noted in eqs. (23) and (28). Next consider the viscoelastic case. We have
and if w denotes the positive jump in aP/at (0,O) we may replace a0 by -o/Q& On integrating (3.33), we find that if ‘J’=(l-dQ;[Qf+Po] -
,T{Q; + ‘yAo@“(Ao)~,
(3.34)
is positive or zero, no breakdown can occur along the front. On the other hand if ?P < 0, breakdown will occur at time ti given by t”, = -7 In{-V/WT[Qi +
~Ao@“(Ao)I~
x(l+PolQ;).
(3.35)
Consider again the simplest case when @(A) = A-A-‘. Then F=y’(l-m){4-(4+2m+G)i + (1 + m)P’},
with Q as in (3.11) and
(3.30) g= LLlJ
*
(3.37)
then Cy2 0. IAs pointed out earlier 0 4 P < 1 with
(3.31)
A0 + co as P + 1. Thus for low initial pressures large jumps can be tolerated in the derivative aP/at at time t =0 without shocks developing on the
wavefront. As PO increases the allowable range of c3 is decreased.
where o=o,exp{
-(y)(l+P,io$)t/T}, (3.32)
along Here initial front.
where we have taken G = UT/-~, p = PO/my. It follows that if GS4/b+(l+m)F-2(2+m),
It now follows that
(3.36)
C”, : x T+jmot,where Q. is given by (3.10). a0 is tl. _Iinitial jump. It follows that any discontinuity decays exponentially along the In addition we have
dx= a{Q: + ‘yAo@‘(Ao)V2Qo, dT
(3.33)
4. Conclusion In the preceding paragraphs we have developed a possible model for the propagation of pressure pulses in nonlinear elastic and viscoelastic fluid filled tubes. In the elastic case, for a MooneyRivlin material, we have displayed a simple exact solution exhibiting waves propagating without distortion. If the quantities involved are infinitesimal
R.J. Tait, T.B. Moodie / Waves in @id filled tubes
this reduces to a simple progressing wave of fixed velocity, whereas the general case shows that the velocity is a function of the initial pressure. For other types of elastic materials one encounters shock formation. One would not expect the model to be appropriate in the neighbourhood of a shock. In the viscoelastic case the model predicts that any initial discontinuities in derivatives will decay in magnitude along the wave front. The condition for shock formation then depends on several parameters, w, 7, m, PO. In the simple generalisation of a Mooney-Rivlin material which we have chosen, if the parameters fall in a reasonable range no shock will form on the wavefront and we conjecture that none will form behind the wavefront. To proceed further will apparently require numerical calculation of the solution. Whether or not for more complicated viscoelastic materials shock formation will occur will apparently have to be answered by computational methods.
203
References 111G. Rudinger,
“Shock waves in mathematical models of the aorta”, J. Appl. Me& 37, 34 (1970). 121J.W. Lambert, “Fluid flow in a nonrigid tube”, Doctoral Dissertafion Series 19,418, University Microfilms Inc., Ann Arbor, Mich. (1956). [31 T.B. Moodie and J.B. Haddow, “Waves in thin walled elastic tubes containing an incompressible inviscid fluid”, Inr. J. Non-Linear Mechanics 12, 223 (1977). possible motions of incomr41 R.L. Fosdick, “Dynamically pressible, isotropic, simple materials”, Arch. Rat. Mech. Anal. 29, 272 (1968). deformations of incompress[51 M.M. Carroll, “Controllable ible simple materials”, In?. J. Eng. SC. 5, 515 (1967). [61 W.G. Frasher, “What is known about the physiology of larger blood vessels”, Biomechanics Symposium, ed. Y.C. Fung, New York, ASME, App. Mech. Div., 1 (1966). [71 F.J. Lockett, Nonlinear Viscoelastic Solids, Academic Press (1972). PI D.A. Spence, Nonlinear Wave Propagation in Viscoelastic Materials in Nonlinear Elasticify, ed. R.W. Dickey, Academic Press (1973). hyperbolic systems and waves”, [91 A. Jeffrey, “Quasilinear Research Notes in Math., Pitman (1976).