INELASTIC WAVE PROPAGATION
UNDER COMBINED STRESS STATES
By
CHARLES DANIEL MYERS
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA IN PARTIAL
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1973
TO PEGGY
2. CKh'NO WLEDC UENT S
I would like to thank Professor Martin A. Eisenberg, Chairman of
the Supervisory Committee, not only for his untiring efforts during
the development and preparation of the material contained in this
manuscript, but also for being a counselor, teacher, and friend during
both my undergraduate and graduate studies. I am also indebted to
Professors L. E. Malvern and E. K. Walsh for their helpful criticism
and encouragement during my doctoral studies. In addition, I would
like to express my appreciation to the other members of my Supervisory
Committee: Professors U. H. Kurzweg, C. A. Ross, and R. C. Fluck.
A special word of thanks is extended to Professor N. Cristescu
for his many helpful discussions during the development of this
dissertation.
I thank my wife, Peggy, for her encouragement, moral support, and
understanding during the course of my studies. I also thank Peggy for
typing and proofreading the drafts of this dissertation. I appreciate
the efforts of Mrs. Edna Larrick for the final typing of the manuscript
and Mrs. Helen Reed for the final preparation of figures.
I acknowledge financial support from the National Defense Education
Act, the National Science Foundation, and the University of Florida
which made myv studies possible.
I also acknowledge the Northeast Florida Regional Computing Center
for the use of its IBM 370 Model 165 digital computer without which
the scope of this work would have been greatly curtailed.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS . . . . ... ...... iii
LIST OF TABLES . . . . ... ... vi
LIST OF FIGURES . . . .. . . vii
ABSTRACT . . . . ... . ... ix
CHAPTER 1. INTRODUCTION . ... . . 1
CHAPTER 2. THEORETICAL DEVELOPMENT . . .. . 12
CHAPTER 3. D]IENSIONLESS EQUATIONS AND NUMERICAL PROCEDURES 29
3.1. Wave Speeds as a Function of the State of Stress .. 29
3.2. Characteristic Solution in Terms of
Dimensionless Variables . . . ... 46
3.3. Numerical Grid for Characteristic Solution . .. 51
3.4. Finite Difference Eouations . . ... 58
3.5. Solution to the Finite Difference Equations . .. 69
3.6. Calculation of the Strains .... . . . 77
CHAPTER 4. RESULTS AND DISCUSSION . . . ... 80
4.1. Introduction . . . ... .. 80
4.2. Grid Size Effects . ... . . .. 81
4.3. Effects of Radial Inertia . . ... .84
4.4. Effects of StrainRate Dependence . . .. 110
CHAPTER 5. SUIIARY . . . . ... . 118
APPENDIX A. CONSTITUTIVE EQUATIONS . . . ... 125
A.1. Comments on the Constitutive Equation . ... 125
A.2. Rate Independent Incremental Plasticity Theory ... 127
A.3. Rate Dependent Plasticity Theory . . .. 135
A.4. Dimensionless Expressions for the Functions
0(s,a) and .(s, . . . . 137
TABLE OF COIM'ENTS (Continued)
Page
APPENDIX B. CHARACTERISTICS AND EQUATIONS
ALONG THE CHARACTERISTICS . . ... 140
B.1. Equations for the Characteristics . . .. 140
B.2. Equations along the Characteristics . ... 141
B.3. Reducing Equations to Simpler Case . . .. 156
B.4. Uncoupled Waves . . . .... .. 159
B.5. Elastic Waves .... . . . . 162
APPENDIX C. PROGRAMS FOR DETERMINING THE PLASTIC
WAVE SPEEDS . . . ... ... 164
APPENDIX D. SOLUTION TO THE FINITE DIFFERENCE EQUATIONS
IN THE CHARACTERISTIC PLANE . . ... 171
D.1. Equations for Fully Coupled Waves . . .. 171
D.2. Equations for Uncoupled Waves . . ... 176
D.3. Solution at a Regular Grid Point
for Fully Coupled Waves . . . ... 180
D.4. Solution at a Regular Grid Point
for Uncoupled Waves ................. 183
D.5. Solution at a Boundary Point (X= 0)
for Fully Coupled Waves . . . ... 185
D.6. Solution at a Boundary Point (X= 0)
for Uncoupled Waves . . . ... 191
APPENDIX E. COMPUTER PROGRAM FOR CHARACTERISTIC PLANE
SOLUTION . . . ... . 194
E.I. General Description of the Program . . ... 194
E.2. Initial Conditions . . . . 196
E.3. Calculation of A . . . .... .. 201
E.4. Input Data . . . . . 202
E.5. Listing of the Program . . . ... 204
LIST OF REFERENCES . . . . . 231
BIOGRAPHICAL SKETCH . . . .... . 236
LIST OF TABLES
Table Page
1 Normalized Longitudinal Stress ) . . 36
s
2 Normalized Hoop Stress ( ............. 38
s
3 Normalized Shear Stress .) . . . 40
S
LIST OF FIGURES
Figure Page
2.1 Coordinate System for the ThinWalled Tube . ... 13
2.2 Stresses on an Element of the Tube . . ... 14
3.1 Yield Surface Representation in Spherical Coordinates 31
3.2 Plastic Wave Speeds as Functions of $ and y for
Poisson's Ratio of 0.30 . . . ... 41
3.3 Values of $ at v= 0 for which c =c =c ...... 45
f s 2
3.4 Numerical Grid in the Characteristic Plane . ... 53
3.5 Regular Element in Numerical Grid . . ... 54
3.6 Boundary Element in Numerical Grid . . ... 55
3.7 Location of the Characteristic Lines Passing Through P 57
3.8 Numerical Representation of the Characteristic Lines
in a Regular Element . . . . ... 59
3.9 Representation of the Characteristic Lines in
a Boundary Element ... . . . .. 60
4.1 Grid Size Effects on the Longitudinal Strain
at X = 3.75 . . . . ... . 83
4.2 Grid Size Effects on the Longitudinal Velocity
at X = 3.75 . . . . ... . 85
4.3 Grid Size Effects on the Stress Trajectories
at X = 3,75 . . . . ... . 86
4.4 Longitudinal Strain Versus Time at X = 3.75
for Data Set 1 . . . . ... . 88
4.5 Change in Shear Strain Versus Time at X = 3.75
for Data Set 1 . . . . . .89
LIST OF FIGURES (Continued)
Figure
4.6 Transverse Velocity Versus Time for Data Set 1
Without Radial Inertia . . . .
4.7 Longitudinal Velocity Versus Time for Data Set 1
Without Radial Inertia . . . .
4.8 Longitudinal Strain Versus Time for Data Set 1 .
4.9 Maximum Radial Velocity Versus X for Data Set 1
With Radial Inertia . . . .
4.10 Change in Shear Strain Versus Time for Data Set 1
Without Radial Inertia . . . .
4.11 Longitudinal Strain Versus X for Data Set 1 .
4.12 Stress Trajectories for Data Set 1
Without Radial Inertia . . . .
4.13 Strain Trajectories for Data Set 1
Without Radial Inertia . . . .
4.14 Shear Stress Versus Longitudinal Stress for Data S
With Radial Inertia . . . .
Page
. 93
. 94
. 96
. 97
. 98
 99
100
et 1
. 102
4.15 Stress Trajectories for Data Set 1 With Radial Inertia
4.16 Hoop Stress Versus Longitudinal Stress for Data Set 1
With Radial Inertia . . . . .
4.17 Stress Trajectories for Data Set 2
Without Radial Inertia . . . . .
4.18 Stress Trajectories for Data Set 2 With Radial Inertia
4.19 Hoop Stress Versus Longitudinal Stress for Data Set 2
With Radial Inertia . . . . .
4.20 Shear Strain Versus Time for Data Set 2 . .
4.21 Change in Longitudinal Strain Versus Time forData Set 2
4.22 Longitudinal Strain Versus Time for Data Set 3 . .
4.23 Stress Trajectory at X= 0 for Data Set 3 . .
4.24 Stress Trajectory at X= .25 for Data Set 3 . .
viii
104
105
107
108
109
111
112
114
116
117
Abstract of Dissertation Presented to the
Graduate Council of the University of Florida in Partial
Fulfillment of the Requirements for the Degree of Doctor of Philosophy
INELASTIC WAVE PROPAGATION UNDER
COMBINED STRESS STATES
By
Charles Daniel Myers
August, 1973
Chairman: Dr. M. A. Eisenberg
Major Department: Engineering Science, Mechanics
and Aerospace Engineering
The purpose of this dissertation was to investigate the effects
of radial inertia and material strainrate dependence on the propa
gation of inelastic waves of combined stress along a thinwalled tube.
A general quasilinear constitutive equation for multiaxial stress (and
strain) states was introduced. The equations of motion and the
straindisplacement equations, along with the constitutive equations,
were written to form a set of nine simultaneous hyperbolic, quasilinear,
partial differential equations. This set of equations was reduced to
a set of six equations which was then used to determine the expres
sions for the characteristic lines and the equations along the char
acteristic lines.
For combined torsional and longitudinal loading, two distinct
wave speeds were found. The values of these two wave speeds were
found as functions of the state of stress. Including radial inertia
effect in the formulation of the problem was shown to significantly
increase the wave speeds for a given stress state. Also certain
critical combinations of Poisson's ratio and the "effective tangent
mocdulus" caused the two wave speeds to be equal when the shear stress
vanished.
The equations for the characteristics and the equations along the
characteristics were written in terms of dimensionless variables.
These equations were then written as first order finite difference
equations. A computer code was written in the Fortran IV language,
and several problems were solved using an IBM 370 model 165 digital
computer. In order to obtain these solutions two particular forms of
the constitutive equation were used; one form represented a strain
rate independent material while the other form represented a strain
rate dependent material.
The strain at the impact end was considerably larger when radial
inertia effects were included than when radial inertia effects were
not included in the problem. However, radial inertia effects were
found to have little influence on the solution more than two diameters
from the impact end. The strain at the impact end was lowered by
including strainrate dependence of the material. For any particular
set of initial conditions and boundary conditions, the stress trajec
tories behaved in the same manner, at least qualitatively, whether or
not radial inertia effects or strainrate dependence were included.
The details of the stress trajectories were more complicated when
radial inertia effects were included since the trajectories were
threedimensional.
CHAPTER 1
INTRODUCTION
Stress wave propagation is the mechanism by which forces and
displacements are transmitted from one part of a structure to another.
Stress waves arise when a transient force is applied to a structure,
and they propagate through the structure reflecting (at least partially)
back into the structure whenever they encounter a boundary. After
several reflections the amplitude of the stress waves diminishes and
the structure reaches a state of equilibrium. In many engineering
problems the time required to reach equilibrium is very short, and for
practical purposes the structure can be assumed to reach equilibrium
instantly. Problems in which the forces are applied slowly or in which
the state of stress is required a long time after the forces are applied
are examples of instances when wave propagation effects may be neglected.
However, in many cases, the forces are applied rapidly (such as during
impact loading or explosive loading), and failure is most likely to
occur in the structure almost immediately after the application of these
forces. In these cases when it is necessary to determine the state of
stress during and immediately following the loading, wave propagation
effects may be significant and should be included in the analysis of
the problem. In order to understand the development of the theory of
stress wave propagation and the application of this theory to modern
engineering problems, it is instructive to review briefly the history
of wave propagation research.
The first serious attempt (at least in this century) to understand
nonlinear wave propagation in solids was made by Donnell (1930). In
this paper, Donnell used energy principles and impulsemomentum expres
sions to find the particle velocity and the elastic wave speed for
longitudinal waves. He also predicted that if a material with a bilin
ear stressstrain curve were impacted at the end by a stress above the
yield stress, two stress waves would propagate with distinct velocities.
However, after the publication of this paper interest in wave propaga
tion subsided until the early 1940's.
A more general theory of longitudinal stress wave propagation was
developed independently by Taylor (1940), von Karman (1942) and
Rakhmatulin (1945) by assuming that the material exhibited a nonlinear
stressstrain curve above the yield point. This stressstrain curve
was assumed to be independent of the rate of straining. Using this
theory the velocity of propagation of the longitudinal waves was found
to be given by
1 da
where c is the wave speed, p is the density of the material, C is the
da
stress, and e is the strain. Thus represents the slope of the stress
de
strain curve or the tangent modulus. This theory also considered the
stressstrain curve of the material which was obtained for the static
case to be valid in the dynamic case. With this assumption, the stress
and strain followed a unique functional relationship as long as no
unloading occurred. Because of this, the tangent modulus could be
written as a function of the stress (or strain) only, so that the
velocity of propagation then became a function of the level of stress
(or strain). This immediately led to the conclusion that a given
level of stress (or strain) propagated at a specific speed, and the
stress wave changed shape as it propagated along a prismatic bar for
stresses in the nonlinear region of the stressstrain curve. For a
bilinear stressstrain curve, the results of Donnell (1930) were again
predicted. However, these theories did not account for the lateral
inertia effects in the bar or the dependence of the stressstrain curve
on the rate of strain, and so more complex theories and constitutive
equations were proposed to account for these phenomena.
By the late 1940's many investigators including Davis (1938),
Manjoine (1940), and Clark and Wood (1950) had experimentally observed
the effect of the rate of strain on the stressstrain curve for several
materials. In order to incorporate this strainrate effect into the
constitutive equations used to study plastic wave propagation,
Sokolovsky (1948a, 1948b) and Malvern (1949, 1951a, 1951b) independ
ently introduced onedimensional constitutive equations in which the
stress was a function of the plastic strain and the plastic strain rate.
By selecting a particular form of this constitutive equation, Malvern
(1951a, 1951b) was able to obtain a numerical solution which predicted
several experimentally observed phenomena. However, his numerical solu
tion did not apparently predict a region of constant strain near the
impact end such as had been observed by Duwez and Clark (1947) and
others. This new strainrate dependent constitutive equation also
predicted that, if a bar were strained statically above the yield stress
and then impacted, the first increment of strain would propagate with
the elastic wave velocity and not the velocity given by the tangent
modulus.in the strainrate independent theory. Since this prediction
was quite different from that of the strainrate independent theory
several investigators tried to verify one or the other. Bell (1951)
published the results of his experiments with aluminum which showed
that, for a bar stressed above the yield point, the initial strain
pulse propagated with the elastic wave velocity. These results were
in accordance with the strainrate dependent model of Malvern (1951a,
1951b) as were the experimental results of Sternglass and Stuart (1953)
which were obtained using copper, Alter and Curtis (1956) which were
obtained using lead, Bell and Stein (1962) which were obtained using
aluminum, and Bianchi (1964) which were obtained using copper.
Encouraged by these experimental results, many investigators
continued the development of more general constitutive equations to
describe material behavior. Perzyna (1963) generalized the semilinear
constitutive equation of Malvern (1951a, 1951b) to multiaxial states of
stress. At about this same time Cristescu (1964) introduced full quasi
linear constitutive relations for a onedimensional problem. This
quasilinear equation was used immediately by Lubliner (1964) to show
that the strainrate independent constitutive equation of Taylor (1940),
von Karman (1942), and Rakhmatulin (1945), and the strainrate depend
ent constitutive equation of Malvern (1951a, 1951b) and Sokolovsky
(1948a, 1948b) were both special cases of this more general constitu
tive equation. Later Cristescu (1967a) gave a generalization for multi
dimensional stress states of the quasilinear constitutive equation as
well as an extensive summary of the developments in dynamic plasticity
until that time. Lindholm (1967) developed a constitutive equation
for combined stress states of aluminum which included strainrate
effects and temperature dependence. He also presented extensive data
for onedimensional loading and combined stress loading at several
strain rates and temperatures which were used in empirically determin
ing the constants used in his generalized constitutive equation.
While these more general constitutive equations were being developed,
it was shown by Malvern (1965), by Wood and Phillips (1967), and by
Efron and Malvern (1969) that the semilinear equation of Malvern"
(1951a, 1951b) did indeed predict a region of constant or nearly con
stant strain near the impact end if the solution was obtained long
enough after impact. Suliciu, Malvern, and Cristescu (1972) have shown
that a region of constant strain is not possible for the semilinear
constitutive equation but may be approached asymptotically. They have
also shown that a region of constant strain is possible when the quasi
linear constitutive equation is used. However, in the interpretation
of experimental results it has been difficult to differentiate between
a region of truly constant strain and a region in which the constant
strain is approached asymptotically.
The experiments of Sternglass and Stuart (1953), Alter and Curtis
(1956), and others were believed by many investigators to be proof of
the strainrate dependence of some materials. This led to the exten
sive development of constitutive equations just discussed. However,
other investigators sought to explain the experimentally observed phe
nomena by including radial inertia in the formulation of the wave prop
agation problem. Plass and Ripperger (1960) introduced radial inertia
effects into the problem of longitudinal impact and used the constitu
tive equation of Malvern (1951a, 1951b). In order to find a character
istic solution, all of the variables were averaged at each cross section
6
and these averaged variables were used. The :results of this work were
given by Tapley and Plass (1961) but were somewv.hat inconclusive. More
work including radial inertia effects was published by Hunter and
Johnson (1964), and a year later DeVault (1965) showed that, at least
qualitatively, many observations formerly attributed to a material
strainrate effect could be explained by including radial inertia effects
in the formulation of the problem of longitudinal impact of a bar.
Shea (1968) obtained good agreement between theory and experiment for the
propagation of longitudinal waves in a lead bar. He used the strain
rate dependent constitutive equation of Malvern (1951b) and the
"correction" for radial inertia proposed by DeVault (1965). Mok (1972)
used the same averaging technique for the variables as Plass and
Ripperger (1960) for the problem of longitudinal impact of a bar with
radial inertia effects included. He used the strainrate independent
constitutive equations and agreed in essence with DeVault (1965) that
radial inertia effects could explain, at least qualitatively, those
experimental results usually attributed to strainrate sensitivity
of the material. Since radial inertia is always present in an experi
ment using longitudinal impact it seemed that the only way to conclu
sively determine strainrate effects in a material would be to perform
the experiments using a torsional wave.
In an effort to determine the strainrate dependence of various
materials, several investigators have recently conducted theoretical
and experimental studies concerning the propagation of torsional waves.
Convery and Pugh (1968) gave the results of their experiments in which
a tube was stressed statically above the yield stress in torsion and
then subjected to a suddenly applied incremental torsional load. The
strain caused by this inciremcntal load was lound to propagate with the
elastic shear wave velocity. This seemed to be proof that the strain
rate dependent theory was correct, but Convery and Pugh (1968) cau
tioned against that conclusion. For Bell (1960, 1963) and Bell and
Stein (1962) had asserted that (based on experimental results with
annealed aluminum), while an increment of strain may propagate with the
elastic wave velocity, the larger amplitude strains propagate with the
wave velocity predicted by the strainrate independent theory.
Nicholas and Garey (1969) tested aluminum samples in torsion at high
strain rates and found very little strainrate dependence. However,
Yew and Richardson (1969) were able to measure some strainrate depen
dence in copper.
Another problem which was encountered in wave propagation studies
r
was that of unloading. The two most common unloading cases were when
the applied load was reduced and when waves were reflected from a bound
ary. Unloading was examined for longitudinal plastic wave propagation
by Lee (1953) using the strainrate independent constitutive equation
and by Cristescu (1965), Lubliner and Valathur (1969), and Cristescu
(1972) using the quasilinear constitutive equation. In all of these
investigations, regions of unloading and boundaries between regions of
unloading and loading in the characteristic plane were predicted but
the results have not been verified experimentally.
Many investigators in recent years have become interested in the
behavior of materials under combined stress and, more specifically, the
propagation of waves of combined stress. One of the first discussions
of combined stress wave propagation was given by Ranhmatulin (1958).
In this paper he developed the equations which must be solved for
elasticplastic wave propagation under combined stress. Strainrate
independent constitutive equations were used and only the problem for
the elastic case was solved. He found that the shear wave did not
affect the longitudinal wave in the elastic case. A similar discussion
of combined stress wave propagation was presented by Cristescu (1959).
Until now nothing has been said about the plasticity theory used.
The two plastic strain theories were the total strain theory proposed
by IIencky (1924) and the incremental strain theory proposed by Prandtl
(1924) and Reuss (1930). These two plasticity theories along with many
other developments in plasticity theory were presented in detail by
Hill (1950). The different plasticity theories were not presented
earlier because in many cases both theories gave the same results.
For instance, when a strainrate independent constitutive equation was
used, the two plasticity theories led to identical results when one
dimensional (either longitudinal or torsional) stress wave propagation
was studied, when combined stresses were used if the loading was pro
portional, or even when unloading occurred in onedimensional problems.
However, when strainrate dependent material behavior of nonpropor
tional loading under combined stresses was considered, most investiga
tors used the incremental strain theory. Shammamy and Sidebottom
(1967) showed that the incremental strain theory more accurately pre
dicted the experimental results when various metal tubes were subjected
to nonproportional static loading in tension (compression) and torsion.
Interest in the propagation of waves of combined stress continued
and Cli ton (1966) presented the results of his study of combined longi
tudinal and torsional plastic wave propagation in a thinwalled tube.
Strainrate independent material behavior and incremental strain
theory were used while radial inertia effects were ignored. The thin
walled tube allowed Clifton to eliminate any dependence on the radial
coordinate so that a solution could be obtained in the characteristic
plane. (Earlier, Plass and Ripperger (1960) had used a rod and averaged
the variables over the cross section in order to eliminate the dependence
on the radial coordinate.) The results of this investigation were
based on a simple wave solution which resulted from applying a step
velocity impact at the end of the tube. Clifton (1966) found that when
the tube was stressed into the plastic range, an impact at the end of
the tube caused waves with two different speeds to propagate. These
waves were called the fast wave and the slow wave, and each wave was
found to carry both longitudinal and torsional stresses. Two special
cases were examined. The first case involved statically prestressing
the tube above the yield stress in torsion and then applying a longi
tudinal velocity at the end. In this case the fast wave caused almost
neutral loading, that is, as the fast wave passed a point on the tube,
the shear stress decreased and the longitudinal stress increased in such
a way that the stress state at that point remained close to the initial
loading surface. Then as the slow wave passed the same point, loading
occurred so that the stress path was normal to the initial loading sur
face. The second case was for a tube with a static longitudinal plastic
prestress impacted by a torsional velocity at the end. In this case the
fast wave caused unloading along the longitudinal stress axis followed
by an increase in shear stress at a constant value of longitudinal
stress and then the slow wave caused loading such lhail the stress path
was normal to the initial loading surface. Clifton (1966) also found
that for a given initial loading surface, the two wave speeds depended
upon the particular stress state on the initial loading surface, and
that for one particular initial loading surface the fast and slow wave
speeds were equal when the shear stress vanished.
This work of Clifton (1966) was a significant step forward in the
investigation of waves of combined stress. An extension of this work
was presented by Clifton (1968) in which the simple wave solution was
used along with unloading at the impact end. In this way certain unload
ing boundaries for combined stress states were determined. Two years
later Lipkin and Clifton (1970) published their experimental results
from combined stress wave propagation tests and compared these results
to the simple wave solution developed earlier. Agreement between the
simple wave theory and the experiments was fair.
Cristescu (1967b) formulated the problem of combined stress wave
propagation in a thinwalled tube using general quasilinear constitu
tive equations but again ignoring radial inertia effects. The equa
tions for the characteristic lines and the equations along these char
axteristic lines were determined. No numerical results were given
but the two waves (fast wave and slow wave) were shown to be coupled
during loading. Again Cristescu (1971) showed that the coupling of the
waves of combined stress depended on the constitutive equations and
yield conditions used.
This concludes a brief survey of the history of the development
of plastic wave propagation theory. No attempt ,vas m3nd1" tc givO
a complete historical background. For more information the reader is
directed to Hopkins (1961), Kolsky (1963), Olszak, Mroz, and Perzyna
(1963), and Cristescu (1967a, 1968).
The remainder of this dissertation will be devoted to solving
the problem of combined stress wave propagation in a thinwalled tube
when radial inertia effects are included. A general quasilinear
constitutive equation for multiple states of stress will be presented,
and it will be shown to be a generalization of the constitutive equa
tions of both Lipkin and Clifton (1970) and Cristescu (1972). But
first the wave propagation problem itself must be developed.
CHAPTER 2
THEORETICAL DEVELOPMENT
The specific problem to be considered here is that of the propagation
of inelastic waves of combined stress along a semiinfinite thinwalled
tube, with the effects of radial inertia included. The material consti
tutive equation used is a generalization for multiple states of stress
of the quasilinear constitutive equation used by Cristescu (1972) for
a single stress component, and is a special case of the very general
quasilinear constitutive equation given by Cristescu (1967a). The coor
dinate system used is shown in Figure 2.1, and the stresses on an ele
ment of the tube are shown in Figure 2.2, where r is the mean radius of
the tube.
The problem is assumed to be axisymmetric so that there is no
dependence on e. Since the tube considered is thinwalled, the stresses
Sr, Tr, and Trx are assumed to be negligibly small as are the strains
e r and e The strain e is not included in the problem. Stability
re rx r
of the tube wall and thermal effects are not included in the formulation
of the problem, and only small strains are used. The strain rate is
assumed separable into elastic, plastic, and viscoplastic parts. The
radial displacement is very small compared to the tube radius, and
plane sections of the tube remain plane. The material is assumed to be
isotropic and homogeneous, to obey the von Mises yield condition, and
to be isotropically workhardening. All unloading is assumed to be
elastic.
Figure 2.1 Coordinate System for the ThinWalled Tube
ax
rx
Figure 2.2 Stresses on an Element of the Tube
15
The equations of motion in the cylindrical coordinates shown in
Figure 2.1 are given by
1 1
C + T + + T =
x,x rx,r r ex,e r rx x,tt
1 1
T + C + + ) = p Ur
rx,x r,r r r@e, r r r,tt
1 2
T + T +  T = 0 U
6x,x re,r r e,e r re9 ,tt
which, under the assumptions given above, become
(xx =p uxtt (2.1)
X,x P xtt 
 = u (2.2)
r r,tt
o
Texx = P utt (2.3)
where the subscripts following the comma represent partial differentia
tion with respect to the variables x (the coordinate along the tube
axis) or t (time). The density of the material is p, and ux, u and
u are the displacements of any point in the x, r, and 9 direction,
respectively.
For the cylindrical coordinates of Figure 2.1, the strain
displacement equations are given by
e =
r r,r
C (1 + (u )
r 9,0 r
e = U
x x,x
r 11 r 
r 2(r r,6 + Ur 
1
S= (u + u )
= (u + u )
Px 2 x r x,
and under the above restrictions, these reduce to the following three
equations
S= u (2.4)
x X,x
1
e = u (2.5)
o
1
S= u (2.6)
ex 2 e,x
Defining the velocities v v, and va as u t, ur,t and ue ,
respectively, equations (2.1) to (2.6) become
xx = x,t (2.7)
r e = r,t v(2.8)
o
Tx,x= p V,t (2.9)
S = v (2.10)
x,t X,X
1
=  v (2.11)
xt = vx (2.12)
ex,t 2 e,x
Under the assumptions used here, the variables no longer depend on r,
so that the problem becomes twodimensional (the independent variables
are x and t) and can be solved by the method of characteristics.
The equations necessary for completion of the set of simultaneous
partial differential equations describing the behavior of the body are
the constitutive equations. Cirstescr. (1972) uses a full quasilinear
constitutive equation for ,a in:.;le sot'itu diinal stress as
Ot =E t+ (, ~ + (~ ) (2.13)
As a generalization of this equation to a constitutive equation
governing multiaxial states of stress and strain, the following equation
is used
+ 3+ 3 ij
S l+v + (s,A)s + (s,) (2.14)
ij E ij E ijkk A)s (2.14)
s
where the dot represents partial differentiation with respect to time,
s.. is the deviatoric stress, 6.. is the Kronecker delta, v is Poisson's
ij ij
ratio, E is Young's modulus, 0(s,iA) and '(s,A) are material response
functions as yet unspecified, and s and A are defined as
s / s..s.. (2.15)
S 2 ij ij
/2 .P "P s
A ij ij dt + (2.16)
v 3 / 13ij E
*P
and e. is the inelastic portion of the strain rate which, using
ij
equation (2.14) can be written as
P 3 "3 sii
P. 3 0(s,A)s + (s,A) (2.17)
3ij 2L I
s
when the elastic, plastic, and viscoplastic portions of the strain
rate are assumed to be separable. The constitutive equation (2.14)
is a special case of the equation
k1
e. = fk. 0 + gi
ij ij kl ij
given by Cristescu (1967a). The form of equation (2.14) was chosen as
Lhe general constitutive equation because it contains terms which may be
considered separately as elastic, plastic, and viscoplastic strain
rate terms, because the inelastic strainrate tensor is proportional
to the corresponding deviatoric stress tensor, and because it reduces
to the form of equation (2.13) when the only stress present is the longi
tudinal stress. This simplification to the form of equation (2.13) is
shown in Appendix A.
The functions 0(s,A) and *(s,A) are functions which depend on the
particular material being studied. The function 0(s,A) is a measure of
the rate insensitive inelastic workhardening, and the function r(s,A)
is a measure of the viscoplastic strain rate due to the strainrate
sensitivity of the material. In the classical rate independent plastic
ity theory, t(s,) vanishes. When s < 0 or when s < a (the current
"yield stress"), 0(s,4) is set equal to zero. The unloading conditions
when ((s,6) = 0 are stated in equation (A.3.1).
Two separate materials are modeled in the numerical work done.
One is a 3003H 14 aluminum alloy used in the experimental work of
Lipkin and Clifton (1970). This material is assumed to be insensitive
to strain rate and the functions 0(s,A) and r(s,A) are obtained using
the classical PrandtlReuss incremental plasticity theory with iso
tropic workhardening and the stressstrain curve for uniaxial tension.
(See Appendix A.) The other material used is a commercially pure
aluminum dead annealed at 1100F. This material is assumed to be
strainrate sensitive, and the functions 0(s,A) and '(s,A) are
obtained from the data given by Cristescu (1972). (See Appendix A.)
Since the stresses a, T and T are assumed to vanish,
r r rx
equation (2.141) as applied to the present problem r reduces to
1 v 1 (s,)s+(sA)
2s
v 1 1
xt E x,t E E ,t 2
2s
e+1 ( 0 ( )s+ (s,A)s + (s.A) (2.18)
+ S
xt E +x,t (sA)s (s, )
2s
where the deviatoric stresses are
1
s = s (2a ()
x 11 3 x
1
s = s (+ a )
r =33 3 x + )
s = S = s1 = 2 Tx
sx 12 21 ex
sr = s = .
re rx
Using these deviatoric stresses, the expression for s becomes
1 1
ds = t L sijij 2 i sklkl s L 'ijsi
1 3 1 3 a F
s L2 sijsi= LssIl + s22s22+ s33s33+ sl2sl2+ s21s21
S 3 o 2 2 2 2
s (Cx + ~ ) + 2T
4s
i 1 r 1
s 1 (2 ) (20 )a + 6T7 (2.19)
x x,t x ,t x x
2s
and equations (2.18) become
(2c % ,7"2 .(s,Z) (2 a7,) (a 2cr) 
L ix ,t lJ i' j et
6(20 a )T (2c c )
42 ex ex (S A)] + t (S,) (2.20)
+ 2 (s,A) Text + ((s,A) (2.20)
4s 2s
(20 a,)(a 2 j F (a 2c
,t E 2 (s A) x,t L+ 2 ( c ,t
4s 4s
(6(2crQ )IxTx j 2crQ cr
6 ( 2 ,7 C Y x ) T e x 2 a 9 C 7x
+ ((s,A) Txt + (s,A) (2.21)
L 4s 2s
r3T e(2x ) 3Te] 3r e (2e 1)
ex, t 2= L0( xt 2(A ,t
4s 4s
18T 3T
P1_+ Ox ex
+  0(sA) T +  (s,6) (2.22)
4s 2s
The equations (2.7) to (2.12) and equations (2.20) to (2.22) form
a set of nine simultaneous hyperbolic quasilinear partial differential
equations for the nine unknowns ax' cre T6 x' V Vx', V, e', e and
e6x. A special case of this system of equations is the set of equations
obtained by neglecting radial inertia effects. When radial inertia
effects are ignored, the variables a', E and vr are not included
directly in the problem formulation. This case can be incorporated
into the more general formulation by multiplying a~, e vr, and their
derivatives by the dummy variable "a," where "a" has the value of 1
when radial inertia effects are included and the value of 0 when radial
inertia effects are neglected. Also the equation of motion, the kine
matic equation and the constitutive equation for motion in the radial
direction must be multiplied by "a." Doing this, and defining the
quantities
2
1 (2 a% )
A = + (s,6)
1 E 2(
4s
2 L
(2j ac,)(c
x x
2
4s
2a~ )
6Te (2y aOe )
A3 2 (s,A)
4s
A=
4 E
(7 2ac )2
+ s(2)
4s
6 T Gex(ax 2aca ) 
A = 02(, A)
4s 2
4s
(2.23)
the nine simultaneous equations (2.7) to (2.12) and equations (2.20) to
(2.22) can be written as
,x=vx,t (2.24)
a
r a = apvt (2.25)
ex,x ,t (2.26)
x,x
a
ae = r v
e,t r r
o
1
ex,t = Ve,x
ex, = Ax, + +aA + A +
x,t 1 x,t 2 e,t 3 ex,t
as ,t= aA 2 + aAa ,t
9,t 2 x,t 4 e,t
2 aaa
x ,
(s,6)
2s
2aa a
'2s
3T
1 a 1 ex
, = A3x + A + T +  (sA)
Ex,t 2 3 ,t 3 A5 ,t 2 56 6xt 2S
2s
(2.27)
(2.28)
(2.29)
(2.30)
(2.31)
(2.32)
x,t
1
3(5.1_. )_1
Eliminating the strain rates from the last six of these equations,
and defining
2a au,
x e
x (,)
2s
2aa a
= x (sA)
2s
(2.33)
the system of
equations for
3T
Ox 4(s,/) J
2s
nine equations reduces to the following system of six
the unknown variables ex, 1O, Tx' Vx, v and v
a x= xt
X,X X,t
a
 a7
r
o
ap vr,t
r,t
T6x,x = Pv,t
v = Aa + aA + A T + X
x,x 1 x,t 2 9,t 3 ,t x
a 
 v = aA +aAo + aAT + a
r r aA2 x,t 4 9,t 5 9x,t a
o
O,x= A3 x,t + aAO 9,t + A5 x,t + 2ex
2 2 2 2 1
s = (c ax 0 + a + 3T Ox
(2.34)
(2.35)
(2.36)
(2.37)
Since the equations (2.24) to (2.26) and (2.34) to (2.36) form
a system of hyperbolic equations, they can be solved by the method of
characteristics. To do this, first the equations for the characteristic
lines must be determined, and then the equations along these characteristic
where
A w, B b (2.38)
where
p 0 0 0 0 0
0 ap 0 0 0 0
0 0 p 0 0 0
A = (2.39)
0 0 0 A! aA A3
0 0 0 aA2 aA aA5
0 0 0 A3 aA5 A6
0 0 0 1 0 0 
0 0 0 0 0 0
0 0 0 0 0 1
B = (2.40)
1 0 0 0 0 0
0 0 0 0 0 0
0 0 1 0 0 0
V
x
v
r
w = (2.41)
CT
x
LT:xJ
a,.
and
0
r
o
0
x
x
av
o
2a
r @Xe
 2te X
and if the slope of the characteristic line (or the wave speed ) is
denoted by c where
dx
dc t
then the equation for the characteristic lines is given by
cA B =0
and from the calculations shown in Appendix B, equation (2.
[ 2 2 Cr 2,2 c
a pc J L(pc ) [a}
(2.43)
(2.44)
44) yields
2 (
 (pc ) {b + A4j = 0 (2.45)
where
S2 2 2
a = A 4 A6 + 2aA2A 3 5 aA 1 A5 aA26 A3A
(2.46)
2 2
= A1 aA2 + A A a. 5
1 4 2 4 6 5
Setting the first factor in (2.45) equal to zero,
dx
ac = a = 0 (twice)
dt
and, setting the second factor in (2.45) equal to zero,
r= I 1 2 1il 2
c = L (b 4aA4) 4
2ap
(2.42)
(2.47)
(2.48)
If the wave speeds in equation (2.48) are denoted by
c = i  b (b (2. 49)
f L9_ fj
2ap
and
c = ( 4aA ) (2.50)
2ap
where cf is the fast wave speed and cs is the slow wave speed, then
the slopes of the characteristic lines are given by
c = 0 (twice)
c = f
c = c
s
Equations (2.47) and (2.48) are the six equations for the character
istic lines for the set of six simultaneous, hyperbolic, quasilinear
partial differential equations of (2.38). When radial inertia effects
are not included (that is, when a= 0), the equations (2.47) vanish
identically and the remaining four simultaneous equations of (2.38) have
the characteristics given by equations (2.48). For this case (a=0) and
when 0(s,A) is obtained from incremental, rate independent plasticity
theory with isotropic workhardening, equations (2.48) reduce to those
given by Clifton (1966) as shown in Appendix B.
The equations along the characteristics can be obtained in two
different ways, both of which are discussed in Appendix B. The result
ing equations along the characteristic lines of equation (2.48) are
0 =
5 A A ]dTex
5 3 4iJ x
2 {r
+ pc2)(A2 6 35 2 1L
o
1
c
c
+ 2(pc )aA2A5A3A4 ( dt
2s
( 2aaG)
2s
(s,A) dt
[(P c) (iA6A5) A dv  (pc 2) (AiA5A2A3) A dvg
c
S 2 6 
S (P c ) (A115A 23) 5 dTgx (p 2 ) (A2A6A3A5 dex
pc pc
[ 3 ] [ (2Jxa a) )]dt
2 x (s,
+ (PC )(2 6AA 5 A (s,) 2dt
2s
1
+ .
2
PC
2 2 2 (2 AA)]
(pc ) (A16A3) (pc )(A1 )+1 1
o
(a 2aa )
x 6
+2
2s
2 3ex3T
+ 2 (PC)(A 1 A 2 A35 (s,6)]dt
2s
23 ( C 2 A jidA aA 2 d
O =p c aA2 A34dvx oc c)( 1 4aA2)4 dv
+ F(pc2)
L
(A A4aA2)Aj dTx
1 4 2 4~
+ ( c2a23
+ (pc )LaA2A53 A4] dx
2 2 ( x c
+ (pc2)2 aA2AAA 3J x (sA) dt
2s
1 2 2
[(p2) (A" 2 (A 4 dv, ~ 6 La 2
+ 2 [(PC (A4A6aA)AX d [(pc ) (A4 
pc
[(2a ac7 )
(s,A) dt
2s
(2.51)
(2.52)
(s,A) dt
S2r 9 _2 r'o7 (, 2ao ) ]
+ a(pc ) (PC ) 1(A.A2A)i I) dt
0
+ 2(pc(p pc) (A1A 4aA2) 4I L ,) dt (2.53)
2s
These three equations each represent four equations, one equation
in differential form along each of the four charaxeristic lines of
equation (2.48). When the waves are coupled, equations (2.51), (2.52),
and (2.53) are identical. That is, by multiplying equation (2.52) by
the quantity
(pC ) (aA25 3A4)
24
(pc )(A1A A2A3) A5
and using equation (2.48), equation (2.51) is obtained; or by multiply
ing equation (2.52) by the quantity
2 2
pc (pc 2)(A a2) A
2
(pc )(A A A2A3) A5
and using equation (2.48), equation (2.53) is found. When the numerator
and denominator of these multiplying quantities do not vanish, equations
(2.51), (2.52), and (2.53) are identical. However, when the waves
become uncoupled, a phenomenon discussed in Appendix B, A3 and A5 vanish.
In this case the multiplying factors used above become undefined and
the equations (2.51), (2.52), and (2.53) are not the same. When
A =A =0 the equations (2.51) and (2.53) reduce to equations (B.4.4)
and (B.4.6), respectively. Under these conditions, equation (2.52)
also reduces to the form of equation (B.4.4).
The equations along the characteristics of equation (2.47) may be
obtained (since :nlonig Lhese characteristics there is no variation in x)
directly from equations (2.25) and (2.35) by multiplying each of these
equations by the increment of time, dt. These equations then yield
aco
S dt = ap dv (2.54)
r r
0
av
r dt = aA 2dxa + aA4do + aA 5dTrx + a4 dt (2.55)
The equations for the six characteristics and the equations along
these characteristics, along with the appropriate initial conditions
and boundary conditions, represent a complete mathematical formulation
of the problem, and the solution to these equations is the solution to
the problem posed here. The solution to these equations will be
obtained by using a finite difference numerical technique which will
be discussed in the next chapter.
CHAPTER 3
DIMENSIONLESS EQUATIONS AND NUMERICAL PROCEDURES
3.1 Wave Speeds as a Function of
the State of Stress
In this chapter the numerical schemes used to find the solution to
the wave propagation problem of Chapter 2 will be presented. In this
first section the dependence of the wave speeds on the stress state
will be shown. The stresses a T, and have already been assumed
r r' rx
negligibly small so that the scalar representation of the stress state
is given by equation (A.1.2) as
(2 9 9x9 2 2 
s = au T + a + 3T 2 (3.1.1)
x x 6 9 6
Next, the new variables a ae and Tx will be defined so that the
x 9 ex
surface s = constant can be represented in terms of these variables as
a sphere, and the stress state on this surface in terms of these new
variables can be described in terms of spherical coordinates. Now
defining,
/ 1
= 2( + aca)
x 2x G
a (ae ) (3.1.2)
*'Bx = /3 T
x ex
equation (3.1.1) can be written as
12 + 2 2T] (3.1.3)
s = / + a/ + T (3.1.3)
x 9 9xJ
and defining the angles y and S6 as shown in Figure 3.1 these new
variables defined in equation (3.1.2) Con Ibe written as
x = s cos y cos 6
e, = s cos y sin 6
(3.1.4)
x = s sin y
The angle y is the complementary angle to the one normally used in
spherical coordinates. It is used here to facilitate comparison of
results obtained later on to already published results.
From equations (3.1.2) and (3.1.4),
9 s cos v sin 6
/ 
S s cos v cos 6
x
so that the cx and ae
x 9
a axis : 0 =
( + a(') 7 + ac
tan =
S(aa a) a x_
2 9 x
axes are located by
and tan 6 = and 6 = 60
and tan 6 and 6 = +60.
In order for the equations (3.1.4) to reduce automatically to the
simpler case when radial inertia is not considered, the angle 6 is
defined as
6 = a6' + (a1)600
(3.1.7)
so that when radial inertia effects are included, a=l, and 6= 6', and
when radial inertia effects are not included, a= 0, and 6 = 60, which
from equation (3.1.6) automatically causes (a to vanish as it should.
e
(3.1.5)
(3.1.6)
x
ex
Figure 3.1 Yield Surface Representation in Spherical Coordinates
Using the uniaxial s;rs strain curve in the form of equation
(A.2.13), the universal str'ess:strain curve can Ie r ..itten as
(3.1.8)
S n
S= B )
y
and letting Et(s) be the tangent modulus of this curve, yields
1 dA 1 nl
S + Bn (s 
Et(s) ds y
(3.1.9)
and from equation (A.2.18) this becomes
1 1
+ O(s)
Et(s) E
() 1 1
E (s)
(3.1.10)
Now g = n(s) is defined so that
where
Et(s) = 1(s) E
0 I $(s) < 1
and when 5=1, the material is elastic, and when P=0, the material
is perfectly plastic. Using equation (3.1.11) in equation (3.1.10),
0(s)can be written as
Inverting equations (3
1 1
(s) = (  1) .
.1.2), the stresses a
x x 3
/ 1 ,
ao = 0 +  Ca
ao8 O Fx
(3.1.12)
re given by
(3.1.13)
1 1
78x :
(3.1.11)
and using equation (3.1.4) thse become
i 1
S CO ' jlC,; S c .
1 1
aoe = s cos ycos 6 +  sin
T  s sin y
Now substituting equations (3.1.12) and (3.1.14)
S(3.1.14)
into equation (2.23),
F 1 ) 2 2 1 21
I 1+ 1) (1)(cos Y) (2 cos 6  sin 6 cos 6  sin 6)
1r 1 1 2 2 1
E [ ( 1) )(cos y(2 cos  sin 6 cos 6  sin 6)
(cos 6  sin 62 cos 6 2 sin 6)
1 1 6 1 1
E  sin v) (cos v)(2 cos 6  sin 6 cos 6  sin 6)
S3 I 1 2 2 2N
= 1 + 14 1) (cos 6  sin 6 2 cos  sin 6) cos
E1 i (1 1 1 2 ]
E (1( sin y)(cos y)(cos 6  sin 6 2 cos 6  sin 6)
1F 36 1 1 2
[2(+1 + (1) ( sin Y)
E 4 B3
1 14 11 2 2 2
S[1 + ( 1) (cos v)(cs6 2 sin 6 cos + 3 sin 6)
1 I 1 1 2 2 2
1 + ( 1)(cos 2y) (cos2 6+3 sin 6)
) () (I 1) (sin y cos y) (cos 6 ,3 sin 6)
E 2 J (3.1.15)
1 1 1 2 2 2
1 +(+ 1) cos y(cos 6 + 2A/3cos 6 sin 6+ 3 sin 26)
E v( 1) (sin y cos y)(cos 6 /3 sin 6)
S2(1+v) + 3(11) sin y
A3
A4
A5
A 6
or
A1
2
A3
A4
A
A5
A6
The elastic wave speeds arc defined from equations (B.5.4) as
c /
o / c
c1 = (3.1.16)
(1 )
c2 ~ p
and the wave speeds from equation (2.48) can be written in dimension
less form as
c = L b Vb 4aj a (G/p)
c2 k2 ap L 4
C2
c b 4A4a
1 Lv j 24Aa]
Ea
or
12 lv F 2 3
C 2 = E E2 b (E2b )2 4(EA )(E3 a) (3.1.17)
3 4
E a
By defining the dimensionless functions from equation (3.1.15) as
Ai = E A i = 1,2,...,6 (3.1.18)
1
and
3 2 2 2
a = Ea = A1A4A + 2aA2A3A5 A3A4 aA1A aA2A6
S(3.1.19)
b' = E = A1A aA2 + A4A aA ,
1 4 2 46 5 ./
then the wave speeds in dimensionless form become
/2 l+v r; / i ~2 ~) / c \2
c = Vb b 4A4a = (3.1.20)
a 2
A computer program was written to solve this equation for the two
positive wave speeds as functions of the angles y and 6 for specified
values of v and 8. This pDrol i i lin ised in Appendix C. This pro
gram also calculated the viuecs of the normalizz:ed stresses a /s,
0 /s, and ex/s as functions of v and 6. Those results are given in
Tables 1, 2, and 3. The wave speeds are shown in Figure 3.2 for
the case when 6 =60o, which corresponds to j = 0. Also plotted in
Figure 3.2 are the results given by Clifton (1966). It is obvious
that the results are not the same and that including radial inertia
effects in the formulation of the problem can have significant effects
on the wave speeds and that, for any given state of stress, the waves
are always faster when radial inertia effects are included. The
results plotted in Figure 3.2 do not correspond to the case when
a= That is, although a = 0 when 6 =60 a does not necessarily
vanish for this case. Ahen a= 0, the results obtained were identical
to those of Clifton (1966) as they should be, since a= 0 corresponds
to the absence of radial inertia effects.
An interesting phenomenon can be observed by remembering that the
physical presence of radial inertia is due to the Poisson effect.
That is, the longitudinal (fast) wave speed would be expected to be
the same when a=0 (no radial inertia effects) as when v=0 (the cause
of the radial inertia effects vanishes). However, in the fomulation
of this problem it is tacitly assumed that Poisson's ratio for the
1
inelastic portion of the material behavior is , or that the material
behavior in the inelastic range is incompressible. When the material
is elastic ( = 1), this Poisson effect can be studied directly.
Comparing equations (3.1.15) with (B.5.1) when 1 = 1, the wave speeds
are given by equation (B.5.4) as
\
TABLE 1 NORMALtZLD ,LONGI'L:.UD.NAL STRESS ()
s
Gamma
Delta
00 100 200 300 400
900 0.57735 0.56858 0.54253 0.50000 0.44228
80 0.74223 0.73095 0.69747 0.64279 0.56858
70 0.88455 0.87111 0.83121 0.76605 0.67761
600 1.00000 0.98481 0.93969 0.86603 0.76604
500 1.08506 1.06858 1.01963 0.93969 0.83121
40 1.13716 1.11988 1.06858 0.98481 0.87111
30 1.15470 1.13716 1.08506 1.00000 0.88455
200 1.13716 1.11988 1.06858 0.98481 0.87111
100 1.08506 1.06858 1.01963 0.93969 0.83121
00 1.00000 0.98481 0.93969 0.86602 0.76604
100 0.88455 0.87111 0.83121 0.76604 0.67761
200 0.74223 0.73095 0.69746 0.64279 0.56858
300 0.57735 0.56858 0.54253 0.50000 0.44228
400 0.39493 0.38893 0.37111 0.34202 0.30253
500 0.20051 0.19746 0.18842 0.17365 0.15360
600 0.00000 0.00000 0.00000 0.00000 0.00000
700 0.20051 0.19747 0.18842 0.17365 0.15360
800 0.39493 0.38893 0.37111 0.34202 0.30254
0.56858 0.54253
900 0.57735
0.50000 0.44228
TABLE 1 (Continued)
Gamma
Delta
500 600 700 800 900
900
800
70
600
500
400
300
200
100
00
100
200
300
400
500
600
70
800
0.37111
0.47710
0.56858
0.64279
0.69747
0.73095
0.74223
0.73095
0.69747
0.64279
0.56858
0.47709
0.37111
0.25386
0.12889
0.00000
0.12889
0.25386
0.28868
0.37111
0.44228
0.50000
0.54253
0.56858
0.57735
0.56858
0.54253
0.50000
0.44228
0.37111
0.28867
0.19747
0.10026
0.00000
0.10026
0.19747
0.19747
0.25386
0.30253
0.34202
0.37111
0.38893
0.39493
0.38893
0.37111
0.34202
0.30253
0.25386
0.19747
0.13507
0.06858
0.00000
0.06858
0.13507
0.10026
0.12889
0.15360
0.17365
0.18842
0.19747
0.20051
0.19747
0.18842
0.17365
0.15360
0.12889
0.10026
0.06858
0.03482
0.00000
0.03482
0.06858
0.00000
0:00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.28868 0.19747
900 0.37111
0.10026 0.00000
sa,
TABLE, 2 NOWRALIZFD ITOOIP S TiESS
S
Gamma
Delta
00 100 200 300 400
900
800
700
600
50c
400
300
200
100
00
100
200
300
40
500
60"
700
800
900
0.57735
0.39493
0.20051
0.00000
0.20051
0.39493
0.57735
0.74223
0.88455
1.00000
1.08506
1.13716
1.15470
1.13716
1.08506
1.00000
0.88455
0.74223
0.57735
0.56858
0.38893
0.19746
0.00000
0.19747
0.38893
0.56858
0.73095
0.87111
0.98481
1.06858
1.11988
1.13716
1.11988
1.06858
0.98481
0.87111
0.73095
0.56858
0.54253
0.37111
0.18842
0.00000
0.18842
0.37111
0.54253
0.69747
0.83121
0.93969
1.01963
1.06858
1.08506
1.06858
1.01963
0.93969
0.83121
0.69746
0.54253
0.50000
0.34202
0.17365
0.00000
0.17365
0.34202
0.50000
0.64279
0.76605
0.86603
0.93969
0.98481
1.00000
0.98481
0.93969
0.86602
0.76604
0.64279
0.50000
0.44227
0.30253
0.15360
0.00000
0.15360
0.30254
0.44228
0.56858
0.67761
0.76604
0.83121
0.87111
0.88455
0.87111
0.83121
0.76604
0.67761
0.56858
0.44228
TABLE 2 (Continued)
Gamma
Delta
500 600 700 800 900
900
800
700
600
500
400
300
200
100
00
100
200
300
400
500
600
700
800
900
0.37111
0.25386
0.12889
0.00000
0.12889
0.25386
0.37111
0.47710
0.56858
0.64279
0.69747
0.73095
0.74223
0.73095
0.69747
0.64279
0.56858
0.47709
0.37111
0.28867
0.19746
0.10026
0.00000
0.10026
0.19747
0.28868
0.37111
0.44228
0.50000
0.54253
0.56858
0.57735
0.56858
0.54253
0.50000
0.44228
0.37111.
0.28867
0.19747
0.13507
0.06858
0.00000
0.06858
0.13507
0.19747
0.25386
0.30253
0.34202
0.37111
0.38893
0.39493
0.38893
0.37111
0.34202
0.30253
0.25386
0.19747
0.10026
0.06858
0.03482
0.00000
0.03482
0.06858
0.10026
0.12889
0.15360
0.17565
0.18842
0.19747
0.20051
0.19747
0.18842
0.17365
0.15360
0.12889
0.10026
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
/T
TABLE 3 NORLALIZED SHEAR STRESS \
T x
Gamma Value of x Gamma Value of Tx
s s
00 0.0 500 0.44228
100 0.10026 600 0.50000
200 0.19747 700 0.54253
300 0.28868 800 0.56858
400 0.37111 900 0.57735
_ E
c =
f2
p (1 av2)
S= E 
Cs p 2p (1 +) 2
It is now obvious that the fast wave speed is the same when a=0 and
v 0 as when a= 1 and v= O. The slow wave speed (and consequently c2)
is the same when a=1 as when a=0, although it does depend on v.
Because of this dependence of c2 on v, the dimensionless fast wave
speed of equation (3.1.20) will have values when a=0 and v O differ
ent from those when a=1 and v =0.
Also in Figure 3.2 it can be seen that when a=0 and y=0 the
fast and slow wave speeds are the same for = .385 and = .30. There
is usually some value of $ for which the two wave speeds are equal
at y=0 for each combination of values of v and 6. The condition for
which this is true can be obtained from equation (3.1.20) and is
Without Radial Inertia
Clifton (1I66)
 With Radial Inertia
and ca = (6 = 600); 6 0
S= .01
y, Degrees
Figure 3.2 Plastic Wave Speeds as Functions of $ and y for Poisson's Ratio of 0.30
0, 8
ti'hen = 0
A1
2
A
A3
A4
A =
5
A6
6
b/2 4a A = 0 .
equations (3.1.18) and (3.1.19) are
1 + 1(1)(cos 6 2v sin 6 cos 6 + 3 sin26)
43
v + (1)(cos 63 sin25)
0
1 1 2 2
1 ( 1)(cos 6 +23 sin 6 cos 6+ 3 sin 6)
1+45
(3.1.21)
(3.1.22)
0
2(
:1 + )
2 2
a = AIAA aAA6
b/ = A1A + AA aA2
and using the same manipulations as in Appendix B, Section 4,
equation (3.1.21) becomes
0 = A(A A ) + aA2
where Al, A2, A4, and A6 are given by equation (3.1.22). Now
fl = () = f cos26 2 3 sin 6 cos 6 + 3 sin2 ]
f2 =f2(6) = cos26 + 2v sin 6 cos 6 + 3 sin26
f= f(6) = 1 cos26 3 sin26]
1
z = z( ) 1
c c
c
(3.1.23)
(3.1.24)
defining
(3.1.25)
the critical value of $ is found from
0= L L c f L 2(1,.,) (I + zcl) I+ a (v cf3
0,= z fl 2+af3 + z ([2(1+l ) 13 f + a2vf3
+ [2(1 +v) 1+ av (3.1.26)
The expressions for the critical value of 3 will now be found
for the two separate cases of a= and a=1. First, when a= 0,
equation (3.1.24) becomes
A6 = A
2(1+v) 1+ z f
c 1
1 + 2v (3.1.27)
z 
c f
1
c z +1
C j
and in this case 6= 60" from equation (3.1.7) and
fl = cos2(60) 273 sin (600) cos (600) +3 sin2(600)
f 24( )( ) + 3( )
1 4 4 2 2 2
f =1
and
1
c (1+ 2v) +1
1
S2 ) (3.1.28)
c 2(1+ v)
When a = 1,
2 1 9
fl 2 af3 6 (cos" 2 ,':3sin cos c + 3 sin )
22 2 2 21]
(cos52 + 2/3 sin cos 6+3 sin6) (3 sin6 cos2)2
ff2 af = 0
1 2 3
and equation (3.1.26) becomes
[2(l +v) 1+ V
c [(2)(1 +) 1] f, f +2vf
(1+v)
z (3.1.29)
S 3 (1 +v) sin 6 cos + 3v sin26
A short computer program was written to calculate the critical
values of (using equation (3.1.28) when a=0 and equations (3.1.29)
and (3.1.27) when a= 1) for various values of v and 6 when y= 0.
This program is shown in Appendix C and the results are plotted in
Figure 3.3. The only values of 3 which are physically possible are
between 0 and 1 and therefore only values of $ in this range are
plotted in Figure 3.3. For all other values of 6, there is no phys
ically possible critical value of S; that is, there is no value of
such that the fast and slow wave speeds are equal at y= 0.
For the case when a = 0 (6 =60), for any value of v the critical
value of $ is smaller when radial inertia effects are included.
v = 0, a 0
v = .25,
=.30, a
Figure 3.3 Values of 3 at y=00 for which cf cs = c
s 2
.6
.4 _
0 30 60 90
6, Degrees
3.2 Characterirsi.c Solution in Terms
ol' uL'i ; loni:iolon ?le s C:: vria bles
In order to make the numerical solution in the characteristic
plane more general, the equations for the characteristics and the
equations along the characteristics given in Chapter 2 will be written
in terms of dimensionless variables. The dimensionless variables
used are
c
x 1 s
X T = t s = 
2r 2r E
0 0
x ex
S S T
x E E E
V VV V
v v v
x e r
V v = V =
x 9 r 
Sc c
c f s
c c c (3.2.1)
s
c c c
r
C1
,(Sf) 2o
0 (S
Y(s,A) = 0
s
where c1 given by equation (3.1.16) is the elastic longitudinal wave
speed when radial inertia effects are included. In terms of these
variables, the functions defined by equations (3.1.18) can be written
from equations (2.23) as
(2S aS )2
x 9
A1 EAI = 1 9(st)
A1 = E2 1 2 O(s,6)
4s
(2S aS(S 2aS)
A2 = EA2 L + 4s (s,
4s
6T(2S aS )
A3 = E 3 2 2(s)
4s
(S 2aS )2
A4 = EA4 2 (s,) +
4s
A = EA =
5 5.
A6 = EA6 =
6 6
6T(S 2aS )
4x
2
9T2
2(1+ v) + 2 @(s,A).
s
Using these and equations (3.1.19), the fast and slow wave speeds can
be written from equations (2.49) and (2.50) as
cl
c
s
s (
c1
S {b + (b2 4aA4 (1 2
2a p
= [E { b (b'2 4a'A4 (1
2a p
/2s
f =  b' + (b'2 4aAA4)2
2a
es b' (b2 4a'A4
L 2a'
(3.2.3)
and the wave speeds in dimensionless form are given by equations
(3.2.3) and (2.47) and (2.48) as
(3.2.2)
48
c = cf
c = c (3.2.4)
c = 0 (twice)
where the wave speeds were written in dimensionless form by dividing
the wave speeds by cl. This was done because cl is the fastest wave
speed possible in the problem considered here, and thus all the
dimensionless wae speeds have values in the range
1 c < 1 .
,When radial inertia effects are not included, the maximum value of the
fast wave speed is
C 
Cf p=
max
so that, for a=0, the maximum value of the dimensionless wave speed is
Cf
I cfI < max E/
1 E
p (1 v2)
Next, the equations along the characteristics for fully coupled
waves will be written in dimensionless form. Along the vertical char
acteristics (c=0), the equations can be written from equations (2.33),
(2.54) and (2.55) as
aS dT = a L 1 dV
2(1 o2)
or
a [2(1 v2)S] dT = adV (3.2.5)
.I2aS *S a
a(2V dT) = aA2dS A!dS +AArC +  (s,. j
r L 2 x 4 + 5 \ s 0 o
(3.2.6)
The equations along the nonvertical characteristics (c= cf, cs)
can be written in three different forms from equations (2.51), (2.52)
and (2.53). In dimensionless form these become, respectively,
O = I (A4A6 aA5) j dx A2A5 A3A dV
1 v v
2 2
1 c 2
S[aA2A5 A3A4] dT + 2 ( 2(A4A6 5 A4dS
c 1 v
[(2 ) r92S aS]
+ C 2(A4A6 _A2) < s "0 o0(s,)dT
K 2 \ n r S 2aS
+ a L( 22(A A(3 A A ) A 2V+
2+ (A2 A6A A3A5) 2J L2Vr s( )dT
22 ]
L+ 6 V2 2 L _A o(s')Jd
(3.2.7)
2
S= c 2 (A A AA5 )A dV 
c 1 \
+ ( 2,) c <2A
c 1 v
+ (. (A2A A3A5)
L \1 v 2A
r 2
L c (A5A2A )A dV
c L5 2v 5
 A5 dT
 A3A5) A dSx
2Sx aS)dT
A2] Lx o(s dT
i 2 "0
(1 2 r( e2 \2 3 2
+ ( ) [ 2) (A1A6 ) (A1 +A 1
c 1 v 1 V
and
F S 2aS 
2V + ( (sA) dT
r s 0o
+ 6 c 2 (A1A5 A2A3) (s,A)]dT (3.2.8)
and
3 2
0 c 2 aA2A A3A4 d c 2) (A1A4aA2) A dV
(1 v 2) 1 v v
+ c2 (A1A4 aA 2)A4 dT+ L_ 2A5A3A 4 dSx
1 v 1 V
2 2 r 2S a
+ ) LaA2AA3A4 L s o(s, l dT
1 i
+ a )L (AlA5A 2A3)A 2Vr x s(s'A)j dT
+ 6 2) (A aA ] )AA tsi AJdT (3.2.9)
1v 1 v
WVhen the equations are uncoupled (A3 = A 0), the equations along the
characteristics are given by equations (B.4.7) and (B.4.8). In dimen
sionless form, the equations (B.4.8) along the vertical characteristics
(c=0) become
a[2(a v )S ] dT = adV (3.2.10)
er
and
2aS S
a[2V dT] = a[AdS + A4dS +  (s,)dT] (3.2.11)
The equations (B.4.7) along the nonvertical characteristics (c= c )
for the fast waves are given by
1 1 2 2S aS
_0 1 dVd + 2 d +  dTA4
f c
r S 2aS
+ 2 L2Vr + ( ) o dT (3.2.12)
and the equations (B.4.7) along the nonvertical characteristics
(c= c ) for the slow waves can be written as
s
2
c 6c2
0 = s dV + dT + () o(s,T)d (3.2.13)
L(1v2) u (lv2)
The equations (3.2.5), (3.2.6), and (3.2.7), or (3.2.8), or
(3.2.9) are the equations along the characteristics for the fully
coupled waves written in differential form in terms of the dimension
less variables. The solution to this set of equations will be obtained
numerically by writing them in finite difference form, and then solving
the resulting set of algebraic equations simultaneously. When the
waves are uncoupled, the equations along the characteristics are given
by equations (3.2.10), (3.2.11), (3.2.12), and (3.2.13). These equa
tions will also be written in finite difference form and solved (when
applicable) in the same manner as described for fully coupled waves.
The procedure for obtaining these finite difference solutions is out
lined in the next sections.
3.3 Numerical Grid for Characteristic Solution
Since the slopes of the characteristic lines at any point in the
characteristic plane depend upon the state of stress at that point and
upon the history of the deformation at the corresponding location along
the axis of the tube, the equations for the characteristic lines cannot
be determined before the solution (in terms of stresses) is known.
Because of this, the slope of the characteristic lines and the solution
to the problem must be determined at each point simultaneously. This
is done by using the iterative numerical technique described below.
The numerical grid shown in Figure 3.4 will be used. There are
two types of elements in this grid: boundary elements and regular
elements. All of the regular elements are alike, and all the boundary
elements are like the righthand side of a regular element. A detailed
picture of a single regular element is shown in Figure 3.5, and a
single boundary element is shown in Figure 3.6. The grid is defined
in terms of the dimensionless variables given in equations (3.2.1) and
(3.2.3). It is diamond shaped with the straight outer lines corre
sponding to the characteristic lines for elastic longitudinal waves
with radial inertia effects included. These outer characteristic lines
have slopes of either c= +1 or c= 1, which can be seen from equations
(3.2.1) and (B.5.4). The vertical straight line corresponds to the
two vertical characteristic lines, and the straight inner nonvertical
lines correspond to the characteristic lines (through the point P) for
the elastic shear waves. For both types of elements, the problem reduces
to that of determining the values of the stresses and velocities at the
point P, when their values at the points B, R, and L are known.
This grid with all the elements constant in size simplifies the
writing of the finite difference equations. The diamond shape allows
the vertical characteristic lines to automatically connect point P (at
which the solution is desired) with point B (at which the solution is
known) and makes the finite difference equations along the vertical
53
T
AFigX 3.4 Nu l Gd in te Chc
Figure 3.4 Numerical Grid in the Characteristic Plane
2AT
2AT
77
Ic = c = 

AI
1c
B
a
Figure 3.5 Regular Element in Numerical Grid
op
55
T
P
c = c = 
C 1 I
c c= c =1
B
 AX  k
LX
Figure 3.6 Boundary Element in Numerical Grid
characteristic lines very easy to obtain. The boundary lines for each
element are c= This is the smallest value of c which insures that
all characteristic lines passing through tne point P will intersect
the line LB between the points L and D if the lines have a positive
slope at P or will intersect the line RB between the points R and B
if they have a negative slope at P. This is true since all of the
waves considered here will propagate with a speed less than or equal
to the speed of an elastic longitudinal wave with radial inertia effects
included. A larger value of c could be used, but the element size
would increase (for a given distance along the T axis), and the solution
would be inherently less accurate.
The straight lines representing the elastic shear wave character
jstic lines are added to the grid elements as a convenience. The
results of Section 3.1 show that the fast wave speed always occurs in
the range
c2 : cf 1 c1
and the slow wave speed always occurs in the range
0 O c s c2
s 2
Therefore, these characteristic lines c=c2 divide each element so
that a characteristic line through P lies in one of the upper triangles
(P L LB or P RRB) if it is for a fast wave and in one of the lower
triangles (P BLB or P BRB) if it is for a slow wave. This is
shown for the regular grid elements in Figure 3.7. These characteristic
lines for the fast and slow waves will not, in general, be straight.
P
e /I \ \=\
c = c2 / \ 0C 2
/ c c = cC
c = 0
B
  ^ x
Figure 3.7 Location of the Characteristic Lines Passing Through P
3.4 Finite Difference Equations
General Discussion
1hile the actual characteristic lines for the fast and slow waves
are seldom straight, they can be represented as straight lines within
each grid element without introducing significant errors if the grid
elements are small. From the discussion in Section 3.3, it is known
that the slope of the characteristic lines at any point cannot be
determined before the solution at that point is known. Because of
this the solution at the point P (Figures 3.5, 3.6, and 3.7) must be
obtained by an iterative technique. Within any grid element, the slope
of each characteristic line will be constant during each iteration
although the slope of each characteristic line will change from one
iteration to the next as the solution at P is approached. These
straight lines are used to represent the characteristic lines for
th
c= c and c = c during each iteration and are shown for the i iter
f s
ation as c=cfi and c= cs in Figure 3.8 for a regular grid element
and in Figure 3.9 for a boundary grid element. The points LLB, LBB,
RBB, and RRB are the intersections of the lines shown in Figure 3.8.
Each element has its own coordinate system X T which is also shown
in Figures 3.8 and 3.9, and the finite difference equations are written
in terms of this local coordinate system so that the finite difference
equations for each element are the same.
First order finite difference equations will be written along each
characteristic line. The coefficients of the dependent variables in
these equations will in general be functions of the stresses and (s,A).
Thus, in order to linearize the equations, the coefficients for each
T 2
bx
Figure 3.8 Numerical Representation of the Characteristic Lines in a Regular Grid Element
60
I
ST
P
c = C1
C i=
C = cf
R
RRBI
c=c )\
s. \
1
AT
RBB
c=0 RBB
B
^  ^X  =
X
Figure 3.9 Numerical Representation of the Characteristic Lines
in a Boundary Grid Element
iteration will be calculated using the solution obtained in the
previous iteration. In this way the coefficients are :llays known
quantities.
One other scheme will be used with the coefficients in the finite
difference equations in order to reduce the time required for compu
tation. Normally each coefficient used is the average value of that
coefficient at the end points of the interval over which the finite
difference equations are written. As an example, consider the charac
teristic line from point LLB to point P, and let one term in the finite
difference equation along this characteristic line be
U(S S ).
xP xLLB
As a rule the value of the coefficient is calculated as
1
U = (Up + ULLB)
If this method is used, the coefficients of each variable in the equa
tions along the characteristics of positive slope will be different
from the coefficients of the corresponding variables in the equations
along the characteristics of negative slope. For instance, one term
in the equation along c=+cf. can be represented as
S(SxP xLLB
and the corresponding term in the equation along c=:cfi as
UR(Sx SxRRB
where
1
U = (U + ULLB)
S= +
R 2 p RRB
and U is the value of the coefficient at point P calculated from the
solution from the previous iteration.
When the coefficients are calculated in this manner, the number of
equations which must be solved simultaneously cannot be conveniently
reduced below five. However, if the coefficients are calculated in
such a way that the coefficient of any variable in the equation along
c=+cfi is equal to the coefficient of that same variable in the
equation along c=cfi (so that UL = UR, etc.), then by adding these
two equations and subtracting one from the other, two different equa
tions can be obtained, each with fewer variables than the two original
equations. If this procedure is applied to the equations along c=Csi,
then the set of five simultaneous equations can be reduced to at most
a set of three simultaneous equations and a set of two simultaneous
equations. This is shown in Appendix D. Since this set (or sets) of
equations must be solved during each iteration, the savings in compu
tation time is significant.
One way to make the coefficients of similar terms equal is to
calculate the coefficients from values of the variables obtained at
point P during the previous iteration such as
UL = R = U *
The coefficients will be calculated in a somewhat more accurate manner
by using a weighted average of the value of each coefficient between
point P and point B, that is
UL = UR = alUP + (1 )UB
This gives the value of the coefficients at a point nearer the center
of each grid element. For this work, the value of ai is chosen arbi
trarily as .625, so that the point at which the coefficients are calcu
lated is at approximately the same location along the T'axis
(Figures 3.8 and 3.9) as the centers of the four characteristic lines
c= cfi, Csi.
The values of all quantities at the points LB, LLB, and LBB will
be obtained by linear interpolation between the points L and B.
Similarly, the values of all quantities at the points RB, RRB, and
RBB will be obtained by linear interpolation between the points R and
B. From Figure 3.8, the times T1, T2, and T3 can be written as
2c 2c 2c
s 2 f
T, = AT T c AT T +c AT
1 1 +c 2 1 + c 3 1 +c
s 2 f
and the interpolation constants for the points LB and RB are
2 2c2
CLRB 
AT 1+ c2
T2 1 2
CLRBI = 1
AT 1 + c2
Using subscripts to denote the grid point, the values of any quantity F
at the points LB and RB are
FLB = CLRB*FL + CLRBI*FB
LB L B
FRB = CLRB*FR + CLRBI*FB
RB R B
The interpolation constants for the points LLB, LBB, RBB, and RRB are
T T 2(c c)
CON1 
AT T2 (1 + cf)( c2)
T1 Cs ( + c 2)
CON2 =  s
T c2(1 + c )
2 2 s
CON3 = 1 CON1
CON4 = 1 CON2
so that the value of F at each of these points is
FLLB = CON1*FL + CON3*FLB
FBB = CON2FLB + CON4*FB
FRRB = CON1*R + CON3*FRB
FRBB = CON2*F + CON4*FB
RBB RB B
For Fully Coupled Waves
When the equations are fully coupled, the equations along the
nonvertical characteristics (c=c cs) are given by either
equation (3.2.7) (3.2.8), or (3.2.9). Equation (3.2.7) will be used,
and the values of AV, A2, A3, A4, A5, and A6 of equation (3.2.2) will
th
be calculated for the i iteration as described for U earlier in this
section and defined as Ali, A2i, A3i A4i A5i, and A6i, respectively.
th
The coefficients will then be defined for the i iteration as
2
IR c f (A .A
If \ 2/ 4i 6i
1 v
2
C
R2f = 2 (AA 6i
1 v
2
c
/ s
R is 2 (A 4iA6i
1 v
2
R2s 2) (2i 6i
1 v
2
 aA.) 
5i
 A3iAi) A2i
3i 5i 2i
2
 aA ) 
5i
 Ai ) A2i
*3i 5i 2i
Rfs= aA iA i Ai A
fs 2i 5i 3i 4i
Thus, the finite difference equations when the waves are fully coupled
can be written directly from equations (3.2.5), (3.2.6) and (3.2.7),
and using the last subscript to represent the point in the numerical
grids of Figures 3.8 and 3.9 the equation
along c = 0 is
2 [SP SB AT = a(V
2a(l v 2AT = a(V2
L 2 rP
rB )
(3.4.2)
along c = 0 is
rP rB A + 2B\
2a VB (2LT) = a 2P )(S S )
2 1 2 ) xP xB
A4P + A4B (S
+ ( 2 (SeP
S()A5P + A5B)( B
 S ) + (A( T )
OB \ 2 P B
1 P xp + B SxB
+ 2  + 2aS (2LT)]
SoP/ s
P
(3.4.3)
(3.4.1)
along c = + c is
Rf
0 f V
Cf xP
+ R f(T p
fs P
c Rfs
xLLB 2 VP
1 v
 V )LLB
eLLB
(1 2)Rf
) + (S S )
LLB 2 xP xLLB
Cf
+ R [lf1 (2Sxp
 aS + 2SLL
GP xLLB
 aS ) 2AT 1
6LLB 1 + cf_
+ FV +V + '1 + (Sxp 2aS +S 2aS
2irP rLLB 1 2 xP xLLB GLLB J
[2 T f
* T
L1 + c j
6c Rfs1 l1
2 2 (Tp + P LLB
1 v
2[rAT
/+f
(3.4.4)
along c = cf is
Rlf
0 = (V
C xP
f
c R
f fs
 V ) + (V V )
xRRB 2 GP RRB
1 v
+ R (TP RR) +
fs P RRB
+ R f 1i
(1 v2)Rf
2 R if(S S )
2 xP xRRB
f
(2Sxp aS p + 2SRRB
xP 6P xRRB
 aS )L cfJ
6RRB 1 + c
+ aR2Vr +VrRRB + (SP 2aSgp +SxRB aSRR)
2
S2T 1 R fs 1 R 2T (3.4.5)
26T* L (T (3.4.5)
1+ c 2 P RRB c
L1 C
67
along c = + c is
s
1s s fs
c xP xLBB 2 eP 6LBB
s lv
(1 V2)R
s
+ Rs 1(2Sxp aS + BB (aS BB+
+ aR V + V + (S 2aS +S 2aSL)
s rP rLBB 2 P P xLBB LBB
6c 2R
2 2T 6csRf 1 T ) j (3.4.6)
+ c 2 L2 TP + LBB Ll+ (3
s 1 s
along c = c is
s
R cR
Is s fs
0 = (V V ) + (V V )
c xP xRBB 2 GP eRBB
s 1v
(1 v2)R
+ (T TRBB) + 2 (S SRBB)
fs P RBB 2 xP xRBB
C
s
+ R11 [+(2Sx aS +2S aS )] [2
1S 1 xP GP xRBB GRBB 1 +c
s
+ aR2s [Vrp VrRBB + 1{(Sx 2aS p+ SxRBB 2aSRBB)}
c 2 R
+C + (T + T ) 1 (3.4.7)
L +c2 2 [ (TP + RBB 1 +c RB J
S 1 v s
where the values of *P / and s / obtained at point P from the previous
iteration are used and
oP roB
S= al + (1 a) s (3.4.8)
sp/ B
with al defined earlier in this section.
For Uncoupled Waves
Then the waves become uncoupled as described ir Section B. l, the
equations along the characteristic lines have a simpler form given by
equations (3.2.10) to (3.2.13). Using the averaging technique already
described in this section for the coefficients, the equations along
the characteristic lines for uncoupled waves simplify. The equation
along c = 0 is
Sp + SB
2 QP GB
2a(1 v )( )(2AT) = a(Vp V ) (3.4.9)
2 rP rB
along c = 0 is
[(A +A A +A
a(V + V )(2T) = a (2 2B)(S S) + 4( )(S s
rP rB 2 xP xB 2 OP GB
S(2aS S )* o (2aS S )P oB
+ a op + 9B B (2LT) (3.4.10)
2L Sp, s B
along c = + c is
2
S(V V ) + (S S
S= x LLB 2 xP xLLB
f c
1 2,T
+ 1 {(2S aS + 2SLLB aS ALLB ] i
frP r P xLLB GLLB 1 C B 4i
+ aA V+V + l(S 2aS +iS 2aS ) 2/T
+ VrPVrLLB+ P G xlLB LLB I
(3.4.11)
along c = c is
0 = A (vPx ) + 12 (Sx S
4 c xP xRRB 2 xP xRRB
f Cf
+ 1{ (2SxP aS e+2SxRRB2aS eRB)} i
+ aA v +V +'{(S 2aS SX 2aS, 1 26T.
2iL rP YrRRB 1 (SxP2aS PxRRB SRRB L
(3.4.12)
along c = + c is
s
c
C
0= (V p V LBB) + (T T LBB
2 P LLBB P LBB
S .(Tp + TLBB) 2 A (3.4.13)
2 L2P
lv s
along c = c is
s
c
s
0 (V e V ) + (T T )
2 eP eRBB P RBB
2
+ 22 (P + RBB) (3.4+
1v s
3.5 Solution to the Finite Difference Equations
The solution to the finite difference equations of Section 3.4 are
given here for any iteration. The solutions consist of expressions for
VxP V p, VrP, Sx, S and Tp in terms of known quantities, including
quantities calculated during a previous iteration. The solutions given
in this section are obtained using Cramer's rule as shown in Appendix D,
and the definitions of the variables used in Appendix D will not be
repeated here.
At a Regular Grid Point for Fully Coupled Waves
The solution to the finite difference equations along the charac
teristic lines at a regular grid point in the case of fully coupled
waves is given here. The longitudinal and transverse velocities from
equations (D.3.5) and (D.3.6) are
1
V = (Ds RHSBA D2 RHSDC) (3.5.1)
1
V (Df RHSDC Ds RHSBA) (3.5.2)
When radial inertia effects are included, the stresses at point P
are given by equations (D.3.13), (D.3.14), and (D.3.15). These stresses
are
T = RHSF(D2 D1D7s) + RHSG(D1D D2D4f)
+ RHSH(D4fD7s D4sD7f)] (3.5.3)
1 r
S RHSF(A D DD ) + RHSG(D2D AD )
xP 2L 5Q 7s 23s 2 3f 5Q 7f
+ RHSH(D3sDf D3fD7s) (3.5.4)
S HSF(D1D A D4s) + RHSG(A5QD D1D)
P =2 53s 5Q 4s Q 4f 3f
+ RHSH(D3fD4s D4fD3s) (3.5.5)
and the radial velocity of equation (D.1.3) is
aVrp = a(D3 Q3S p). (3.5.6)
When radial inertia effects are not included, the hoop stress,
S p, and the radial velocity, VrP, automatically vanish, and the shear
veocty V,
71
stress and the longitudinal stress are given by equations (D.3.18) and
(D.3.19), respectively, as
S= (D RHSF D4f RHSG) (3.5.7)
P 4s 4f
1
S = (D RHSG D RHSF). (3.5.8)
xP 3 3f 3s
3
At a Regular Grid Point for Uncoupled Waves
When the waves are uncoupled, the solution to the finite differ
ence equations has a much simpler form. In this case, the shear stress,
the transverse velocity, and the longitudinal velocity of equations
(D.4.1), (D.4.2), and (D.4.3), respectively, are
T = (RHSCE + RHSDE) (3.5.9)
P 2F2
2s
1
V = f (RHSBEM RHSAEM) .(3.5.11)
xP 2F
When radial inertia effects are included, the longitudinal and
hoop stresses from equations (D.4.8) and (D .4.9) are
1
SP = (D2 RHS3 F5f RHSEEM) (3.5.12)
5
1
S = (F RHSEEM D RHS3) (3.5.13)
OP 6 2f2 1
5
and the radial velocity is again given by equation (3.5.6).
When radial inertia effects are not considered, both the hoop stress
and radial velocity vanish, and the longitudinal stress of equation (D.4.10)
is
RHS3
S S (3.5.14)
xP F22
2f2
At a Boundary Point (X=0) for Fully Coupled Waves
In a boundary element, there are only four characteristic lines
(cO, c= O, c =c c= c ) and consequently only four equations along
these characteristic lines. Since the equations along the character
istic lines are written in terms of six unknown variables at point P,
the solution at a boundary point can be obtained only if two of these
variables are prescribed at each boundary point. The hoop stress and
the radial velocity do not enter the formulation of the problem when
radial inertia effects are omitted, and therefore these variables are
not specified at the boundary. Thus, the four remaining variables, two
of which may be specified at any boundary point, are the longitudinal
stress, the longitudinal velocity, the shear stress, and the transverse
velocity. From a purely physical standpoint, it is also reasonable to
specify the longitudinal and transverse variables at the boundary since
these are the quantities which are normally associated with the impact
at the end of the tube and which can be measured more readily than
radial velocity and hoop stress. Only two of the four variables Sx,
Vx, T and V can be specified at any one boundary point. Furthermore,
at a given boundary point V and S cannot both be specified since they
x x
are not independent. Also, both T and V6 cannot be given at the same
boundary point. Therefore, four combinations of variables to be speci
fied on the boundary will be considered: for Case I, S and T will be
x
given at the boundary, for Case II, Vx and V will be given, for Case III,
S and Ve will be given, and for Case IV, Vx and T will be given. The
solution to the finite difference equations at a boundary point for each
of these four cases when the waves are fully coupled is given below.
Case I: Traction boundary conditions
ihen Sp amd are known, then from equations (D.5.1), (D.5.5)
and (D.5.6) the solution to the finite difference equations at P is
a
Sp D (RHSH DS A5QT p) (3.5.15)
Vp = (B2s RHS1 Bf RHS2) (3.5.16)
4
V = 1 (B RHS2 B RHS1) (3.5.17)
OP A if is
and Vrp is given by (3.5.6).
Case II: Kinematic boundary conditions
When Vxp and Vp are given, the solution at P is given by
equations (D.5.11), (D.5.12), and (D.5.13) when radial inertia effects
are included as
T RHS4(D42 D2 + B7s D1 RHS5(D D2+BfD
+ RHSH(D4s2B7f D4f2B7)] (3.5.18)
Sx [RHS4(D3s D + AD) + RS5(D D + B A)
xP L6 3s2 2 5Q 7s 3f2 2 7f5Q
RHSH(D3s2Bf D3f2B7)] (3.5.19)
S RHS4(D D A D ) RHS5(D D A D
P = RS 3s2 1 A5Q4s2 3f2D1 5Q4f2
+ RHSH(D3f2D4 D4f2D3s2)] (3.5.20)
and Vrp is given by equation (3.5.6). When radial inertia effects are
not included, Vrp and Sp are zero and the solution given by equations
(D.5.16) and (D.5.17) is
T = A (D sRHS4 D RHS5) (3.5.21)
P 4 4s2 4f2
7
1
Sxp = (D3 RHS5 D s2RHS) (3.5.22)
xP A 3f2 3s2
Case III: Mixed boundary conditions
When Sxp and Vp are known, the solution when radial inertia
effects are included is given by equations (D.5.22), (D.5.23), and
(D.5.24) as
S RHS6(D D + AB) RS7(D D + A5Q B 7f)
xP L 3s2 2 5Q7s 3f22 5Q7f
8
+ RHS8(D3s2B7f D3f2B7s) (3.5.23)
Tp F (D2RHS7 + B7sRHS8) B s(D2RHS6 + B RHS8) (3.5.24)
P A Li~f 8 2 2
S 1 Bf(D3s2RHS8 A5QRHS7) B (D3f2RHS8A5RHS6)] (3.5.25)
and Vrp again is found from equation (3.5.6). When radial inertia
effects are not included, VrP= S= 0 and from equations (D.5.27)
and (D.5.28), the solution at P becomes
V = (D 3s2RHS6 D 3f2RHS7) (3.5.26)
xP A 3s2 3f2
p = 1 (B RHS7 B RHS6) (3.5.27)
9 if Is
9
Case IV: Mixed boundary conditions
When Vp and T are known at the boundary, the solution at P is
found from equations (D.5.33), (D.5.34), and (D.5.35) when radial
inertia effects are included to be
i F
LP 10 RHS9(Ds2D + D1B7s) RHS10(D4f2D2 + D1B7f)
7
+ RHSll(D4s2B7s D4f2B7s)J (3.5.28)
SxP 0 B2f(D2RHS10 + B7sRHS11) B2s(D2RHS9 + BsRHS11) (3.5.29)
Sp = LB2f(D4s2RHS11 D1RHS10) B2s(D4f2RHS11 D1RHS9) (3.5.30)
where again Vrp is given by equation (3.5.6). When radial inertia
effects are not included, Vr and S vanish, and from equations (D.5.38)
and (D.5.39), the solution at P is found to be
11
SP = (B2fRHS10 B2sRHS9) (3.5.32)
11
At a Boundary Point (X= 0) for Uncoupled Waves
When the waves are uncoupled, the solutionsto the finite differ
ence equations are obtained at the boundary points for the same four
cases outlined above. When radial inertia terms are included in the
formulation of the problem, the expression for VrP is given by
equation (3.5.6), and in all cases when radial inertia terms are not
included both Vrp and S vanish. In all four cases the solutions can
be found in Appendix D.
Case I: Traction boundary conditions
When T and Sxp are known at a boundary point, then from equations
(D.6.1), (D.6.2), and (D.6.3) at that point
V = (RISDE F Tp) (3.5.33)
GP z9 2s P
a
S = (RHSEEM D S ) (3.5.34)
2
1
VP (RHSBEM aFfSp F2f Sx) (3.5.35)
XP F1 5f P 2f xP
Case II: Kinematic boundary conditions
When Vxp and Vp are prescribed at a boundary point, then the
solution at that boundary point is given by equations (D.6.4), (D.6.9),
(D.6.10), and (D.6.11). When radial inertia is included the solution is
S 1 (D RHS12 Ff RHSEEM) (3.5.36)
12
S P 1 (F RHSEEM D RHS12) (3.5.37)
P 2f 1
12
and when no radial inertia effects are included the solution becomes
RHS12
S (3.5.38)
xP F2f
The shear stress in both cases is
1
T (RHSDE Z V ) (3.5.39)
P F2 2 OP
Case III: Mixed boundary conditions
When Sp and Vp are known at a boundary point, then from
equations (D.6.12), (D.6.13), and (D.6.14), the solution at that point
is
Tp (RHSDE Z2 V ) (3.5.40)
F2s
a
S = D2 (RHSEEM DSpS ) (3.5.41)
1
Vxp (RHSBEM F S aF Sep) (
P 2f xP f P). (3.5.42)
if
Case IV: Mixed boundary conditions
Vihen VXP and are given at a boundary point, then the solution
at that point is given by equations (D.6.9), (D.6.10), (D.6.11) and
(D.6.15) i.e.,
1
p = (RHSDE F T ) (3.5.43)
oP 2 2s P
and Sxp and Sp are given by equations (3.5.36), (3.5.37), and
(3.5.38).
3.6 Calculation of the Strains
At any grid point P, the solution is obtained by an iterative
technique. Once this is done, the values of Sx, Se, T Vx, Vr, and
V are known at P as well as at points L, B, and R (see Figures 3.5
and 3.6). The strains at point P can be computed very easily from
equations (2.27), (2.28), and (2.29). These equations can be written
in dimensionless form using equation (3.2.1) as
x x
x x (3.6.1)
T 2 (3.6.2)
= 2Vr (3.6.3)
For a regular grid element, these equations can be written in
finite difference form as
e V V V
xP xB xR xL
2ZT 2t X
C. C V V
exP OxB 1 (R L
2AT 2 2AX
OP B Vrp VrB
9P eB rP rB
= 2 ( )
2AT 2
where the final subscript on each variable denotes the point in the
grid element where that variable is evaluated. Now, since the outer
dX
grid lines defined in Section 3.3 have slopes of c = 1, AX and
AT are equal so that the expressions for the strains at point P are
C = x+ V V (3.6.4)
xP xB xR XL
1
CxP = exB + (V V ) (3.6.5)
e = CeB + 2(Vrp + VrB) AT (3.6.6)
For a boundary grid element, equations (3.6.1), (3.6.2), and
(3.6.3) can be written in finite difference form as
1
e C V 1 (V + V )
xP xB xR 2 xP xB
2LT 4x
GxP OxB 1 R 2 (Vp + B)
2AT L AX
eP V + V
eP BA2 [frP 2 rB
2 T 2
and again since AX = AT, the strains at the boundary point P are given
by
79
S= xB + 2V V VxB (3.6.7)
xP xB xR xP xB
1
xP = eB + VR 2 (Ve + ) (3.6.8)
e" = 6eB + 2(VrP + VrB) AT (3.6.9)
where equations (3.6.6) and (3.6.9) are the same expression.
CHAPTER 4
RESULTS AND DISCUSSION
4.1 Introduction
In Chapter 2 the problem of inelastic wave propagation was
formulated and the equations for this problem were found. In Chapter 3
these equations were written in finite difference form and from them
expressions for the stresses and the velocities at the points in the
numerical grid (Figure 3.4) were determined. Next a computer code
(shown in Appendix E) was written to facilitate the calculation of
the stresses, velocities, and strains at the grid points in the charac
teristic plane. Now, in this chapter the results obtained by using
this computer code will be discussed for several different combinations
of initial conditions and boundary conditions.
The computer code is written so that the boundary conditions are
specified by reading in values of two variables at each grid point
along the boundary (X=0). By specifying the boundary conditions in
this manner, any variable given as one of the boundary conditions can
have any functional shape. All of the data presented in this chapter
were obtained using the kinematic boundary conditions (Case II), that
is, by assigning values to the longitudinal velocity (V ) and the trans
verse velocity (V ) at the impact end of the tube. Furthermore, the
same functional form was chosen for the two velocities in each case.
This form consists of assuming that each of the velocities at the
boundary increases linearly up to its final value (denoted by Vxf
or VG ) during a period of time called the rise time (T ) and then
remains constant. That is
T
TV Vqf if 0 T < T
T R R
V@ (X = 0) =
Vef if T > TR
T Vxf if 0 < T < T
Tr R R
R
V (X = ) =
Vxf if T > TR
Now that the computer code is set up, it would be advantageous to
compare the results from it to data which have already been published.
This is done in the following section by using the data of Lipkin and
Clifton (1970), and some interesting effects of the size of the numer
ical grid are noted. Then, finally, the effects of radial inertia and
strainrate dependence on the propagation of inelastic stress waves are
discussed.
4.2 Effects of Numerical Grid Size
Lipkin and Clifton (1970) published the results of three different
experiments where a thinwalled tube was given an initial static shear
stress and then impacted longitudinally. In this section the initial
conditions and boundary conditions from one of these experiments will
be used and the results obtained from the computer code will be compared
with the experimental and theoretical results of Lipkin and Clifton
(1970). The data vhich will be used are
0
T O
0
x
= initial static shear stress = 3480 psi
= initial static longitudinal stress = 0
Xf = final longitudinal boundary velocity = 500 ips
vf = final transverse boundary velocity = 23 ips
= rise time = 9.6 p sec
which can be written in terms of the dimensionless quantities for input
to the computer code as
0
T
o 6x
E
= .0003480
0
o0
S x 0
S O 0
x E
Cl
T
R 2r
No radial inertia effects or rate
this section.
v
xf
Xf
V .002404
x
f c1
ef
V .0001106
f c1
t = 4.00
dependence will be considered in
The results from three different computer runs will now be made.
Each computer run used these initial conditions and boundary conditions
but had different grid sizes. The three grid sizes used were
AX=AT = .25, AX=AT= .125, and AX=AT= .05. The longitudinal strain
versus time obtained by using the computer code in Appendix E is shown
in Figure 4.1 along with the experimental results and the simple wave
solution of Lipkin and Clifton (1970). From this it can be seen that
Simple Wave Solution, Lipkin and Clifton (1970)
Experimental Results, Lipkin and Clifton (1970)
 X .05
 X .125
 X= .25
Figure 4.1
u G60
Time, T
Grid Size Effects on the Longitudinal Strain at X = 3.75
.012
.008
.00,4
for the small grid size the strain follows closely the strain obtained
by Lipkin and Clifton (1970) for a simple ;'.ave with ;n instantaneously
applied velocity at the boundary. For the larger grid sizes the strain
versustime curve is smoother and follows more closely the experimental
results of Lipkin and Clifton (1970). Apparently, the larger grid
sizes tend to smooth out the data and eliminate the distinction between
the fast and slow wave speeds. For instance, in Figure 4.1, the' simple
wave solution of Lipkin and Clifton (1970) exhibits a region where the
longitudinal strain has the constant value of 0.00085. The strain
remains at this constant value from just after the fast wave passes
until the arrival of the slow wave.
From these computer runs other quantities of interest can also be
plotted and the same grid size effect can be observed. This is shown
in Figure 4.2 for the longitudinal velocity versus time. The grid size
has a much smaller effect on the stress trajectory than on the time
history curves. The stress trajectory is shown in Figure 4.3.
Because the details of the solution depend on the size of the
numerical grid, all subsequent computer runs will be made using a small
grid. This small grid size necessitates a large amount of computer
time to obtain a solution more than 1.0 diameter from the impact end,
and most of the results given below are obtained near the end of the
tube.
4.3 Effects of Radial Inertia
In order to determine the effects of radial inertia, four separate
computer runs were made using the computer code in Section E.5. The
generalization of the uniaxial stressstrain curve of Lipkin and Clifton
4/'
// X= .05
,/I * X = .125
/ X = .25
J/
20
10
0
Time, T
Figure 4.2 Grid Size Effects on the Longitudinal Velocity at X = 3.75
C
V
x
r4
o
X
o
>,
4(
o
'0
0
11
C3
C,
*H
4 
. X =.05
* X = .125
X = .25
..... 5
.rLr'
Yield Surface
after
Static Preload
6 9
Longitudinal Stress, S x '04
A
Figure 4.3 Grid Size Effects on the Stress Trajectories at X = 3.75
4
2
0
(1970) was used. This constitutive equation (shown in Appendix A) was
for strainrate independent material behavior.
The first two computer runs (one including and one not including
radial inertia effects) were made using the initial conditions and the
boundary conditions which Lipkin and Clifton (1970) used in one of
their experiments. These input data used were
T = 4.00
R
AX = AT = .050
S= 0
x
Data Set 1
T = .0003480
V = .002404
x
f
V = .0001106
1 f
These data represent a tube with an applied static pretorque
(above the yield stress) impacted longitudinally at one end. The time
history curves of the longitudinal strain and the change in shear
strain are shown in Figures 4.4 and 4.5, respectively, for the section
of the tube 3.75 diameters from the impact end. The simple wave solu
tion and the experimental results of Lipkin and Clifton (1970) are
also shown in these figures. It can be seen in Figure 4.4 that the
longitudinal strain obtained in this work follows the experimental
results more closely than does the simple wave solution. Most of the
improvement over the simple wave solution is the result of using
a finite rise time (T = 4.0) for the impact velocity. The fast wave
has passed the point X= 3.75 at the time when the longitudinal strain
I 
i/
/ /
/
 =sc.
Solution Without Radial Inertia
 Solution With Radial Incrtia
Simple Wave Solution Lipkin and Cliflton (i1 70)
 Experimental Results Lipkin and Clifton (1.'70)
Figure 4.4 Longitudinal Strain Versus Time at X = 3.75 for Data Set 1
. 012 
.008
.004_
.7
 I I I
0 20 40 60
Time, T
~
.003
Solution Without Radial Inertia
Solution With Radial Inertia
 Simple Wave Solution Lipkin and Clifton (1970)
 Experimental Results Lipkin and Clifton (1970)
I
 I
4
/ Time, T
. I 1
. 001
Figure 4.5 Change in Shear Strain Versus Time at X = 3.75 for Data Set 1
has reached the value of 0.00085. For the simple wave solution this
time is approximately T 6 Where for the finite rise time (T 41) this
time is approximately T=10. The difference in time when the fast wave
has passed can thus be accounted forby the finite rise time.
As the slow wave passes a point on the tube, the longitudinal com
pressive strain begins to increase to values larger than 0.00085. The
higher levels of strain (e x .008) occur later (in the results given
x
here) than in the simple wave solution. Again this can be accounted for
by the finite rise time.
The inclusion of a finite rise time in the theoretical solution
gives results which resemble the experimental data more closely than
the simple wave solution. It can also be seen that including radial
inertia effects in the formulation of the problem gives longitudinal
strains which are somewhat closer to the experimental data than the
corresponding strains when radial inertia effects are ignored.
The change in shear strain versus time curve in Figure 4.5
exhibits the same rise time effect as the longitudinal strain. The
results obtained here are much closer to the experimental data than
the results for the simple wave solution. The final value of the
shear strain appears to be low. Since this shear strain is calculated
from the values of the transverse velocity, it may be that the final
value of the transverse velocity should be larger.
This can be seen more easily by examining the transverse velocity
at several distance from the impact end as shown in Figure 4.6.
A transverse velocity is induced when the tube is impacted with a
longitudinal velocity, if the tube is statically preloaded in torsion.
From Figure 4.6 it can be seen that the transverse velocity induced
