- Email: [email protected]
Inflation waves in non-linear distensible tubes
hr. J. Engng Sri. Vol. 22, No. Printed in Great Britain.
INFLATION
2, pp.
101-118,
1984
0020-7225/84 Pergamon
WAVES IN NON-LINEAR
DISTENSIBLE
$3.00 + .OO Press Ltd.
TUBESt
D. ELAD, A. FOUX, Y. KIVITY and Y. LANIR Department of Biomedical Engineering, Technion-Israel Institute of Technology, Haifa 32000, Israel Abstract-A mathematical model for large amplitude wave propagation in a thin walled distensible tube is developed. The tube wall is considered as a membranic shell made of an incompressible, non-linear viscoelastic material with cylindrical orthotropy. The fluid is regarded as incompressible and inviscid and the flow is quasi-one-dimensional. The case of a pressure step applied at one end of a uniform straight tube is solved as an example. The system of partial differential equations, describing the motions of the fluid and the wall, are integrated numerically by using a two-step explicit scheme. Flow and deformation variables as well as the wave velocity are determined in time and space.
1. INTRODUCTION
rupture of the aorta in vehicular collision, which is responsible for significant percentage of traffic fatalities[l-31 was the motivation to this work. The problem is one of dynamics of fluid-filled distensible tubes. The object of this work is to study one phase of the general problem, namely the wave phenomena. Waves in distensible tubes filled with incompressible fluid may be sustained by inflation of the distensible tube wall and are sometimes called inflation waves. Such waves are developed through the interaction of the fluid and the tube wall and may be of large deformations. This interaction depends strongly upon the mechanical properties of the wall material as well as on the mechanical constraints imposed on the tube by its surrounding. There have been many studies of fluid flow through distensible tubes in recent years. Many of them were motivated by problems of blood flow in the vascular system. An often used method for analysis of such problems is based on the assumption that the fluid-wall interaction can be described by the “tube law” relating the local cross-sectional area of the tube to the local transmural pressure and to time[4-91. Beyond first approximation, it is unlikely that a simple tube law can account for the biaxial behaviour of a material which may be non-linear viscoelastic and orthotropic. In the present work, it is intended to present a method for explicit solution of the problem of wave propagation in fluid-filled distensible tubes, with particular emphasis on local stresses and stretches. An attempt is made to formulate the problem in a general frame through the use of non-linear shell theory. This allows for such factors as tapering and curvature of the tube as well as different boundary and loading conditions to be considered. For illustration purpose a specific solution is presented for the case of a cylindrical straight tube subjected to a pressure step at one end. TRAUMATIC
2. ANALYSIS
2.1 Tube geometry and deformation The tube, regarded as a thin membranic shell, is defined by a reference surface. The inner surface of the tube is conveniently chosen to be the reference surface. As it is also the outer surface of the fluid, the interaction will be defined on it. In order to define this reference surface a fixed Cartesian frame x, y, z is chosen in such a way that, under loading of internal pressure and axial tension, the curved axis of the tube is assumed to remain in the x-y plane (see Fig. 1). The radius vector to a point of the tube axis is
~(5, t) = ~(5, tY + y(5, t)j tThis research was supported by a grant from the United States-Israel (BSF), Jerusalem, Israel. 103
(1) Binational Science Foundation
D. ELAD et al.
104
where 5 is a non-dimensional Lagrangian parameter along the tube axis, t is the time and i, j are unit vectors in the x and y directions respectively. Real lengths along the tube axis, denoted by Z, are given by AZ = aA<
(2)
u2 = (x,J2 + 01,J2
(3)
where TVis a scaling variable given by
The inner surface of the she11is formed by closed planar curves C, enclosing the axis. A point on the curve C is described by the radius vector r(& 0, t) emanating from the point t on the axis, The angle 8, defining the direction of this vector in the plane of C, is another Lagrangian parameter. The following is assumed with respect to the deformation of the tube: (a) The curve C remains planar during defo~ation and the normal to the plane enclosed by C at the intersection with the tube axis, remains tangent to the axis. (b) The cross-sectional curves C are circular before and during deformation, hence, rrO= 0. (c) The deformations are in the principle directions without shear, hence there is no torsion of the tube and circumferential displa~ment are negligible.
The first assumption is analogous to Kirchhoff s hypothesis for thin plates [lo]. The second assumption ensures orthogonality of the shell coordinates, It is valid in straight tubes and is a first approximation in curved tubes. Similar considerations and the loading conditions lead to the third assumption. Under these assumptions the radius vector r at any time is given by r(r, B, t) = Ir(& 8, t)l[ - cos 8n + sin Ok]
(4)
where k is a unit vector in the z direction. The variable unit vectors 1, II, k (Fig. 1) represent a local orthogonal frame of reference that move with the curved axis as the tube deforms. They may be conventionally used to calculate the covariant base vectors at a point on the shell. The vectors 1 and II are unit vectors in the x-y plane, tangent and perpendicular to the curved axis respectively. The reference surface of the tube shell at any given time is defined in the fixed coordinate system by the radius vector to a point on it, namely
This vector enables to calculate the deformations and can be used in the equilibrium equations of the shell [I 1, 121.The parameters r and 8 may serve as a system of curvilinear coordinates on the shell surface (see Fig. 1) and are variably named “shell coordinates”, “material coordinates” or “Lagrangian coordinates”. For simplicity of notation the time variable t will not be indicated explicitly in all cases. The angle d, (Fig. 1) in the x-y plane is the angle between the tangent to the axis (defined by 1) and the x-axis, hence:
cosc#=;x,”
I
1
sin &,= ;y,(.
The unit vectors I and n are given by 1 = p,s = x,,i + y,d = cos 4i + sin (bj (7) n = -sin #i + cos f$j.
fnfiatiou
waves in
non-linear distensible tubes
in the
x-y
105
plane
Parallel to the
Diameter in the
Of the tube x-y
plane
4
n
k
fixed system
- unit vectors along
local Cartesian coordinates .on the curved axis R,{ R,,N-
Vectors alOng
curvilinear Cartesian coordinates on the shell
Fig. 1. The tube model and
It can be shown that the first derivatives
itsreference
systems.
of 1 and n are
These variables simplify the mathematical representation of the covariant base vectors. The radius vector to a point on the shell surface, after imerting eqn (4) into eqn (5) iS
R(t, 8, t) = p([, 1) - r cos On -t r sin Ok. The covariant
(9)
base vectors will now be
R, = R, = (a + r& cm S)l - (r,ccm f3)n+ fr,esinB)k Rr2 = R,.
= (r sin O)n +
The notation of the shell coot4inates f , 2 respectively.
(r cm O)k.
<,19 will be replaced in the following
equations
by
106
D.
F&AD
et
al.
The covariant base vectors are used to calculate other geometrical components af the metric tensor AaB will be
A,, = R,2- R,2 =
r2
A,, = R,, qR,2=
T,~T,~ =
A = AI,A22-
The base vectors are ~r~ndic~~ar
variables. The
0
(11)
A:, = r”[(a + F#,j cos ey+ (r*,)“l. to each ather since AI2 = 0,
The unit vector normal to the surface is
- cos B(a + r&l
cas @n-I- sin@(a+ r& cos @}k].
w
Hence, the vectors R,,; R,,; N form a local arthogonaf system that move with the shell surface. From N and the derivatives of the base vectors (R,aB), the components of the curvature tensor can be calculated. Hence,
thus,
Inflation waves in non-linear distensible tubes
The components
of length elements along the shell coordinates
107
are given by [1 1, 121
4 =&T& As,
For a given deformation be written as
(16)
= ,,&A&
that occurred during a time period t, the extension ratios may
(As,),
A1=o,= (17)
(AsJr
12=(d=
II--(&),
G422h
= 0’
2.2 Constitutive relations Inflation waves are associated with large deformations of the tube wall, hence constitutive relations of the wall material must be valid for finite deformations. The wall material is considered as an incompressible, non-linear viscoelastic material with cyiindrical orthotropy. To facilitate the numerical calculations, it is necessary to express the viscoelastic response in a very simple mathematical form which accounts for the effect of time. We chose the approach of Collins and Kivity [13, 141 which state that:
(a
o,(t) = D&l,, A,)
not summed).
(18)
Here, D,--is the quasi-static non-linear, biaxial elastic response and b&/A, is a uniaxial rate dependent function where & is the stretch rate and b, is a material parameter. Although eqn (18) ignores the fading memory aspect of the viscoelastic behavior it does account for the effect of rate. This may be justified as an ad hoc constitutive modeling for analysis of a single abrupt event (a large amplitude inflation wave) in which the deformation history plays a small role, but the effect of stretching-rate may be significant. The case under consideration is one where material axes coincide with geometric axes. With load conditions of internal pressure and axial tension the case may be reduced to one of biaxial normal stresses if the small radial stress is neglected. A model for the non-linear elastic response (D,) of such material was proposed by Foux et al. [15], and adopted in this work. Non-linearity is introduced, in this model, through the definitions of generalized non-linear strain and stress in such a way that the stress-strain relations are quasi-linear. This enables inversion of the stress-strain equations, an important feature in the solution of mixed boundary value problems. The rational and justifications of the following choice for biaxial nonlinear orthotropic stress-strain law is given in Foux et al. [15]. Here we shall outline the major points and results. The strain is defined by the generalized strain measure as proposed by Seth [ 161
c,+n:1)
(19)
where n is the non-linearity material parameter. The associated generalized stress is obtained by differentiation of the strain energy density function with respect to the generalized strain. The generalized stress is given by s, ==$ 01 where ca is the true (Eulerian) stress.
(a not summed)
(20)
108
D. ELAD er al.
The strain-stress relations are expressed, in analogy with the generalized Hooke’s law for an orthotropic material [17], as
(21)
where E8 are material parameters analogous to Young’s moduli and vE6are lateral-strainparameters analogous to Poisson’s ratios. Following Foux et al. [IS] we will take E, to be constants. The lateral-strain-parameters must therefore be strain dependent to allow for material incompressibility which requires that L&3 = 1.
(22)
Equation (22) imposes a restriction on the strain dependent which are thus functionally related. From symmetry we have: V21 ---_
Vtz
Et
-5
lateral-strain-parameters
(231
It is also assumed 1151that V3t
=
(A;&)-”
-w2t
V32
=
(A~Jp,)-(’
-831)Vt2
(24) where pPla= E,/EB and K is another material parameter. The value of the lateral strain parameters may thus be calculated for a given state of deformation from eqs (19)-(24). The final form of the constitutive equations expressed in terms of the extensions and their time derivatives is given by [15]
(25)
with n, E,, fizt, &, K, 6, and 6, being material parameters. 2.‘3 Dynamic ~~~i~ibriurnequations for the tube wall A membranic shell element of the tube wall (see Fig. 2) is subjected to inte~al forces per unit length (T”), due to the membrane stresses, and to external forces per unit area (q). The external forces are due to transmural pressure (P) and inertia forces. The equation of equilibrium may be written as 111, 121
(T’fiL, + W’JA,),,+ qA = 0 where T’ = a,h
(26)
Inflation waves in non-linear distensible tubes
109
Reference surface
Fig. 2. Membranic shell element of the tube wall.
and q = PN - pwha.
(28)
Here, pW is the wall density, h is the thickness, a is the acceleration components in the directions of the shell coordinates are
vector and its
(29) Inserting eqns (27j(29) into (26) yields a vectorial equation that includes derivatives of the base vectors. These are easily expressed in terms of the base vectors themselves as follows R,,, = C$R,, + BgBN. Hence, the scalar equilib~um equations in the directions of the shell coordinates, and N, are given respectively by
(30) R,,; R,,
(3lc) The assumption of circular cross sections implies that eqn (31b) is identically zero. 2.4 Equations ofmotion of the fluid Inflation waves in distensible tubes are associated with relatively high fluid velocities and consequently high Reynolds numbers. Also radial flows are negligible with respect to axial flows. For these reasons it will be assumed that the fluid is homogenous incompressible and inviscid and that the flow is quasi-one-dimensional. The Navier-Stokes equations for this case can be written as
s,,+ (SU),, u,,+
=0
(;,l+;+r >,j^ -0
(continuity)
(32)
(momentum).
(331
110
D. ELAD et at.
Here S is the internal cross-sectional area of the tube, U the fluid velocity (assumed to be uniform over the cross section), P the dynamic fluid pressure, pp the constant density of the fluid, 4” the length of the curved axis of the tube and Y’is potential of the body forces. Here, Y,= = G,(c) where G,(t) is the axial body force per unit mass. 2.5 Fluid-wall interaction Fluid inviscidity allows relative movement between the fluid and the wall in the longitudinal direction. However, continuity between the fluid and the wall demands that at any point along the tube axis the internal cross-sectional area of the tube be equal to the external cross-sectional area of the fluid. In addition the fluid pressure equals the pressure exerted on the internal surface of the tube.
3. NUMERICAL
SCHEME
3.1 Numerical method
The motion of the fluid and the wall are described by eqs (18) (31a), (31c), (32) and (33). These equations together with the initial and boundary condition are solved by numerical integration. The numerical solution is obtained by means of an explicit scheme based on MacCormack’s two-step method [19,20]. This method is of second order accuracy in both time and space and was originally developed to analyse shock waves in a compressible fluid. The coupling at the interface between the wall and the fluid is done by the “master and slave” procedure. The fluid is taken as the “master” and the lower mass wall is regarded as the “slave”. The interaction at the interface is solved by using two explicit phases: First the tube wall is regarded as rigid and fixed in space while the fluid body moves and assumes new cross-sectional areas. Then the fluid is kept still and its dynamic pressure is exerted on the shell which subsequently deforms to fit the new fluid geometry. This method deviates from the conventional methods for problems of gas-metal interaction, where equilibrium of the shell is considered first, after which the compressible fluid is adjusted to the new geometry [18]. The reason for not following the same procedure here is the absence of equation of state for the fluid, as in our case it is incompressible. However, under the present reverse procedure it is very difficult to account for the inertia forces on the wall material. Ignoring of the wall mass is expected to have a small effect on the solution since it is lower than that of the fluid within the tube. The equilibrium equations of the tube wall are described in the Lagrangian coordinates of the shell while those of the fluid are described in fixed Eulerian coordinates. The matching between variables described in the Lagrangian grid and those described in the Eulerian one is done by linear interpolation. To define the shell geometry at any time, under the given loading conditions, it is enough to define the lines of intersection of the reference surface and the x-y plane. Along one such line, a series of grid points are defined with mid-points halfway between them. Shell cells of small width are defined between two mid-points (see Fig. 3). This choice of representation enables natural placing of variables obtained by space di~erentiation and equilib~um equations attain physical meaning and will be of first order
Shell element
Fig. 3. Representation
,Grid point
of the numerical grid.
Inflation waves in non-linear distensible tubes
111
consistency. The following variables are defined at the grid points (i): x; y; r; r,,,; a,,; +,r; A,,; B,,; ri,; A,; &; cr2and P. The following variables are defined at the mid-points (i + 4): r,l; a; 4; +,,I; A,,; Bz; ri2; 1,; A,; c1 and h. The governing finite-difference equations are written for two steps according to MacCormack’s scheme. The fluid variables are defined in a fixed Eulerian grid indicated by j. The tube wall variables are defined in Lagrangian grid and mid points as described above and indicated by i and i + f respectively. The time is advanced by At in two steps. The first step yields approximate values and is indicated by n + 1. The second step which gives the final values is indicated by n + 1. Space step in the Eulerian coordinates is Ax and that in the Lagrangian is AZ. The governing equations in the two steps are First step Sin+’ =
(34)
(35)
(36)
(37)
Second Step [(Sqf
-
(38)
(39)
(40)
(41)
The value of variables defined at grid points but needed at mid points is calculated by linear interpolation and vice-versa. Variables that are not yet known at n + 1 will be taken from n. The numerical scheme is shown in the flow chart (Fig. 4). The solution is explicit and each calculation is carried out explicitly along the tube. 3.2 Initial deformation The tube may be initially loaded by static internal pressure or axial tension. The
112
D. ELAD rl al.
_
---
t Initial
t and boundary conditions
Initial
geotnetry
t
t
Calculation of static deformation of the tube
Calculation of At (Eq. 45)
t
Evaluation of geometric parameters: I h ,*,+.A ,L' UP ' Trt ' RP
Load condition: Pressure step at one end
(Lqs. 3,6,11.13.15)
(Eq. 45)
r,
c r---
STEP 11 ---- .- -
II
,#
IIt1 and r
Calculation of $+'
j
of the fluid at Eulerian grid
i ,
-
points (Eqs. 38.39)
1
-
-:
t
I
/
x and rj
Calculation of $7
1
/ 1
STEP I
I
of the fluid at Eulerian grid , 1 points (Eqs. 34.35)
1 I
t
/
x Determination of r,
'
of the
tube at Lagrangfan grid
tube at lagrangian grid
points by interpolation
,
I
/-
in+l Calculation of (oL)i+,
1
(Eq. 41)
Evaluation of Itcratlon
“+I (a~)~+,
by
(Es. 25)
,
I
I I
I
Rtermination
geometry of
of the new
I ’
the tube
of
p:+'
Lktennination of
,
x Pj
of the
fluid at Eulerian grid points
e
I Determin8tion
I
of the
I /
Fig. 4. Ftow chart of numerical solution.
I
Inflation waves in non-linear distensible tubes
113
calculation of the initial geometry and stresses in the tube is performed by using the computer code in a quasi-static mode, by letting the load attain its value ass~~toti~alIy. 3.3 Stability consideration The governing equation has both hyperbolic and parabolic aspects because the tube wall is made of a viscoelastic material. However, it was shown by Kivity (7, 143 that the stability criterion associates with the parabolic form is more restrictive. An approximate parabolic stability criterion is developed as follows: The Lagrangian grid of the shell is regarded as the Eulerian grid of the fluid. Also the inertia forces on a shell element are neglected and the tube is assumed to be straight. These approximations are introduced to the equilibrium eqn (31a) and to the Auid eqns (32) and (33) to yield
where Q is the discharge, W(Q; f&) is a function of lower order derivatives with respect to z and p is (431 The form of eqn (42) is similar to the non-linear heat equation. Richtmeyer and Morton [21] have shown that the lower order derivatives, W(Q; Q,& do not influence the stability of the solution. Stability will be insured if
Hence, the time step must be (45) where STAB is a stability factor (I 1.O) used in the computer code. 4. RESULTS
AND DISCUSSION
We shall #us&ate the predictive capacity of the present model by solving the problem of a pressure step applied at one end of a straight circular tube while there is a continuous supply of fluid thereafter. A wave front is thus developed moving away from that end. At some point along the tube, the wave reaches a steady sta.te and the velocity of the wave front as well as its profile become stationary. We shall analyze this stationary wave and evaluate the effect of various parameters on the solution. As this research was motivated by problems in the circulatory system the present illustration will be carried out with parameters found in blood vessels and under initial loading conditions which are in the physiological range.
The tube is assumed to be semi-jn~nite and cylindrical. The unstressed radius and thickness ratio were taken as r, = 1.25 cm,
h/r = 0.08.
The material parameters: n, E,, &,, a,, and ICwere determined from experimental data of biaxial tests on the carotid artery [22] by multiple regression (151. These parameters are: n = 3, E, = 0.2230 x IO6 dyn/cm2; /I$, = 5.5, j& = 34.5, IC= 9.0. To compare between orthotropic and isotropic materials we assumed the parameters for the later case to be n = 3, E, = 1.226 x IO6dyn/cm*, &, = #& = K = 1 .O. ES Vol. 22. No. 2-B
114
D. ELAD
ef nl.
The viscosity parameter 6, in eqn (18) was taken as b, = 0.16 set [12]. It was assumed that the wall density equals that of the fluid, namely pf = pw= 1.05 gr/cm3. 4.2 Numerical parameters It was found that stability can be obtained with a cell length (Fig. 3) of 1 cm and the stability factor STAB = 0.8. Increasing the mesh density or decreasing the time step considerably increased the calculation time with no significant change in results. 4.3 Initial conditions The fluid in the tube is assumed to be initially at rest hence, U,”= 0. The initial internal pressure in the tube was taken in some cases as zero (P, = 0) and in other cases as the mean physiologic blood pressure, P,, = 100 mmHg (100 mmHg = 0.133 x lo6 dyn/cm2). The initial static deformation and wall stresses were calculated when non-zero initial pressure had been assumed. 4.4 Boundary conditions
One open end of the tube (x = 0) is restrained with respect to axial motion but free to change its diameter. The other open end (x = L) is free. A pressure step of AP = 150 mmHg is applied at x = 0. Its time course is given by
P(t) = P, + AP sin’ P(t) = PO+ AP
OstsT t>t
(46)
where T is the rise time of the pressure, taken as 20 msec. We assume that Bernoulli’s law holds hence P(t) = P + ;pJJ2
(47)
4.5 Results
The evaluated variables P, U, A,, ,I,, (T,, 02, ii,, & are presented graphically along the tube axis at 10 msec intervals for the case of stationary wave, in an orthotropic cylindrial tube with P,, = 0 (Fig. 5a). The same results but with P,, = 100 mmHg are shown in Fig. 5b. It can be seen that the wave front, in the later case is more gradual, the wave speed is higher and the fluid velocity is lower. This is due to the increased stiffness of the non-linear wall material. It can also be seen that the longitudinal stretch (1,) decreases as the wave front passes while the corresponding stress (r, increases. In order to follow the transition from increasing i,, with the wave front, at PO= 0 to decreasing 1, at P,, = 100 mmHg a series of runs were conducted at P,, = 20 and 50 mmHg with the same Ap = 150 mmHg. These are shown in Fig. 6. It is seen that the trend of decreasing ;1, at the wave front starts at P,, = 50 mmHg. This can be explained in the following way: From the numerical results (Fig. 6) it is seen that the increase of ~7~at the wave front is higher than the corresponding increase of cl. The nonlinear biaxial orthotropic stress-strain relation can be represented by equi-stress diagram of cr, and (TV in the A,, AZplane (Fig. 7). If at the range of high I, and 4 one follows a path of large step of rr2 and smaller one of o,, it becomes apparent that L, may decrease while A, increases. A stationary wave propagation in an isotropic tube was also simulated. Results are shown in Fig. 8. The overall stiffness of the isotropic tube in this case was higher than that of the orthotropic one. This affected the wave speed which was higher. It is also seen that the longitudinal stretch (A,) at the leading edge of the wave front decreases in the case of PO= 0, but increases at higher pressures. This can also be explained in a similar manner, by analyzing the isotropic equi-stress diagram.
0.28 0.13
1.699 1.b82
P
D. ELAD
116
et al.
Al 0 Fig. 7. Equi-stress
lines of u, and o2 in the i,, ;1> plane.
“.
5. 0. -3.
7.
2 0.
a
Fig.
8. Propagation
b
in an isotropic cylindrical tube. (a) P,, = 0; (17) of a stationary wave P,, = 100 mmHg (unstressed length of tube L = 70 cm).
Inflation waves in non-linear 4.6
distensible
117
tubes
Numerical accuracy
The accuracy of the numerical solution is verified in two ways: First the wave velocity, calculated numerically, is compared with analytic results obtained by assuming steady state. In this method the wave speed is calculated from the known fluid velocity and the tube geometry far enough from the wave front [23]. Comparison of results obtained with PO= 0 after 90 msec shows that the speed obtained numerically is 572 cm/set while that obtained analytically is 581 cm/set. The deviation is only 1.5% even though complete steady state has not been reached at that time. Another check is done by calculating the balance of mass, momentum and energy of the whole system, fluid and wall. Calculations show that deviation in mass is about O.l%, that of the momentum is 1.3% and that of energy is 2%. These results prove that the numerical solution is adequate. 4.7 Influence of viscoelastic parameters To check the influence of the viscoelastic parameter 6, in eqn (25) solutions for stationary waves at P,, = 0 were made with various values of bJO.08 and 0.32 set, Fig. 9) and compared to the case of b, = 0.16 set (Fig. 5a). It can be seen that an increase in the material viscosity increases the wave speed while the slope of the wave profile becomes more moderate. The stresses decrease because the wall stiffness is increased. Also there is a noticeable increase in computational time because of the reduction in the time step (dt), (eqn 45). 5. CONCLUSIONS
The theory and numerical results presented in the preceding sections show that large amplitude wave phenomena in distensible tubes subjected to loading conditions of internal pressure and axial tension can be analysed by using non-linear shell theory. This approach has many advantages over the methods in which the constitutive relations for the wall material are introduced through the “tube law” which relates local cross sectional areas
b Fig. 9. Influence
of viscoelastic parameter (b) on stationary wave propagation in an orthotropic cylindrical tube with P,, = 0. (a) b = 0.08 set; (b) b = 0.32 sec.
D. ELAD et al.
I18
of the tube to the transmural pressure. This is especially so when the tube material is incompressible and orthotropic. The numerical scheme developed in the present work seems adequate for solving problems of large waves in distensible tubes since the error is within 1 to 2%. Although particular constitutive relations [t 51were chosen to represent the material in the present work, the method is of general validity. Any constitutive relations, suitable for representing a particular wall material, may be used provided that the strain-stress equations can be easily inverted. authors wish to express their gratitude to Frof, A, Libai, Aeronautical Engineering. for his generous help in applying the theory of thin shell in this work.
~r~~~~~~~~~~~~~~--~e
Techion,
REFERENCES [I] R. M. GREENDYKE, J. Am. Med. Assoc. 195, 527 (1966). [2] G. STAFFORD, AUG. N. Z. J. Surg. 47, 175, (1977). [3] C. W. AKINS, M. J. BUCKLEY, W. DAGGETT, J. B. MGLLDUFF and W. G. GERALD, Ann. Thor. Surg. 31, 30.5 (1981). [4] J. LURENZ and H. ZELLER, I?% Cu& es Pressure Surgeq Bedford, En&zrzd, paper 35 (1972). [5] [6] [7] [8] [9] [IO]
G. RUDINGER,
.7. Appi. Me&. 37, 34 (1970).
M. ANLIKER, R. L. ROCKWELL and E. OGDEN, S. Appl. Math. Phys. (ZAMP) 22,217 and 563 (1971). Y. KIVITY and R. COLLINS, J. Bjomechunics 7, 67 (1974). A. H. SHAPIRO, J. Biomech. EnKny 99, 126 (1977). S. J. COWLEY, J. F&d Mech. 116, 459 (1982). Y. C. FUNG, S&d Mechaprics, p. 4%. Prentice-Hal& New York (1965). 11I] A. E. GREEN and W. ZERNA, ~~e~~~~~~~~ Elasticity, OxFord University Press, London (1954). [I4 A. LIBAI, Private ~o~~~~j~t~~~. [13] R. COLLINS, and W. C. L. HUIJ,J. Bomechanics 5, 333 (1972). [14] Y. KIVITY and R. COLLINS, Lecture Notes in Camp. Sci. 11, 213 (1974). [l5] A. FOUX, D. ELAD, Y. LANIR and Y. KIVITY, To be published. [I61 B. R. SETH, Second order effects in elasticity. Plasficity and Fluid Dynamics (Edited by M Reiner and D. Abir), p. 162. McMillan, New York (1964). [I7] S. G. LEKHNITSKI, Thwry qf Elasticity Qf A~~sufrop~c Elastic Body. Holden-Day, San Francisco (1963). [lg] Y. KIVITY and D. TZUR, Proc. 4th Int. Symp. B&tics, Mwttrey, ~~~rnj~ f1978). 1191 R. W. MACCORMACK, AIAA paper No. 69-354 (1949). [20] R. W. MACCQRMACK, Lecture Notes in Phys. 8, X51 (1970). [21] R. D. RICHTMYER and K. W. MORTON, Difference Metho& for Initial-Value Problems, Sect. 12.7. Interscience, New York (1967). [22f R. H. COX, J. Appi. Physiol. 36, 381 (1975). [23] Y. KIVITY and R. COLLINS, Arch. Mech. (Warsaw) 26, 921 (1974). (Rece&w&25 February 1983)
INFLATION
2, pp.
101-118,
1984
0020-7225/84 Pergamon
WAVES IN NON-LINEAR
DISTENSIBLE
$3.00 + .OO Press Ltd.
TUBESt
D. ELAD, A. FOUX, Y. KIVITY and Y. LANIR Department of Biomedical Engineering, Technion-Israel Institute of Technology, Haifa 32000, Israel Abstract-A mathematical model for large amplitude wave propagation in a thin walled distensible tube is developed. The tube wall is considered as a membranic shell made of an incompressible, non-linear viscoelastic material with cylindrical orthotropy. The fluid is regarded as incompressible and inviscid and the flow is quasi-one-dimensional. The case of a pressure step applied at one end of a uniform straight tube is solved as an example. The system of partial differential equations, describing the motions of the fluid and the wall, are integrated numerically by using a two-step explicit scheme. Flow and deformation variables as well as the wave velocity are determined in time and space.
1. INTRODUCTION
rupture of the aorta in vehicular collision, which is responsible for significant percentage of traffic fatalities[l-31 was the motivation to this work. The problem is one of dynamics of fluid-filled distensible tubes. The object of this work is to study one phase of the general problem, namely the wave phenomena. Waves in distensible tubes filled with incompressible fluid may be sustained by inflation of the distensible tube wall and are sometimes called inflation waves. Such waves are developed through the interaction of the fluid and the tube wall and may be of large deformations. This interaction depends strongly upon the mechanical properties of the wall material as well as on the mechanical constraints imposed on the tube by its surrounding. There have been many studies of fluid flow through distensible tubes in recent years. Many of them were motivated by problems of blood flow in the vascular system. An often used method for analysis of such problems is based on the assumption that the fluid-wall interaction can be described by the “tube law” relating the local cross-sectional area of the tube to the local transmural pressure and to time[4-91. Beyond first approximation, it is unlikely that a simple tube law can account for the biaxial behaviour of a material which may be non-linear viscoelastic and orthotropic. In the present work, it is intended to present a method for explicit solution of the problem of wave propagation in fluid-filled distensible tubes, with particular emphasis on local stresses and stretches. An attempt is made to formulate the problem in a general frame through the use of non-linear shell theory. This allows for such factors as tapering and curvature of the tube as well as different boundary and loading conditions to be considered. For illustration purpose a specific solution is presented for the case of a cylindrical straight tube subjected to a pressure step at one end. TRAUMATIC
2. ANALYSIS
2.1 Tube geometry and deformation The tube, regarded as a thin membranic shell, is defined by a reference surface. The inner surface of the tube is conveniently chosen to be the reference surface. As it is also the outer surface of the fluid, the interaction will be defined on it. In order to define this reference surface a fixed Cartesian frame x, y, z is chosen in such a way that, under loading of internal pressure and axial tension, the curved axis of the tube is assumed to remain in the x-y plane (see Fig. 1). The radius vector to a point of the tube axis is
~(5, t) = ~(5, tY + y(5, t)j tThis research was supported by a grant from the United States-Israel (BSF), Jerusalem, Israel. 103
(1) Binational Science Foundation
D. ELAD et al.
104
where 5 is a non-dimensional Lagrangian parameter along the tube axis, t is the time and i, j are unit vectors in the x and y directions respectively. Real lengths along the tube axis, denoted by Z, are given by AZ = aA<
(2)
u2 = (x,J2 + 01,J2
(3)
where TVis a scaling variable given by
The inner surface of the she11is formed by closed planar curves C, enclosing the axis. A point on the curve C is described by the radius vector r(& 0, t) emanating from the point t on the axis, The angle 8, defining the direction of this vector in the plane of C, is another Lagrangian parameter. The following is assumed with respect to the deformation of the tube: (a) The curve C remains planar during defo~ation and the normal to the plane enclosed by C at the intersection with the tube axis, remains tangent to the axis. (b) The cross-sectional curves C are circular before and during deformation, hence, rrO= 0. (c) The deformations are in the principle directions without shear, hence there is no torsion of the tube and circumferential displa~ment are negligible.
The first assumption is analogous to Kirchhoff s hypothesis for thin plates [lo]. The second assumption ensures orthogonality of the shell coordinates, It is valid in straight tubes and is a first approximation in curved tubes. Similar considerations and the loading conditions lead to the third assumption. Under these assumptions the radius vector r at any time is given by r(r, B, t) = Ir(& 8, t)l[ - cos 8n + sin Ok]
(4)
where k is a unit vector in the z direction. The variable unit vectors 1, II, k (Fig. 1) represent a local orthogonal frame of reference that move with the curved axis as the tube deforms. They may be conventionally used to calculate the covariant base vectors at a point on the shell. The vectors 1 and II are unit vectors in the x-y plane, tangent and perpendicular to the curved axis respectively. The reference surface of the tube shell at any given time is defined in the fixed coordinate system by the radius vector to a point on it, namely
This vector enables to calculate the deformations and can be used in the equilibrium equations of the shell [I 1, 121.The parameters r and 8 may serve as a system of curvilinear coordinates on the shell surface (see Fig. 1) and are variably named “shell coordinates”, “material coordinates” or “Lagrangian coordinates”. For simplicity of notation the time variable t will not be indicated explicitly in all cases. The angle d, (Fig. 1) in the x-y plane is the angle between the tangent to the axis (defined by 1) and the x-axis, hence:
cosc#=;x,”
I
1
sin &,= ;y,(.
The unit vectors I and n are given by 1 = p,s = x,,i + y,d = cos 4i + sin (bj (7) n = -sin #i + cos f$j.
fnfiatiou
waves in
non-linear distensible tubes
in the
x-y
105
plane
Parallel to the
Diameter in the
Of the tube x-y
plane
4
n
k
fixed system
- unit vectors along
local Cartesian coordinates .on the curved axis R,{ R,,N-
Vectors alOng
curvilinear Cartesian coordinates on the shell
Fig. 1. The tube model and
It can be shown that the first derivatives
itsreference
systems.
of 1 and n are
These variables simplify the mathematical representation of the covariant base vectors. The radius vector to a point on the shell surface, after imerting eqn (4) into eqn (5) iS
R(t, 8, t) = p([, 1) - r cos On -t r sin Ok. The covariant
(9)
base vectors will now be
R, = R, = (a + r& cm S)l - (r,ccm f3)n+ fr,esinB)k Rr2 = R,.
= (r sin O)n +
The notation of the shell coot4inates f , 2 respectively.
(r cm O)k.
<,19 will be replaced in the following
equations
by
106
D.
F&AD
et
al.
The covariant base vectors are used to calculate other geometrical components af the metric tensor AaB will be
A,, = R,2- R,2 =
r2
A,, = R,, qR,2=
T,~T,~ =
A = AI,A22-
The base vectors are ~r~ndic~~ar
variables. The
0
(11)
A:, = r”[(a + F#,j cos ey+ (r*,)“l. to each ather since AI2 = 0,
The unit vector normal to the surface is
- cos B(a + r&l
cas @n-I- sin@(a+ r& cos @}k].
w
Hence, the vectors R,,; R,,; N form a local arthogonaf system that move with the shell surface. From N and the derivatives of the base vectors (R,aB), the components of the curvature tensor can be calculated. Hence,
thus,
Inflation waves in non-linear distensible tubes
The components
of length elements along the shell coordinates
107
are given by [1 1, 121
4 =&T& As,
For a given deformation be written as
(16)
= ,,&A&
that occurred during a time period t, the extension ratios may
(As,),
A1=o,= (17)
(AsJr
12=(d=
II--(&),
G422h
= 0’
2.2 Constitutive relations Inflation waves are associated with large deformations of the tube wall, hence constitutive relations of the wall material must be valid for finite deformations. The wall material is considered as an incompressible, non-linear viscoelastic material with cyiindrical orthotropy. To facilitate the numerical calculations, it is necessary to express the viscoelastic response in a very simple mathematical form which accounts for the effect of time. We chose the approach of Collins and Kivity [13, 141 which state that:
(a
o,(t) = D&l,, A,)
not summed).
(18)
Here, D,--is the quasi-static non-linear, biaxial elastic response and b&/A, is a uniaxial rate dependent function where & is the stretch rate and b, is a material parameter. Although eqn (18) ignores the fading memory aspect of the viscoelastic behavior it does account for the effect of rate. This may be justified as an ad hoc constitutive modeling for analysis of a single abrupt event (a large amplitude inflation wave) in which the deformation history plays a small role, but the effect of stretching-rate may be significant. The case under consideration is one where material axes coincide with geometric axes. With load conditions of internal pressure and axial tension the case may be reduced to one of biaxial normal stresses if the small radial stress is neglected. A model for the non-linear elastic response (D,) of such material was proposed by Foux et al. [15], and adopted in this work. Non-linearity is introduced, in this model, through the definitions of generalized non-linear strain and stress in such a way that the stress-strain relations are quasi-linear. This enables inversion of the stress-strain equations, an important feature in the solution of mixed boundary value problems. The rational and justifications of the following choice for biaxial nonlinear orthotropic stress-strain law is given in Foux et al. [15]. Here we shall outline the major points and results. The strain is defined by the generalized strain measure as proposed by Seth [ 161
c,+n:1)
(19)
where n is the non-linearity material parameter. The associated generalized stress is obtained by differentiation of the strain energy density function with respect to the generalized strain. The generalized stress is given by s, ==$ 01 where ca is the true (Eulerian) stress.
(a not summed)
(20)
108
D. ELAD er al.
The strain-stress relations are expressed, in analogy with the generalized Hooke’s law for an orthotropic material [17], as
(21)
where E8 are material parameters analogous to Young’s moduli and vE6are lateral-strainparameters analogous to Poisson’s ratios. Following Foux et al. [IS] we will take E, to be constants. The lateral-strain-parameters must therefore be strain dependent to allow for material incompressibility which requires that L&3 = 1.
(22)
Equation (22) imposes a restriction on the strain dependent which are thus functionally related. From symmetry we have: V21 ---_
Vtz
Et
-5
lateral-strain-parameters
(231
It is also assumed 1151that V3t
=
(A;&)-”
-w2t
V32
=
(A~Jp,)-(’
-831)Vt2
(24) where pPla= E,/EB and K is another material parameter. The value of the lateral strain parameters may thus be calculated for a given state of deformation from eqs (19)-(24). The final form of the constitutive equations expressed in terms of the extensions and their time derivatives is given by [15]
(25)
with n, E,, fizt, &, K, 6, and 6, being material parameters. 2.‘3 Dynamic ~~~i~ibriurnequations for the tube wall A membranic shell element of the tube wall (see Fig. 2) is subjected to inte~al forces per unit length (T”), due to the membrane stresses, and to external forces per unit area (q). The external forces are due to transmural pressure (P) and inertia forces. The equation of equilibrium may be written as 111, 121
(T’fiL, + W’JA,),,+ qA = 0 where T’ = a,h
(26)
Inflation waves in non-linear distensible tubes
109
Reference surface
Fig. 2. Membranic shell element of the tube wall.
and q = PN - pwha.
(28)
Here, pW is the wall density, h is the thickness, a is the acceleration components in the directions of the shell coordinates are
vector and its
(29) Inserting eqns (27j(29) into (26) yields a vectorial equation that includes derivatives of the base vectors. These are easily expressed in terms of the base vectors themselves as follows R,,, = C$R,, + BgBN. Hence, the scalar equilib~um equations in the directions of the shell coordinates, and N, are given respectively by
(30) R,,; R,,
(3lc) The assumption of circular cross sections implies that eqn (31b) is identically zero. 2.4 Equations ofmotion of the fluid Inflation waves in distensible tubes are associated with relatively high fluid velocities and consequently high Reynolds numbers. Also radial flows are negligible with respect to axial flows. For these reasons it will be assumed that the fluid is homogenous incompressible and inviscid and that the flow is quasi-one-dimensional. The Navier-Stokes equations for this case can be written as
s,,+ (SU),, u,,+
=0
(;,l+;+r >,j^ -0
(continuity)
(32)
(momentum).
(331
110
D. ELAD et at.
Here S is the internal cross-sectional area of the tube, U the fluid velocity (assumed to be uniform over the cross section), P the dynamic fluid pressure, pp the constant density of the fluid, 4” the length of the curved axis of the tube and Y’is potential of the body forces. Here, Y,= = G,(c) where G,(t) is the axial body force per unit mass. 2.5 Fluid-wall interaction Fluid inviscidity allows relative movement between the fluid and the wall in the longitudinal direction. However, continuity between the fluid and the wall demands that at any point along the tube axis the internal cross-sectional area of the tube be equal to the external cross-sectional area of the fluid. In addition the fluid pressure equals the pressure exerted on the internal surface of the tube.
3. NUMERICAL
SCHEME
3.1 Numerical method
The motion of the fluid and the wall are described by eqs (18) (31a), (31c), (32) and (33). These equations together with the initial and boundary condition are solved by numerical integration. The numerical solution is obtained by means of an explicit scheme based on MacCormack’s two-step method [19,20]. This method is of second order accuracy in both time and space and was originally developed to analyse shock waves in a compressible fluid. The coupling at the interface between the wall and the fluid is done by the “master and slave” procedure. The fluid is taken as the “master” and the lower mass wall is regarded as the “slave”. The interaction at the interface is solved by using two explicit phases: First the tube wall is regarded as rigid and fixed in space while the fluid body moves and assumes new cross-sectional areas. Then the fluid is kept still and its dynamic pressure is exerted on the shell which subsequently deforms to fit the new fluid geometry. This method deviates from the conventional methods for problems of gas-metal interaction, where equilibrium of the shell is considered first, after which the compressible fluid is adjusted to the new geometry [18]. The reason for not following the same procedure here is the absence of equation of state for the fluid, as in our case it is incompressible. However, under the present reverse procedure it is very difficult to account for the inertia forces on the wall material. Ignoring of the wall mass is expected to have a small effect on the solution since it is lower than that of the fluid within the tube. The equilibrium equations of the tube wall are described in the Lagrangian coordinates of the shell while those of the fluid are described in fixed Eulerian coordinates. The matching between variables described in the Lagrangian grid and those described in the Eulerian one is done by linear interpolation. To define the shell geometry at any time, under the given loading conditions, it is enough to define the lines of intersection of the reference surface and the x-y plane. Along one such line, a series of grid points are defined with mid-points halfway between them. Shell cells of small width are defined between two mid-points (see Fig. 3). This choice of representation enables natural placing of variables obtained by space di~erentiation and equilib~um equations attain physical meaning and will be of first order
Shell element
Fig. 3. Representation
,Grid point
of the numerical grid.
Inflation waves in non-linear distensible tubes
111
consistency. The following variables are defined at the grid points (i): x; y; r; r,,,; a,,; +,r; A,,; B,,; ri,; A,; &; cr2and P. The following variables are defined at the mid-points (i + 4): r,l; a; 4; +,,I; A,,; Bz; ri2; 1,; A,; c1 and h. The governing finite-difference equations are written for two steps according to MacCormack’s scheme. The fluid variables are defined in a fixed Eulerian grid indicated by j. The tube wall variables are defined in Lagrangian grid and mid points as described above and indicated by i and i + f respectively. The time is advanced by At in two steps. The first step yields approximate values and is indicated by n + 1. The second step which gives the final values is indicated by n + 1. Space step in the Eulerian coordinates is Ax and that in the Lagrangian is AZ. The governing equations in the two steps are First step Sin+’ =
(34)
(35)
(36)
(37)
Second Step [(Sqf
-
(38)
(39)
(40)
(41)
The value of variables defined at grid points but needed at mid points is calculated by linear interpolation and vice-versa. Variables that are not yet known at n + 1 will be taken from n. The numerical scheme is shown in the flow chart (Fig. 4). The solution is explicit and each calculation is carried out explicitly along the tube. 3.2 Initial deformation The tube may be initially loaded by static internal pressure or axial tension. The
112
D. ELAD rl al.
_
---
t Initial
t and boundary conditions
Initial
geotnetry
t
t
Calculation of static deformation of the tube
Calculation of At (Eq. 45)
t
Evaluation of geometric parameters: I h ,*,+.A ,L' UP ' Trt ' RP
Load condition: Pressure step at one end
(Lqs. 3,6,11.13.15)
(Eq. 45)
r,
c r---
STEP 11 ---- .- -
II
,#
IIt1 and r
Calculation of $+'
j
of the fluid at Eulerian grid
i ,
-
points (Eqs. 38.39)
1
-
-:
t
I
/
x and rj
Calculation of $7
1
/ 1
STEP I
I
of the fluid at Eulerian grid , 1 points (Eqs. 34.35)
1 I
t
/
x Determination of r,
'
of the
tube at Lagrangfan grid
tube at lagrangian grid
points by interpolation
,
I
/-
in+l Calculation of (oL)i+,
1
(Eq. 41)
Evaluation of Itcratlon
“+I (a~)~+,
by
(Es. 25)
,
I
I I
I
Rtermination
geometry of
of the new
I ’
the tube
of
p:+'
Lktennination of
,
x Pj
of the
fluid at Eulerian grid points
e
I Determin8tion
I
of the
I /
Fig. 4. Ftow chart of numerical solution.
I
Inflation waves in non-linear distensible tubes
113
calculation of the initial geometry and stresses in the tube is performed by using the computer code in a quasi-static mode, by letting the load attain its value ass~~toti~alIy. 3.3 Stability consideration The governing equation has both hyperbolic and parabolic aspects because the tube wall is made of a viscoelastic material. However, it was shown by Kivity (7, 143 that the stability criterion associates with the parabolic form is more restrictive. An approximate parabolic stability criterion is developed as follows: The Lagrangian grid of the shell is regarded as the Eulerian grid of the fluid. Also the inertia forces on a shell element are neglected and the tube is assumed to be straight. These approximations are introduced to the equilibrium eqn (31a) and to the Auid eqns (32) and (33) to yield
where Q is the discharge, W(Q; f&) is a function of lower order derivatives with respect to z and p is (431 The form of eqn (42) is similar to the non-linear heat equation. Richtmeyer and Morton [21] have shown that the lower order derivatives, W(Q; Q,& do not influence the stability of the solution. Stability will be insured if
Hence, the time step must be (45) where STAB is a stability factor (I 1.O) used in the computer code. 4. RESULTS
AND DISCUSSION
We shall #us&ate the predictive capacity of the present model by solving the problem of a pressure step applied at one end of a straight circular tube while there is a continuous supply of fluid thereafter. A wave front is thus developed moving away from that end. At some point along the tube, the wave reaches a steady sta.te and the velocity of the wave front as well as its profile become stationary. We shall analyze this stationary wave and evaluate the effect of various parameters on the solution. As this research was motivated by problems in the circulatory system the present illustration will be carried out with parameters found in blood vessels and under initial loading conditions which are in the physiological range.
The tube is assumed to be semi-jn~nite and cylindrical. The unstressed radius and thickness ratio were taken as r, = 1.25 cm,
h/r = 0.08.
The material parameters: n, E,, &,, a,, and ICwere determined from experimental data of biaxial tests on the carotid artery [22] by multiple regression (151. These parameters are: n = 3, E, = 0.2230 x IO6 dyn/cm2; /I$, = 5.5, j& = 34.5, IC= 9.0. To compare between orthotropic and isotropic materials we assumed the parameters for the later case to be n = 3, E, = 1.226 x IO6dyn/cm*, &, = #& = K = 1 .O. ES Vol. 22. No. 2-B
114
D. ELAD
ef nl.
The viscosity parameter 6, in eqn (18) was taken as b, = 0.16 set [12]. It was assumed that the wall density equals that of the fluid, namely pf = pw= 1.05 gr/cm3. 4.2 Numerical parameters It was found that stability can be obtained with a cell length (Fig. 3) of 1 cm and the stability factor STAB = 0.8. Increasing the mesh density or decreasing the time step considerably increased the calculation time with no significant change in results. 4.3 Initial conditions The fluid in the tube is assumed to be initially at rest hence, U,”= 0. The initial internal pressure in the tube was taken in some cases as zero (P, = 0) and in other cases as the mean physiologic blood pressure, P,, = 100 mmHg (100 mmHg = 0.133 x lo6 dyn/cm2). The initial static deformation and wall stresses were calculated when non-zero initial pressure had been assumed. 4.4 Boundary conditions
One open end of the tube (x = 0) is restrained with respect to axial motion but free to change its diameter. The other open end (x = L) is free. A pressure step of AP = 150 mmHg is applied at x = 0. Its time course is given by
P(t) = P, + AP sin’ P(t) = PO+ AP
OstsT t>t
(46)
where T is the rise time of the pressure, taken as 20 msec. We assume that Bernoulli’s law holds hence P(t) = P + ;pJJ2
(47)
4.5 Results
The evaluated variables P, U, A,, ,I,, (T,, 02, ii,, & are presented graphically along the tube axis at 10 msec intervals for the case of stationary wave, in an orthotropic cylindrial tube with P,, = 0 (Fig. 5a). The same results but with P,, = 100 mmHg are shown in Fig. 5b. It can be seen that the wave front, in the later case is more gradual, the wave speed is higher and the fluid velocity is lower. This is due to the increased stiffness of the non-linear wall material. It can also be seen that the longitudinal stretch (1,) decreases as the wave front passes while the corresponding stress (r, increases. In order to follow the transition from increasing i,, with the wave front, at PO= 0 to decreasing 1, at P,, = 100 mmHg a series of runs were conducted at P,, = 20 and 50 mmHg with the same Ap = 150 mmHg. These are shown in Fig. 6. It is seen that the trend of decreasing ;1, at the wave front starts at P,, = 50 mmHg. This can be explained in the following way: From the numerical results (Fig. 6) it is seen that the increase of ~7~at the wave front is higher than the corresponding increase of cl. The nonlinear biaxial orthotropic stress-strain relation can be represented by equi-stress diagram of cr, and (TV in the A,, AZplane (Fig. 7). If at the range of high I, and 4 one follows a path of large step of rr2 and smaller one of o,, it becomes apparent that L, may decrease while A, increases. A stationary wave propagation in an isotropic tube was also simulated. Results are shown in Fig. 8. The overall stiffness of the isotropic tube in this case was higher than that of the orthotropic one. This affected the wave speed which was higher. It is also seen that the longitudinal stretch (A,) at the leading edge of the wave front decreases in the case of PO= 0, but increases at higher pressures. This can also be explained in a similar manner, by analyzing the isotropic equi-stress diagram.
0.28 0.13
1.699 1.b82
P
D. ELAD
116
et al.
Al 0 Fig. 7. Equi-stress
lines of u, and o2 in the i,, ;1> plane.
“.
5. 0. -3.
7.
2 0.
a
Fig.
8. Propagation
b
in an isotropic cylindrical tube. (a) P,, = 0; (17) of a stationary wave P,, = 100 mmHg (unstressed length of tube L = 70 cm).
Inflation waves in non-linear 4.6
distensible
117
tubes
Numerical accuracy
The accuracy of the numerical solution is verified in two ways: First the wave velocity, calculated numerically, is compared with analytic results obtained by assuming steady state. In this method the wave speed is calculated from the known fluid velocity and the tube geometry far enough from the wave front [23]. Comparison of results obtained with PO= 0 after 90 msec shows that the speed obtained numerically is 572 cm/set while that obtained analytically is 581 cm/set. The deviation is only 1.5% even though complete steady state has not been reached at that time. Another check is done by calculating the balance of mass, momentum and energy of the whole system, fluid and wall. Calculations show that deviation in mass is about O.l%, that of the momentum is 1.3% and that of energy is 2%. These results prove that the numerical solution is adequate. 4.7 Influence of viscoelastic parameters To check the influence of the viscoelastic parameter 6, in eqn (25) solutions for stationary waves at P,, = 0 were made with various values of bJO.08 and 0.32 set, Fig. 9) and compared to the case of b, = 0.16 set (Fig. 5a). It can be seen that an increase in the material viscosity increases the wave speed while the slope of the wave profile becomes more moderate. The stresses decrease because the wall stiffness is increased. Also there is a noticeable increase in computational time because of the reduction in the time step (dt), (eqn 45). 5. CONCLUSIONS
The theory and numerical results presented in the preceding sections show that large amplitude wave phenomena in distensible tubes subjected to loading conditions of internal pressure and axial tension can be analysed by using non-linear shell theory. This approach has many advantages over the methods in which the constitutive relations for the wall material are introduced through the “tube law” which relates local cross sectional areas
b Fig. 9. Influence
of viscoelastic parameter (b) on stationary wave propagation in an orthotropic cylindrical tube with P,, = 0. (a) b = 0.08 set; (b) b = 0.32 sec.
D. ELAD et al.
I18
of the tube to the transmural pressure. This is especially so when the tube material is incompressible and orthotropic. The numerical scheme developed in the present work seems adequate for solving problems of large waves in distensible tubes since the error is within 1 to 2%. Although particular constitutive relations [t 51were chosen to represent the material in the present work, the method is of general validity. Any constitutive relations, suitable for representing a particular wall material, may be used provided that the strain-stress equations can be easily inverted. authors wish to express their gratitude to Frof, A, Libai, Aeronautical Engineering. for his generous help in applying the theory of thin shell in this work.
~r~~~~~~~~~~~~~~--~e
Techion,
REFERENCES [I] R. M. GREENDYKE, J. Am. Med. Assoc. 195, 527 (1966). [2] G. STAFFORD, AUG. N. Z. J. Surg. 47, 175, (1977). [3] C. W. AKINS, M. J. BUCKLEY, W. DAGGETT, J. B. MGLLDUFF and W. G. GERALD, Ann. Thor. Surg. 31, 30.5 (1981). [4] J. LURENZ and H. ZELLER, I?% Cu& es Pressure Surgeq Bedford, En&zrzd, paper 35 (1972). [5] [6] [7] [8] [9] [IO]
G. RUDINGER,
.7. Appi. Me&. 37, 34 (1970).
M. ANLIKER, R. L. ROCKWELL and E. OGDEN, S. Appl. Math. Phys. (ZAMP) 22,217 and 563 (1971). Y. KIVITY and R. COLLINS, J. Bjomechunics 7, 67 (1974). A. H. SHAPIRO, J. Biomech. EnKny 99, 126 (1977). S. J. COWLEY, J. F&d Mech. 116, 459 (1982). Y. C. FUNG, S&d Mechaprics, p. 4%. Prentice-Hal& New York (1965). 11I] A. E. GREEN and W. ZERNA, ~~e~~~~~~~~ Elasticity, OxFord University Press, London (1954). [I4 A. LIBAI, Private ~o~~~~j~t~~~. [13] R. COLLINS, and W. C. L. HUIJ,J. Bomechanics 5, 333 (1972). [14] Y. KIVITY and R. COLLINS, Lecture Notes in Camp. Sci. 11, 213 (1974). [l5] A. FOUX, D. ELAD, Y. LANIR and Y. KIVITY, To be published. [I61 B. R. SETH, Second order effects in elasticity. Plasficity and Fluid Dynamics (Edited by M Reiner and D. Abir), p. 162. McMillan, New York (1964). [I7] S. G. LEKHNITSKI, Thwry qf Elasticity Qf A~~sufrop~c Elastic Body. Holden-Day, San Francisco (1963). [lg] Y. KIVITY and D. TZUR, Proc. 4th Int. Symp. B&tics, Mwttrey, ~~~rnj~ f1978). 1191 R. W. MACCORMACK, AIAA paper No. 69-354 (1949). [20] R. W. MACCQRMACK, Lecture Notes in Phys. 8, X51 (1970). [21] R. D. RICHTMYER and K. W. MORTON, Difference Metho& for Initial-Value Problems, Sect. 12.7. Interscience, New York (1967). [22f R. H. COX, J. Appi. Physiol. 36, 381 (1975). [23] Y. KIVITY and R. COLLINS, Arch. Mech. (Warsaw) 26, 921 (1974). (Rece&w&25 February 1983)