D
4 6 5
A,
dax+aA2do0 + A3dTe
J
2 i r 
+ I AA_ AA + 2 j aAdcf_ + aA^dcr^ + aA,_dT,
3 5 2 6
pc
2 x 4 0 5 0x
J[]
[p 2J [ai2I5 V j [Vx + aA5doe + A6dT0
j
[
+ I (pc2)(A4A6aA2) A4
'][
p c2 I I aA2A5 A3A4
1] +xdt {
J [2,eJ
A3A5A2A6 +
pc o
147
and regrouping the terms in this equation yields
r
2. 
= rL(PC)(A4A6aA) A
c
dv^ pc
LaV5vJdv
r 2 t t T t2 7 t 7 t2t t r r r2:
(pc)'
+ _(pC) (A1A4A6aA1A+aA2A3A5aA2A6+aA2A3A5;A4)(14a2)jdaj
+ [(Pi2)(W6^5+W5Â¥4V^5'V4V (VVVV]adCTe
[(p;
+ Â¡ (pc ) (A3A4A6aA3A5+aA3A5aA2A5A6+aA2A5A6A3A4A6) (AgA4 aA2A5J d 0X
[<
+  (p c2) (A4AgaA2) aJ 4x dt a j^(p c) (A3A5A2Ag)+A2  lr I dt
f (pc2)[aA2A5A3A4][2^0J dt .
Using the characteristic equation (2.48), the above equation along these
nonvertical characteristics can be written in differential form as
0 r[j
dvxpo_aA2A53jdv0
h [(pH2) (4 a2)jdax+ [a2534JciTex
PC
[(pc2)(4A6aX2)4] [tlrjdt a[(Pc2)(3526)+2] [^Jdt
o
J[2*eJdt
2
+ (p c )
aA2A5_A3A'
(B.2.11)
The second method for obtaining the equations along the nonvertical
characteristics may be shown as follows. Consider a curve in the x t
(characteristic) plane described parametrically by
x = x(Â§)
t = t(Â§)
(B.2.12)
148
and assume that the six variables CT ct, T. v v and v. are known
x 8 0x x r 8
along this curve as functions of the parameter Â§. Then the derivatives
of these functions can be computed as
dv
"Hr = vx,x xÂ§ + vx,t t$
dve
dg ve,x xg + ve,t
dv.
r
= v x,_ + v t, _
r,x Â§ r,t Â§
d a
x
jr* = ct x + CT t,
dg x,x I x,t Â§
da.
"5T CTe,x x,g + CTe,t ti
dT
6x
dc~ = T6x,x XÂ§ + T8x,t tfg
(B.2.13)
where the subscript following the comma denotes differentiation with
respect to that variable, i.e.,
v
x, x
9v
x
=
V
X,
dv
X
t ~ ~5F
dx
Â§ dg
The six equations (B.2.13) and the equations (2.24), (2.25), (2.26),
(2.34), (2.35), and (2.36) form a set of twelve simultaneous partial
differential equations which can be written in matrix form as
Op OOOOIOOOO
OOOapOOOOOOO
OOOOOp 0000 1
1 000000 A O aA 0
X Cj
0 0 0 0 0 0 0 aA2 0 aA4 0
0 0 0 0 1 0 0 A 0 aA^ 0
d 5
x,g t,^ 0 0 0 0 0 0 0 0 0
0 0ax,^at,^0 0 0 0 0 0 0
0 0 0 0 x,_ t,_ 0 0 0 0 0
Â§ ?
0 0 0 0 0 0 x,g t,g 0 0 0
0 0 0 0 0 0 0 0 ax,^ at,g. 0
0000000000 x,,
0
V
X,X
0
0
Vx,t
a cr0
0
r
o
0
V
r, x
0
>1
CO
Vr ,t
Ilf
T X
a5
V0, x
a(r+e)
o
6
V0,t
2^0x
0
a
X X
dv
X
db
0
x, t
dv
r
a dÂ§
0
CTe,x
dv0
dÂ£
0
CTe,t
da
X
dÂ§
0
T0X,X
da0
a di
tÂ§
T0x,t
dT0x
dT
(B.2.14)
149
150
or as
where
MZ = d
(B.2.15)
M =
Op 0000 1 00000
OOOapOOOOOOOO
OOOOOp 0000 1 0
1 000000A OaAOA
12 3
0 0 0 0 0 0 0 aA 0 aA 0 aA_
2 4 5
0 0 0 0 1 0 0 A 0 aA_ 0 A
3 5 6
x,Â£ t,g 0000000000
0 0 ax, a^Â§ 0 0 0 0 0 0 0
0 0 0 0 x,g t,g 0 0 0 0 0 0
0 0 0 0 0 0 x,g t,g 0 0 0 0
0 0 0 0 0 0 0 0 ax,^ at, 0 0
0000000000
XS *Â§
(B.2.16)
z = 1v v 4.t v v V. V.  a a *, aQ *, t .
x,x vx,t vr,x 'r,t v0,x v0,t ^x,x wx,t u0,x w0,t 0x,x T0x,1j
(B.2.17)
and
L aCTe r r N\ dVx dVr dV0 dax da0 dT0x
d = o, o, V HV> _2^9x "Hr a^r 1 "dT dT adT^T
o o
(B.2.18)
151
Now, since all of the components of M and d are known at a given point
(x,t), the solution to equation (B.2.15) can be found directly if the
determinant of M does not vanish. However, if the determinant of M
vanishes, the solution for Z cannot be found. This is represented by
 M j = 0 (B.2.19)
which is the equation for the characteristic lines along which discon
tinuities in the variables may propagate. Evaluation of this determi
nant results in the characteristic equation (2.45). When equation
(B.2.19) is satisfied, the system of equations (B.2.15) will be shown
below to reduce to a set of ordinary differential equations along the
characteristic lines, that is, along each characteristic line one ordi
nary differential equation involving the six variables v^, Vq, v^,
, and t must be satisfied. If M is defined as the matrix M with
0' 0x n
t h
its n column replaced by d, then the solution to equation (B.2.15) can
be written in terms of Cramer's rule as
M
Z = (B.2.20)
11 I M I
t h
where Z is the n element of Z. However, by equation (B.2.19), the
n
denominal;or vanishes and for a solution to equation (B.2.15) to exist,
it is necessary that
 M 1=0 (B.2.21)
i ~n 1
Equations (B.2.21) are twelve equations, one for each column in M which
is replaced. Next, the calculations involved in equation (B.2.21) will
be carried out as an example for the case where n=l. This yields
152
r o p
ac
1 0
M, 
alr*e
o
2^9x 0
dv
df
dv
I
"dT
dv
s
e
d9
da
y.
dÂ§
da
0
dT
0X
d?
0 0 0
0 ap 0
0 0 0
0 0 0
0 0 0
0 01
0 0 0
ax,g at, 0
0 0 x, g.
0 0 0
0 0 0
0 0 0
0
0
p
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
xÂ§
0
0
0 0 0 0
0 0 0 0
0 0 0 1
A 0 aX2 0
a 0 aA. 0
2 4
A 0 a A 0
o o
0 0 0 0
0 0 0 0
0 0 0 0
t, 0 0 0
0 ax,^ at,^ 0
0 0 0 x, ^
A,
0
0
0
0
0
0
or, after some manipulation
153
^ 3 5 r ______ o __ i) __ r)
_M^ ) = a p x,_ p ^A^A,A+2aAAA^AAaAH aA.
46 235 34 15 26
]
4 2 2
3
dVQ
P ^ > p
Xr
dÂ£
4 2
4
da
X
a p t,.
XÂ£
dÂ§
4 2 2
a p t,g
4
x, _
s
[*:
4 2 2
3
dv
X
a p
XF
a
dÂ§
4 ^ 3
2
da
X
p i
XÂ§
dÂ§
aA2A5_A3A4
J
r.
de I AlA4A6+2ail2A3A5_A3A4_aA2A6 &A
3
)(A2A6A3A5)2^0x(A2A5aA3A4)
]
4 3 2 dT 0x
1 ptÂ§ xÂ§ dT L^Wa
J
5V^]+ a4pt,Â§ x,Â§ [^(r^eX] = 
Along the characteristic lines of equation (2.47), x,_ = 0 and this
equation is identically satisfied. Along the characteristic lines of
equation (2.48), t,F is never equal to zero since this would correspond
5
to an infinite wave speed. Therefore, dividing by t,P and noting that
dx/dÂ£ dx
t, ~~ dt/dÂ§ ~ dt
S
= c
(B.2.22)
the above equation becomes
dv
0 = aV2
iV^]
dv
d
]+ if hAvJ
pc(aA2A5A3A4) j + I aAA_AA,
da
x
"dT
[p c2 (^46+2a2352X4a26a142) (A^^]
4/igdrt3^ 4j c r a \ ~ y g jj3 c ^2A6_A3A5^'
+ 2iojpÂ£2(a52s5y4)]t,5} .
154
Using equation (2.48) in this equation yields
4 2f 1 T 2 2
0 = a pc [i[pc (A4A6_aA5)_A4
%
db
*[ 
VrvJar=i [<
pc
[P;2(i446^)iJ 1*. + a[,
pe2(A4A6^)Sj
do
x
"dT
pc (AgAg AgAg) A2j
dt
~~ V df
o
+ 2 [pc2(a2A534)j i6x
This is an ordinary differential equation along the characteristic lines
of equation (2.48). This equation can be written in differential form
by multiplying by the increment dÂ§ along the characteristic line. Doing
this, and noting that the equation is satisfied if the term in the
bracket vanishes, this equation finally reduces to
0 =  dvx_i5 [aV5_vJdve+ [fvvvj^ex
[p^'hV^sh] +xdt
pc
+ a[p2(A2A6A3A5)~2] (/"O^ 2 [pc2(aV5A3A4)]'t'exdt
o
(B.2.23)
which is identical to equation (B.2.11) above. Note that equation (2.48)
represents four characteristic lines and that equation (B.2.23) repre
sents different equations along each of these four lines since along
each line c has a different value. The process just described can be
carried out for values of n from 2 to 12. Doing this, two additional
forms of the equations along the characteristics are found to be
155
  [(pE2)(26A35)2]dvj!4 [<1VV3>5]dve
> A,A !
2
P c
i r ~2 
l''V'lf3' 5j u 9x + ~=2 L(p C ) A2A6" VV _A2J dCTx
(pe2) (1 A AJ.4 I d
[(PE2)(A2A6A3A5)A2]?xdt+ 2 [(DE2)(A1A5A2A3)A5]^0xdt
~2 j^(p c')2(AiA62) (pe2) (A. T
pe
i+V+1] (/^e)dt (B224)
and,
23
O = pc~ [j^2A5A3A4J dvx pe ^(pc2)(A1
Jdve
[^(p C2) (A1A4a2)  J dTfa + 2
5vJ
r
(PC >LaA2A3A3A4J^x
J
2 2
+ (p c )
,aA2A5
t dt + 2(p c)
[(pc2)(A1A4aA2)A4]1J/exdt
+ ap^^p^XA^^j (^*e)dt
(B.2.25)
Equations (B.2.23), (B.2.24), and (B.2.25) are all forms of the same
equation when the waves are fully coupled. Uncoupled waves will be
discussed in Section B.4.
The equations along the vertical characteristics (2.47) can be
written directly from equations (2.25) and (2.35), since along these
characteristics only t varies (x is constant) and partial differential
equations with all partial differentiation with respect to t become
ordinary differential equations. Thus, along the vertical character
istics of equation (2.47), equations (2.25) and (2.35) may be written as
a ov
and
'a dv
0 r
= St
aVr dCTx dCTfi dTPx
= aA + aA  + aA
r 2 dt 4 dt
o
s dt + a*e
156
of in differential form these become
e
dt = ap dv
r
r
o
(B.2.26)
B.3. Reducing Equations to Simpler Case
When radial inertia effects are not included and when the functions
0(s,A) and \f(s,) are obtained from the rate independent incremental
plasticity theory, then the equation (B.1.2) for the characteristics
and the equations (B.2.23) and (B.2.25) along the characteristics
reduce to those given by Clifton (1966) for the von Mises yield condi
tion as will be shown below.
Under these conditions, equation (A.1.2) becomes
2
(B.3.1)
s
and Clifton's expression for k can be written as
Now Clifton defined the function H(k) for the von Mises yield condition
as
H (k)
1
3
157
and this function H(s) can be written in terms of 0(s,A) by using
equations (A.2.16) and (A.2.18) to get
H(s)
H(s) =
3Bn(sa )
y
n1
2
s
30(s)
or
1 2
0 (s ,A) = 0(s) = S H(s) .
(B.3.2)
Using equation (B.3.2) in equation (2.23), and remembering that
a = 0,
where
and
A5 2 CTxTGx H
V 5 + 3Te* H
(B.3.3)
2(1+y) = 1
E G
H = H (s) .
158
Using equation (B.3.3) in equations (2.46)
a = A.
rv1 ax
3Tt
Sx 2 2 2 \
(22 2\1
i
6T
at
+
S!
a
+
 H + a T H 
E x 0x /
(Ve* H )J
a = A.
r~
LEG
a2 3tB
x 0x
+ 3GH + ~E
= ]
 1 x 1 2
b A4 ^ + gH + + 3T H
]
G 0x
and the equation (2.45) for the characteristics becomes after dividing
by A^ (which is always greater than zero)
(P^2)2 ih + 3^H+ ^H) (PS2)(i+^H+ ^+3T0xH) + 1 = O (B34)
which is the characteristic equation given by Clifton (1966).
The equations along the characteristics can be reduced to the two
forms given by Clifton (1966). When a = 0, and ijj = Q=fQ =0,
x 0 0x
equation (B.2.23) becomes
0= \ [(pc2)X4X64]dvx+pc[3jdve]34]dTex+ 4(pc2)(/Ic)jda,
4 6 4J x
or dividing by A
0 = [(p52)6i][ i dvx + A_ dCTJ + a3 [P;civedT6x]
c pc
and using equation (B.3.3), this equation along the characteristic
becomes
[(I
+ 3T H
P c
J + [V0XH][P
ex/~r2j[dc7x~DcdvT H npdvedT
6 9x
J (B.3.
5)
0
159
Also, under the conditions given here, equation (B.2.25) becomes
(p o (A.AJA, dv. +
23 r~
i r
0 = P c _~A3A4
JdvxpcL
[(p52XA1A4)A4JdTexpi2 [vjda,
2
and dividing by (p c ), this becomes
o=[h^][d'
pc
P c dvft J + A3 ^p c dvx doj
6x ^ 0
and again using equation (B.3.3) this is
2
0 = [^2 ^iTH}][pSdv8dTexJ+[v9xH][p3dvxdCTJ
Equations (B.3.5) and (B.3.6) are the equations along the character
istics given by Clifton (1966) for the von Mises yield condition.
B.4. Uncoupled Waves
When the stress waves are fully coupled, the equations along the
characteristics given by equations (B.2.23), (B.2.24), and (B.2.25) are
equivalent. Each of these equations represents one equation along each
of the four nonvertical characteristics. The waves are fully coupled
when none of the coefficients of the variables in these equations vanish.
However, it can be seen that some of the coefficients vanish in these
equations when A =A =0. This condition can occur when loading is
o o
along the TQ axis (a =cto = 0), or when loading is perpendicular to the
OX X W
T axis in the plane T. =0, or when loading occurs at stress levels
0x 0X
within the yield surface (0(s,A) =0), or when unloading occurs (0(s,A) = O)
160
Under these conditions (A = A =0), the equation for the nonvertical
o o
characteristics of equation (2.45) becomes
(pc)2(14A6 aA2A6) (pc2) (^ + a2) + 4 = 0
and as in equation (2.48)
2 2 WA 
C =
2p(X146 aX2X6)
2 ^A1A4 + A4A6 ~ ^2^ ^ iA4^A6 Al^ + aA2^
2pA6(AiA4 aA2)
and using equations (2.49) and (2.50),
2 (A1A4 + A4A6 aA2) + (A4A6 ~ A1A4 + aA2)
2PA6(A1A4 aA2}
2
Cf
A.
p (V4 ~ ^
(B.4.1)
2 (A1A4 + A4A6 aA2) ~ (A4A6 A144 + aA2)
2PA6(A1A4 ^
2 1
c =
PAC
(B.4.2)
and the four nonvertical characteristics are
c = cf, cs.
The equations along the characteristics will be obtained from
equation (B.2.23). For A =A =0, this equation becomes
o O
o = [(P;2)(6)i][a4][
dv + ~ da + i dt
x 2 x Tx J
c pc
[;2)V1][Â¡2]Lr*Jdt
+ a
(B.4.3)
161
However, along the characteristic lines for the slow wave speeds
(c = cg) equation (B.4.3) is identically satisfied by using equa
tion (B.4.2) and no information can be obtained from equation (B.4.3)
along the characteristic lines c = c Along the characteristic lines
for the fast waves (c= c^) given by equation (B.4.1), the equations
(B.4.3) become
0 =
(pcp6lJ [ J [t idvx+ A dax+ dt]
Pc.
+ aj^(p c)Agl
(B.4.4)
In this case when the waves are uncoupled, equation (B.2.23) yields
no information about the variables vA and T. and by inspection, neither
9 9x
does equation (B.2.24). Therefore, we will next investigate the equa
tion along the characteristics of equation (B.2.25). This equation for
A_ =A_= 0 reduces to
O
0 = [(pc2)^^ aA2) T
<Â¡[
2
^Cdve+ dTex + 2p hxdt
*]
(B.4.5)
However, along the characteristic lines for the fast wave speeds
(c = c^.) of equation (B.4.1), equation (B.4.5) is identically satisfied
Along the characteristic lines for the slow wave speeds (c = c ) given
by equation (B.4.2), equation (B.4.5) becomes
0 =
^2' '"'ij r
f sdve + dV + 2fcs't'exdt
]
(B.4.6)
Thus, in the case of uncoupled waves the nonvertical characteristic
lines are given by equations (B.4.1) and (B.4.2) and the equations along
these lines are given by (B.4.4) and (B.4.6) which may be simplified.
162
The equation along c = c is
0 =
dv +
x
~2 dcrx + x dt 1A4
rv.
PC4
+a(A)
L_r
\
dt
and along c = cg is
2
0 = + pcg dv0 + dtQv + 2pco if dt
9x
's 0X
(B.4.7)
J
For the uncoupled waves, the equations along the vertical character
istics (c=0) can be written directly from equation (B.2.26) as
acr.
6 ^ j
dt = an dv
r r
o
v
a( Tlr_ Idt = aA0 dcr + aA. dan
\r Y 0/ 2 x 40.
o J
(B.4.8)
B.5. Elastic Waves
A special case where the waves are uncoupled occurs when 0(s,A) =0.
In this case the waves are elastic and may be either loading at stress
levels within the yield surface or unloading. Since 0(s,A) =0,
equations (2.23) become
(B.5.1)
163
and the equations (B.4.1) and (B.4.2) for the characteristic lines
become
2
Cf
(B.5.2)
and
2 G
c =
S p
(B. 5.3)
which are
(B.5.4)
J
Also, for the elastic case we have iji = ij/ = i(r = 0 and the
x 0 0x
equations (B.4.7) and (B.4.8) along the characteristics further reduce.
The equation along c = c^ is
0 = I + dv + i da j i
I x 2 xl E
c_ P c
,v \
(L)dt
E Vr /
o
(B.5.5)
along c
cs is 0 = + pcs dv0 + dT0x
(B.5.6)
along c = 0 is
aa
dt = ap dv
r r
(B.5.7)
and along c = 0 is
av
r ,, av a
dt = da + da
r E x E 8
where c^ and CÂ£ are given by equation (B.5.4).
(B.5.8)
APPENDIX C
PROGRAMS FOR DETERMINING THE PLASTIC WAVE SPEEDS
In this appendix the listings are given for the two computer
programs used in the first section of Chapter 3. The first program
listed is the one which calculates the plastic wave speeds as functions
of v, 3, y, and 6 for the case when radial inertia is included or for
the case when radial inertia effects are not included. It is also used
to calculate the values of the normalized stresses as functions of y
and 6. The second program listed is the one which calculates the crit
ical value of 8 (as a function of v) at y = 0 for which c =c =c .
I s z
164
c
c
c
c
c
c
c
c
c
c
c
PROGRAM T CALCULATE; THE PLASTIC WAVE SPEEDS
THE EQUATIONS USED ARE (3.1.7), (3.1.18), (3.1.19), AND (3.1.20)
A = ALPHA = 0 FOR CLIFTON'S CASE
A = ALPHA = 1 FOR RADIAL INERTIA EFFECTS
0 E T A = SLOPE OF THE STRESSSTRA IN CURVE IN TENSION DIVIDED BY
YOUNGS MODULUS
NU = POISSONS RATIO
DIMENSION DELTA(20), GAM(IC), A1(10,2C), A2(10,20), A3(1C20)
DIMENSION A4(10,73), A5( 10,20), A6( 10,20), ABAR(10,20)
DIMENSION B B A R(ID 2 0) RAD(10,2?), CFC2(13.20) CSC2(10,20)
DIMENSION SX( 10 23), ST110.20), TAU(ll)
REAL NU
INTEGER A
1 REAC(StlOl) BETA,DEL GAM0,DGAM,DDFL,NU, A
I F ( NU. GT.C.5)GO TO 10
E = SORT(3.0 )
PI = 3.141593
1 = (l./BETA) 1.0
1 = 1
2 GAM(I) = GAMS
C = COS(GAM( I ) )
S = S I N ( G AM ( I ) )
C2 = C * 2
S 2 = S * 2
CS ^ C*S
TAU(I) = S/t
J = 1
CELL = DEL
3 DELTA(J) = ADELL+ (A l)*(PI/3.0)
165
CD = CS(CELTAlJ))
SD = S IN ( DEL T A(J ) )
CD2 = CD *2
SD2 = S D 2
CSC = CD* SD
S X( I ,J) = C *(C D SD/b)
ST( I,J ) = C *(CD + SD/b)
A1 ( I J ) = 1. + .25*Z*C2*(CD2 2.*E*CS0 + 3.C*SD2)
A2(I,J) = NU .25*Z*C2*(3.*S02 CD2 )
A 3 ( I J ) = .5G*E7*CS*(CC E*SD)
A 4 ( I J ) = 1. .2WC2MCD R 2.Â£*CSD + 3.0*SD2)
A5(I,J)=.50*E*Z*CS(CD+E*SD)
A6 ( I J ) = 2.0 ( 1 * iMU) + 3 0 Z S 2
ABAR(I,J)=A1(I,J)*A4(I,J)*A6(I,J)A3(I,J)*A3U,J)*A4
I(I,J)*A3( I,J)* A 5( I,J )Al( I,J )*A5(I,J)*A5(1,J)A2(I ,J)*A2( I ,J)*A6(I
2, J ) )
RBAR(I,J >A4l I,J)*(A6(I,J)+Al(I,J))A*{A2(I J)*A2(I ,J)+A5(I,J)*A5(
1 I J ) )
WAD2 = BBAR( I,J ) *BBARlI,J ) 4.* ABAR( I,J ) *A4(I,J )
IE(RAB2.GT.0.'')G0 TU 13
RAD 3 = RAD2
RAD4 = 110.*SQRT(RAD3)
I F(ORAR( I ,J) .LT.RAD4JG TO 13
RAD { I J ) = .i.O
GO TC 14
13 CONTINUE
RAD( I J ) = SQRT(RAD?)
14 CONTINUE
CFC2 ( I J ) =SORT( ( ( 1. 4UU ) /ABAR( I J ) ) ( BBAR ( I J ) 4RAD { I J) ) )
CSC2(I J)= S Q R T( ((l.+NU)/ABAR(I,J) )*(BBAR( I,J)RAD(I J)) )
DELL= DELL+ DDEL
J = J + 1
IFIJ20) 3,4,4
4 CONTINUE
166
GAM1 = GAMO + UGAM
1 = 1 + 1
I F( I 1 1 ) 2,5,5
5 CONTINU
JJ = J 1
11=11
CO 6 1=1,11
GAM I I ) = GAM(I)* 1 BO./PI
6 CONTINUO
CO 7 J=lfJJ
CELTA!J) = DELTA(J)*1BG./PI
7 C U N T I N U _
WRIT 2!6 102) A,BET A NU
W NIT C(6,103)
WRITE!6,1C4) (GAM( I), 1 = 1, II )
CO d J=1,JJ
WRITE(* 105) DELTA!J),(CFC2!I,J ),1 = 1, I I )
3 CONTINUE
KM IT (6, 106)
HR 1T c{6, 104) !G A M( I ), 1 = 1, I I)
CU 9 J=1,JJ
WRITC(,105 ) DELTA!J),(CSC21I,J ) ,I = 1, I I )
9 CONTINUE
GO TC 1
19 CONTINUE
KRITE!6,107)
WRIF E(6, 1^4 ) (GAM! I ), 1 = 1, I I )
CU 11 J=1,JJ
KRITt(6,105) DELTA(J),(5X(I,J),1=1, II)
11 CONTINUE
WRITE(6, 108 )
WRITtlto,104)!GAM( I), 1 = 1, II )
CU 12 J=i,JJ
WRITE(6, 1C5 ) DEL TA(J),(ST( I J ) 1=1, II)
167
12 CONTINUE
WRITE(6109)
WRITE(6,110)(GAM( I), 1 = 1, I I)
WRITEI6, 111 ) (TAUl I ), 1 = 1, II )
101 F0RMAT(6F1G.4,15)
10? FORMAT('I',5X,'ALPHA = I 1 /6X 'BETA = ,F6.4/6X,NU = ',Fb.4)
103 FORMAT( / / / 5 O X ,FAST PLASTIC WAVE SPEEDS, CF/C2'/)
104 FORMAT(2X,'DELTA*,3X,'GAMMA = 10(F4.1,7X)/)
105 FORMAT(F8.24Xf10FI 1.5)
106 F0KMATI///5JX,SLOW PLASTIC WAVE SPEEDS, CS/C2'/)
107 FORMAT (' 1',50XNORMAL IZED AXIAL STRESS'/)
108 FORMAT (///5jX,'NORMALIZED HOOP STRESS'/)
109 FORMAT(///50X,NORMALIZED SHEAR STRESS'/)
llu FORMAT(10X, 'GAMMA =,10 IF4.1,7X )/)
111 FORMAT(12X,10 FI 1.5)
STOP
END
168
o o n c~. o
THIS PROGRAM CALCULATES THE CRITICAL VALUES OF BETA FuR WHICH
CF/C2 = CS/C2 = l.C AT GAMMA = 0.0
THIS PROGRAM USES EQUATIONS (3.1.7), (3.1.27), (3.1.28), AND
(3.1.29)
DIMENSION DELTA(20) BETA(2u)
REAL NU
INTEGER A, AA(20)
1 REALMS.101) NU,CUTOFF
I F ( CUT Or F G E 1.0 ) GO TO 4
A = I
E = SORT(3.0)
PI = 3.141593
J = 1
CEL = PI/2.
2 CFLTA(J) = A*DEL (A l)*(PI/3.0)
C = CCS( EL T A(J ) )
S = SINIDFLTA(J))
S 2 = S 2
CS = CS
CENOM = (l.+NU)*E*CS + 3.*NU*S?
IF(A.EQ.I)GO TO 5
Z = 1. + 2.Nil
CO TC 6
5 CONTINUE
C = ABS(DENOM)
IF(C.GT.1.0 E4)CO TO 8
BETA(J) = 0.0
GO TO 7
8 CONTINUE
Z = ((1. NU ) *2 )/DENOM
6 CONTINUE
BETA!J) = l./( 1 .*Z )
169
7 CONTINUE
A A ( i ) = A
J J + I
CEL = DEL + PI/1H.0
IF(J.LT.2D)GO TU 2
A = u
IF(J.ED.20)GO TO 2
WRIT E(6,102)
CO 3 J=1,20
CELT A(J) = DELTA!J)*180./PI
WRITE(6103) DELTAJ ) ,AA(J), NU.BETAIJ)
3 CONTINUE
CC TC I
4 CONTINUE
Id FOkMAT( 2 F 10.3 )
102 FORMAT(' IFOR EACH VALUE F DELTA, THE APPROPRIATE VALUE OF Bn TA IS
1 THE ONE BETWEEN 0 AND +1'//' THE VALUc UF BETA GIVEN IS TH~ ON
AT WHICH CF/C2 = CS/C2 = 1.0 AT GAMMA = C//' A = ALPHA = 0 FO
3CLIFT0NS CASE'/' A = ALPHA 1 FOR RADIAL INERTIA EFFECTS////1
4, DELTA',7X,ALPHA*15X, 'NU13X, CRITICAL VALUE OF BETA /)
103 FORMAT IF 20.2,I 13,F2 > 3.F30.6)
STOP
END
o
uj aÂ£ in
APPENDIX D
SOLUTION TO THE FINITE DIFFERENCE EQUATIONS
IN THE CHARACTERISTIC PLANE
In this appendix, the expressions for the stresses and velocities
at a grid point P will be determined. In doing this, many new terms
will be defined. In order to simplify the interpretation of the
computer code in Appendix E, these new terms will be defined in exactly
the same manner in which they are used in this computer code.
D.l. Equations for Fully Coupled Waves
When the waves are fully coupled, the finite difference equations
can be simplified by grouping known quantities together. Thus, letting
Q3 = 2AT(1 v2) ^
D = V Q S
3 rB j 9b
A2Q 2(A2P + A2B^
A4Q 2(A4P + A4B)
A5Q 2(A5P + A5B'>
and
RHSE = >Tj^2V + (S 2aSrsJ ^l + A_S_+A_SftT,+ A_T
I
rB xB
0B s
oB~]
3
B
2Q xB 4Q 0B 5Q B
(D.1.1)
(D.l.2)
the equations (3.4.2) and (3.4.3) along the characteristics (c = 0) can
be written as
aVrP = a(D3"Q3S6P)
(D.l.3)
171
172
and
ir / 2\Jf /
ToP YoP
aRHSE = 2aATV n+ a(Ao^ AT)S +a(A, + a AT)S
rP 2Q Sp/ xP 4Q s /
GP
+ aA T .
5Q P
Now defining the quantities
(D.1.4)
Q4 = 2AT Q3
op'
D, = A AT
1 2Q sp/
lit /
D = A + 2a 5
2 4Q s /
A
AT + Q,
(D.1.5)
and
D4 = 2AT D3
RHSH = RHSE + D,
4 J
the two equations along the vertical characteristics (c = 0) can be
reduced to a single equation by substituting equation (D.1.3) into
equation (D.1.4). This becomes
aDlSxP + aD2S0P + aA5QTP = aRHSH (D.1.6)
The equations along the nonvertical characteristic lines are
given by equations (3.4.4), (3.4.5), (3.4.6), and (3.4.7). These can
be written in a much simpler form by defining the following quantities
as A
Z1 =
Z2 ^
Z3 =
1 v
c
s
2
Z4
1 v
AT
1 
Z5 1 + c,
7 AT
O 1 + C
(D.1.7)
J
173
and
Z7 = VZ5
Z8 VZ6
Z =
6*VZ7
Z10 = 6 Z4'Z8
Z11 = V(2,Rif + aR2f} }
z.o = zc* (2 R + aR0 )
12 8 Is 2s y
and then
If
If
B2f = VRfs
B3f = Rfs*(1V
R,
If
 Z*
4f Z3 11
B5f = V(Rlf+2aB2f> '
B61 = 2'VR2f
y
and
Is
Is c
B2s Z2*Rfs
B3s = Efs(1Z10>
lS
34s = ~Z~ ~ Z12
4
"N
B5s = Vis+2aR2s) (
B = 2*ZC*ROc
6s 6 2s
y
and finally
D3f2 = BfS,(1 + Z9)
D3s2 = BfS<1 + Z10>
Rlf
4f2 ~Z~+ Z11
R
D = + Z
4s2 Z 12
4
The equations along the four nonvertical characteristic lines
(c = c^., cg) can now be written as given below. The equation
along c = + c^ is
BlfVxP B2fV0P + D3f2TP + D4f2SxP aB5fS0P + aB6fVrP
_BlfVxLLB "" B2fV0LLB + B3fTLLB+ B4fSxLLB
+ aB5fS0LLB aB6fVrLLB
(D.1.8)
(D.1.9)
(D.1.10)
(D.1.11)
(D.1.12)
174
along c = c is
BlfVxP + B2fV9P+ D3f2TP + D4f2SxP aB5fS0P + aB6fVrP
BlfVxRRB+ B2fVGRRB+ B3fTRRB+ B4fSxRRB
+ aB5fS0RRB aB6fVrRRB
(D.1.13)
along c = + cs is
" BlsVxp B2sV6P + D3s2TP + D4s2SxP aB5sS0P + aB6sVrP
"B1sVxLBB ~ B2sV6LBB+ B3sTLBB^ B4sSxLBB
+ aB5sS9LBB aB6s^rLBB
(D. 1.14)
and along c=c is
s
B V + B V +D T + D S aB_ S. + aB V
Is xP 2s 0P 3s2 P 4s2 xP 5s 0P 6s rP
BlsVxRBB+ B2sV0RBB+ B3sTRBB+ B4sSxRBB
+ aB5sS0RBB_aB6sVrRBB '
(D.1.15)
Next, let
RHSA BlfVxLLB B2fV0LLB+ B3fTLLB+ B4fSxLLB+ aB5fS0LLB
aB6fVrLLB
RHSB BlfvxRRB+ B2fVeRRB+ B3fTRRB+ B4fSxRRB+ aB5fS0RRB
aB6fVrRRB
RHSC b1svxLBB B2sVeLBB+ B3sTLBB+ B4sSxLBB+ aB5sS0LBB
aB6sVrLBB
RHSD B2sV0RBB+ B3sTRBB+ B4sSxRBB+ B5sS0RBB
HB6s^ rRBB
J
(D.1.16)
175
Now substituting equations (D.1.3) and (D.1.16) into equations (D.1.12)
to (D.1.15) and letting
B7f B5f + Q3B6f
B = B + Q B
7s 5s *3 6s
RHSA1 = RHSA aD B
3 6f
RHSBl = RHSB aD^B
3 6i
RHSCl = RHSC aD B_
3 6s
RHSD1 = RHSD aD B. ,
3 6s s
(D.1.17)
the four equations along the nonvertical characteristics can finally
be written in the form given below. The equation
along c = + c^ is
BlfV B2fV + D312TP + D4f2SxP aB7Â£S8P EHSA1
(D.1.18)
along c = c^ is
Blf'xp + B2f V + D3f2Tp + D4f2SXP aB7fS0P = RHSB1
(D.1.19)
along c = + c^ is
 B V B V + D_ T + D S aB = RHSCl
Is xP 2s 0P 3s2 P 4s2 xP 7s 0P
(D.1.20)
and along c = c^ is
BlsVxP + B2sV0P + D3s2TP + D4s2SxP aB7sS0P
RHSDl .
(D.1.21)
176
The four equations (D.1.18), (D.1.19), (D.1.20), and (D.1.21) along
with equation (D.1.6) now become a set of five simultaneous algebraic
equations which must be solved during each iteration for the unknowns
VxP V0P TP SxP and Sep point P of a regular grid element. The
reason for calculating the coefficients by the method outlined in
Chapter 3 is now apparent. If the coefficients had been calculated in
the normal manner, these five equations would have to be solved simul
taneously for each iteration at a regular grid point when the waves are
fully coupled. However, by calculating these coefficients by the method
used here, these equations reduce to two separate sets of simultaneous
equations, one set with two equations and one set with three equations
(this will be shown in Section D.3), and the computation time required
to solve these two sets of equations is significantly less than the
time to solve one set of five equations. Once the solution is found
at a regular grid point, the radial velocity is obtained from equation
(D.1.3). This will be done in Section D.3.
At a boundary point the solution is obtained by specifying two of
the variables (as described in Chapter 3) and then solving equations
(D.1.19), (D.1.21), and (D.1.6) simultaneously. Again equation (D.1.3)
is used to get V This will be done in Section D.5.
D.2. Equations for Uncoupled Waves
When the waves are uncoupled, the finite difference equations
along the characteristic lines are given in equations (3.4.9) to
(3.4.14). These expressions will now be simplified by grouping known
177
terms together. Using equation (D.1.1), equation (3.4.9) along the
vertical characteristic (c=0), becomes
aVrP = a(D3W
which is the same as equation (D.1.3). Next, using equations
(D.1.1) and (D.1.5), defining
\
S 2aS. i
L2 '
S 2aS
RHSEE = AT 2V + ill
rB s_ toB.
B
J
A2QSxB + A4QS0B
and
(D. 2.2)
RHSEEM= RHSEE + 2T D
J J
and using equation (D.2.1) in equation (3.4.10), the other equation along
the vertical characteristic (c = 0) is
aD S + aD S = aRHSEEM (D.2.3)
1 xP 2 0P .
The finite difference equations along the nonvertical character
istic lines (c = c^, c^) will be simplified by using equations
(D.1.7) and (D.1.8) in equations (3.4.11) to (3.4.14) and noting that
A. is not zero so that the four equations reduce to
4i
v*p+[
1 + 2Z + aZ n
Z3 7 7 A4i'
xP
aZv[1 + 2aSÂ¡i]
A
S. + 2aZ_  V =
GP 5 A . rP
4i
V +
c xLLB
f z (2 + a S + aZ j 1 + 2a Hi 1
LZ3 7 \ A4. /J xLLB 7L A4iJ
0LLB
 2aZ i V
5 A . rLLB
4i
(D.2.4)
178
v + + Z (2 + a ~ )ls aZ j 1
XP LZq 7 V A4i /J xP 7 L
A 1 A .
2aAiJsep+2aZ^v
+ 2a
'4i'
5 A, rP
4i
V + i Z ^2 + a ~ ') I S + aZ
cf xRRB \_Z3 7 V A4i /J xRRB 7
1 S + aZ 1 + 2a =il
yj xRRB 7 L A, J
4i
0RRB
2 aZ  V
5 A rRRB
4i
 z2 Vep + <1 + Z10,TP = Z2Ve.,BB+ (1 zio>tlbb
Z2VSP + (1+Z10)TP= VeRBB+ (1Z10)TRBB
By letting
Flf c
2f
3f
Z7 P + 2a A~
4i
F4f ^ 22
5 A
2i
4i
F = F + Q F
5f 3f ^3 4f
F = 7 (2 + a
6f z3 7l a4.
>
and
F2s = 1 + Z10
F3s = 1 Z10
(D. 2. 5)
(D. 2.6)
(D.2.7)
(D. 2. 8)
(D.2.9)
179
and
RHSAE Fif^xLLB+ F6fSxLLB aF3fS0LLB + aF4fVrLLB
RHSBE = F V + F _S aF S + aF V
If xRRB 6f xRRB 3f 0RRB 4f rRRB
RHSCE = ZV + F T
2 0LBB 3s LBB
RHSDE = Z V. + F T
2 0RBB 3s RBB
and
RHSAEM = RHSAE + aD F
O I
RHSBEM = RHSBE + aD F
o ^xl
> (D.2.10)
y
(D.2.11)
and by substituting equation (D.2.1) into equations (D. 2.4) and (D. 2.5),
the equations (D.2.4), (D.2.5), (D.2.6) and (D.2.7) along the nonvertical
characteristics become as follows. The equation
along c = + c
is
FlfVxP+F2fSxP+aF5fS0P=RHSAEM
(D.2.12)
along c = cp
is
FlfVxP+F2fSxP+aF5fS0P=RHSBEM
(D.2.13)
along c = + c
s
is
 Z2V0P + F2sTP = RHSCE
(D.2.14)
and along c =
 c is
s
Z2V0P + F2sTP=RHSDE *
(D.2.15)
Equations (D.2.3) and (D.2.12) to (D.2.15) form a set of five
simultaneous algebraic equations for the unknowns V^, Vgp, Sxp, sgp
and T Tlieir solution will be obtained at a regular grid point and
a boundary grid point in Sections D.4 and D.6, respectively. In all
cases the radial velocity, V will be found from equation (D.2.1)
after the hoop stress, Sa at the grid point is known.
0P
180
D.3. Solution at a Regular Grid Point
for Fully Coupled Waves
For this case the solution will be obtained by solving equations
(D.1.6) and (D.1.18) to (D.1.21) simultaneously. These equations can
be reduced to a set of two simultaneous equations and a set of three
simultaneous equations as follows. Sutracting equation (D.1.18) from
equation (D.1.19) and subtracting equation (D.1.20) from equation.
(D.1.21) results in two equations involving only the unknowns V and
Vgp. They are
2BlfVxP + 2B2fV6P = RHSB1 RHSA1 (D.3.1)
2EL V + 2B = RHSD1 RHSCl (D.3.2)
Is xP 2s 0P
and if
II
tH
Q
2Blf
D2f "
2B2f
Dls =
2Bls
2s =
2B2s
RHSBA
= RHSBl
 RHSA1
= RHSB
 RHSA
RHSDC
= RHSD1
 RHSCl
= RHSD
 RHSC
(D. 3.3)
then these equations can be written in matrix form as
r .
Dlf
D2f
VxP
RHSBA
^Dls
D2s _
_V9P_
RHSDC
Solving this, the expressions for the longitudinal
velocities at a regular grid point are given by
(D.3.4)
and transverse
V = (D RHSBA D RHSDC) (D.3.5)
xP 2s 2f
V.n = i (D.. RHSDC D RHSBA)
6P A^ If Is
(D.3.6)
181
where
A1 = DlfD2s DlsD2f
(D.3.7)
Two equations involving only stresses can be obtained by adding
equations (D.1.18) and (D.1.19) and by adding equations (D.1.20) and
(D.1.21). Doing this and letting
D3f
2D3f2
D4f =
2D4f2
D7f =
2B7f
3s 
2D3s2
4s =
2D,
4s2
7s =
2B,
7s
RHSF = RHSA1 + RHSBl
RHSG = RHSC1 + RHSD1
these equations become
D3fTF + D4fSxP + aD7fS0P = RHSF
(D. 3.8)
(D. 3. 9)
D T + D S + aD = RHSG
3s P 4s xP 7s 0P
(D.3.10)
The third equation required to solve for the stresses is equation (D.1.6)
which is
aD S + aD S + aA T = aRHSH .
1 xP 2 0P 5Q P
When radial inertia effects are included, that is, when a=1,
these three equations can be written in the matrix form
D3f
D4f
D7f
T
P
rhsf
D3s
4s
D7s
SxP
=
RHSG
_A5Q
D1
D2_
_SeF_
RHSH
(D.3.11)
182
and with
42 = 3f + D7f
the stresses at a regular gird point are given by
TP = ^[RHSF(D2D46D1D7S) + RHSG(DlD7fD2D4f)
+ raSH(D4fDjsD4sD7f)]
'xP h,
LrhSf
+ BiSE<3sD7fV7S)]
6P A2 L
RHSF(JL D A D ) + RHSG(A D D, D )
1 3s 5Q 4s 5Q 4f 1 3f
+ RHSH(D3ffD4sD4fD3s)
]
and the radial velocity, V^p, is given by equation (D.1.3).
When radial inertia effects are not included, that is, when
equation (D.1.6) vanishes, and the two remaining equations (D.3.
(D.3.10) for the stresses, can be written in matrix form as
D3f
I
Q
T
P
RHSF
D3s
I
w
Q
RHSG
and letting
a = D D D D
3 3f 4s 3s 4f
(D.3.12)
(D.3.13)
(D.3.14)
(D.3.15)
a = 0,
9) and
(D.3.16)
(D.3.17)
183
the stresses at a regular grid point are given by
> = 4Â¡(D4sBHSFD4fRHSG)
SxP = iÂ¡
(D.3.18)
(D.3.19)
and the hoop stress, S ,and the radial velocity, V are both zero.
H 0P J rP'
D.4. Solution at a Regular Grid Point
for Uncoupled Waves
For this case equations (D.2.3) and (D.2.12) to (D.2.15) must be
solved simultaneously. Once this is done V can be calculated from
equation (D.2.1). An expression for the shear stress is found by
adding equations (D.2.14) and (D.2.15) and dividing by to get
1 P
RHSCE + RHSDE
2F
(D.4.1)
2s
The tangential velocity is obtained by subtracting equation
(D.2.14) from equation (D.2.15) and dividing by 2Z^ to get
ep
RHSDE RHSCE
2Z
(D. 4.2)
By subtracting equation (D.2.12) from equation (D.2.13) and
dividing this by the expression for the longitudinal velocity
is found to be
V
RHSBEM RHSAEM
xP
2F
(D. 4. 3)
If
184
Adding equations (D.2.12) and (D.2.13) and letting
F = 2F
2f 2 2f
F
5f 2
2F
5f
RHS3 = RHSAEM + RHSBEM J
results in the equation
F2f2SxP + aF5f2SGP RHS3
(D. 4.4)
(D. 4.5)
When radial inertia effects are included (a = l), equations
(D.2.3) and (D.4.5) can be written as
F2f2
F5f2
"SxP~
RHS3
_ D1
D2
_ Sqp _
RHSEEM
(D.4.6)
and by defining
D F
2 2f2
 DlF5f2
(D.4.7)
the expressions for the hoop stress and the longitudinal stress can
be written as
S = j(D RHS3 F RHSEEM) (D.4.8)
xP A 2 5f2
5
S0p = r'(F2f2RHSEEM D!RHS3) (D.4.9)
5
When radial inertia effects are not included (a = 0), equation
(D.2.3) vanishes, the hoop stress and the radial velocity are zero,
and the expression for S is found directly from equation (D.4.5) to be
xP
S
xP
RHS3
F
2f2
(D.4.10)
185
D.5. Solution at a Boundary Point (X=0)
for Fully Coupled Waves
At a boundary point, the solution when the waves are fully coupled
is obtained by specifying two of the variables at the point P and then
using equations (D.1.6), (D.1.19), and (D.1.21). Once again the
radial velocity is calculated using equation (D.1.3) after SQri has been
0P
determined for all cases considered. The solution will now be devel
oped for each of the four sets of boundary conditions discussed in
Chapter 3.
Case I: Traction boundary conditions
When S and T are specified, then the hoop stress can be deter
xP P
mined directly from equation (D.1.6) and is
sep = 5dfflSH  a5(jV .
(D.5.1)
Once this is known, two new quantities can be defined as
RHS1 SB1 D3f2Tp D4f2Sxp + aB7fSep
RHS2 = RHSDl D T D. S + aB S_
3s2 P 4s2 xP 7s 0P
(D. 5.2)
so that equations (D.1.19) and (D.1.21) can be written in matrix form as
Blf
B2f
VxP
RHS1
_ Bls
B2s
1
<
CD
'P
l
RHS2
(D.5.2)
With
*4 =
(D.5.4)
186
the velocities at point P are given by
Vxp = (B2sRHS1 B2fRHS2)
4
V6P = A
4
Case II: Kinematic boundary conditions
When V^p and V are known at a boundary point, by letting
and
RHS4 = RHSB1 
RHS5 = RHSD1 EL V Bn V.
Is xP 2s 0P
the three equations for the stresses from equations (D.1.6), (D.
and (D.1.21) can be written as
D3f2TP + D4f2S*P B7fS6P = RHS4
+ D S aB S. = RHS5 \
3s2 P 4s2 xP 7s 0P /
aA T + aD S n + aD = aRHSH
5Q P 1 xP 2 0P J
When radial inertia effects are included (a=1) these equations
be written in matrix form as
D3f2
D4f2
B7f
r t n
p
" RHS4
D3s2
D4s2
"B7s
SXP
=
RHS5
1
>
Oi
D1
D2
_ s6p _
RHSH
_ _
so that with
A6 D3f2(D2D4s2 + D1B7s) D4f 2 (D2D3s2 + A5QB7s)
B7f(DlD3s2 A5QD4s2)
(D.5.5)
(D.5.6)
(D.5.7)
1.19),
(D. 5.8)
can
(D. 5. 9)
(D.5.10)
187
the expressions for the stresses become
Tp = i [rHS4(D2D4s2+ DiB7s) RHS5(D2D4f2+ W
6
+ RHSH(D4s2B7fD4f2B,
*>]
(D.5.11)
xP _
b
 RI,S4
 RHSH
]
(D.5.12)
Sep = rr[RHS4 RHS5(DlD3f2 A5QD4f2^
4 RHSH(D3f2D4s2 04f2D3s2)j
(D.5.13)
When radial inertia effects are not included (a=0) the radial
velocity and hoop stress vanish as does the third equation of (D.5.8).
The remaining two equations of (D.5.8) become
D3f2
D4f2
1
o,
t
L
RHS4
D3s2
4s2
1
CO
$
L
RHS5
(D.5.14)
and letting
A7 ~ D3f2D4s2 3s2D4f2
(D.5.15)
the two remaining stresses are given by
TP = ;
(D.5.16)
SxP = A7(D3f2RHS5"D3s2RHS4)
(D.5.17)
188
Case III: Mixed boundary conditions
When S^p and VQp are known at a boundary point
RHS6 = RHSB1 B2fVep D4f2Sxp
RHS7 = RHSD1 B V D S
2s 0P 4s2 xP
RHS8 = RHSH D S
1 xP
(D.5.18)
so that equations (D.1.19), (D.1.21), and (D.1.6) can be written as
VxP + D3f2TP aB7fS0P = EHS6
B> V + D T aB S_ = RHS7
Is xP 3s2 P 7s 6P
Vp + aD2sep = aRHSB y
CD. 5.19)
When radial inertia effects are included (a=l) these equations become
[hi
D3f2
_B7f
V
xP
RHS6
:
Bls
D3s2
"B7s
T
P
RHS7
L
A5Q
3 j
r
w
CD
"0
L
RHS8
(D.5.20)
and if
A8 Blf (D2D3s2 + A5QB7s) Bls('D2D3f2 + A5QB7f}
(D.5.21)
the expressions for the unknown variables become
V.. = j^RHS6(DD0_0 + A^BJ RHS7 (DnD0J,0 + A^B,,^)
xP A
8
2 3s2 5Q 7s
2 3f2 5Q 7i'
+ KHS8(D3s2B7f D312B,
*>]
(D.5.22)
189
] r n
Tp = [Blf(D2RHS7 + B7sRHS8) Bls(D2RHS6 + B7fRHS8) J
8 ~
S6P " [Blf
o
(D.5.23)
(D.5.24)
When radial inertia effects are not included (a = 0), the last of
equations (D.5.19) vanishes, Vrp and Sgp are zero, and the first two
equations of (D.5.19) become
B_ ^
D
V
RHS6 j
If
3f 2
xP
B,
Do o
T
RHS7
Is
3s2
P
1
J
so that with
BlfD3s2 BlsD3f2
the solution is
V
xP
T
P
= Â¡r
9
= t^(B, RHS7 B, RHS6)
If Is
Case IV: Mixed boundary conditions
When V^p and Tp are specified at a boundary point,
RHS9 = BHSB1 BlfVxp D3f2Tp >
RHS10 RHSD1 BlsVxp D3s2Tp l
(D.5.25)
(D. 5.26)
(D. 5.27)
(D. 5. 28)
(D. 5. 29)
RHS11 = RHSH AT
o y p
190
so that equations (D.1.19), (D.1.21), and (D.1.6) can be written as
B2fV6P +
D4f2SxP
aB7fS6P =
RHS9
B2sV6P +
D4s2SxP
aB7sS0P =
RHS10
(D. 5. 30)
aDlSxP +
aD2S9P =
aRHSll .
When radial inertia effects are included (a=l) these can be
written as
B2f
D
4f 2
"B7f "
vep
RHS9
B2s
4s2
B7s
SxP
=
RHS10
0
D1
2 _
_S0P_
RHS11
(D.5.31)
and when
S10 =B2f B2s(D2D4f2 + DlB7f>
(D.5.32)
the solution is given by
1
V
9P A. LBHS9(D2D4S2 + D1B7s> 11,5101D2D4f2 + DlB7f>
10
+ RHSll(B7fD4s2 B7S
4f2)J
(D. 5. 33)
'xP ~~ A _B2f(D2RHS1+ b7sRHS11) B2s(d2RHS9 + B?fRHSll) I (D.3.34)
10
= ~ [B2f (D4s0RHS11 D1RHS10) B2s(D4f2RHSll D1RHS9)J (D. 5. 35)
8P
10
When radial inertia effects are not included (a=0), the last of
equations (D.5.30) vanishes as do the hoop stress and the radial
velocity, and the first two equations of (D.5.30) become
191
B2f
4f2
" V
RHS9
_ B2s
D4s2
s
xP
RHS10
so that with
11 B2fD4s2
B2sD4f2
the solution is given by
V
S
0P
xP
4s2
RHS9
2f
RHS10
D4f2RHS10)
B RHS9)
2s
(D. 5.36)
(D. 5.37)
(D. 5. 38)
(D. 5.39)
D.6. Solution at a Boundary Point (X=0)
for Uncoupled Waves
When the waves are uncoupled, the solution at a boundary grid point
is obtained by specifying two of the variables at the point and then
solving simultaneously the three equations (D.2.3), (D.2.13), and
(D.2.15). This will be done below for the four cases discussed in
Chapter 3. In all cases the radial velocity is calculated from equa
tion (D.2.1) after the hoop stress is found.
Case I: Traction boundary conditions
When S^p and are known, equations (D.2.3) and (D.2.15) yield
the expressions
S6P = ^(RHSEEM ~ Vxp) (D.6.1)
V0p = ^(RHSDE F2sV (D.6.2)
2
and once is known, the longitudinal velocity is found from equa
tion (D.2.13) to be
V = (RHSBEM 2F S ^ aF S0,_) .
xP F 2 xP 5f 9P
(D.6.3)
192
Case II: Kinematic boundary conditions
When V and V. are known, equation (D.2.15) becomes
xP SP
Tp = ACRHSDE Z2V0p) (D.6.4)
If
RHS12 = RHSBEM .V _
If xP
then equation (D.2.13) is
F2fSxP + aF5fSGP = RHS12
(D.6.5)
(D.6.6)
and when radial inertia effects are included (a=l), equations
(D.6.6) and (D.2.3) can be written as
F
F_
S
RHS12
2f
5f
xP
D.
D
S.
RHSEEM
1
2
0P
(D.6.7)
By letting
J12
D F
2 2f
DlF5f
(D. 6. 8)
the expressions for the longitudinal and the hoop stresses can be
written as
S = i(D RHS12 F RHSEEM) (D.6.9)
xP A12 2 5f
S = ^(FRHSEEM D..RHS12) (D.6.10)
6P 2f 1
When radial inertia effects are not included (a=0), equation
(D.2.3) vanishes along with S0p and V^p, and the expression for
is found from equation (D.6.6) to be
RHS12
(D.6.11)
193
Case III: Mixed boundary conditions
When S and V. are known at a boundary point, the expressions
xP 9P
for the hoop stress and the shear stress are found from equations
(D.2.3) and (D.2.15) to be
s0p = ^(RHSEEM D1Sxp) (D.6.12)
Tp = (RHSDE Z2V ) (D.6.13)
2s
and once SQ is known, the longitudinal velocity is seen from equation
UP
(D.2.13) to be given by
V = F (D'6'14)
Case IV: Mixed boundary conditions
When V and T are known at a boundary point, using equation
xP P
(D.2.15), the tangential velocity is given by
V0P = ZÂ¡(RHSDE ~ F2STP>
Using equation (D.6.5), equation (D.2.13) can be written as
(D.6.15)
F_S + aF S. = RHS12 .
2f xP 5f 0P
However, this equation and equation (D.2.3) are the same two equations
solved for Case II above. Therefore, the longitudinal stress and the
hoop stress are given by equations (D.6.9) and (D.6.10) when radial
inertia effects are included (a = l), and when radial inertia effects
are not included (a = 0) the hoop stress vanishes and the longitudinal
stress is given by equation (D.6.11).
APPENDIX E
COMPUTER PROGRAM FOR CHARACTERISTIC PLANE SOLUTION
E.l. General Description of the Program
This program is written to solve the finite difference equations of
Chapter 3. Two types of grid elements are used: the regular grid ele
ment shown in Figure 3.5 and the boundary grid element shown in Fig
ure 3.6. These two types of elements are fitted together to form the
complete grid network of Figure 3.1, and the solutionsto the finite
difference equations are obtained at each point in this grid network
using equations developed in Appendix D and rewritten in Chapter 3.
To obtain the solution at a grid point P, however, the values of each
variable must be known at the points L, B, and R for a regular element
and at the points B and R for a boundary element. For this reason,
some initial conditions must be specified.
The initial conditions are specified in the region
X S T
since the fastest elastic wave can propagate with a speed no greater
than
For the problem considered here, this amounts to specifying all of the
variables at each grid point along the line
X = T,
that is, at the points (0,0), (AX,AT) (2AX,2AT) (a^AX.a^AT).
The initial conditions used are discussed in the next section.
194
195
Once the initial conditions are known, the solution can be found
at each grid point along the line
T = X + 2AT
beginning with the point (0,2AT). Once this solution is known, the
solution can be found at the point (AX,3AT) and then at (2AX,4AT),
(3AX, 5AT) ,...,( (a^l)AX, (aVl)AT) The solution is then found along
the line
T = X + 4AT
in the same manner, and then along succeeding lines until the solution
is found at the point (0,2a^AT).
At each grid point (except along the line X = T) the solution
must be obtained by an iterative technique using the equations from
Chapter 3 (or Appendix D). Since some of the coefficients in these
equations are determined from the values of some of the functions cal
culated at P for the previous iteration, these coefficients must be
specified for the first iteration. This is done in every case by assum
ing that the waves are elastic, that is, by assuming on the first iter
ation (i=1) that
$(s,A) = (s A) = 0
and
Ali = 1
A.. = 1
4i
A2i 
A = 0
3i
A = 0
5i
A = 2(1+ v) .
6i
Using these values, the solution can be found at P. The values
variables obtained at P in this manner are then used to compute
of the
the
values of the coefficients for the second iteration to obtain a second
196
solution at P. This is done until two successive solutions at P differ
by some value which is specified in the input. Once this convergence
is reached, this iterative technique is begun at the next grid point.
E.2. Initial Conditions
Initially the tube is assumed to be at rest, loaded by a state of
constant stress (and, therefore, constant strain). Designating the
initial state by a superscript o, the initial conditions are given in
terms of the dimensionless variables as
V = V = V = 0
x 0 r
and since the body is at rest
if = 0
o
and from equation (2.8),
s.o
It is assumed that initially the tube is in a state of static
prestress (T = T, S = S) and that this stress state is reached with
x x
out any unloading so that the strains can be uniquely determined.
Furthermore, to simplify matters, it will be assumed that this initial
state of stress (and strain) is reached by proportional loading so that
if the initial stress state is given by a?, (or s. for the deviatoric
ij iJ
stress), then the stress state at any time during the static preloading
is given by
a. = C(\)o. .
iJ iJ
s = C(X)s.
ij ij
\
J
(E.2.1)
197
where \ 'is a parameter along the static loading path and
ij ij iJ
o
and s. are stresses which have not been made dimensionless. The incre
ij
ment of plastic strain during this static loading is given by equation
(A.2.4) as
de. = s d\
ij ij
or
deP. = C(\)s?.d\
ij ij
and the plastic strain along the loading path becomes
P or1
e. = s. C(\)d\
ij iJ ^
where \ is a dummy variable, or
p = _H.
'ij C(\)
C(\)d\
and letting
U =
C(\)
C(\)d\
this plastic strain becomes
Now,
P
e. = cus .
13 ij
P P 2
e e = (ju s..s..
ij ij id ij
and from equation (2.15)
P P 2 22
e. e. = 77 ul) s
ij ij 3
or solving for (ju
P P
e. e. .
u = I
2 2
s
(E. 2.2)
(E. 2.3)
198
Equation (2.16) defines the plastic part of the strain as
or
and from equation
ing becomes
s
(E.2.2), the plastic strain during the static preload
P
e .
ij
s .
ij
(E. 2.4)
199
For the static case, A can be obtained from equations (A.2.13)
and (A.2.21) directly from the two uniaxial stressstrain curves.
For the generalization of the Lipkin and Clifton (1970) curve,
P / vn
A = B (s s^> (E. 2.5)
and for the generalization of the Cristescu (1972) curve,
1
2es
z
 (s s ) + A s
p' y y
sJ
6)
If the only stresses present initially are ct and T then
X 0x
o2 o2
X + 31 ex
and
o / o2 o2
s = / S + 3t
AJ X
(E.2.7)
and A is calculated from either equation (E.2.5) or equation (E.2.6)
as
or
and
AP
A = X
,Po / o vn
A B (s s >
y
i
p2e2 ~i [ / oxz i
1 [ir (s y + \ s]*2 [(fry SJ
P k
> (E.2.8)
o Po o
A = A + s
J
and the initial plastic strains are given by equation (E.2.4) as
_Po 3 AP /Sij
ij 2 sO V E
200
so that
.Po
3
AP
(2
0
o
<3
S
X
2
o
s
\3
o
s
X
Po
3
AP
/
(
1 s)
1
AP
0 "
2
o
s
V
3 x/
2
o
s
Po
:0x ~
3
2
APO
o
o
T
A
(E.2. 9)
J
and the total prestrains are then given by
o o Po
e = S + e
xxx
o o Po
e9 = vSx + Ee
o ... o Po
e0x <1 + V)T + e0x
>1
J
(E.2.10)
For the two stressstrain curves considered here, the expressions
for $(s,A) are given in equations (A.4.3) and (A.4.7). Therefore, the
initial values of Â§(s,A) are given by Â§(s) for these two cases as
x. o. / o vnl
s(s) = Bn(s s )
and
(E.2.11)
i
(~2 i r O 
I**,
respectively.
201
E.3. Calculation of A
In order to find the solution (stresses and velocities) at a grid
point P, a value of A must be calculated during each iteration. This
can be done by calculating the strains at point P from the velocities
(as shown in Section 3.6) and then calculating the value of A from the
strains at points P, R, B, and L (see Figures 3.5 and 3.6). The value
of A can also be found from the functions Â§(s,A) and \[r^(s,A) by eval
uating these functions using the values of s and A obtained during the
previous iteration. The second method will be used so the strains will
not have to be calculated during each iteration, and the finite differ
ence expression for the value of A at the point P will now be found.
From equation (2.14), the inelastic strain rate is
s
and then
2
J
From equation (2.16)
and using equation (3.2.1), the dimensionless form becomes
J^$(s,A) H+2i)ro(s,A)JdT+ s
(E.3.1)
202
or
dA ,, ds ds
dT Â§(s,A) dT + 2ij.'o(s,A) + dT
and omitting the functional dependence, this is
dA _. ds
dT = <+1) di+ 2V
(E.3.2)
Along the vertical characteristics from point B to point P, this
equation can be written in finite difference form as
 rs s
AP AB
2AT
il
($p + 1) + ($B +
1} J [W5]* I [2*op+2U
or,
, A r, ?p + b
ap = ,sb + L1 +
0[SpSB] + [*oP+toB] W '
(E. 3. 3)
This expression is used to calculate the value of A at point P during
each iteration.
If A had been calculated from the plastic strain rates, it would
P *p
have been necessary to include expressions for e and e .
r e
E.4. Input Data
In each computer run, the grid size (AX and AT) must be specified
as well as Poisson's ratio (v), the rise time for the pulse at the
boundary (T = (JRISE 1)AT), the maximum number of grid points in the
X direction (MXMAX), and whether or not radial inertia effects and/or
strainrate dependence is included. Other data which are used are
obtained directly from Lipkin and Clifton (1970) and Case XII of
Cristescu (1972).
When using the data of Lipkin and Clifton (1970), the generalized
stressstrain curv becomes
203
A = s+ EÂ¡(ss)n
y
(E.4.1)
in terms of the dimensionless variables, where
n = 1.923 ^
v= 33^si = 3 39 x1Q4
y 10 psi
 n 10 7 1 923 4
B = BE =(6.1313 xlO ) (10 ) = 1.77236 x10 )
(E.4.2)
When using the data of Cristescu (1972), the generalized stress
strain curve is found from equation (A.2.21) to be
i
2ez nr
7 s \2 1
(s s ) + A s + xo
L q' y y J 2 L
P
& SJ
(E.4.3)
where
and
s
y
a
ys
E
1068.24 psi
6
10.2 X 10 psi
1.047294 x 10
When the data of Cristescu (1972) are used, several dimensionless
quantities can be defined for Case XII. This case is for dynamic load
ing only, and the equations for Â§(s,A) and ^(SjA) are given in
Section A.4. Other dimensionless quantities for this case are
204
O
S = #
y e
1500 psi
10.2 X 10 psi
4
1.47059 x 10
4
A = e = 1.47059 x 10
y y
A = e_
.002575
P =
5.6 X 10 psi
10.2 x 106 psi
, 0054902
A = 6 = .0004
m 3.25 x 10 psi
xm = =
E 10.2 X 10 psi
.0031863
(E.4.4)
xn
4
n 9 x 10 psi
J7 0
10.2 X 10 psi
= .0088235
h = .000004
X = .00005
r k 3 1.
xko = 4^ = (2.5 1n)(A. ec. 2= .0012
c 2.08 x 10 in/sec
where the dimensionless coefficient for the coefficient in the rate
dependent term was obtained for a .5inch diameter aluminum tube.
E.5. Listing of the Program
In this section is the listing of the computer program which was
used to obtain the solution in the characteristic plane to the wave
propagation problem presented in this paper.
o o r> o n n r o o rÂ¡ o r. r< n o o
DIMENSION SX ( 601 ) ST( 601 ) ,TAU( 6.' 1 ) VX( 61 ) VTI601 ) VR (601 ) ,trX(60l )
1 F;T ( CU ) FTX ( 6'U ) 0KLTA ( 601 ) HPLAS( 6'Jl ) ,CF ( oOl ) ,CS ( 6M ) t SB (6n > ,
2E X P(6? 1 ) ,FTP(6 3 1 ) ,uTXP(601 ), PSI' (601),PHIU(601) ,1 TeK(601) ,Al (601 ) f
3 A2 ( 60 1 ) f A 3 I 6C l ) A4( 60 1 ) A5( 60 1 ) A6 ( 601 ) X ( 601 ) T ( 60 1) SMAX(6rl )
DIMENSION DELLHB(120),DELGAM(6'>1),DEX(601)
OI MANSION STRES(601)
DIMENSION SiP( 5? ),STRAIN( 50),VXP( 50) VTP( 50
PEAL N,riU
INTEGER A,UNCOUP
C
C UNCOUP = 1 WHEN THE TRANSVERSE ANO THE LONGITUDINAL WAVES.ARE UNCOUPLED
UNCOUP = O WHEN THE TRANSVERSE AND THE LONGITUDINAL WAVES ARE CUUPLED
A O POR NO RADIAL INERTIA EFFECTS
A = 1 FOR RADIAL INERTIA EFFECTS
JRATE = 0 FOP RATE INUEPENOENT CASE
JRATE = 1 FOR RATE DEPENDENT CASE
JDLL = C FOR CLIFTONS DATA
JBELL = 1 FOR BELLS DATA
IN NCNDIMENS 10NALIZINC THE PARAMETERS THE VALUES OF YOUNGS
MODULUS WHICH WERE USED WERE 
>z 1", 200, roo PSI FOR BELLS DATA
C = lO,v0O,''O0 PSI FUR CLIFTONS DATA
kL Afi (5,1C02) UNCOUP
Rf A0(5 1 Or1) A,JR AT 6,MXM AX,NU,H 1,H2,M3,U4,DEL X,PELT,N XKO
R t A D ( 5,1 0 C 2 ) JR ISE,XVFIN,TVFIN ,SXO,TAU j *SMALLBBARSY
REAL' ( 5,10f 2 ) KASE,SYS
PE All (5,10 02) JBELL,DZ,BETA, CHAT,XM,XN,H,XLAM
KEAD(5,lCn3) INCRX, I NCRT, I 0,MI,MJ,I PUNCH,I PUNI,IPUN2,INCRXP,INCRTP
W R I T E ( 6 ,1005) A,JRATE,DELX,CELT,MXMAX,SMALL
205
WRITE!6, 1030 ) JGELLtDZ fOHAT ,XM,XNXLAM
WRI T E ( 6 100 6 ) NU.SY,BDAR,N,XKC t RETA
WRI r E ( 6 100 7 ) JRISEfXVFINt TVFIN .SXJ.TAUO.H
WRITE(6,1006) H1.H2.H3.F4,I.SYS
WRI r::(6, 1' 3 1 ) MI.MJ, I PUNCH KASE
IF(UNCUUP.GT.O) WRITE(6,1C39)
IF! JNCGUP.LE.O) WRITE!6,1040)
I F( I PUNCH.EQ.O)GO TO 15
WRITZ ( 7, 1034) A, MXMAX,I PUNI,I PUN 2,INCRXPINCRTP,DELX.OELT
15 CONTINUE
C THIS BEGINS THE CALCULATION OF CONSTANTS USED
M M A X = 2 *MX MAX
NMAX = MMAX 1
1M A X = NMAX 2
Q = 1. Nil * 2
C2 = SCR T(.3*(1 .NU ) )
uA = 1. + NU
CELT2 = 2.* D E L T
I 2 A = 2 A
Cl = 2 C A
C2 = C/2.
C3 = Q*D ElT 2
C4 = 03*CELT 2
FV1 1. + C 2
CV2 = 1. C2
CL R B = 2 *C2/l)V 1
CLRPI = 1. CLRB
C Y = ABS(SY)
CO 30 L= 1 MX MAX
3J XIL) = (L 1 ) *CELX
C THIS BEGINS THE CALCULATION OF INITIAL CONDITIONS
LINE =
LINER = I NCRTP 1
LINEO = I NCR T 1
206
S BO = SORT (SXO* 2 + 3.*TAUC**2)
CALL PH II PHIP,PHIQtOJBELL, SBOt DCLTAO SYS,BBAR,N,H,XLAM,BETA,DY,XM
1fXNtDZ,0)
IF(S BO.GT .1. E0 8 ) GO TO 60
FI = 0.J
F 2 = 0.0
GO TO 65
60 CUNUNUL
Fl = sxo/suo
F 2 = TAUO/SBC
66 CONTINUE
G1 = PHIP*F1**2
G2 PHIP *F 2 * 2
C 3 = PHIPF1*F2
A10 = I. + GL
A 2 0 = (NO 5C*G1 )
A 30 = 3 G3
A 40 = I. + 2 6 G1
A50 = 1.5 G 3
A6'j = Q l + 9 G 2
A 2 0 2 = A 2 C * 2
A 5 0 2 = A 5 ~ 2
IHUNCOOP.GT .0) GO TO 67
AO = A40(A10*Ao''A30**2)+A*(2.*A20*A30*A5CAin*A502A6,,*A2 0?)
P. = A 4 3 ( A 1 0 + A 60 ) A A 5 C 2 + A2G2)
COE F = 02/AO
PCCr = BC * 2 4. A 4 G A 0
RAO = SORT!ROOT)
CFO = SORT{CEF *(BO + RAO))
CSJ = SORT! COEFM BO RAD))
GO TO 6B
67 CONTINUE
CFO S0RT(0A4?/(A10*A40 A*A202)>
CSO = SORT!0/A60 )
207
60 CONTINUE
EXPO = SXO*PHIQ
ETP3 = .5*A*EXP1
C PL A SO = PEL TAG SBO
ETXPO = 1.5*TAUP*PHIQ
E Xu = SXO + EXPO
ETu = A*lET PO NU*SXO)
ETX1 = Q A T A U C * ETXPO
CG 7 0 MX = 1MXMAX
T(MX ) = (MX 1 ) *OELT
TAU(MX) = TAUO
SX(MX ) = SXO
ST(MX ) = 0.0
VX(MX) = 0.0
VT(MX) = 0.0
VR(MX) = 0.0
SB(MX) = SB.)
CF(MX) = CFG
CS(VX) = CSu
A 1(MX) = All
AP(MX) = A2 0
A3(MX) = A30
A 4(MX) = A41
AS(MX) = A 50
A6(MX ) = A6 0
EX(MX) = HXO
cT(MX) = ETC
ETX(MX) = FJXO
EXP(MX) = EXPO
ETP(MX) = E T PO
ETXP(MX) = ETXPO
ITER(MX) 1
PHIO(MX) = PH IP
PSIO(MX) = 0.0
208
CEL r A(MX) = CELTAO
CPLAS(MX) = DPL AS')
SMAX(MX) = snc
CEL GAM(MX) = 0.0
C EX(MX) = 0.0
STKES(MX) = SORT(SX0**2)
70 CONTINUE
CELLHR(l) = CELTAO
C THIS ENOS THE CALCULATION OF INITIAL CONDITIONS
I W R I T E = 1
GO TU A01
9^ CONTINUE
KUM = 1
100 CONTINUE
MX = 1
JERROR 0
CO TO 300
200 CONTINUE
MX = MX + 1
'30 3 CONTINUE
MT = MX + KLJUNT
T(MX) = M T D E L T
SXk = SX(MXfl)
STR = ST(MX + 1 )
V X R = VX(MX 1 )
VTR = VT(MX+1)
VRR = VRIMX+l)
TAUR = TAU(MX + 1 )
T AUG = TAU(MX)
SXC = S X(MX)
STB = ST I MX)
VXD = VX{MX)
VT B = VT(MX)
VRB = VR(MX)
to
o
IB = Al (MX)
3n 1
302
A2B = A2(MX)
A 3 B A3(MX)
A 4 B A4(MX)
A 5 B = A5 (MX)
A6B = A6(MX)
SBB = SB(MX)
P S I 3 = PSIO(MX)
PMIB = PHIO(MX)
SMAXB = SMAX(MX)
LEL TAB = DELT A(MX )
1F(SD3.LT.l.UE08)G0 TU 301
PSI3B = PSIB/SBC
GO TG 332
CONTINUE
PSIBB =0.0
CONTINUE
PSIBC = 375 *PSIBB
C3 = VR13 U3*STB
C4 = CELT2*D3
A 1 0 C :
A 2 B C =
A3BC =
A4BC
A 5 3 C =
A6BC =
T AUKR
c X k 3
STR3 =
VXRQ =
VTRB :
VRKB
I F ( MX
IFIMT
. 375*A1B
. 375 A2B
. 3 7 5 A 3 B
: 3 7 5 A 4 B
. 375 A53
. 3 75 A6D
= CLRB*TAUR
: CLRB*SXR +
: ULRB*STR +
: CLRB*VXR +
= CLRB*VTR +
: CLRB*VRR +
.GT.DGU TU
+ CLRR I *TAUB
CLREI*SXB
CLRBI*STB
CLRBIVXC
CLRBI* V T B
CLRBI*VRB
310
GT.JR ISE)GO TO 303
210
30 3
3 0 4
305
3 1C
IFCKASE.EQ.1) READ!5,10 C 4 ) SXI.TAUI
IFIKASE.EQ.2)
IF(K AS L .EQ.3)
IF(KASE.EQ.4)
GO TO 334
CONFINUE
IF(KASE.EQ.I)
IF(KASE.EQ.1)
IF(KASE.EQ.2)
IF(KASE.EQ.2)
IF(KASE.EQ.3)
IF(KASE.EQ.3)
IF(KASE.EQ.4)
IFIKASE.EQ.4)
R F AD( S, 100 4)
READ! 5 1004)
READ(5, 1004)
SXI = XVFIN
TAUI = TVFIN
VXI = XVFIN
VII = TVF1N
SXI = XVFIN
V r I = TVFIN
VXI = XVFIN
TAUI = TVFIN
CGNT I NUE
IF(MT.GT .2)GO TO 30 5
DEL TAP DELTAB
GO TO 32^
VXI,VTI
SXI,VTI
VXI,TAUI
DELTAP = 2.* DEL T AB DELLHB(MT3)
CG 10 320
CONTINUE
sxl = sxirxi)
STL = ST(MX1)
VXL = VX(MX 1)
VTL = VT(MXl)
VRL = VR(MXl)
TAUL = TAU(MXl)
CELTAP = DELTAIMX+1) + DELTA(MXl)
TAULB = CLRB*TA'JL + CLRB I *T AUB
SXLB = CLRB*SXL + CLRBI*SXB
STLO = CLRB*STL + CLRBI*STB
VXLB = CLRB*VXL + CLRCI*VXR
VTL B = CLRB*VTL CLRBI *VTB
VRL B = CLRO*VRL + CLRBI*VRB
DELTAB
320 CONTINUE
I = 1
A1P = 1.0
A2 P = NU
A 3 P = 0.0
A 4 P = 1.0
A5P = 0.0
A6P = Q1
PHIP = 0.0
CO TC 355
34t CONTINUO
1 = 1 + 1
I F { I.LE.MDGO TG 344
IF(JERROR.GT.0)GO TO 341
JERRCR = JERROR + 1
WRIT E(6 1033)
341 WRIT E(6110 3 2) X(MX ) ,T(MX ),DIFF,DENOM,ERROR
GU TC 430
344 CONTINUE
CALL PHI(PHIP,PH IQ,JRATE,JQuLL, 50PC,DEL TAP, SY,BtfAR,N, H XLAM, 3ETA, D
1Y,XM,XN,DZ,10)
I E ( S BPC. GT. SMAX3.UR I .GE .MJ )GU TO 351
PHIP = 0.0
351 CONTINUE
IE(SBPC.GT.1.CE08)GO TO 353
FI = 0.0
F 2 = 0
F 3 = 0.0
GO TC 354
353 CONTINUE
FI = (2.*SX I A S T I J/SCPC
F 2 = (SXI I 2A*ST I ) /SBPC
F 3 = TAU1/S3PC
354 CONTINUE
212
Y 1 = 2 5 PH I P
Y2 = Yl* F1
Y 3 = PH IP *F 3
YA = 1.5 Y3
A 1P = 1.0 + Y 2 F1
A 2 P = { NU + Y 2 F 2 )
A3P = Y A F1
AA P = 1.0 + Y1*F2**2
A5 P = Y A *F 2
A6P = Q1 + 9.*Y3*F3
355 CONTINUE
!F( PHIFJ.LT. L.0EO5.OR.PHIP.LT.1.0E05JG TO 356
1 F( I.CT.2 ) GO TO 360
356 CONTINUE
All

Al P
A 2 I
A2P
A3 I
=
A3 P
AA I
A'+P
A 5 I
A 5 P
A 6 I

A6P
A2C

A2P
m A (.
=
AtP
A 5 C
=
A5P
CO
TO
365
CONTINUE
All
=
.625*A1P
A IRC
A2 I
=
.625*A2P
+
A2BC
A 3 I

.625*A3P
+
A3 OC
AA I
=
625*AAP
+
AA RC
A 5 I
,625*A5P
f
A5BC
A 6 I
=
.625*A6P
+
A6BC
A2C
=
. 5 ( A 2 P
+
A2B )
AA 0

. 5 ( A A P
+
AAB )
A5G
=
. 5 ( A 5 P
+
A5B )
213
365 CONTINUE
366
367
RHSc = DLLTM2. *VRBPSIR0*( 12A*5TBSXB) ) + A2U SXB + A4Q* S TB + A50 TAU6
RHSH = RHSE 04
A 2 I 2 = A 2 I * 2
A512 = A 5 I *2
IF(UNCOUP.GT.0) GO TO 366
A3 = A4I(A1I*A6IA3I**2)+A*(2.*A2I*A3I*A5IA1I*A5I2A6*A213)
BR = A4IMA11 + A61 } A*(A2I2 + A 51 2)
C02F = Q2/A3
ROOT = 03**2 4.*A4l*AB
RAC = SURTI ROOT)
CFI = SURTI COFFM BB + R A D ) )
CSI = SURTICUEFMBB RAO))
GO TC 367
CONTINUE
CFI = SURTIQ*A4I/I AlIA4I A A 2I2) )
CSI = SURTIQ/A6I)
CUNTINUc
CV 3 = 1. + CFI
CV4 = I. + CSI
CUN I = 2 IC FI C2 ) /IDV2*0V3)
C0N2 = 0V1*CSI/I C2*0V4)
CON 3 = L. CON 1
CON4 = 1. CON 2
TAURRO = CON I T AUR + C0N3*TAUR0
SXRRB = C ON I SX R + CN3*SXR0
STKRB = CON I*ST R + C0N3STRB
VXkRC = C ON1 V X R + CN3*VXR0
VTRRB = CON 1*V T R + C0N3*VTRB
VRRRB = CONI*VRR + C0N3*VRKB
TAURCB = C0N2*TAURB + C0N4*TAUB
SXRBB = C0N2*SXRB + C0N4*SXB
STREG = C0N2*STR8 + C0N4*STB
VXRBB = C CN2 *V X R B + C0N4*VXB
214
VTRBB = CON2 *VTRB CON4*VTB
VRKD3 = CON2*VRRB + CON4*VR8
Z1 = CFI/O
l = CSI/Q
Z 3 = Z 1 C F I
Z 4 = Z 2 C S I
Z 5 = CELT/DV3
Z = DELT/DV4
RFS = A*A2I*A5I A 3 I *A4 I
RFSA = A BS(RFS)
IF(RFSA.LT. 1 .OF^6.OR.UNCOUP.GT.O) GO TO 360
CVS = A4 I *A6 I A* A 5 I 2
PV6 = A2 I A6 I A 3 I A 5 I
P1F = Z3nV5 A4I
RIS = Z4 *DV 1 A4 I
P2F = Z 3 CV6 A 2 I
R2S = Z4DV6 A? I
BIS = R1S/CSI
81F = RIF/C FI
B2S = Z 2 *RFS
B2F = Z1RFS
RAS = 2. Z6 R2S
6F = 2.*Z5*R2F
368 CONTINUE
IF(I.GT.1)GO TO 370
PS IP = 0.0
RATIO = 0 .0
R A T I 0 P = 0 0
GO TC 330
370 CONTINUO
CALL PS I(PS IP,JRATE,SBPCDEL TAPDY.SY.BLTAXKO,DHAT,DZ)
SBUIF = SBPC SBB
OELTAP = DELTAB + SBUIF + DELT2MPSIP + PSIB)
lFISGPC.LT.SMAXB.ANU.I.LT.MJ)GO TO 371
215
C fc L T A P = CELTAP + .5MPHIP + PHIO)*SBDIF
371 CUM I NUE
1F(SUPC.GT.1.0E08)GO TU 372
RATICP 0.0
GU TO 373
372 CONTINUE
FAT ICP = PSIP/SBPC
3 73 CONTI.NUt
I F( I.GT.2 )G0 TO 375
RATIC = RAT I OP
GO TO 300
375 CONTINUE
RATIC = .625 *RATIOP + PSIBC
380 CONTINUE
77 = Z5*RATI0
ZR = Z6 K AT IU
79 6.*73*77
7 lu = 6.* Z4 Z 8
Cl = A20 OELT *RATIOP
L2 = AMJ + A*DELT2*RATI0P + 04
1F(RFSA.LT. 1 .0E06.UR.UNCOUP.GT.O) GO TO 382
711
=
17* ( 2.* R1F
+ A R 2 F )
712
=
Z8*(2.*R1S
* A *R2 S )
CV7
=
UF/Z3
CV8
=
RIS/Z4
B3S
=
RF S *(1 
Z 10 )
P3F
=
RF S*( 1 ~
79)
04 S
=
DVH Z 12
B 4 F
=
CV7 Zll
35 S
=
Z8*(RIS +
I 2 A R 2 S )
05 F
=
Z 7 ( R 1 F +
I2A*R2F )
0 7 S
=
05S + 03* B6S
B7F
=
B6F + G 3* B6F
D3S2
= R F S* ( 1. +
Z10)
216
H3F2 = R F S* ( 1. + 7.9)
G4S2 = 0V8 + 712
C4F2 = [)V 7 + 711
RHSB = 1 F V X RR B + B2F*VTRRG + B 3F* TAURK 4 B4F SXRRB +
1 A*(S5F* STRRB B6F *VRRRB )
RHS D = B 1SVXRBB 4 B2S* VTRBB 4 B 3S T AURBB 4 B4S*SXRBB +
1 A *(B5S*STRBB B6S VRR BB )
GO TC 333
382 CONTINUE
ARATIO = A2I/A4I
7 13 = 1 7Z3
714 = 7 7 ( 2 4 AARA TIU )
F2S = 1. 4 7 10
F3S = 1. 710
F1F =1./CFI
F2F = 713 4 714
F6F = 713 714
F3F = 7 7 *( I2A*ARATIU 4 1.)
F4F = 2. *75*ARATIO
F 5 F = F 3 F 4 Q 3 F 4 F
PHSBE = F1F*VXRRB 4 F6F SXRRB 4 A*{F3F*STRRB 4 F4F *VRRRB)
RHSUE = 72*VTRBR 4 F 3S* T AURBB
KHSBEM = RHSBE 4 A*F4F*D3
PHSEt = (2.*VR8 PS I BB* l I 2A*STB SXB))*DELT 4 A2QSXB 4
RHSEeM = RHSEE 4 D4
383 CONTINUE
IF(MX.GT.1)CO TC 385
IF(RFSA.LT. 1.0E06.OR.UNCOUP.GT.Q) GO TO 384
IFUASE.NE. 1 ) GU TU 500
RHSBM = RHSB C3F2*TAUI D4F2SXI
R FIS CM = RHSD C3S2*TAUI C4S2*SXI
R FIS EM = RHSH A5Q*TAUI D1*SXI
STI A*RHSEM/D2
RHS1 = RHSBM 4 A*(B7F*STI B6F*D3)
A4Q*S TB
217
RHS2 = RHSDM + A*(B7S*STI R6S*D3)
CEL 4 = B 1 F B 2 S B1SB2F
VXI = (RHS1*82S RHS2*C2F)/DEL4
VTI = (RHS2*D1F RHS 1*BIS)/DEL4
GO TC 390
580 I F(A.EO.O )GU TO 585
IFIKASE.NE.2)GO TO 581
ZA = D4S2*D2 + C1*B7S
ZB = D3S 2 *02 + A5q* B7S
ZC = C35 2 *01 A50* D4S2
Cr L6 = D3F2*ZA D4F2*Z3 B7F*ZC
RHS4 = RHSB AD3*B6F B1F*VXI B2F*VTI
RHS5 = RHSD A D 3* B6S B1S*VXI B2S*VT1
TAUI = ( RHS4*Z ARFS5* ( D4F 2 02 + 1* 3 7F ) RHSH* {D4S2 3 7FD4F2 *B7S ) ) /DL:L6
SXI=lRHS4*ZB+RFS5*(03F2*D2+37F*A5)RHSH*(03S2*B7FD3F2*B7S))/DEL
16
ST1 = (RHS4*ZCRHS5*(O3F2 *D14F2 *A5U)+RHSH*(D3F2U4S2D4F2*D3S2) )/
1CEL6
GO TC 39 J
581 IFKASE.NE.3 ) GO TO 582
RHS6 = RHSB A*D3*B6F B?F*VTI D4F2*SXI
RHS7 = RHSD A*D3B6S B2S*VTI D4S2*SXI
RHS 8 = R H SH D1*SX1
ZA = C2*D3S2 + A 5Q B 7 S
ZB = D 2 D 3 F 2 + A 5 Q H 7 F
CtLB = 81F* Z A B 1 S Z B
VXI = {R H S6 Z A R H S 7 Z B + RHS8 *{D3S2*B7F 3F2*B7S) )/DEL8
TAUI = ( B1F*(RHS7*D2 + RHS8*B7S)B1S*(RHS6*02 + RGS8*B7F) J/DEL8
ST I = (B1F*(D3S2*RHS8A5Q*RHS7)BlS*(D3F2*RHS8A5Q*KHS6) )/EL8
GO TC 390
582 CONTINUE
RHS9 = RHSB A*D3*B6F B1F*VXI D3F2*TAUI
RHS10= RHSD A*D3*B6S B1S*VXI D3S2TAUI
RHS11= RHSH A 5Q*T AUI
to
oo
7 A = C 2 D 4 S 2 + C1*B7S
ZB = 02 C4F 2 + C1*B7F
DELI") = B 2F Z A B2S*ZB
VT I = ( RHS 9 *7 A RHS1CWB + RHS 11 ( D4S2* B 7FD4F 2 *B 7S ) )/DE L 1 0
SXI = (B2FM 02RHS107S*RHS11)B2S*(D2*RHS9+37F*kHSl1) )/UEL10
ST I = (B2F *(D4S2*RHS11D1*RHS10)R2S*(D4F2*RHS11D1*RHS9))/0EL10
GU TC 390
585 IFIKASE.NE.2)GO TO 536
RHS4 = RHSB A*D3*H6F D1F*VXI B2F*VTI
RHS 5 = RFSD A*D3*36S B1S*VXI B2S*VTI
CEL7 = U3F2*U4S2 D4F2*D3S2
T AUI = (RHS4 *04S2 RHS5*D4F2)/DEL7
SXI = (RHS5*D3F2 RHS4 D 3S 2 ) /DF L 7
ST l ^ 0.0
GU TU 390
586 IE(KASE.NE.i)GO TO 587
RHS6 = RHSB AD3*86F 32F*VTI 04F2*SXI
R HS 7 = RHSD A*D3*B6S B2S*VTI 04S2*SXI
GE 19 = B1F*D3S2 BLSD3F2
VXI = IRHS6*D3S2 RHS7 *D3F2 ) 70EE9
TAU = (RHS7*B1F RHS6*B13 )/OEL9
ST1 = 0.0
GU T 390
587 CONTINUE
RHS9 ^ RHSB AC3*R6F B1F*VXI D3F2*TAUI
RHS10 RHSO AD3*6S B1S*VXI D3S2*TAUI
CECIL = B2F*C4S2 R2S*04F2
VT I = ( RHS9 *04S 2 RHS 10D4F 2 )/PEL 11
SXI = (RHS19*B2F RHS9*B2S)/DEL 11
ST I = 0.0
GO TO 39 U
384 CONTINUE
IF(KASE.E0.2.OR.KAS EEQ4)GO TO 594
STi = A*(RHS E EM D1*SXI)/D2
219
VTL8B = C0N2 *VT LB + CCN4*VTB
VKLBB = C0N2 *VRLB + C0N4*VRB
IF(RFSA.LT.1.0E06.UR.UNC0UP.GT.0) GO TO 386
CIS = 2.*B1S
Cm
\J1
U1
\ji
cn
CD
vO
vC
vC
sO
'J1
CD
O
a
f'
IV)
<
to
CO
*
<
<
<
co
co
1
r>
X
H
X
>
X
H
X
H
X
>
a
a
r~
r~
d
d
r~
r~
r~
r
f~
c
*>
10
cc
CO
r~
r*
r"
r
r*
r*
r
i
i
CD
CD
CD
CD
CD
CD
CD
CD
CD
r*
*
o
C3
CD
d
II
II
It
II
ii
!l
ii
11
c
CJ
II
II
rn
o
O
n
o
o
o
o
o
O
o
o
n
o
o
o
o
CL
n
o
d
r
*
d
'
d
d
d
7
d
o
ro
ro
ro
V
r *
*
P*
*
*
ro
*
*
*
*
r*
<
CO
CO
*
<
<
<
co
co
X
X
i
X
1
X
H
X
H
r*
r*
r~
>
r~
r
r~
r*
r~
>
CD
X
CD
d
d
r*
4
4
4
4
4
r~
4
4
4
CD
o
n
O
o
ro
4
o
o
n
4
c
t
CL
o
d
o
r
7
d
T
.
d
7
o
V
d
o
CO
CM
cm
co
CO
o
T
CD
*
*
7]
*
*
*
d
c
<
<
co
CO
CO
<
Co
co
73
H
X
H
X
X
i
X
r*
r
r~
r~
r*
H
CD
CD
CD
i
CD
CJ
r
3:
CD
>
>
CL
c
r
CD
CD
c
u)
l
*
co
LO
cn
co
CO
X
Cl
1
o
!
t
>
T1
5
X
m
o
X
X
X
CL
>
d
d
r*
Co
c
1
7C
H
i
H
1
<
1
II
O
>
II
ii
ro
c
ti
II
ro
CL
o
II
00
m
II
00
m
ii
'J'
O
X
Cj
ii
CO
CO
X
xJ
T1
o
d
sC
JL
T
o
X
m
f\J
ro
d
g
o
co
t
X
o
X
o
CO
X
c
X
v
ro
~r
T
O
co
Â¥
7D
ro
o
CO
co
fT;
d
73
X
X
CL
X
d
m

X
CO
ro
X
m
m
1
o
VI
XI
ro
H
3
1
o
m
ro
XI
a
1
xi
rr.
1
i
ro
rx
1
X
1
CD
rx
co
ro
CD
a
O
T1
ro
*
*
1
n
*
O'
*
<
on
sjl
*
n
<
>
i
v0
o
X
XI
i
CD
*
u*.
<
*
73
XI
X
*
73
X
**
n
X
co
X
rx
ro
CO
m
ro
ro
co
*
m
co
r^
3:
X
a
o
x
rn
r*
r"
ro ro
ozz
VXI = (RHSBEM F 2 F SX I AF6FSTI )/FlF
IF(KASE.LQ.3)G0 TO 692
VTI = (RHSDE F 2 S T AUI )/Z2
02 S
=
2.* B2 S
06 S
=
2. B6S
DIF
=
2.*B1F
C2F

2 B2 F
C6F
=
2. B6F
C3S
2.*D3S2
C 3 F
=
2.D3F2
04 S
=
2. 04 S 2
C4F
=
2 C4 F 2
DBS
=
2.BBS
DBF
RHSA
=
2 BBF
= BIF*VXLLB B2F*VTLLB + B3F TAULLB
1 A*(B5F*STLLB B6F*VRLL3)
RHSC = B1S*VXLBB B2S VTL BB B3S*TAULBB
L A*(B5S*STLBB B6S*VRLB6)
OFL1 = O 1F*D2S D15*D2F
RHSBA = RHSB RHSA
RHi.CC = RHS 0 RHSC
VXI = (D 2 S R H S B A 02F*RHSDC)/DEL 1
VTI = (01 F*RFSDC D 1S*RHSBA)/DEL 1
C7F = DBF Q3*C6F
C7S = DBS Q3C6S
RHSF = RHSA + RFSB A D 3 D 6 F
RHSC = RHSC + RFSD A*D3*D6S
GO TO 3B 8
64F*SXLLB +
+ 34 S* SXLBB +
386 CONTINUE
RHS C E = Z2*VTLSB + F3S*TAULBB
RHSAE = F1F*VXLLB + F6F*SXLLB + A*(F3F* STL LB + F4F*VRLLB)
FHSAEM = RHSAE + A*C3*F4F
TAUI = (RHSCE + RHSDE )/(2.*F2S)
VTI = (RHSDE RHSCF )/(2.*Z2)
VXI = (KHSBEM RHS A EM ) /(2.*FIF)
F2F2 = 2.*F2F
FBF2 = 2.F5F
221
RHSi = RHSAEM + RHSBEM
380 CONTINUE
IF(A.E0.0 ) GO TO 337
IFIRFSA.LT.1.QEG6.OR.UNCOUP.GT.O) GO T 389
ZA = D2*4S DI D7 S
ZB = ASQ C7S C2 l) 3S
ZC = 0 1 C 3S A5Q*D4S
CEL2 = U3F#ZA + D4F*ZB + l)7F*ZC
TAUI = (RHSF*ZA*RHSG*
SXI = (RHSF*ZB+RHSG*(02*03FA5Q*D7F)+RHSH*(D3S*07FD3F* 07S) )/0LL2
ST I = (RHSF*ZC + RHSG*(A5Q*D4FD1 *D3F)+RHSH*<3F*D4S3S*4F) )/DE L2
GO TO 390
389 CONTINUE
CELS = D 2 F 2 F 2 D1*F5F2
SXI ( R H S3 02 RHSEEMF5F2)/DCL5
STI = (RHSEEM*F2F2 RHS3*01)/DFL5
GO TO 39 J
387 CONTINUE
IF(RFSA.LT.1 .OE96.OR.UNCOUP.GT.O) GO TU 381
DEL 3 = 0 3 F* D4S D3S*D4F
STI = 0.0
TAUI = (RHSF*C4S RHSG*D4F )/DEL 3
SXI = (RHSGD3F RHSF*D3S)/DEL 3
GO TO 390
361 CONTINUE
SXI = RHS3/F2F2
STI J.O
394 CONTINUE
SBP(l) = SGRT(SXI**2 + 3.*TAUI**2 + A*STI*(STI SXI))
STRAIN!I ) = DELTAP
VXP{ I ) = VXI
V T P ( I ) = VT I
SBPC = SBP(I)
IF(I.EG.1)GO TO 340
222
riFFl = ADS(SOP{ I ) SB P( I 1 ) )
DIFF2 = ABS(STRAIN! I ) STRAIN! II) )
CIFF3 = ADS ( VXP ( I ) VXP(ID)
CIFF4 = A BS(VT P( I ) VTP(Il))
CIFF = U1*DIFF1 + H2 C 1 FF2 + H3*DIFF3 * H4*DIFF4
DENOM = A BS ( VX P ( I ) ) + ABS(VTP(I))+ ABSISBPUJ) +
FRRCR = CIFF/DENOM
IF(ERROR.GT.SMALL)GO TO 340
430 CONTINUE
IFIMX.GT.1)00 TO 391
Â£X(MX) = Â£X(MX) + 2*VXR VXI VXB
ET X(MX) = ETX(MX) + VTR .5*(V TI + VT3)
GC TO 392
391 CONTINUE
FX(MX) = EX(MX) + VXR VXL
ETX(MX) = ETX(MX) + .5*(VTR V TL)
392 CONTINUE
VR{MX ) = A*(D3 03*STI )
ET(MX) = ET(MX) + DELT2*(VR(MX) + VRB)
TAU(MX) = TAUI
VX(MX ) = VXI
VT f M X ) = VT I
SX(MX) = SXI
ST(MX) = ST I
Al(MX) = A1P
A 2(M X) = A2P
A 3 ( M X ) = A 3 P
A4 ( MX ) = A4P
A^tMX) = A5P
A6(MX ) = A6P
CF(MX) = CFI
CS(MX) = CSI
S 3(MX) = SBPC
EXP(MX ) = EX(MX) (SXI NU*STI)
AB S(STRAIN(I))
to
to
FTP(MX) = ET(MX) {STI A*NU*SXI)
LTXP(MX) = E T X ( M X ) QA* TAUI
ITER(MX) = I
PS 10(MX) = PS IP
PHI. (MX) = PH IP
UELTA(MX) = CELTAP
PL AS(MX ) = CELTAP SBPC
CELGAM(MX) = 2.*(E T X(MX) ETXO)
CtX(MX) = EX(MX ) cXC
STKES(MX) = SQRT(SXI**2 + A*STI*(STI
IF(SBPC.LT.SMAXO)GO TO 393
SMAX(MX) = SBPC
393 CONTI CUE
IF(MX.GT.1)GO TO 399
CELLFB(MTU) = CELTAP
399 CONTINUE
IF(KCUNT.EQ. IMAX ) GO TU 400
MTT = MX + MT
IF(MTOT.LT.NMAX )G0 TO 2C0
LINE = LINF + 1
IWRITE = 2
GU TO 4)1
397 CONTINUE
KUUNT = KOUNT + 2
IF(KCUNT.LT.NMAX)GO TO 100
400 CONTINUE
LINE = LIME + 1
IWR I TE = 3
491 CONTINUE
LINEO = LINEO + 1
IFLINEO.LT.INCRTJGO TO 402
LINEO = O
W R I T t ( 6 1QG0) LINE
WRITE(6,1010) (X(L).LltMX, INCRX)
SXI ) )
224
402
WRITEl6, 1009) (T(L)*L=1,MX, INCRX)
WklTE(6,lCll)
WRITE(6, 1012)
WRITE(6, 1038 )
WRITE(6,1013)
W R I r E ( 6 10 1 4 )
W R I T F ( 6 ,1015)
WRIT6(6, 1016)
WRIT E{6 i 1Cl7)
W K1 F c(6 1^18 )
WRITEI6, 1019)
WRITE(6, 1020 )
WRITE6,1021 )
WRITE(6, 102 2 )
WRITE(6, in?3 )
WR I TE (6, K'24 )
WRITc(6,102 9)
WRITE16, 1^26)
W R 1 T E(o 1027)
WRITE(6, 102 8 )
WRITE(6,1029)
WRIT F(6, 1036)
WRITlId,1037)
(SX(L ) ,1=1,MX,INCRX)
(ST(L ) ,L= 1,MX,INCRX)
(STRES(L),L=1,MX, I NCR X)
(TAU(L ) ,L = ltMX, INCRX )
(VX(L ) L= 1 .MX,INCRX)
(VT(L ) L = 1 MX.INCRX)
(V R(L ) L = 1 M Xi INCRX)
{CF(L ) L = 1,MX,INCRX)
(CSC L ) ,L l ,MX, INCRX)
C SB(L ) i L = 1 *MX, INCRX)
(EX(L),L=1,MX, INCRX )
( E T(L ) L=i,MX, INCRX)
( ETXtL )L=1MX, INCRX)
(EX P(L ) ,L = 1,MX,INCRX)
(E T P{L ),L = 1,MX, INCRX)
(ETXP(L ) t1 = 1,MX,INCRX)
( ITER(L),L = 1,MX,INCRX)
(CFL T A (l_ ) L = 1 MX, INCRX )
(CPLAS(L ) ,L = i,MX, INCRX)
(PH 10(L),L = l,MX, INCRX)
(LELCAMCL),L=1,MX,INCRX)
l f: F X ( L ),L = 1,MX, INCRX)
CONTINUE
IF( 1 PUNCH.E.0)C0 T0 404
LIN E P = LINIP + I
IFCLINEP.LT. INCRTPJGO TO 4U4
LINEP =
WRITt t 7, 1^3 5) (X(L ) ,L = 1,MX, INCRXP )
WRITEC7,1C 3 5 ) (T(L),L=1,MX,INCRXP)
WRITE(7, 1036) (SB(L ) L= 1,MX,I NCRXP)
WRITEC 7, 103 5) (CELTA(L),L = 1,MX, INCRXP)
IF( I PUNI.EQ.C)GO TO 40 3
WRITEC7, 1C35) {SXIL ) ,L= 1,MX,INCRXP)
225
WRITE 7* 1C35) (EX(L ),L = 1,MX,INCRXP)
WRITE<7, 103 5 ) (VX(L),L=1, MX,INCRXP)
W R I T E ( 7 1335) (CEX(L ), L = 1,MXINCRXP )
IF ( A.EC1.C )G0 T 403
WRI TCl7, 103 5) (ST(L ) L= 1 ,MX,INCRXP)
WPITb(7, 1035) (ET(L ) ,L=lfMX,INCRXP)
WRI TE(7, 1035 ) (VR(L ) ,L=1 ,MX,INCRXP)
403 CUNTI l\IUL
I F ( IPUN2.EQ.OGO T0 404
WRITE!7,1035) (TAU(L ),L = 1MXvINCRXP)
WRI TE(7, 10 35 ) (VT(L ) ,L=l,MX INCRXP)
WRITE(7, 1035 ) (ETX(L),L=1,MX,INCRXP)
WRITEI 7, 103 5 ) (CELGAML),L = 1,MX,INCRXP)
404 CONTINUE
IF! IWRITE 2) 90, 397,405
4C5 CONTINUE
LODO FORMAT(' 1 ,2 X,'LINE =',I5///)
100 1 FORMAT!3 15,5F5.2,4F10.2)
1002 FORMAT!I5.7F10.3,15)
ICC 3 FORMAT(1015)
1004 FORMAT(2 E207)
1005 FORMAT( 1 ,9X,A = , I 1, 15X, 'JR ATE = ,I 1,11X,ELX = *,F9.6,4X,
I'OlLT = ,F9.6,4X, MXMA X = ,I 4,BX, SMALL = ',F8.6/)
1V.6 FORMAT ( 19X, NU = F 5.3 1CX SY = F12.9,3X 1 BOAR F 9. 1,4 X N
1= ,F7.4,9X, 'XK:: = F 8.6,6 X BETA = ',F8.6/)
1007 FORMAT!1 O X,J RIS F
i,*SXO = ,F9.6,5X
iC C 8 FORMAT(1TX, 'HI =
1= F4.2, 11X, n
10 0 9 FORMAT(/3 X, 'TIME
1010 FORMAT(/ IX,1X =
1011 FORMAT!/3X,SX =
1012 FORMAT(/3X,ST =
1013 FUKMATI/3X,'TAU =
= I 4,BX, 'XVFIN = ',F 9.6,3 X, 'TVFIN = *,F8.6,3X
TAUO = ,F8.6,5X, H = ',F8.6/)
,F4.2,11X,H2 = ,F4.2, 11X, *H3 = ,F4.2,11X,H4
12, 13X, SYS = ,F 1 2.9/)
',4X,71F9.56X)F9.5/(8X8(6X,F9.5)))
',4X,7!F9.5,6X),F9.5/!8X,0(6X,F9.5)))
,8 E15.5/( 10X.8E15.5) )
' t 8E155/(10X,6E15.5))
.8E15.5/10X.8E15.5))
226
1314
FORMAT(/3X,
VX =
6cl5.5/( 10X,8E15.5) )
1015
FORMAT(/3X,
' VT =
',8E15.5/1 10X.8E15.5) )
1016
FORMAT(/ 3X,
VR =
*,8E15.5/1 10X,8E15.5> )
1017
FORMAT(/3X,
CF =
',8E15.5/( IDX,RE15.5) )
l 0 1 ri
FORMAT(/3X,
' cs =
'8E 15.5/( 1DX,815.5 ) )
1019
FORMAT(/3X,
SB =
,8E15.5 / ( inx be 155 ) )
1C 2 >j
FORMAT(/ 3X T
EX =
,8E15.5/t 11X,BE 15.6) )
1021
FORMAT(/3 X,
' ET =
,8E15.5/( 10X,E15.5 ) )
102 2
FORMAT(/3X,
FTX =
',8E15.5/ 10X,8E15.5 ) )
1023
FORMAT{/3X,
! EXP =
'v0E15.5/( 10X.8E15.5) )
102 4
FORMAT(/3X,
' ETP a
',8E i 5.5/( 10X.8E15.5 ) )
102 5
FORMAT(/3X,
ETXP =
'8E15.5/<10X.8E15.5M
10 2 6
FOkMAT(/3X,
' ITER =
,01 15/( 10 X 8115) )
102 7
FORMAT(/3 X,
'DELTA
= '(8E15.5/( 10X t 8E15.5 ) )
1 0 2 H
FORMAT(/JX,
CPLAS
= ,8tl5.5/( 10X,8E15.5) )
1020
FORMAT(/3X,
PHI =
*8E15.5/( 10X.8E15.5) )
13 30
FORMAT(10X,
' JBEIL
= 'tlltllXf'D/ = ',rfi.6,7X, 'DHAT = F8.6,5X XM
i
= ',FK.8,
5 X,'XN
= *F10.85Xt XLAM = ,F8.6/)
1031
FORMAT(1GX,
MI = '
, I 3,12X MJ = ,I 3,12X, I PUNCH = ,I 38X 'KAS2 =
J.
I2/ )
10 3 2
FORMAT(10 X,
X = ,
F 9.6 5 X, T = F9.6,5X,'DIFF = ,L13.5,5X,'DNOM
1
= 'ttl3.5.
5 X,ERRUR = ,E13.5//)
103 3
FORMAT( 1 )
10 34
FORMAT(615,
2F10.6)
10 3 5
FORMAT(7 E 11
. 6 )
10 36
FOKMAT(/3X,
'OGAM =
',6E15.5/( 10X,8E15.5) )
1037
FORMAT(/3X,
DEX =
,8E15.5/(10X.8E15.5) )
1038
FORMAT{/3X,
'STRES
= ',PE15.5/( 10X.8E15.6) )
IO 3 > FORMAT(//////50X '
1040 rORMAT(//////50X, '
STOP
ENL
WAVES ARE UNCOUPLED ******#**#)
WAVES ARE COUPLED ********#)
227
SUlii (JUT I NE PHI ( P PQ K, J 81 S U YS B t Y H XLAM, OE TA Y XM XN DZ I 0 )
IF(K.GT.O) GO TO 4
YS = ABS(YS)
IF(S.GT.YS)GO TO 1
SZ = BETA*SuRT( CZ)
IF( IO.EQ.O) C = S
P = .. j
PQ O.J
no tc o
1 CON T I NU 
IF(JB.GT.OIGO TC 2
/ = S YS
ZY = Z *(Y1 )
C = 8 *ZY
P = C Y
IF( Ij.GT.OIGO TO 10
CP c*z
c = dp s
PQ = OP/S
G TC n
2 CONTINUE
I F ( m.GT .0 ) GO T 3
SZ = BETA*SQRT(CZ )
C = (S/GLTA)**2
IF(C.LT.CZ) C = (2.0*SQRT(DZ)/B E T A ) *(S YS) DY
PQ  /S 1.0
3 CONTINUE
IF(S.LT.YS) P = 0.0
IF(S.Gfc.YS.AND.S.LE.SZ) P = 2.O SORT(DZ)/BETA 1.0
IF(S.GT.SZ) P = 2 O S/(B E T A * 2 ) 1.0
GO TC lu
4 CONTINUE
C THE DATA USED HERE ARE FROM CRISTESCU (1972) CASE XII
CZ1 = SORT(DZ)
228
SYH = YS + H
DC = (SYH/BETA)*2
I F(CC.GT.DZ) GO TO 5
DC = 2.*H*DZ1/BETA + DY
5 CONTINUE
I F ( DGTCY)GO TG 6
F = YS
GO TO 8
b CONTINUE
IF(D.GT.DZ)GO TO 7
F = YS + EET A*(C 0Y)/(2.*DZI)
GO TO 8
7 CONTINUE
F = 01T A S O R T ( )
3 CONTINUE
CSTAR = H + XLAM*(D DC)
FM = F * OS TAR
IFIS.GT.FMJGO TC 9
P = 0.0
00 TC 10
0 CONTINUE
XA = XM + XN*SvjRT(D)
P = (3.*(C DY DSTAR + (XA/3.)**1.5)**(2./3.))/XA 1.0
10 CONTINUE
RETURN
END
229
SUCAGUTINE PS I(PK,S*D,DY,YS.BETA,XKtHAT,DZ)
IFIK.GT.O )GO TO 20
P = 0.0
GO TC 25
2 j CONTINUE
C THE CATA USED HERE ARE FROM CRISTESCU (1972) CASE XII
071 = SwRT(UZ)
IF(D.GT.CY)GO TO 21
YS = ABS(YS)
F = YS
GO TC 23
21 CONTINUE
I F(C.GT.DZ)GO TO 2 2
F = YS + .5*BETAMD DYJ/0Z1
CO TO 23
22 CONTINUE
F = BETA*SORT(D )
23 CONTINUE
IF(S.GT.F.ANC.C.GE.CY)G0 TO 24
F = 0.0
GO TC 25
24 CONTINUE
XK = XK:*(1.0 EXP(D/DHAT ) )
P = XK ( S F )
25 CONTINUE
RETURN
END
230
LIST OF REFERENCES
Alter, B. E. K. and Curtis, C. W. (1956), "Effect of Strain Rate on the
Propagation of a Plastic Strain Pulse along a Lead Bar," Journal
of Applied Physics, Vol. 27, pp. 10791085.
Bell, J. F. (1951), "Propagation of Plastic Waves in Prestressed Bars,"
Technical Report No. 5, Johns Hopkins University, Navy Contract
N6ONR243VIII.
Bell, J. F. (1960), "Propagation of Large Amplitude Waves in Annealed
Aluminum," Journal of Applied Physics, Vol. 31, pp. 277282.
Bell, J. F. (1963), "Single, TemperatureDependent StressStrain Law
for the Dynamic Plastic Deformation of Annealed FaceCentered
Cubic Metals," Journal of Applied Physics, Vol. 34, pp. 134141.
Bell, J. F. and Stein, A. (1962), "The Incremental Loading Wave in the
Prestressed Plastic Field," Journal de Mecanique, Vol. 1,
pp. 395412.
Bianchi, G. (1964), "Some Experimental and Theoretical Studies on the
Propagation of Longitudinal Plastic Waves in a StrainRate
Dependent Material," Stress Waves in Anelastic Solids, I.U.T.A.M.
Symposium, Brown University (H. Kolsky and W. Prager, Eds.),
Berlin: SpringerVerlag, pp. 101117.
Clark, D. S. and Wood, D. S. (1950), "The Tensile Impact Properties of
Some Metals and Alloys," Transactions of the American Society of
Metals, Vol. 42, pp. 4574.
Clifton, R. J. (1966), "An Analysis of Combined Longitudinal and
Torsional Plastic Waves in a ThinWalled Tube," Proceedings of
the Fifth U.S. National Congress of Applied Mechanics, ASME,
N.Y., pp. 465480.
Clifton, R. J. (1968), "ElasticPlastic Boundaries in Combined
Longitudinal and Torsional Plastic Wave Propagation, Journal of
Applied Mechanics, Vol. 35, pp. 782786.
Convery, E. and Pugh, H. L. D. (1968), "Velocity of Torsional Waves in
Metals Stressed Statically into the Plastic Range, Journal of
Mechanical Engineering Science, Vol. 10, pp. 153164.
Cristescu, N. (1959) "On the Propagation of ElasticPlastic. Waves in
the Case of Combined Stress," Prikladnaia Matematika i Mekhanika,
Vol. 23, pp. 11241128.
Cristescu, N. (1964), "Some Problems in the Mechanics of Extensible
Strings," Stress Waves in Anelastic Solids, I.U.T.A.M. Ss^mposium,
Brown University (H. Kolsky and W. Prager, Eds.), Berlin:
SpringerVerlag, pp. 118132.
231
232
Cristescu, N. (1965), "Loading/Unloading Criteria for Rate Sensitive
Materials," Archiwum Mechaniki Stosowanej, Vol. 17, pp. 291305.
Cristescu, N. (1967a), Dynamic Plasticity, New York: John Wiley and
Sons, Inc.
Cristescu, N. (1967b), "Dynamic Plasticity under Combined Stress,"
Mechanical Behavior of Materials under Dynamic Loads (U.S.
Lindholm, Ed.), New York: SpringerVerlag, pp. 329342.
Cristescu, N. (1968), "Dynamic Plasticity," Applied Mechanics Reviews,
Vol. 21, pp. 659668.
Cristescu, N. (1971), "On the Coupling of Plastic W'aves as Related to
the Yield Condition," Romanian Journal of Technical Science,
Applied Mathematics, Vol. 16, pp. 797809.
Cristescu, N. (1972), "A Procedure for Determining the Constitutive
Equations for Materials Exhibiting Both TimeDependent and Time
Independent Plasticity," International Journal of Solids and
Structures, Vol. 8, pp. 511531.
Davis, E. A. (1938), "The Effect of Speed of Stretching and the Rate of
Loading on the Yielding of Mild Steel," Transactions of ASME,
Vol. 60, pp. A137A140.
DeVault, G. P. (1965), "The Effect of Lateral Inertia on the Propa
gation of Plastic Strain in a Cylindrical Rod," Journal of the
Mechanics and Physics of Solids, Vol. 13, pp. 5568.
Donnell, L. H. (1930), "Longitudinal Wave Transmission and Impact,"
Transactions of ASME, Vol. 52, pp. 153167.
Duwez, P. E. and Clark, D. S. (1947), "An Experimental Study of the
Propagation of Plastic Deformation under Conditions of Longitu
dinal Impact," Proceedings of ASTM, Vol. 47, pp. 502532.
Efron, L. and Malvern, L. E. (1969), "Electromagnetic Velocity
Transducer Studies of Plastic Waves in Aluminum Bars,"
Experimental Mechanics, Vol. 9, pp. 255262.
Hencky, H. (1924), "Zur Theorie Plastischer Deformationen und der
Hierdurch im Material Hervorgerufenen Nebenspannungen,"
Proceedings of the First International Congress for Applied
Mechanics, Delft, pp. 312317.
Hill, R. (1950), The Mathematical Theory of Plasticity, England:
Oxford.
Hopkins, H. G. (1961), "Dynamic Anelastic Deformation of Metals,"
Applied Mechanics Reviews, Vol. 14, pp. 417431.
233
Hunter, S. C, and Johnson, I. A. (1964), "The Propagation of Small
Amplitude ElasticPlastic Waves in Prestressed Cylindrical Bars,"
Stress Waves in Anelastic Solids, I.U.T.A.M. Symposium, Brown
University (II. Kolsky and W. Prager, Eds.), Berlin: Springer
Verlag, pp. 149165.
Karman, T. von (1942), "On the Propagation of Plastic Deformation in
Solids, National Defense Research Commission Report Number A29
(OSRD No. 365).
Kolsky, H. (1963), Stress Waves in Solids, New York: Dover Publications,
Inc.
Lee, E. H. (1953), "A Boundary Value Problem in the Theory of Plastic
Wave Propagation," Quarterly of Applied Mathematics, Vol. 10,
pp. 335346.
Lindholm, U. S. (1967), "Some Experiments in Dynamic Plasticity under
Combined Stress," Mechanical Behavior of Materials under Dynamic
Loads (U.S. Lindholm, Ed.), New York: SpringerVerlag, pp. 7795.
Lipkin, J. and Clifton, R. J. (1970), "Plastic Waves of Combined Stress
Due to Longitudinal Impact of a Pretorqued Tube," Journal of
Applied Mechanics, Vol. 37, pp. 11071120.
Lubliner, J. (1964), "A Generalized Theory of StrainRateDependent
Plastic Wave Propagation in Bars, Journal of the Mechanics and
Physics of Solids, Vol. 12, pp. 5965.
Lubliner, J. and Valathur, M. (1969), "Some Wave Propagation Problems
in PlasticViscoplastic Materials, International Journal of
Solids and Structures, Vol. 5, pp. 12751298.
Malvern, L. E. (1949), "The Propagation of Longitudinal Waves of Plastic
Deformation in a Material Exhibiting a StrainRate Effect,
Brown University Dissertation.
Malvern, L. E. (1951a), "Plastic Wave Propagation in a Bar of Material
Exhibiting a Strain Rate Effect," Quarterly of Applied Mathematics,
Mathematics, Vol. 8, pp. 405411.
Malvern, L. E. (1951b), "The Propagation of Longitudinal Waves of
Plastic Deformation in a Bar of Material Exhibiting a StrainRate
Effect," Journal of Applied Mechanics, Vol. 18, pp. 203208.
Malvern, L. E. (1965), "Experimental Studies of StrainRate Effects and
Plastic Wave Propagation in Annealed Aluminum," Behavior of
Materials under Dynamic Loading, ASME, pp. 8192.
Malvern, L. E. (1969), Introduction to the Mechanics of a Continuous
Medium, Englewood Cliffs, N.J.: PrenticeHall, Inc.
234
Manjoine, M. J. (1944), "The Influence of Rate of Strain and Temper
ature on Yield Stresses in Mild Steel," Journal of Applied
Mechanics, Yol. 11, A211A218.
Mok, C. H. (1972), "Effects of Lateral Inertia on the Propagation of
ElasticPlastic Waves," BRLCR70, Contract No. DAAD0570C0224,
Aberdeen Proving Ground, Md.
Nicholas, T. and Garey, G. F. (1969), "Torsion Testing of Aluminum at
High Rates of Strain," AFMLTR69172.
Olszak, W., Mroz, A. and Perzyna, P. (1963), Recent Trends in the
Development of the Theory of Plasticity, New York: Pergamon .Press.
Perzyna, P. (1963), "The Constitutive Equations for Rate Sensitive
Plastic Materials," Quarterly of Applied Mathematics, Vol. 20,
pp. 321332.
Plass, H. J., Jr., and Ripperger, E. A. (1960), "Current Research on
Plastic Wave Propagation at the University of Texas," Plasticity
(E. H. Lee and P. S. Symonds, Eds.), New York: Pergamon Press, Inc.,
pp. 453487.
Prandtl, L. (1924), "Spannungsverteilung in Plastischen Koerpern,"
Proceedings of the First International Congress for Applied
Mechanics, Delft, pp. 4354.
Rakhmatulin, Kh. A. (1945), "Propagation of an Unloading Wave,"
Prikladnaia Matematika i Mekhanika, Vol. 9, pp. 91100.
Rakhmatulin, Kh. A. (1958), "On the Propagation of ElasticPlastic
Waves Owing to Combined Loading," Prikladnaia Matematika i
Mekhanika, Vol. 22, pp. 10791088.
Reuss, A. (1930), Beruecksichtigung der Elastischen Formaenderunger
in der Plastizitaetstheorie," Zeitschrift fur Angwandte Mathematik
und Mechanik, Vol. 10, pp. 266274.
Shammamy, M. R. and Sidebottom, O. M. (1967), "incremental Versus Total
Strain Theories for Proportionate and Nonproportionate Loading of
TorsionTension Members," Experimental Mechanics, Vol. 7,
pp. 497505.
Shea, J. H. (1968), "Propagation of Plastic Strain Pulses in Cylindrical
Lead Bars," Journal of Applied Physics, Vol. 39, pp. 40044011.
Sokolovsky, V. V. (1948a), "Propagation of ElasticViscoPlastic Waves
in Bars," Doklady Akad. Nauk SSSR, Vol. 60, pp. 775778.
Sokolovsky, V. V. (1948b), "Propagation of ElasticViscoPlastic Waves
in Bars," Prikladnaia Matematika i Mekhanika, Vol. 12, pp. 261280.
235
Sternglass, E. J. and Stuart, E. J. (1953), '"An Experimental Study of
the Propagation of Transient Longitudinal Deformations in Elasto
plastic Media," Journal of Applied Mechanics, Vol. 20, pp. 427434.
Suliciu, I., Malvern, L. E. and Cristescu, N. (1972), "Remarks Concern
ing the 'Plateau' in Dynamic Plasticity," Archiwum Mechaniki
Stosowanej, Vol. 24, pp. 9991011.
Tapley, B. D. and Plass, H. J., Jr. (1961), "The Propagation of Plastic
Waves in a SemiInfinite Cylinder of a StrainRateDependent
Material," Developments in Mechanics, Vol. 1 (J. E. Lay and
L. E. Malvern, Eds.), Amsterdam: NorthHolland, pp. 256267.
Taylor, G. I, (1940), "Propagation of Earth Waves from an Explosion,"
British Official Report, R. C. 70.
Wood, E. R. and Phillips, A. (1967), "On the Theory of Plastic Wave
Propagation in a Bar," Journal of the Mechanics and Physics of
Solids, Vol. 15, pp. 241254.
Yew, C. H. and Richardson, H. A., Jr. (1969), "The StrainRate Effect
in the Incremental Plastic Wave in Copper," Experimental Mechanics,
Vol. 9, pp. 366373.
BIOGRAPHICAL SKETCH
Charles Daniel Myers was born July 27, 1945, in Brevard County,
Florida. He was graduated from Terry Parker High School in Jacksonville,
Florida, in June, 1963, and received his Bachelor of Science degree in
Engineering Science from the University of Florida in June, 1968.
He was employed in the Structural Analysis Section of the Martin Marietta
Corporation in Orlando, Florida. While on a leave of absence from the
Martin Marietta Corporation, he received his Master of Engineering
Degree in Engineering Mechanics from the University of Florida in
December, 1970. He returned to the Martin Marietta Corporation in
September, 1972. He is a member of Tau Beta Pi (having served as presi
dent of the Florida Alpha Chapter during the year 19671968), Sigma
Tau, Omicron Delta Kappa, Phi Kappa Phi, and the Society of Engineer
ing Science. He was recently selected for Outstanding Young Men in
America 1973.
Charles Daniel Myers is married to the former Peggy Lee Hufham
of Jacksonville, Florida. They have one son, Daniel Lee.
236
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
MS A. Eisenberg, Chairman
^Associate Professor of Engineering Science,
Mechanics and Aerospace Engineering
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
L. E. Malvern
Professor of Engineering Science, Mechanics
and Aerospace Engineering
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
E. K. Walsh
Associate Professor of Engineering Science,
Mechanics and Aerospace Engineering
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
T4 Pfnv^w^cr
U. H. Kurzweg
Associate Professor of Engineering Science,
Mechanics and Aerospace Engineering
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
R. C. Fluck
Associate Professor of Agricultural
Engineering
This dissertation was submitted to the Dean of the College of Engineering
and to the Graduate Council, and was accepted as partial fulfillment
of the requirements for the degree of Doctor of Philosophy.
August, 1973
Dean', College of Engineering
Dean, Graduate School
Internet Distribution Consent Agreement
In reference to the following dissertation:
AUTHOR: Myers, Charles
TITLE: Inelastic wave propagation under combined stress states, (record
number: 577578
PUBLICATION DATE: 1973
I, Charles D. Myers as copyright holder for the
aforementioned dissertation, hereby grant specific and limited archive and distribution rights to
the Board of Trustees of the University of Florida and its agents. I authorize the University of
Florida to digitize and distribute the dissertation described above for nonprofit, educational
purposes via the Internet or successive technologies.
This is a nonexclusive grant of permissions for specific offline and online uses for an
indefinite term. Offline uses shall be limited to those specifically allowed by "Fair Use" as
prescribed by the terms of United States copyright legislation (cf, Title 17, U.S. Code) as well as
to the maintenance and preservation of a digital archive copy. Digitization allows the University
of Florida to generate image and textbased versions as appropriate and to provide and enhance
access using search software.
This grant of permissions prohibits use of the digitized versions for commercial use or profit.
Signature of Copyright Holder
Charles D. Myers
Printed or Typed Name of Copyright Holder/Licensee
Personal Information Blurred
Printed or Typed Phone Number and Email Address of Copyright Holder/Licensee
5~ Z7 &&
Date of Signature
Please print, sign and return to:
Cathleen Martyniak
UF Dissertation Project
Preservation Department
190
so that equations (D.1.19), (D.1.21), and (D.1.6) can be written as
B2fV6P +
D4f2SxP
aB7fS6P =
RHS9
B2sV6P +
D4s2SxP
aB7sS0P =
RHS10
(D. 5. 30)
aDlSxP +
aD2S9P =
aRHSll .
When radial inertia effects are included (a=l) these can be
written as
B2f
D
4f 2
"B7f "
vep
RHS9
B2s
4s2
B7s
SxP
=
RHS10
0
D1
2 _
_S0P_
RHS11
(D.5.31)
and when
S10 =B2f B2s(D2D4f2 + DlB7f>
(D.5.32)
the solution is given by
1
V
9P A. LBHS9(D2D4S2 + D1B7s> 11,5101D2D4f2 + DlB7f>
10
+ RHSll(B7fD4s2 B7S
4f2)J
(D. 5. 33)
'xP ~~ A _B2f(D2RHS1+ b7sRHS11) B2s(d2RHS9 + B?fRHSll) I (D.3.34)
10
= ~ [B2f (D4s0RHS11 D1RHS10) B2s(D4f2RHSll D1RHS9)J (D. 5. 35)
8P
10
When radial inertia effects are not included (a=0), the last of
equations (D.5.30) vanishes as do the hoop stress and the radial
velocity, and the first two equations of (D.5.30) become
6
and these averaged variables were used. The results of this work were
given by Tapley and Plass (1961) but were somewhat 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
184
Adding equations (D.2.12) and (D.2.13) and letting
F = 2F
2f 2 2f
F
5f 2
2F
5f
RHS3 = RHSAEM + RHSBEM J
results in the equation
F2f2SxP + aF5f2SGP RHS3
(D. 4.4)
(D. 4.5)
When radial inertia effects are included (a = l), equations
(D.2.3) and (D.4.5) can be written as
F2f2
F5f2
"SxP~
RHS3
_ D1
D2
_ Sqp _
RHSEEM
(D.4.6)
and by defining
D F
2 2f2
 DlF5f2
(D.4.7)
the expressions for the hoop stress and the longitudinal stress can
be written as
S = j(D RHS3 F RHSEEM) (D.4.8)
xP A 2 5f2
5
S0p = r'(F2f2RHSEEM D!RHS3) (D.4.9)
5
When radial inertia effects are not included (a = 0), equation
(D.2.3) vanishes, the hoop stress and the radial velocity are zero,
and the expression for S is found directly from equation (D.4.5) to be
xP
S
xP
RHS3
F
2f2
(D.4.10)
19
Since the stresses cr f T and T are assumed to vanish,
r r8 rx
equation (2.14) as applied to the present problem reduces to
\
I I
e :,t = 1 CTx,t e,tV<2ffx'V L0<5'i,5 + *<:4)J
2 s
X,
ee,t = iCTx(t + iCTe,t + ^(2aeCTx)L
1 + V
3
+ T,
'9x,t E 9x,t 8x
2s
0 (s A) s + if(s.A)J
0(s,A)s + ty(s,A)l
(2.IS)
where the deviatoric stresses are
s = s = (2ct ctq)
x 11 3 x 9
Se = S22 = 3<29 ax1
Sr = S33 = 3(x+V
S0x S12 S21 T0x
s = s = 0 .
r 0 rx
Using these deviatoric stresses, the expression for s becomes
9s 9 3
S ~~ dt 9t L2 SijSij_
"Â¡e i rs ~Â¡ssr3 1
jj ~~ 2 _2 SklSklJ ~St L2 SijSijJ
1 3 P s s 1 3 3 r
s =
2s
_3_ 5
4s
S11S11 + S22S22 + S33S33 + S12S12 + S21S21
]
2 2 2 2
 q(o +
 ot _3 x o x 0 0
J
L [<2ax V4x,t+ <% x)4e,t+ 6TexTex,t]
(2.19)
SUlii (JUT I NE PHI ( P PQ K, J 81 S U YS B t Y H XLAM, OE TA Y XM XN DZ I 0 )
IF(K.GT.O) GO TO 4
YS = ABS(YS)
IF(S.GT.YS)GO TO 1
SZ = BETA*SuRT( CZ)
IF( IO.EQ.O) C = S
P = .. j
PQ O.J
no tc o
1 CON T I NU 
IF(JB.GT.OIGO TC 2
/ = S YS
ZY = Z *(Y1 )
C = 8 *ZY
P = C Y
IF( Ij.GT.OIGO TO 10
CP c*z
c = dp s
PQ = OP/S
G TC n
2 CONTINUE
I F ( m.GT .0 ) GO T 3
SZ = BETA*SQRT(CZ )
C = (S/GLTA)**2
IF(C.LT.CZ) C = (2.0*SQRT(DZ)/B E T A ) *(S YS) DY
PQ  /S 1.0
3 CONTINUE
IF(S.LT.YS) P = 0.0
IF(S.Gfc.YS.AND.S.LE.SZ) P = 2.O SORT(DZ)/BETA 1.0
IF(S.GT.SZ) P = 2 O S/(B E T A * 2 ) 1.0
GO TC lu
4 CONTINUE
C THE DATA USED HERE ARE FROM CRISTESCU (1972) CASE XII
CZ1 = SORT(DZ)
228
139
and
F(A) = f(A)
E
r K(A)
o
XK
XKo [l exp (4)]}
r K
o o
J
$(s,A) = E0 (s A) = E0(A)
EX
/ a \3/2 2/3
3/21 2/3
() }
 1
a/E
$(s,A) = X
r / v 3/2, 2/3
4vi+t(Â¥) 1
 1
xa
(A. 4.10)
(A.4.11)
where
a
xa =
(A.4.12)
Shear Stress, T x10
4.0 _
Final Stress State
Figure 4.17 Stress Trajectories for Data Set 2 Without Radial Inertia
107
FTP(MX) = ET(MX) {STI A*NU*SXI)
LTXP(MX) = E T X ( M X ) QA* TAUI
ITER(MX) = I
PS 10(MX) = PS IP
PHI. (MX) = PH IP
UELTA(MX) = CELTAP
PL AS(MX ) = CELTAP SBPC
CELGAM(MX) = 2.*(E T X(MX) ETXO)
CtX(MX) = EX(MX ) cXC
STKES(MX) = SQRT(SXI**2 + A*STI*(STI
IF(SBPC.LT.SMAXO)GO TO 393
SMAX(MX) = SBPC
393 CONTI CUE
IF(MX.GT.1)GO TO 399
CELLFB(MTU) = CELTAP
399 CONTINUE
IF(KCUNT.EQ. IMAX ) GO TU 400
MTT = MX + MT
IF(MTOT.LT.NMAX )G0 TO 2C0
LINE = LINF + 1
IWRITE = 2
GU TO 4)1
397 CONTINUE
KUUNT = KOUNT + 2
IF(KCUNT.LT.NMAX)GO TO 100
400 CONTINUE
LINE = LIME + 1
IWR I TE = 3
491 CONTINUE
LINEO = LINEO + 1
IFLINEO.LT.INCRTJGO TO 402
LINEO = O
W R I T t ( 6 1QG0) LINE
WRITE(6,1010) (X(L).LltMX, INCRX)
SXI ) )
224
T
Figure 3.7 Location of the Characteristic Lines Passing Through P
m
i
152
r o p
ac
1 0
M, 
alr*e
o
2^9x 0
dv
df
dv
I
"dT
dv
s
e
d9
da
y.
dÂ§
da
0
dT
0X
d?
0 0 0
0 ap 0
0 0 0
0 0 0
0 0 0
0 01
0 0 0
ax,g at, 0
0 0 x, g.
0 0 0
0 0 0
0 0 0
0
0
p
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
xÂ§
0
0
0 0 0 0
0 0 0 0
0 0 0 1
A 0 aX2 0
a 0 aA. 0
2 4
A 0 a A 0
o o
0 0 0 0
0 0 0 0
0 0 0 0
t, 0 0 0
0 ax,^ at,^ 0
0 0 0 x, ^
A,
0
0
0
0
0
0
or, after some manipulation
WRITE!6, 1030 ) JGELLtDZ fOHAT ,XM,XNXLAM
WRI T E ( 6 100 6 ) NU.SY,BDAR,N,XKC t RETA
WRI r E ( 6 100 7 ) JRISEfXVFINt TVFIN .SXJ.TAUO.H
WRITE(6,1006) H1.H2.H3.F4,I.SYS
WRI r::(6, 1' 3 1 ) MI.MJ, I PUNCH KASE
IF(UNCUUP.GT.O) WRITE(6,1C39)
IF! JNCGUP.LE.O) WRITE!6,1040)
I F( I PUNCH.EQ.O)GO TO 15
WRITZ ( 7, 1034) A, MXMAX,I PUNI,I PUN 2,INCRXPINCRTP,DELX.OELT
15 CONTINUE
C THIS BEGINS THE CALCULATION OF CONSTANTS USED
M M A X = 2 *MX MAX
NMAX = MMAX 1
1M A X = NMAX 2
Q = 1. Nil * 2
C2 = SCR T(.3*(1 .NU ) )
uA = 1. + NU
CELT2 = 2.* D E L T
I 2 A = 2 A
Cl = 2 C A
C2 = C/2.
C3 = Q*D ElT 2
C4 = 03*CELT 2
FV1 1. + C 2
CV2 = 1. C2
CL R B = 2 *C2/l)V 1
CLRPI = 1. CLRB
C Y = ABS(SY)
CO 30 L= 1 MX MAX
3J XIL) = (L 1 ) *CELX
C THIS BEGINS THE CALCULATION OF INITIAL CONDITIONS
LINE =
LINER = I NCRTP 1
LINEO = I NCR T 1
206
RHSi = RHSAEM + RHSBEM
380 CONTINUE
IF(A.E0.0 ) GO TO 337
IFIRFSA.LT.1.QEG6.OR.UNCOUP.GT.O) GO T 389
ZA = D2*4S DI D7 S
ZB = ASQ C7S C2 l) 3S
ZC = 0 1 C 3S A5Q*D4S
CEL2 = U3F#ZA + D4F*ZB + l)7F*ZC
TAUI = (RHSF*ZA*RHSG*
SXI = (RHSF*ZB+RHSG*(02*03FA5Q*D7F)+RHSH*(D3S*07FD3F* 07S) )/0LL2
ST I = (RHSF*ZC + RHSG*(A5Q*D4FD1 *D3F)+RHSH*<3F*D4S3S*4F) )/DE L2
GO TO 390
389 CONTINUE
CELS = D 2 F 2 F 2 D1*F5F2
SXI ( R H S3 02 RHSEEMF5F2)/DCL5
STI = (RHSEEM*F2F2 RHS3*01)/DFL5
GO TO 39 J
387 CONTINUE
IF(RFSA.LT.1 .OE96.OR.UNCOUP.GT.O) GO TU 381
DEL 3 = 0 3 F* D4S D3S*D4F
STI = 0.0
TAUI = (RHSF*C4S RHSG*D4F )/DEL 3
SXI = (RHSGD3F RHSF*D3S)/DEL 3
GO TO 390
361 CONTINUE
SXI = RHS3/F2F2
STI J.O
394 CONTINUE
SBP(l) = SGRT(SXI**2 + 3.*TAUI**2 + A*STI*(STI SXI))
STRAIN!I ) = DELTAP
VXP{ I ) = VXI
V T P ( I ) = VT I
SBPC = SBP(I)
IF(I.EG.1)GO TO 340
222
172
and
ir / 2\Jf /
ToP YoP
aRHSE = 2aATV n+ a(Ao^ AT)S +a(A, + a AT)S
rP 2Q Sp/ xP 4Q s /
GP
+ aA T .
5Q P
Now defining the quantities
(D.1.4)
Q4 = 2AT Q3
op'
D, = A AT
1 2Q sp/
lit /
D = A + 2a 5
2 4Q s /
A
AT + Q,
(D.1.5)
and
D4 = 2AT D3
RHSH = RHSE + D,
4 J
the two equations along the vertical characteristics (c = 0) can be
reduced to a single equation by substituting equation (D.1.3) into
equation (D.1.4). This becomes
aDlSxP + aD2S0P + aA5QTP = aRHSH (D.1.6)
The equations along the nonvertical characteristic lines are
given by equations (3.4.4), (3.4.5), (3.4.6), and (3.4.7). These can
be written in a much simpler form by defining the following quantities
as A
Z1 =
Z2 ^
Z3 =
1 v
c
s
2
Z4
1 v
AT
1 
Z5 1 + c,
7 AT
O 1 + C
(D.1.7)
J
129
From eqation (A.2.5),
Bf
3a 3a
pq pq
Bf 3
r# t
L2 1J ljJ
Ba 2f pq
pq
(A. 2. 8)
Now from equation (A.2.8) it can be seen that equations (A.2.4) and
(A.2.6) are identical when
d\ = ~ d\ .
The first invariant of the deviatoric stress is
1 ,
s = a o6 cr, = a a,, = 0
pp pp 3 pp kk pp kk
so that
Bf
a K = a
pq ^ pq
Bf
[J
pq 3a
= f .
(A. 2. 9)
pq
Using equations (A.2.8) and (A.2.9), the expression for the plastic
strain increment, equation (A.2.7), becomes
[h si J [ %i] daki
or
de .
ij
f' (w1
9s s
deP .
ij
ij kl
4f3 f' (/)
9s s
P
ij kl
ij
4f3 f' (/)
[f]
da
kl
kl
(A.2.10)
where
3a,
3eP.
P ij
ij 3t kl 3t
kl
APPENDIX A
CONSTITUTIVE EQUATIONS
A.l. Comments on the Constitutive Equation
The constitutive equation used in the general formulation is given
by equation (2.14) as
3
1 + v
e, =  a. .
 6
.a,, +
ij E ij E ij kk 2
j~0(s,A)s + \r(s,A)J (A.1.1)
where s is the scalar representation of the stress state given by
Â¡3
/ s s
*J 2 ij ij
and s is the deviatoric stress given by
ij
s. a. 6 .a .
xj ij 3 ij kk
Under the assumption that a = T = T = 0, s reduces to
r rx r0
/ 3 r 2 2 2
s = J 2 LS11 + =22 + s33 + 2S
i!
and with
S11 = S,
3(2cx aCTe}
s22 se
3(2acre CJx)
33
= s
3(ax + raB>
S12 S21 S0x T0x
this becomes
2 2 2
cT aa oG + aac + 3T
x x 0 6 0x
(A.1.2)
125
VTRBB = CON2 *VTRB CON4*VTB
VRKD3 = CON2*VRRB + CON4*VR8
Z1 = CFI/O
l = CSI/Q
Z 3 = Z 1 C F I
Z 4 = Z 2 C S I
Z 5 = CELT/DV3
Z = DELT/DV4
RFS = A*A2I*A5I A 3 I *A4 I
RFSA = A BS(RFS)
IF(RFSA.LT. 1 .OF^6.OR.UNCOUP.GT.O) GO TO 360
CVS = A4 I *A6 I A* A 5 I 2
PV6 = A2 I A6 I A 3 I A 5 I
P1F = Z3nV5 A4I
RIS = Z4 *DV 1 A4 I
P2F = Z 3 CV6 A 2 I
R2S = Z4DV6 A? I
BIS = R1S/CSI
81F = RIF/C FI
B2S = Z 2 *RFS
B2F = Z1RFS
RAS = 2. Z6 R2S
6F = 2.*Z5*R2F
368 CONTINUE
IF(I.GT.1)GO TO 370
PS IP = 0.0
RATIO = 0 .0
R A T I 0 P = 0 0
GO TC 330
370 CONTINUO
CALL PS I(PS IP,JRATE,SBPCDEL TAPDY.SY.BLTAXKO,DHAT,DZ)
SBUIF = SBPC SBB
OELTAP = DELTAB + SBUIF + DELT2MPSIP + PSIB)
lFISGPC.LT.SMAXB.ANU.I.LT.MJ)GO TO 371
215
204
O
S = #
y e
1500 psi
10.2 X 10 psi
4
1.47059 x 10
4
A = e = 1.47059 x 10
y y
A = e_
.002575
P =
5.6 X 10 psi
10.2 x 106 psi
, 0054902
A = 6 = .0004
m 3.25 x 10 psi
xm = =
E 10.2 X 10 psi
.0031863
(E.4.4)
xn
4
n 9 x 10 psi
J7 0
10.2 X 10 psi
= .0088235
h = .000004
X = .00005
r k 3 1.
xko = 4^ = (2.5 1n)(A. ec. 2= .0012
c 2.08 x 10 in/sec
where the dimensionless coefficient for the coefficient in the rate
dependent term was obtained for a .5inch diameter aluminum tube.
E.5. Listing of the Program
In this section is the listing of the computer program which was
used to obtain the solution in the characteristic plane to the wave
propagation problem presented in this paper.
103
where now S represents the component of the stress vector (s) in the
N
SrS^ plane. By plotting a new stress traj ectory of T versus S we
x b N
can effectively represent this threedimensional stress trajectory as
a twodimensional one; this is shown in Figure 4.15. Again from this
figure it can be seen that at distances greater than 1.0 diameter from
the impact end the tube exhibits neutral loading followed by loading.
This is the same general behavior as the case when radial inertia is
not included, except that now the stress trajectory is not planar.
This is easier to visualize if these same stress trajectories are pro
jected onto the plane T = 0; the stress trajectory at X=1.0 is shown
in Figure 4.16. Since in any experimental wave propagation problem of
this type radial inertia effects are present, interpretation of the
experimental data without considering the hoop stress could lead to the
false conclusion (for a set of initial conditions and boundary condi
tions similar to these) that unloading occurred. Therefore, extreme
care must always be exercised when interpreting experimental results.
The other two computer runs made used artificial data where the
initial conditions and the boundary conditions were described by the
input data
" tr = 40
AX = AT .025
S = .0006040
J X
Data Set 2 <
o
T =0
V =0
Xf
V0f = 00100 .
42
b' 4a'A = 0 .
4
When y = > equations (3.1.18) and (3.1.19) are
11 2 r 2
A = 1 + (pr 1) (cos 62 a/3 sin 6 cos 6 + 3 sin 6)
1 4 P
A^ = v + ~(tt 1) (cos"6 3 sin26)
A
a3 = 0
1 i ' 2
A = 1 + (g 1) (cos 6 + 2 a/3 sin 6 cos 6+3 sin 6)
T l)
4V3
A5 =
Ag = 2(1+v)
J
and
a = A A .A aA A
14 6 2 6
A1A4 + A4A6
aA,
and using the same manipulations as in Appendix B, Section 4,
equation (3.1.21) becomes
0 = WV + ^2
where A^, A^, A^, and Ag are given by equation (3.1.22). Now
f ^ = f^(6) = i ^cos^6 2 f/3 sin 6 cos 6 + 3 sin26^
(6)
i j^cos26 + 2 ^3 sin 6 Cos 6+3 sin26^
f = f (6) = i r
3 3V 4 L
z = z(P ) = g
c c 6
cos26 3 sin6
]
(3.1.21)
(3.1.22)
(3.1.23)
(3.1.24)
defining
/ (3.1.25)
7
then subjected to a suddenly applied incremental torsional load. The
strain caused by this incremental load was found 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
f
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
63
This gives the value of the coefficients at a point nearer the center
of each grid element. For this work, the value of 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 Taxis
(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 T T and T can be written as
L Â£ *3
2c
2c,
2c.
T, =
1 1 + c
AT
T =
2 1 + c.
AT
= r
3 1 + c
AT
and the interpolation constants for the points LB and RB are
2c,
CLRB
AT 1 + c.
Cl,RBI = 1 
T2 ^
AT ~ 1 + c.
Using subscripts to denote the grid point, the values of any quantity F
at the points LB and RB are
FT = CLRBFt + CLRBI F
LB L B
F = CLRB*F_ + CLRBIF
RB R B
Op OOOOIOOOO
OOOapOOOOOOO
OOOOOp 0000 1
1 000000 A O aA 0
X Cj
0 0 0 0 0 0 0 aA2 0 aA4 0
0 0 0 0 1 0 0 A 0 aA^ 0
d 5
x,g t,^ 0 0 0 0 0 0 0 0 0
0 0ax,^at,^0 0 0 0 0 0 0
0 0 0 0 x,_ t,_ 0 0 0 0 0
Â§ ?
0 0 0 0 0 0 x,g t,g 0 0 0
0 0 0 0 0 0 0 0 ax,^ at,g. 0
0000000000 x,,
0
V
X,X
0
0
Vx,t
a cr0
0
r
o
0
V
r, x
0
>1
CO
Vr ,t
Ilf
T X
a5
V0, x
a(r+e)
o
6
V0,t
2^0x
0
a
X X
dv
X
db
0
x, t
dv
r
a dÂ§
0
CTe,x
dv0
dÂ£
0
CTe,t
da
X
dÂ§
0
T0X,X
da0
a di
tÂ§
T0x,t
dT0x
dT
(B.2.14)
149
10
by an increase in shear stress at a constant value of longitudinal
stress and then the slow wave caused loading such that 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.
BIOGRAPHICAL SKETCH
Charles Daniel Myers was born July 27, 1945, in Brevard County,
Florida. He was graduated from Terry Parker High School in Jacksonville,
Florida, in June, 1963, and received his Bachelor of Science degree in
Engineering Science from the University of Florida in June, 1968.
He was employed in the Structural Analysis Section of the Martin Marietta
Corporation in Orlando, Florida. While on a leave of absence from the
Martin Marietta Corporation, he received his Master of Engineering
Degree in Engineering Mechanics from the University of Florida in
December, 1970. He returned to the Martin Marietta Corporation in
September, 1972. He is a member of Tau Beta Pi (having served as presi
dent of the Florida Alpha Chapter during the year 19671968), Sigma
Tau, Omicron Delta Kappa, Phi Kappa Phi, and the Society of Engineer
ing Science. He was recently selected for Outstanding Young Men in
America 1973.
Charles Daniel Myers is married to the former Peggy Lee Hufham
of Jacksonville, Florida. They have one son, Daniel Lee.
236
142
methodthe system of equations is written in the form of equation (2.38)
and the equations along the characteristics are given by
1T A dw = ir b dt
(B.2.1)
where 1 is the left eigenvector and is determined from
1T (cA B) =0
(B.2.2)
which becomes
L1! *2 x3 x4 x5 X6J
Cp
0
0
1
0
0
0
acp
0
0
0
0
0
0
cp
0
0
1
1
0
0
CA1
acA
2
CA3
0
0
0
acA
2
acA^
4
acA
5
0
0
1
CA3
acA_
5
CA6
This is the matrix form of the following six simultaneous algebraic
equations for the six components of the left eigenvector
pl1 + l4 = 0 ^
acpl.
= 0
P13 + 16 =
^1 + CA114 + acA215 + CA316 =
acAl + acA 1 + acA_l_ = 0
2 4 4 5 5 6
1 + cA 1 + acA 1 + cA 1 = 0
3 3 4 5 5 6 6
(B.2.3)
However, from equation (2.38) the equation for the characteristics is
cA B
0
75
V6P = iT~ LRHS9 RHS10
+ KHS11(D4s2B7s D4f2B7s)]
SxP = 4^ [B2fS11> B2s(D2RHS9 + B^RHSll)]
S6P = 4^ LB2Â£
(3. 5.28)
(3.5.29)
(3. 5. 30)
where again is given by equation (3.5.6). When radial inertia
effects are not included, and vanish, and from equations (D.5.38)
and (D.5.39), the solution at P is found to be
V
S = (B RHS10 B rhS9) (3.5.32)
xP A.,., 2f 2s
At a Boundary Point (X=0) for Uncoupled Wraves
When the waves are uncoupled, the solutions to 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 is given by
equation (3.5.6), and in all cases when radial inertia terms are not
included both V and vanish. In all four cases the solutions can
rP 0P
be found in Appendix D.
Case I: Traction boundary conditions
When Tp and S^p are known at a boundary point, then from equations
(D.6.1), (D.6.2), and (D.6.3) at that point
79
e = e + 2V V V
xP xB xR xP xB
:9xP eexB + V0R 2
E0P = e0B + 2 iT
where equations (3.6.6) and (3.6.9) are the same expression.
(3.6.7)
(3. 6. 8)
(3.6.9)
APPENDIX B
CHARACTERISTICS AND EQUATIONS ALONG THE CHARACTERISTICS
B.l. Equations for the Characteristics
The equations necessary for determining the characteristics are
given by equations (2.38) to (2.44). Thus, the characteristics c are
given by the equation
or
acp
Co
0
0
1
0
0
0
acp
0
0
0
0
0
0
cp
0
0
1
1
0
0
O 1
>1
H1
acA2
CA3
0
0
0
acA2
acA.
4
acA_
5
0
0
1
O
>
CO
acA,_ cA^,
5 6
cp
0
1
0
0
0
cp
0
0
1
1
0
CA1
acA2
CA3
= 0
0
0
acA2
acA
4
acA
0
0
1
CA3
acA
o
CA6
0
(B.1.1)
140
81
boundary increases linearly up to its final value (denoted by Vx^
or Vgf) during a period of time called the rise time (T ) and then
remains constant. That is
'f
if 0 < T < T.
R
R
Ve (X = 0, =/
if T > T.
R
r
if 0 < T < T.
R
V (X = 0) = /
if T > T.
R '
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
Plastic Wave Speeds,
Figure 3.2 Plastic Wave Speeds as Functions of Â¡3 and y for Poisson's Ratio of 0.30
/
7 ex
X
Figure 3.1
Yield Surface Representation in Spherical Coordinates
co
k*
Hoop Stress, SQxlO
109
Figure 4.19 Hoop Stress Versus Longitudinal Stress for Data Set 2
With Radial Inertia
127
If 0(s,A) is assumed to vanish in equation (A.1.1), and if
t (s ,A) = yf (y l) ,
\''3 K
then the constitutive equation (A.1.1) reduces to the generalized semi
linear constitutive equation of Perzyna (1963).
The constitutive equation (A.1.1) can also be reduced to the equa
tion obtained from the rate independent equations of incremental plastic
ity as v/ill be shown in the next section.
A.2. Rate Independent Incremental Plasticity Theory
The constitutive equation derived in this section is obtained when
the material is loading in the manner outlined by Malvern (1969).
If the material is strainrate independent, if it exhibits isotropic
workhardening, and if its yield function, f, depends only on the state
of stress, then
f(C. .) = F (A.2.1)
ij
where F is the workhardening parameter. If F, and, consequently, the
size of the yield surface, is assumed to depend only on the total
P
plastic work, w then
F = F(vF) = a* (A. 2.2)
where vF = f a.. de^. (A.2.3)
J ij iJ
*
and ct is the magnitude of the stress vector.
During loading the material is assumed to obey the PrandtlReuss
flow law
ds.
ij
dX s..
ij
(A.2.4)
where d\ is a scalar function, and the von Mises yield condition
f = a
rr
y 2 s s
ij ij
s
(A.2.5)
144
The fourth equation of equations (B.2.3) has not been used. However,
when equations (B.2.5), (B.2.7), and (B.2.8) are substituted into it,
the resulting equation is the characteristic equation (B.1.2) and no
new information is obtained. From these results, the left eigenvector
can be written as
1 =
1
pc
1
pc
0
aA A A A
Ad o 4
 2 4
A A aA
4 o 5 Z '
Pc
^3^5 ^2^6 + 2
V6 ^5
A4
2
pc
aA A A A
Ad
 2 A4
A4A6 ^5
(B.2.10)
A dw
(2.
39) and
(2
41),
p
0
0
0
0
o ~
dv
X
0
ap
0
0
0
0
dv
r
0
0
P
0
0
0
dvQ
0
0
0
A1
^2
>1
CO
da
X
0
0
0
aA2
^4
^5
dCTe
0
0
0
A3
aA
D
>i
02
1
K1
CD
ir
1
TO PEGGY
Longitudinal Strain,
Figure 4.1 Grid Size Effects on the Longitudinal Strain at X = 3.75
oo
CO
TABLE OF CONTENTS (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 152
APPENDIX C. PROGRAMS FOR DETERMINING THE PLASTIC
WAVE SPEEDS 164
APPENDIX D. SOLUTION TO THE FINITE DIFFERENCE EQUATIONS
IN THE CHARACTERISTIC PLANE 171
D.l. Equations for Fully Coupled Waves 171
D.2. Equations for Uncoupled Waves 17&
D.3. Solution at a Regular Grid Point
for Fully Coupled Waves 130
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.l. General Description of the Program I94
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
v
128
so that'the assumed flow law is derivable from the plasticpotential
equation
o 3f(o\ .)
deP = dX y 1J
(A. 2.6)
The form of the von Mises yield condition, equation (A.2.5), is chosen
so that under uniaxial tension, f = a .
x
From equations (A.2.3) and (A.2.6) an expression for d\ can be
obtained as
dW*3 = a. ,deP = a. ,d\
and
1J 1J ij oaTj
dX
dW1"
df
ij So..
ij
de
p
P iJ
ij df
^mn da
mn
Now, from equations (A.2.1) and (A.2.2)
so that
and
da, = ^ dW1" = dvF = f' (vP) dvF
dCTki kl 9W
F
dr
dvF
1 Sf
33 CTki
de
F'GV13) ~wkl
df df
p v
1J f'(W)P a
mn oa
mn
(A.2.7)
87
(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
f T =4.00
R
AX = AT = .050
s = 0
Data Set 1 /
T = .0003480
V = .002404
Xf
V. = .0001106
l Gf
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
R
has passed the point X=3.75 at the time when the longitudinal strain
Shear Stress, T x10
108
106
These data represent the impacting of the end of a tube with a torsional
(transverse) velocity pulse when the tube is preloaded in compression
above the yield stress. First the stress trajectories obtained by not
including radial inertia effects will be examined. The stress trajec
tories for this case are shown in Figure 4.17 at various distances from
the impact end. Here the stress trajectory develops (at distances from
the end of the tube greater than 1.0 diameter) into the type given qual
itatively by Clifton (1966), that is, the tube undergoes unloading
along the S^axis followed by an increase in T at a constant value of
S^_ until the yield surface is reached followed by loading normal to the
yield surface. Again the stress trajectory does not exhibit this
behavior near the impact end, but rather approaches this behavior as
the distance from the impact end of the tube increases. As explained
above, the stress trajectory near the end of the tube does not behave
in the manner given by Clifton (1966) using the simple wave solution
because the fast and slow waves are not distinct near the end of the
tube and because the input boundary conditions force the solution near
the impact end.
The stress trajectories (T versus S ) are shown in Figure 4.18
N
for the case when radial inertia effects are present. The stress
trajectories are plotted using the combined hoop stress and longitu
dinal stress of equation (4.3.1) since this gives an effective three
dimensional representation of the stress trajectories. The stress
trajectories exhibit a "ringing" effect after loading has taken place.
This "ringing" is caused by the hoop stress as can be seen from
Figure 4.19.
113
inertia effects). For all four cases the initial conditions and the
boundary conditions were specified by
. 40
T
R
AX = AT = .025
o
0
x
o
.0000900
T
.000250
x
f
These input data represent a tube with a static pretorque which
is impacted longitudinally at the end. The longitudinal strain versus
time curves obtained by using these data are shown in Figure 4.22.
When no strainrate dependence is included, the results are the same
qualitatively as those obtained earlier for the Lipkin and Clifton
(1970) data. However, when strainrate dependence is included (and
the tube is subjected to constant velocity impact), the maximum longi
tudinal strain is reduced.
This lower value of total longitudinal strain may be better under
stood by considering a longitudinal impact of an unstressed tube of
a strainrate dependent material. The stressstrain curve for such
a material is steeper at higher strain rates; and therefore, if the
maximum applied stress is the same for two different loading conditions
the total strain will be larger for the condition when the rate of load
ing (or straining) is smaller. For this combined stress case, the
44
When a = 1,
2 2 ~ 2 2
(cosJ5 + 2 */3 sin 6 cos $ + 3 sin 6) (3 sin 6 cos
and equation (3.1.26) becomes
[2(1 + v) 1 iv2]
z
c [(2) (1 + v) 1] f2 f1 + 2vf3
 (1+v)2
(3.1.29)
z
2
(1 + v) sin 6 cos 6 + 3v sin 6
A short computer program was written to calculate the critical
values of 3 (using equation (3.1.28) when a = 0 and equations (3.1.29)
and (3.1.27) when a = l) 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 3 in this range are
plotted in Figure 3.3. For all other values of 6, there is no phys
ically possible critical value of 3; that is, there is no value of
3 such that the fast and slow wave speeds are equal at y =0.
For the case when cQ = 0 (6 = 60), for any value of v the critical
W
value of 3 is smaller when radial inertia effects are included.
137
A. 4. Dimensionless Expressions for the Functions
0(s,A) and ijr(s,A)
Using the definitions for the dimensionless variables in equation
(3.2.1), the expressions for 0(s,A) and if(s,A) given in Sections A. 2
and A.3 can be converted into dimensionless form. First, when the
uniaxial stressstrain curve in equation (A.2.13) is used, the functions
0(s,A) and ir(s,A) are given by equations (A.2.18) as
(A.4.1)
and using equation (3.2.1) the dimensionless form of these variables is
r
(A.4.2)
and
5(s,A) = E0(s,A)
EBn(s ct )n
y
n / \nl
 E Bn(s s )
y
$(s,A) = Bn(s Sy)n 1
(A. 4.3)
where
n
(A.4.4)
B = BE
When the uniaxial stressstrain curve of equation (A.2.21) is
used, the functions 0(s,A) and ir(s,A) are given in equation (A.2.25) as
77
Case IV: Mixed boundary conditions
When V^p 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.,
v9p = 4 <3'5'43)
and S^p and Sgp 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 S S., T V V and
Xu X P
Vg 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
9e Sv
x x
T ~ X
(3.6.1)
Se
i Bve
9x
dT 2 5x
(3. 6. 2)
de6
~W = 2Vr
(3.6. 3)
For a regular grid element, these equations can be written in
finite difference form as
232
Cristescu, N. (1965), "Loading/Unloading Criteria for Rate Sensitive
Materials," Archiwum Mechaniki Stosowanej, Vol. 17, pp. 291305.
Cristescu, N. (1967a), Dynamic Plasticity, New York: John Wiley and
Sons, Inc.
Cristescu, N. (1967b), "Dynamic Plasticity under Combined Stress,"
Mechanical Behavior of Materials under Dynamic Loads (U.S.
Lindholm, Ed.), New York: SpringerVerlag, pp. 329342.
Cristescu, N. (1968), "Dynamic Plasticity," Applied Mechanics Reviews,
Vol. 21, pp. 659668.
Cristescu, N. (1971), "On the Coupling of Plastic W'aves as Related to
the Yield Condition," Romanian Journal of Technical Science,
Applied Mathematics, Vol. 16, pp. 797809.
Cristescu, N. (1972), "A Procedure for Determining the Constitutive
Equations for Materials Exhibiting Both TimeDependent and Time
Independent Plasticity," International Journal of Solids and
Structures, Vol. 8, pp. 511531.
Davis, E. A. (1938), "The Effect of Speed of Stretching and the Rate of
Loading on the Yielding of Mild Steel," Transactions of ASME,
Vol. 60, pp. A137A140.
DeVault, G. P. (1965), "The Effect of Lateral Inertia on the Propa
gation of Plastic Strain in a Cylindrical Rod," Journal of the
Mechanics and Physics of Solids, Vol. 13, pp. 5568.
Donnell, L. H. (1930), "Longitudinal Wave Transmission and Impact,"
Transactions of ASME, Vol. 52, pp. 153167.
Duwez, P. E. and Clark, D. S. (1947), "An Experimental Study of the
Propagation of Plastic Deformation under Conditions of Longitu
dinal Impact," Proceedings of ASTM, Vol. 47, pp. 502532.
Efron, L. and Malvern, L. E. (1969), "Electromagnetic Velocity
Transducer Studies of Plastic Waves in Aluminum Bars,"
Experimental Mechanics, Vol. 9, pp. 255262.
Hencky, H. (1924), "Zur Theorie Plastischer Deformationen und der
Hierdurch im Material Hervorgerufenen Nebenspannungen,"
Proceedings of the First International Congress for Applied
Mechanics, Delft, pp. 312317.
Hill, R. (1950), The Mathematical Theory of Plasticity, England:
Oxford.
Hopkins, H. G. (1961), "Dynamic Anelastic Deformation of Metals,"
Applied Mechanics Reviews, Vol. 14, pp. 417431.
Shear Stress, 7 x10
1. 2
rr
Classical
With Radial Inertia
With Rate Dependence
With Radial Inertia and Rate Dependence
Figure 4.24 Stress Trajectory at X = .25 for Data Set 3
117
27
_2 r 2   ry (d 2aa )
+ a(pc ) (pc ) (A1AÂ¡593)A I j ^ +  iHs,A) ldt
^ o 2s
s,A)J'
2 r 2 21 r1 sx i
+ 2(p Oj^Cp O (14aA^)A4J i)r(s,A)Jdt .
(2.53)
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
(p2)(ay5 y4>
(pc") (A A. A A ) A
1 3 O
and using equation (2.48), equation (2.51) is obtained; or by multiply
ing equation (2.52) by the quantity
pc_(pc ) (A1A4 aA2) A4J
(p c")(xa5 23) 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, A and A vanish.
o
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)
3 5
and (B.4.6), respectively. Under these conditions, equation (2.52)
also reduces to the form of equation (B.4.4).
132
r a* a i
/ = B(a* a )n I a* Zl
y L n +1 J
W = (a*a )(na*+a )
n+1 y y
and
dW13
(vf)
dF
dvf
da da
dW
i = Bn(a*) (a* a )n 1
f'(/) y
(A.2.16)
Using equations (A.2.5) and (A.2.16) in (A.2.12), the expression for
the plastic strain rate for the static stressstrain curve of equation
(A.2.13) becomes
oD / \nl 3a*
.P 3Bn ~5t
e. =
ij
2a
ij
(A. 2.17)
Comparing equation (A.2.17) with equation (A.1.1), the functions
0(s,A) and ir(s,A) for the rate independent plasticity theory using
equation (A.2.13) are
t(s,A) = 0
0 (s,A) = Bn(s ay)n
(A.2.18)
where s = a The constitutive equation (A.1.1) during loading for
this case is
1+v v c 3 / \n1 ij
e. = a.. 6 . a, + Bn( s a ) s 
ij E ij E ij kk 2 y 
(A.2.19)
The second form of the uniaxial stressstrain curve is the one
used by Cristescu (1972) as the relaxation boundary for dynamic load
ing. Although Cristescu (1972) uses a quasistatic curve which is
ACKNOWLEDGMENTS
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 my 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.
iii
WRITE 7* 1C35) (EX(L ),L = 1,MX,INCRXP)
WRITE<7, 103 5 ) (VX(L),L=1, MX,INCRXP)
W R I T E ( 7 1335) (CEX(L ), L = 1,MXINCRXP )
IF ( A.EC1.C )G0 T 403
WRI TCl7, 103 5) (ST(L ) L= 1 ,MX,INCRXP)
WPITb(7, 1035) (ET(L ) ,L=lfMX,INCRXP)
WRI TE(7, 1035 ) (VR(L ) ,L=1 ,MX,INCRXP)
403 CUNTI l\IUL
I F ( IPUN2.EQ.OGO T0 404
WRITE!7,1035) (TAU(L ),L = 1MXvINCRXP)
WRI TE(7, 10 35 ) (VT(L ) ,L=l,MX INCRXP)
WRITE(7, 1035 ) (ETX(L),L=1,MX,INCRXP)
WRITEI 7, 103 5 ) (CELGAML),L = 1,MX,INCRXP)
404 CONTINUE
IF! IWRITE 2) 90, 397,405
4C5 CONTINUE
LODO FORMAT(' 1 ,2 X,'LINE =',I5///)
100 1 FORMAT!3 15,5F5.2,4F10.2)
1002 FORMAT!I5.7F10.3,15)
ICC 3 FORMAT(1015)
1004 FORMAT(2 E207)
1005 FORMAT( 1 ,9X,A = , I 1, 15X, 'JR ATE = ,I 1,11X,ELX = *,F9.6,4X,
I'OlLT = ,F9.6,4X, MXMA X = ,I 4,BX, SMALL = ',F8.6/)
1V.6 FORMAT ( 19X, NU = F 5.3 1CX SY = F12.9,3X 1 BOAR F 9. 1,4 X N
1= ,F7.4,9X, 'XK:: = F 8.6,6 X BETA = ',F8.6/)
1007 FORMAT!1 O X,J RIS F
i,*SXO = ,F9.6,5X
iC C 8 FORMAT(1TX, 'HI =
1= F4.2, 11X, n
10 0 9 FORMAT(/3 X, 'TIME
1010 FORMAT(/ IX,1X =
1011 FORMAT!/3X,SX =
1012 FORMAT(/3X,ST =
1013 FUKMATI/3X,'TAU =
= I 4,BX, 'XVFIN = ',F 9.6,3 X, 'TVFIN = *,F8.6,3X
TAUO = ,F8.6,5X, H = ',F8.6/)
,F4.2,11X,H2 = ,F4.2, 11X, *H3 = ,F4.2,11X,H4
12, 13X, SYS = ,F 1 2.9/)
',4X,71F9.56X)F9.5/(8X8(6X,F9.5)))
',4X,7!F9.5,6X),F9.5/!8X,0(6X,F9.5)))
,8 E15.5/( 10X.8E15.5) )
' t 8E155/(10X,6E15.5))
.8E15.5/10X.8E15.5))
226
Shear Stress, t x10
Figure 4.14 Shear Stress Versus Longitudinal Stress for Data Set 1 With Radial Inertia
102
176
The four equations (D.1.18), (D.1.19), (D.1.20), and (D.1.21) along
with equation (D.1.6) now become a set of five simultaneous algebraic
equations which must be solved during each iteration for the unknowns
VxP V0P TP SxP and Sep point P of a regular grid element. The
reason for calculating the coefficients by the method outlined in
Chapter 3 is now apparent. If the coefficients had been calculated in
the normal manner, these five equations would have to be solved simul
taneously for each iteration at a regular grid point when the waves are
fully coupled. However, by calculating these coefficients by the method
used here, these equations reduce to two separate sets of simultaneous
equations, one set with two equations and one set with three equations
(this will be shown in Section D.3), and the computation time required
to solve these two sets of equations is significantly less than the
time to solve one set of five equations. Once the solution is found
at a regular grid point, the radial velocity is obtained from equation
(D.1.3). This will be done in Section D.3.
At a boundary point the solution is obtained by specifying two of
the variables (as described in Chapter 3) and then solving equations
(D.1.19), (D.1.21), and (D.1.6) simultaneously. Again equation (D.1.3)
is used to get V This will be done in Section D.5.
D.2. Equations for Uncoupled Waves
When the waves are uncoupled, the finite difference equations
along the characteristic lines are given in equations (3.4.9) to
(3.4.14). These expressions will now be simplified by grouping known
126
If the only stress present is a then
x
S = CT
and from equation (A.1.1) the plastic portions of the strain rates are
P *P *P
e~ = = = 0 .
'6x 9r
rx
P P 1P
e. = e = e
9 r 2 x
If, in addition, each plastic strain component is monotonically increas
ing (or monotonically decreasing), then the expression for A in equa
tion (2.16) becomes
CT
._ .'P.2 P 2 P.2 x
A = / I / CO + (eQ) + (er) dt +
P , X
e dt +
x E
A 
P e
e + e
x x
A = e
and equation (A.1.1) reduces to
1 3 f
e = ct + Â¡
x E x 2 L
0(ct e )ct + A(ct
X, X X X
, X J
3 CTx
c = ci+0(CT,e)cr+t(CT,e)
x Ex xxx xx
(A.1.3)
which is the same as equation (2.13) used by Cristescu (1972).
Again, when the only stress present is ct the function 0(cr e )
X XX
can be set equal to zero and
'Kct e ) = i g(o e )
xx E xx
to obtain the semilinear constitutive equation of Malvern (1951a,1951b) as
1 1 ,
e = ct +g(CT,e)
x E x E xx
52
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
Shear Stress, T x10
6 _
4
Figure 4.12
Stress Trajectories for Data Set 1 Without Radial Inertia
to
o
26
O =
(p c2) (A
A .X aA2) A. dv pd aA A A A dvQ
46 5' dj x M l. 2 o 3 4j 9
3A A A A
Z D O <Â£
]dV
pe
(p c2) (A4A6aA2) A4J dav +
(p c2) (A A aA2)
x L 4 6 5
d]
(2o acr )
x 0
2s
']
i>(s,A) dt
v (o 2actq)
o ] rv (.0 j
+ al (pe ) (A A.A A )A + t(s,A) dt
2 6 o 5 2J [_r _l
o 2s
2
2(pe )
~A3A4J [^(s,)Jdt
aA2A5
2s
(2.51)
0 =
(P2)(A2A6A3A5)A2]dVx^
(pe ) (A1A5A2A3)A5 Â¡dvf
J
4Â£(Ps2)(A1A523)xjdTex+
pe pe
ex' 2Llfc(Vei3M%
ti
d][
(2a aa )
x 0
2s
iKs>A)Jdt
1j[(pS2)2(V6
pe
+ 2
1 FV 'u j i
a2) (pc )(A +A ) + l j ZL + t(s,A) dt
x D L_r0 2s J
i r3t
.(PC )(AiA5A2A3)A5J
0X
^ 2
, A)J <
t(s,A) dt
(2.52)
23
0 = p c
n
aAA AA I dv o c
2 5 3 4J x
(p c2) (A1A4aA2) A J dv
e
[(pc2)(A1A4aX2)Xjd
VsvJ
9x
, 2.
+ (p c )
aA A A A. da
. 2 5 3 4J x
a 2,2
+ (p c )
(2a aa.)
x 6
2s
,A)J'
f(s,) dt
CEL r A(MX) = CELTAO
CPLAS(MX) = DPL AS')
SMAX(MX) = snc
CEL GAM(MX) = 0.0
C EX(MX) = 0.0
STKES(MX) = SORT(SX0**2)
70 CONTINUE
CELLHR(l) = CELTAO
C THIS ENOS THE CALCULATION OF INITIAL CONDITIONS
I W R I T E = 1
GO TU A01
9^ CONTINUE
KUM = 1
100 CONTINUE
MX = 1
JERROR 0
CO TO 300
200 CONTINUE
MX = MX + 1
'30 3 CONTINUE
MT = MX + KLJUNT
T(MX) = M T D E L T
SXk = SX(MXfl)
STR = ST(MX + 1 )
V X R = VX(MX 1 )
VTR = VT(MX+1)
VRR = VRIMX+l)
TAUR = TAU(MX + 1 )
T AUG = TAU(MX)
SXC = S X(MX)
STB = ST I MX)
VXD = VX{MX)
VT B = VT(MX)
VRB = VR(MX)
to
o
121
along the tube exhibited neutral loading as the fast wave passed,
followed by loading normal to the yield surface as the slow wave
passed. However, at distances closer than 1.0 diameter from the
impact end the stress trajectory at a point along the tube followed
a path in stress space which represented a combination of neutral
loading and loading. This behavior was caused by the finite rise time;
since the velocity pulse was not a step function, it took about 1.0
diameter for the two waves to separate. In other words, the fast
wave starting at T = T overtook the slow wave starting at T = 0 in
R
about 1.0 diameter.
When radial inertia effects were included, the stress trajectories
became more complicated as they were no longer planar. For this case
the longitudinal stress was plotted against the shear stress and it
appeared that at distances greater than 1.0 diameter from the impact
end unloading occurred. However, when the stress trajectories were
plotted in the threedimensional stress space the behavior was found
to exhibit neutral loading followed by loading.
The strains and velocities were shown at several tube locations
and constant state regions were observed at points along the tube
after the fast wave had passed and before the slow wave had arrived.
Also, radial inertia was shown to have little effect more than two
diameters from the impact end.
For the second case when the impact velocity was torsional, the
same types of results were obtained. Here the stress trajectory
obtained was like the one predicted by Clifton (1966) with his simple
wave solution.
Critical Value of
Figure 3.3 Values of 3 at Y = 0 for which c = c = c
' f s 2
160
Under these conditions (A = A =0), the equation for the nonvertical
o o
characteristics of equation (2.45) becomes
(pc)2(14A6 aA2A6) (pc2) (^ + a2) + 4 = 0
and as in equation (2.48)
2 2 WA 
C =
2p(X146 aX2X6)
2 ^A1A4 + A4A6 ~ ^2^ ^ iA4^A6 Al^ + aA2^
2pA6(AiA4 aA2)
and using equations (2.49) and (2.50),
2 (A1A4 + A4A6 aA2) + (A4A6 ~ A1A4 + aA2)
2PA6(A1A4 aA2}
2
Cf
A.
p (V4 ~ ^
(B.4.1)
2 (A1A4 + A4A6 aA2) ~ (A4A6 A144 + aA2)
2PA6(A1A4 ^
2 1
c =
PAC
(B.4.2)
and the four nonvertical characteristics are
c = cf, cs.
The equations along the characteristics will be obtained from
equation (B.2.23). For A =A =0, this equation becomes
o O
o = [(P;2)(6)i][a4][
dv + ~ da + i dt
x 2 x Tx J
c pc
[;2)V1][Â¡2]Lr*Jdt
+ a
(B.4.3)
riFFl = ADS(SOP{ I ) SB P( I 1 ) )
DIFF2 = ABS(STRAIN! I ) STRAIN! II) )
CIFF3 = ADS ( VXP ( I ) VXP(ID)
CIFF4 = A BS(VT P( I ) VTP(Il))
CIFF = U1*DIFF1 + H2 C 1 FF2 + H3*DIFF3 * H4*DIFF4
DENOM = A BS ( VX P ( I ) ) + ABS(VTP(I))+ ABSISBPUJ) +
FRRCR = CIFF/DENOM
IF(ERROR.GT.SMALL)GO TO 340
430 CONTINUE
IFIMX.GT.1)00 TO 391
Â£X(MX) = Â£X(MX) + 2*VXR VXI VXB
ET X(MX) = ETX(MX) + VTR .5*(V TI + VT3)
GC TO 392
391 CONTINUE
FX(MX) = EX(MX) + VXR VXL
ETX(MX) = ETX(MX) + .5*(VTR V TL)
392 CONTINUE
VR{MX ) = A*(D3 03*STI )
ET(MX) = ET(MX) + DELT2*(VR(MX) + VRB)
TAU(MX) = TAUI
VX(MX ) = VXI
VT f M X ) = VT I
SX(MX) = SXI
ST(MX) = ST I
Al(MX) = A1P
A 2(M X) = A2P
A 3 ( M X ) = A 3 P
A4 ( MX ) = A4P
A^tMX) = A5P
A6(MX ) = A6P
CF(MX) = CFI
CS(MX) = CSI
S 3(MX) = SBPC
EXP(MX ) = EX(MX) (SXI NU*STI)
AB S(STRAIN(I))
to
to
141
or
a(cp )'
or
a(cp )
or
,
8
O
O
1
0 Cp 0
1
0 cA^ acA2 cA2
1 0 acA,
2 Zl3
+ acp
0 acA acA acA
i nr D
0 0 acA acA
4 5
1 cA acA cA
0 1 acA
cA
5
6
5 6
cA. acA cA
0 cAw acA
12 3
1 2
ac.4 acA, acA
2 4 5

 2
a (cp )
0 acA acA,
2 4
 acp
cA acA,. cA
1 cA acA
3 5 6 J
3 5
3d~ 
2
 2 2 "1
2 24
cl A A,A + 2aA^
Â¡ir
 A A 
aA A aAA
an c
L 1 4 6 2
3 5
3 4
1 5 2 6J
M 1
2 lT 2*1
_ r
1
c Â¡ A A aA 1
L 4 6 5J
acp_
 ao\J
= 0
cp 0 1
0 acA^ acA^_
1 acA cA
o o
V4 all\
and using the definitions (2.46), this becomes
2 36 2 24 2 2
apc(a)apc(b)+apcA4 =0
^a2oc2J j (p c2) 2 (a) (pc2) (b) + aJ = 0
(B.1.2)
which is the equation for the characteristics given in equation (2.45)
B.2. Equations along the Characteristics
The equations (2.48) along the nonvertical characteristics can be
obtained in two different ways. One form of the equations along the
characteristics will be found below by each method as an example; then
the other forms of these equations will be given. Using the first