é, firs".

.- 0‘.

 

-..u..u........‘...;',;:.x. .V

 

NUMERICAL ANALYSISOF comammiamoma AND; ~ '

TORSION or- A ‘PRISMATIC BAR USING RICtDL-PLASTIC-

AND WORK-HARDENING PLASTICITY THEME-$7

. Thesis {of the Define of .PhQD.‘  
MICHIGAN. STATE UNIVERSITY
Patrick M. Milier
19.6.6

 

aways

I LIBRARY

Michigan 5: to
University

 

 

 

 

 

 

 

ABSTRACT

NUMERICAL ANALYSIS OF COMBINED BENDING AND TORSION
OF A PRISMATIC BAR USING RIGID-PLASTIC AND
WORK-HARDENING PLASTICITY THEORIES

by Patrick M. Miller

The problem formulation assumes an incompressible
plastic material, a prismatic bar having at least one
cross-sectional axis of symmetry, and the stresses to be
independent of the longitudinal coordinate.

The objectives achieved in the study are threefold:

l) to better establish the torque-moment interaction
curve for the perfectly-plastic material,

2) to provide a solution to the combined bending and
torsion problem for both the work—hardening
generalized J2 deformation theory and the counter—
part generalized J2 flow theory, and

3) to compare the predictions of these work-hardening
theories for different loading parameters.

The rigid-perfectly plastic numerical analysis is
confined to a solution to an equation derived by S.
Piechnik, a non-linear second—order partial differential
equation. Using a finite-difference approach, the dis-

crete analogue of this equation is solved by the

Patrick M. Miller
Gauss-Seidel over-relaxation procedure. The results of
the rigid—perfectly plastic analysis confirm the more
limited results obtained by other investigators through a
numerical solution of the Hill-Handelman stress-function
equation. But, because of the better convergence pro—
perties of the Piechnik equation, the solution easily
gives any point on the torque—moment interaction curve.

Equations, considering the effect of work-hardening,
are developed for both the generalized J2 deformation
theory and the counterpart flow theory. With each theory,
a stress-function and a warping-function formulation is
given. This development, in all cases, leads to a system
of two simultaneous partial differential equations. The
first is the governing equation, corresponding to a com-
patibility equation with the stress function and to an
equilibrium equation with the warping function, while the
second equation represents the non-linear behavior of the
constitutive relation.

In determining the discrete analogue for the flow-
theory equations, a backward-difference quotient is
applied to derivatives taken with respect to the time-
like variable. This approach leads to a system of finite-
difference equations, whose solution is comparable to the
solutions of the discrete deformation—theory system of

equations.

Patrick M. Miller

The numerical algorithms apply the Gauss-Seidel
over-relaxation iteration to the discrete analogue of the
governing equation in a manner analogous to the corres-
ponding linear elasticity problem. Concurrently, the
Newton-Raphson iteration is applied to the discrete repre-
sentation of the constitutive relation. During a sweep
through the cross-sectional mesh, these iterations are
performed successively on each unknown discrete variable
at all points. The sweeps are continued until the
greatest change in any of the unknown discrete values of
the governing equation due to the last iteration is less
than a preassigned convergence parameter.

The results of work-hardening analysis lead to the
following conclusions: (1) In both the deformation and
flow-theory calculations, the stress-function equations
and the warping-function equations give comparable results
with approximately the same amount of computation.

(2) The results for the deformation theory agree with

those of flow theory when a constant ratio is maintained
between the curvature and unit angle of twist variables.

(3) The two work-hardening theories predict different
results for a given deformation when this ratio is allowed
to vary. (4) Because the backward-difference quotient was
used with flow-theory formulation, the numerical solution

of the flow-theory equation is of the same degree of numeri-

cal difficulty as that for the deformation-theory equations.

NUMERICAL ANALYSIS OF COMBINED BENDING AND TORSION
OF A PRISMATIC BAR USING RIGID-PLASTIC AND

WORK-HARDENING PLASTICITY THEORIES

BY

Patrick M. Miller

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Metallurgy, Mechanics and Materials Science

1966

ACKNOWLEDGMENTS

The author expresses his appreciation to
Professor L. E. Malvern for help in formulation
and solution of this problem, to Professors w. A.
Bradley and C. S. Duris for suggestions which
contributed to the solution, and to Professors
G. E. Mase and T. Triffet for serving on the

guidance committee.

ii

 

-—-. — —- --———-_— _. m *—

 

TABLE OF CONTENTS

ACKNOWLEDGMENTS O O O O O O O O O O O O O O O O O 0

LIST OF
LIST OF
LIST OF
Chapter

I.

II.

III.

 

TABLE S O O O O O O O O O 0 O O 0 O O O O O
ILLUSTRATIONS. . . . . . . . . . . . . . .

APPENDICES O O O O O O O O O O O 0 O O O 0

INTRODUCTION . . . . . . . . . . . . . . .

1.1 Review of Related Investigations. . .
1.2 Objectives and General Discussion.of
the Present Study . . . . . . . . . .

DERIVATION OF EQUATIONS. . . . . . . . . .

2.1 General Principles of Mechanics . . .

2.2 Rigid-Perfectly Plastic Material
Governed by Hencky Deformation Theory
2.2.1 Warping-Function Formulation .
2.2.2 Stress-Function Formulation. .

2.3 Rigid-Perfectly Plastic Material
Governed by Levy-Mises Theory . . .

2.3.1 Warping-Rate Function Formulati

2.3.2 Stress-Function Formulation.
2.4 Work-Hardening Material Governed by
Generalized J Deformation Law. . .
2.4.1 Warping-Function Formulation
2.4.2 Stress-Function Formulation.
2.5 Work-Hardening Material Governed by
Generalized J Flow Law . . . . . . .

mecca.

2.5.1 Warping-Rate Function Formulation. 36

2.5.2 Stress-Function Formulation. .

NUMERICAL SOLUTIONS FOR A UNIT SQUARE
CYLINDER O O O O O O O O O O O O O O C O O

3.1 Symmetry Properties of the Solution .

3.2 Discretization and Finite-Difference
Operators . . . . . . . . . . . . . .

iii

Page
0 0 ii
. . vi
. .viii
. . x
O O l
. . 1
O O 5
. . 10
O O 10
. . 15
. . 18
O O 20
. . 23
on. 25
. . 27
. . 29
. . 31
. . 33
. 34
. . 37
. . 4O
. . 4O
. . 41

Chapter

IV.

 

CONTENTS (continued)

3.3 Piechnik Equation for Rigid-Perfectly
Plastic Material. . . . . . . . . . .
3.4 Warping-Function Equations for Work—
Hardening Generalized J2 Deformation
Theory. . . . . . . . . . . . . . . .
3.5 Stress-Function Equations for Work-
Hardening Generalized J2 Deformation
Theory. 0 O O O O O O O O O O O O 0 O
3.6 Warping-Function Equations for Work-
Hardening Generalized J Flow Theory.
3.7 Stress-Function Equatiogs for Work-
Hardening Generalized J2 Flow Theory.

NUMERICAL SOLUTION FOR A CIRCULAR BAR OF
UNIT RADIUS O O O O O O O O O O O O O O O O

4.1 Transformation of Equations to Polar
Coordinates . . . . . . . . . . . . .
4.2 Polar-Coordinate Discretization and
Finite-Difference Representation. . .
4.3 Piechnik Equation for Rigid-Perfectly
Plastic Material. . . . . . . . . .
4.4 Warping-Function Equation for Work-
Hardening Generalized J2 Deformation
Theory. 0 I 0 O O O O O O O O O O O 0

RESULTS AND DISCUSSION . . . . . . . . . .

5.1 Effect of Mesh Refinement . . . . . .
5.2 Effect of Over-Relaxation Factor and
Choice of Convergence Parameter . . .
5.3 Results of the Piechnik Equation for
Rigid-Perfectly Plastic Material. . .
5.4 Generalized J2 Deformation-Theory
Solutions . . . . . . . . . . . .
5.4.1 Warping-Function Solution for
Unit Square Bar. . . . . . . .
5.4.2 Stress-Function Solution for
Unit Square Bar. 0 o o e o 0

5.4.3 Warping-Function Solution for a

Unit Circular Bar. . . . .
5.4.4 Approximation of the Elastic-

Plastic Boundary for an Elastic-

Perfectly Plastic Material . .

5.4.5 Comparisons of the Deformation-

Theory Solutions . . . . . . .

iv

Page

44

49

S4
59
68

74

74
78
80

84
88
88
90
94
98
99
100
101

101
102

CONTENTS (continued)

Chapter Page
5.5 Generalized J Flow-Theory Solutions. . . 107
5.5.1 Warpin -Function Solutions . . . . 107
5.5.2 Stress-Function Solutions. . . . . 107
5.5.3 Comparisons for Flow Theory. . . . 107

5.6 Comparison Between Flow and Deformation
Theory. . . . . . . . . . . 108

5.6.1 Conditions Where the Theories

Agree. . . . . 108

5.6.2 Examples of Load Paths Where the
Theories Do Not Agree. . . . . . . 115
5.7 Conclusions . . . . . . . . . . . . . . . 118

VI. POSSIBLE EXTENSIONS OF THE RESEARCH. . . . . . 121

6.1 Related Mechanics Problems. . . . . . . . 121
6.1.1 More General Loading Conditions. . 121
6.1.2 General Boundaries . . . . . . . . 122
6.1.3 Other Theories . . . . . . . . . . 123
6.1.4 Work-Hardening Equations Not

Related to Combined Bending and
Torsion. . . . . . . 125
6.1.5 Experimental Verification. . . . . 126

6.2 Mathematical Investigations . . . . . . 128
6.2.1 Bound on the Theta Increment . . . 128
6 2 2 Effect of Newton-Raphson

Iteration. . . . . . . . . . . . . 129
6 2 3 General Convergence. . . . . . . . 129

APPENDICES O O O O O O O 0 O O O O O O O O O O O O O O 132

BIBLIOGRAPHY O O O O O O O O O O O O O O O O O O I O O 20 1

 

Table

l.

2.

8.

10.

11.

LIST OF TABLES

Page

Effect of mesh refinement on normalized
torque and moment . . . . . . . . . . . . . . . 89

Comparison between numerical and perturbation
solutions for unit radial cross section . . . . 98

Comparison between deformation-theory solution
approximation to perfectly-plastic material
and the Piechnik solution . . . . . . . . . . . 106

Comparison of torques and moments for
different load paths. . . . . . . . . . . . . . 116

Effect of grid dimension on Piechnik equation
warping function and stresses for u = l/‘VS . . 173

Comparison between stresses computed from
Piechnik equation and Hill-Handelman equation
for p = 1m: 0 O O O O O O O O O O O O O O O O l 75

Normalized stress resultants for warping-
function deformation-theory analysis of a
square section. . . . . . . . . . . . . . . . . 176

Typical stress output for warping-function
deformation-theory with n = 2, u = 1.00, and
(9 = 3000 o o o o o e o o o o o o e o o o o o o 180

Normalized stress resultants for stress-
function deformation-theory analysis of a
square section. . . . . . . . . . . . . . . . . 181

Typical stress output for stress-function
deformation-theory with n = 2, u = 1.00, and
® '-'-" 3000 o o o e o o o o o o o o o o o o o o e 183

Normalized stress resultants for warping-

function deformation-theory analysis of a
circular section. . . . . . . . . . . . . . . . 184

vi

TABLES (continued)
Table Page

12. Normalized stress resultants for warping-
function flow-theory analysis of a square
section . . . . . . . . . . . . . . . . . . . . 186

13. Typical stress output for warping-function
flow-theory with n = 2, p = 1.00, and
® = 3000 o o e e o o o o o o o e o o o o o o o 190

14. Normalized stress resultants for stress—
function flow-theory analysis of a square
section 0 O O 0 I O O O O O O O O O O O O O O O 191

15. Typical stress output for stress-function

flow-theory with n = 2, p = 1.00, and
C9 = 3000 o o o o o o o o o o o o e o o o o o o 193

vii

Figure

l.
2.

4.

5.
6.

7.

8.
9.
10.

11.

12.

13.
14.

15.

16.

LIST OF ILLUSTRATIONS

Coordinate system . . . . . . . . . . . . . . .

Relation between unit normal vector and
direction cosines . . . . . . . . . . . . . . .

Ramberg-Osgood stress-strain curves . . . . . .

Coordinate system for unit square cross
seCtion O O O O O O O O O O O O O O O O O O O O

Lattice spacings for square mesh. . . . . . . .

Arbitrary integration paths for a given
deformation(K,6)..............

Relationship between unit normal direction
cosines and coordinate angle. . . . . . . . . .

Lattice spacings for polar-coordinate mesh. . .
Effect of relaxation factor . . . . . . . . . .

Normalized torque vs. bending moment inter-
action curve for a unit square bar. . . . . . .

Comparison of the normalized torque vs. bend-
ing moment interaction curves for the square
and circular bar. . . . . . . . . . . . . . . .

Approximation to elastic-plastic boundaries . .
Normalized second invariant J2 vs. y. . . . . .
Functions CB vs..§EL showing non-radial

stresses. . . . . . . . . . . . . . . . . . . .

Paths chosen for comparing the deformation and
flow-theory calculations at the point
(3000,3000) O O O 0 O O O O I O O O O O O O O 0

Estimation of new approximation by tangent. . .

viii

Page
10

14

29

40
41

67

77
79
91

95

97

103
105

112

116
138

Figure

17.
18.
19.
20.
21.
22.
23.
24.

25.

Basic
Basic
Basic
Basic
Basic
Basic

Basic

ILLUSTRATIONS (continued)

flow
flow
flow
flow
flow
flow

flow

chart
chart
chart
chart
chart
chart

chart

pattern
pattern
pattern
pattern
pattern
pattern

pattern

for
for
for
for
for
for

for

program
program
program
program
program
program

program

WARPI.
WARPO.
STRO .
WARPLV
STPLDT
POLARR

POLARO

Torque-moment vs. unit angle of twist for

warping-function deformation-theory solutions

of a square section . .

Torque-moment vs. unit angle of twist for

warping-function deformation—theory solutions

of a circular section .

ix

Page
140
141
142
143
144
145

146

195

199

LIST OF APPENDICES

Appendix

I. SOLUTION TO FINITE-DIFFERENCE EQUATIONS BY
GAUSS-SEIDEL OVER-RELAXATION METHOD . . . .
II. NEWTON-RAPHSON ITERATIVE METHOD FOR SOLVING
A NON-LINEAR EQUATION . . . . . . . . . . .
III. BASIC FLOW CHART PATTERNS . . . . . . . . .
IV. FORTRAN PROGRAMS. . . . . . . . . . . . . .
v. TORQUE, MOMENT, AND STRESS OUTPUTS. . . . .

VI. GRAPHICAL REPRESENTATION OF TORQUES AND

MOMENTS O O O O O O O O O O O O O O I O O O

Page

132

137
139
147
173

194

I. INTRODUCTION

In the elastic analysis of combined bending and
torsion, the solutions may be determined separately for
pure bending and pure torsion, and then superimposed to
give the solution for the combined problem. In a plastic
medium, this principle of superposition no longer applies,
and the problem must be formulated so that these two phe-
nomena and the interaction between them are considered
simultaneously. This interaction gives rise to a system
of non-linear equations. In general, closed form solu-
tions are not available for these equations; however, the
development of high speed digital computers permits the

problem to be treated by approximate numerical methods.

1.1 Review of Related Investigations

 

The problem of combined bending and torsion of a
rigid-perfectly-plastic bar has been considered by several
investigators. In 1944, Handelman [7], using a stress-
function approach, formulated the governing differential
equation, assuming a material obeying Levy-Mises flow
theory. The work of Handelman was extended to a more
general loading situation by Hill [11]. The stress—
function equation, governing the combined bending and tor-
sion of a prismatic bar, known as the Hill-Handelman

1

2
equation, is highly non-linear and, in addition, is
singular along the cross-sectional neutral axis of
bending.

Because of the non-linear nature of the equation,
early efforts to establish the interaction curve between
the applied moment and torque were confined to the appli-
cation of energy methods. Hill and Siebel [13], in 1953,
established energy techniques which permit the calculation
of upper and lower bounds for the torque and moment inter-
action curve. Using these energy methods, Steele [30]
determined the upper and lower bounds for a square cross
section. The maximum deviation between these upper and
lower bounds was found not to exceed 14 per cent. In
addition, Steele used Southwell's relaxation method to
determine two points on the interaction curve.

Employing an approach similar to that of Steele,
Imegwu [15] obtained additional points on the interaction
curve. It appears that the line singularity mentioned
above gave the most difficulty in obtaining the numerical
solutions to the Hill-Handelman equation.

In 1961, Piechnik [24] formulated the problem using
the warping function, as opposed to the stress-function
formulation of Hill and Handelman. Piechnik determined
the interaction curve in the region where the bending may
be treated as a perturbation of the problem of pure

torsion. Later in 1961, Piechnik and Zyczkowski [25]

3
solved the problem where the torsion may be considered as
a perturbation of the situation of pure bending. Using
both the perturbation solutions and the conditions of
isotropy, these authors fitted a smooth fourth order
interaction curve whose constants were evaluated so that
the location, slope, and curvature at both intercepts
agree with the two perturbation solutions. This curve was
found to lie almost exactly between the upper and lower
bounds of Hill and Siebel.

Concurrently, other researchers were considering
the problem of work-hardening from both the standpoints of
a deformation theory and of a flow theory.1 Ramberg and
Osgood [27] had introduced an empirical three-parameter
uniaxial stress-strain law including the effect of work-
hardening. This law was extended by Nadai to what is

called the generalized J deformation theory. Prager [26]

2
formulated a procedure giving the counterpart flow theory

to a given deformation theory. Using Prager's method, a

counterpart generalized J flow theory has been developed

2
by Prager and Laning [6].
Greenberg, Dorn, and Wetherell [6] applied both
generalized J2 deformation theory and the generalized J2
flow theory to the problem of a unit square bar subjected

 

1Some authors designate the deformation and flow
type theories as "total strain" and "incremental"
theories, respectively.

4
to pure torsion. In the deformation theory, their
governing equation is a non—linear second order partial
differential equation. For flow theory, the governing
equation is treated as a non-linear second order partial
differential equation in the space variables at each level
of the time-like variable. The finite difference approxi-
mations to these equations were solved by methods of
numerical analysis.

The numerical solution of the deformation-theory
governing equation tended to converge at a reasonable
rate, but with the flow-theory equation, the rate was much
slower; and, consequently, comparisons between the theories
were made only for a very limited range of the Ramberg-
Osgood parameters. In all of the cases compared, the two
theories were found to agree, even though the stress paths
were non-radial.

A point in the body is said to be subjected to a
radial stress path if the stress components at that point
maintain the same ratios during the loading, that is, if

d = Cdij’ where the Gij's are constants and C is a mono-

iJ
tonically increasing parameter. For such a radial path,
the flow-theory equations are integrable to those of the
deformation theory [12], [8]. With non-radial stress
paths the two theories are generally expected to give

different results, but the torsion solution of the square

bar seems to indicate that while the radial stress path is

5
a sufficient condition for agreement of the two theories,
it may not be a necessary condition.

Based on Drucker's assumptions [3], Budiansky [1]
proposed a criterion under which a deformation theory
would be an "acceptable competitor" to its counterpart
flow theory even under the conditions of a non-radial
stress path. Thus, on the basis of the "physical sound-
ness" of a plasticity theory, one might expect the defor—
mation and flow theories to yield comparable results even
though the stresses deviate from a radial path during the
loading. However, it should be emphasized that this cri-
terion deals only with the physical acceptability of the
deformation theory and not with the question of integra-
bility of the flow-theory equations to those of the defor-
mation theory.

2 Objectives and General

1.
DIScussion of the PresentStudy

 

The principal objectives achieved in this study are
threefold: (l) to better establish the torque-moment
interaction curve for the perfectly plastic material,

(2) to provide a solution to the combined bending and tor-

sion problem for both generalized J deformation theory

2

and generalized J flow theory, and (3) to compare these

2
solutions for different loadings.
In the establishment of the interaction curve for a

rigid-perfectly plastic material, Piechnik's warping-

6
function equation is solved numerically for a wide range
of combined bending and torsion. Solutions are obtained
for prismatic bars having either a unit square or a unit
radial cross section. Since the warping-function formula—
tion avoids the singularity of the Hill-Handelman equa-
tion, the numerical solution appears to converge through-
out the range of the load parameters. In Section 5.3,
these interaction curves are given, and it is found that
they agree with those points calculated by previous
investigators.

The application of generalized J2 deformation or
flow theory implicitly assumes the entire bar to be
plastic. Thus, upon initial loading, yielding takes
place immediately, followed by isotropic work-hardening,
so that the material properties are governed by an
expanding yield surface in a stress-representation space.
In this study no unloading is considered, i.e., the
second invariant of stress always increases.

Since governing equations including work-hardening
in the combined bending and torsion problem were not
available for either deformation or flow theory, systems
of equations were developed for each theory. The develop-
ment of these governing equations departs somewhat from
the standard approach since it was not possible to deter-
mine a single governing equation. Instead, in both

instances, the governing equations are presented as a

A
u
b

 

(D
0‘ )

I 1.
(D

a
‘1‘,
.h‘

7
system of implicit non—linear partial differential
equations for the unknown functions. These equations are
solved simultaneously by numerical methods.

In the usual processes of determining work-
hardening solutions, the flow-theory solution to a problem
is generally much more difficult to obtain than its
counterpart deformation-theory solution. The application
of the backward finite-difference representation with
respect to the time-like variable in the flow-theory
equations (Section 3.6) leads to a flow-theory system of
equations in which the difficulty in getting the numerical
solution through iteration is comparable to that for
deformation theory. The significance of this result is
further analyzed in Section 5.7.

This approach, in the case of flow theory, seems to
offer a distinct advantage for numerical solution over the
method which Greenberg, Dorn, and Wetherell [6] applied to
the case of pure torsion. Thus, numerous comparisons
between the two theories may be made without difficulty,
including cases of combined bending and torsion not pre—
viously treated.

Numerical solutionsikx‘both the generalized J2
deformation and flow theories are obtained for a prismatic
bar of a unit square cross section. With each theory, the
solutions are obtained for both a stress-function and a

‘warping-function formulation. The stress—function and

8
warping-function formulations appear to be equally
competitive, with neither offering any computational
advantage, and their results were found to agree within
the assumed numerical accuracy.

Differences appear between the flow and deformation-
theory predictions for combined bending and torsion when
the bending curvature ;< and the torsional angle of twist
per unit length 6) do not maintain a constant ratio during
the deformation. When-13 was constant throughout the
loading, the numerical solutions for the two theories
agreed. This agreement would seem to indicate that a
radial straining in the K -(9 generalized-strain space may
be a sufficient condition for integrability of the flow
theory for combined bending and torsion. The results of
Greenberg, Born, and Wetherell are a special case of this
since for pure torsion-ii = 0 throughout the loading. At
the present time a proof has not been given that the con-
dition of a constant ratio is in fact a sufficient condi—
tion for integrability, but the numerical solutions
obtained for the square bar support this hypothesis.

A numerical solution for the generalized J2 defor-
mation theory was also developed for the circular bar,
using polar coordinates. This solution is given only to
illustrate the feasibility of the method in instances

where the coordinates are chosen so as to correspond to

the cross-sectional geometry.

9

The basic equations are developed in Chapter II.
The rigid-perfectly plastic formulation essentially
repeats the work of Piechnik [24]. This formulation is
presented in order to acquaint the reader with the basic
differences between the rigid-perfectly plastic equations
and the work-hardening formulations. The numerical
algorithms for the square bar and the circular bar are
given in Chapters III and IV, respectively. Chapter V is
devoted to the presentation of the results and general
conclusions of the study. Chapter VI suggests possibi-

lities for additional investigations.

II. DERIVATION OF EQUATIONS

2.1 General Principles of Mechanics

A long prismatic bar with a simply—connected cross
section having at least one axis of symmetry is acted on
by a combined bending moment and torque at each end

(Figure 1).

 

Figure 1. Coordinate system

The bending moment M acts in the yz-plane of
symmetry about the Ox axis and the torque T acts in the
xy-plane. The z-axis is parallel to the generators of the
cylinder and is also the axis of twist.

The equilibrium equations for no body forces and
the boundary conditions are

10

11

8 OX aOXZ

XX
“affix-‘3‘“

50x 30 ad 2

 

+ + = 0
OX BY 82
.2332 + ad Z + 8022 = 0
OX ay OZ
(2.1)
an = joxx + moXy + noxz
P = /[O + m0 + no
nY XY YY YZ
Pnz = )zdxz + rnOyz + nUzz

where (an, P , Pn ) are the traction vector components

ny 2

at a boundary point with outer unit normal 3.
The relationship between the displacements and the

small strain components is

 

 

 

 

 

 

au 1 av au
SXXLSY 8w“? [ax+ BY]
_ iii _ 1 63w 63v
E:yy " ay eyz ‘ 7 [3y * 52]. (2.2)
8w 1 8w au
Szz = 3? sz = '2 [ax + 52] ,°

The total displacement of the prismatic bar is
assumed to be the sum of the displacement due to pure

bending and that due to pure torsion. For an

12
incompressible bar, assuming bending parallel to the

yz-plane, this sum is [12], for 692 << 1,

u=-%ny- eyz
v = --2i-K(222 - x2 + y2) + 6x2 (2.3)
w = KYZ + f(x:Y9@):

where I< is curvature, 69 the angle of twist per unit
length and f(x,y,(9) the warping function.

The terms u' = -6yz and v' = (9xz represent a
rigid rotation of the cross section around the z-axis
through the angle 92, provided 6’2 is small compared to
one radian. Since finite rotations are considered, these
displacement expressions are accurate only for small
values of 2. But the strain components derived from them
are independent of z and are assumed to apply along the
whole length of the bar.

The strain components for the prismatic bar are

determined from Equations (2.2) and (2.3). Thus,

 

1
8xx = - §4<y 8xy = O
8 = - le E: = 1 8f + 6x (2.4)
W ’2 yz '2 ay

l
822 = KY 5x2 = ? [—a—3—(v - QY] o

13
This strain tensor satisfies the incompressibility

assumption of the material since

1 1
38m—Exx+ syy+ ezz--'2Ky—§Ky+ Ky-O. (2.5)

Therefore, the strain tensor of (2.4) is also the strain
deviator tensor.

The stresses Oxx’ O , and Oxy are identically zero

YY
for the combined bending and torsion problem in elasti-
city. By analogy witflm the elastic problem, these stresses
are assumed to be zero for the corresponding problem in
plasticity [11]. In addition, the bar is taken to be of
sufficient length so that the stresses are independent

of z. As a result of these assumptions, the equilibrium

equations and lateral boundary condition (2.1) reduce to

 

80 O'
(2.6)
Pnz = 10x2 + mOyz. (b)

The components of the unit normal vector 3 at any
point on the lateral boundary are 3 ( I,m,0) and for the
load vector they are—P (0,0,0). Therefore, the boundary

condition (2.6b) takes the form

0x21 '1' Gyzm = O. (207)

14

From Figure 2 it is observed that the componentsj7

and m of the unit normal vector are given by [34]

j: COS nx=cai%,m=cos ny=--a—§. (2.8)

Y
city—c\
Y
ds \ \

n
d

YA

 

 

 

 

V
x

 

25>

 

Figure 2. Relation between unit normal
vector and direction cosines

Using Equations (2.8), Equation (2.7) is expressed

Oxzdy - oyzdx = O. (2.9)

For the complete specification of a mechanics
problem, the above equations must be supplemented by
appropriate constitutive equations. These physical equa-
tions define the material relationship between stress and
strain.

In the following sections the problem will be

formulated for the constitutive relations of:

15
l. Rigid-perfectly plastic material governed by Hencky
deformation theory.
2. Rigid-perfectly plastic material governed by Levy-

Mises flow theory.

3. Work-hardening material governed by generalized J2
deformation theory.

4. Work-hardening material governed by generalized J2
flow theory.

Each theory will lead to a different formulation of
the problem. In all cases the entire bar will be assumed
to be plastic. Surfaces across which the interior normal
component of stress is discontinuous are permitted, e.g.,
the neutral surface in pure bending.

The development of the rigid-perfectly plastic
equations (Sections 2.2 and 2.3) essentially repeats the
derivation given by Piechnik [24]. The work-hardening
equations (Sections 2.4 and 2.5) represent a formulation

which is original with this study.

2 2 Rigid-Perfectl Plastic Material

_,_‘

Governed'byfiHenéEy eformatiOn Theory

 

For a rigid-perfectly plastic material obeying
Hencky deformation theory the strain deviator tensor is
proportional to the stress deviator tensor [9], [12].

Thus,

(2.10)

r'

(’1

m

16
where T is an unspecified proportionality function, sij
the strain deviator tensor, and Gij the stress deviator
tensor. For this material, the function T is usually
determined through the specification of a yield condition
[12]. The Mises yield condition [19] is used in the

following development, namely

(O -O')2+(O' —O‘)2+(O -Or)2
xx yy yy zz 22 xx
2 2 2 2
+ 6(0’xz + dyz + Oxy) = 6k , (2.11)

where k is the yield stress in pure torsion. .According to

the elastic analogy assumptions, Oxx = Oyy = Oxy = 0, and

this yield condition becomes

02 + 302 + 302 = 3k2. (2.12)
22 xz yz
From the assumption that Oxx = ny = 0 it follows

that the deviatoric components of normal stress are Oix

l 2
— O§y = — 3022 and 0&2 _ 3022. Hence, the following
independent equations follow from the application of
Equation (2.10) and the elastic analogy assumptions on the

stress tensor:

17

3
022 “a?“
__ 1 5f
CYz “ET [7337 " 9"] (2'13)

1 (3f
Oxz=2T[-5'x'- 6y] '

The solution of the problem requires the determina—
tion of functions f and T such that the equilibrium equa-
tion (2.6a) and the yield condition (2.12) are satisfied
identically throughout the cross section, and the boundary
condition (2.9) is satisfied at each boundary point.

These functions may be determined by two different
approaches [24], corresponding respectively to the warping-
function approach or to the stress-function approach for
the pure torsion problem.

The two approaches are summarized here, and the

details carried out in the following sections.

(a) Warping-function method. The stress components
(2.13) are substituted into the yield condition (2.12),
giving T as a function of the warping-function f. The
parameter T found in this manner is substituted along with
Equation (2.13) into the equilibrium Equation (2.6a).
These substitutions result in a single second-order par-

tial differential equation for the unknown function f.

18

(b) Stress-function method. A stress-function Nb,

 

defined by Equations (2.14), is introduced into the system

as a new unknown:

 

Eg=§f£,3{3=- 8T (2.14)
aY'
Stresses defined in this manner in terms of the function
¢Isatisfy the equilibrium Equation (2.6a) identically.
The warping-function f is now eliminated from the shear
strain Equations (2.4) to obtain the following compatibi-

lity equation for the shear strains:

as Z asXZ

5" ay

 

= 9. (2.15)

The stress components (2.13), represented in terms
of NU, are substituted into the yield condition (2.12),
giving T as a function of ND. Substitution of (P along
with these same stress component expressions into (2.15)
produces a single second-order partial differential equa-

tion for the unknown function #1.

2.2.1 Warping-Function Formulation
As a consequence of the stress-strain relations

(2.13) the equilibrium Equation (2.6a) is

 

 

 

 

447% i- 9.)] ”‘SYWS? ex

1
(2.16)

19

The yield condition (2.12) becomes

2 2

9 K’ 2
17y

f
3y + (9x

 

2
] = 3k2o

(2.17)

 

f
"gs-91’

 

 

 

3
+25%

Equation (2.17) is solved for T as a function of f. Thus,

2

2 + f: - 2fx6y + 62y2 + fy

1
T ='§E[3I<2Y
+ 2fy9x + 92x211/2, (2.18)

where the subscript represents differentiation with
respect to that variable.

Equation (2.16), after multiplication by 2w2 yields

(fox + fyy) - cpxu?x - 9y) - «>wa + 9x) = o. (2.19)

Upon substitution for T and its respective derivatives

from (2.18), Equation (2.19) takes the following form:

2 2 2 2 2 2
3[I< y (fxx + fyy) - Kyfy — If faxy] + fyfxx + fxfyy
2 2 ,.
+ 6 y f + 2 xfyfxx - zeyfxfyy _ 2fxfyfxy + zgyfyfxy
-26xff + zezxf + 62x2f =0 (2 20)
x xy Y xy xx ' °

Equation (2.20) was first developed by Piechnik [24] and

is referred to as the Piechnik equation.

f)

u:

p;

(9'

'4-

I“

20

By introducing into (2.20) the dimensionless

coordinates
x=a§,y=a77 (2.21)
and the function
t =7?— (2.22)
a 9
and denoting
u =—é5-, (2.23)

the Piechnik equation is transformed to

2

302772(t;€ + 1:777? ) + tgg (t77 + g) (:7777 (tg

~21: (t 77)(t77 + g) ~3p277t -3p,2€77=
”'7 5 77

(2.24)

Expressed in terms of the dimensionless coordinates

g and 'n and the function t, the boundary condition (2.9)

(tg - 77mm - (t77 + flag = o. (2.25)

2.2.2 Stress-Function Formulation
As a result of the stress-strain relations (2.10)
and the Equations (2.14) defining'¢/, the shear strains

become

21

 

 

e =kT a¢is =-«T a¢j (229

Substitution of these strains into the compatibility

Equation (2.15) gives

.aw

53¢

 

 

+ gy(._5_a‘f).1@=o. (2.27)

From the yield condition (2.12), the normal stress

a
'5')?

 

 

022 is, if lpg and l#& denote partial derivatives,

022 = ifiku - 41$ - (ijl/Z, (2.28)

where choice of sign corresponds to the sign of the
coordinate y. Since symmetry conditions permit restricting
an analysis to the region where y is positive, the posi-
tive sign will be taken for the following development.

According to (2.13) the stress Gzz is

3

Equating Equations (2.28) and (2.29) and solving for T
yields

¢=V3— K

1 . (2.30)
2): 2 2 1/2
[1 - (LIX - l/Jy]

22
Substitution of T from (2.30) into Equation (2.27)
produces the following partial differential equation for

the unknown function #1:

7%.: o, (2.31)

where

 

X = 2y 2 1/2'
[1 - Aux - wy]

Through the use of energy principles, Equation
(2.31) was first developed by Handelman [7] and later
extended to more general loading situations by Hill [11].
This equation is referred to as the Hill-Handelman
equation.

From Equations (2.9) and (2.14) the boundary condi-

tion on 1P for the Hill-Handelman equation is expressed as

23¢ a
dK/J = Tidy + 7:de = O, (2.32)
or

(p = constant A (2.33)

on the boundary. However, since only derivatives of Vb
are required, the constant is taken as zero, making the

boundary condition

([1 = o. (2.34)

23
2.3

Rigid-Perfectlnylastic Material
EBverned'By Leyy-Mises Theory

For a material obeying Levy-Mises theory, the

rate-of-deformation deviator tensor is proportional to the
stress deviator tensor [12].

That is

8' .

ij "" “£1

(2.35)
where the dot represents differentiation with respect to

an appropriate time-like variable and >\ is an unspecified
proportionality function.

The kinematical relation between the velocities and

the rate of deformation tensor, which for small displace-

ments may be identified with the time derivative of the
strain tensor, is

 

 

 

 

_a ' _1[a{r+ 63']
xx ax XY-E ax W
' __a_; ' =1 av} ax}
Syy" ay syz 7[8Y+ OZ] (2.36)
. _a;, .
zz az

 

 

The velocities are determined by differentiating
the displacement (2.3), with respect to time.

For small-
displacement theory, these velocities are (for small 2)

24

- ékxy - éyz

2

<
u

--%I;(222 - x + y2) + (éxy (2.37)

Kyz + F(x,y, 9) ,

2
II

where F(x,y,(9) is the derivative of the warping function
with respect to time. For small-displacement theory the
coordinates may be considered as material coordinates
giving the initial position of the particle.

The rate of deformation tensor is assumed to be
independent of z and is determined from Equations (2.37)

through the kinematical relations (2.36). Thus,

 

. - - l . . — 0
xx ' -§;<y xy -

E =- lpéy 6.: = 1[3F+ 6.x] (2.38)
yy 7 yz 7 81/

O - O . - 1 BF 0

zz"Ky xz_2[_§)—c-9y]'

As a consequence of relation (2.35) the non-zero

components of the stress tensor are

 

3 .
XX 2X

1 BF '
yz 2X [ ay + 9x]

By a procedure analogous to Section 2.2, a warping-
rate function equation and a stress-function equation are

derived.

2.3.1 Warping-Rate Function Formulation
The stresses from (2.39) are substituted into the
equilibrium Equation (2.6a). Differentiation and multi-

plication of this equation by 2 X? result in

>\(Fxx + Fyy) - >\x(Fx — 9y) - )( (FY + 9x) = o. (2.40)

Y

The proportionality factor >\ is obtained through
substitution of (2.39) into the yield condition (2.12).

After simplification, :X is expressed as

)\ = 7%[3 RZyZ

+ F: — 29x93) + 92% + F;

+ 2pyex + 922311”. (2.41)

a“...

Gs

26

The proportionality factor' A is eliminated from

(2.40) through (2.41) to give

'2 2 ° 2 ‘ 2 ' 2 2
3[K y (Fxx + Fyy) - K yFy - K exy] + F Fxx + FxFyy
O 2 2 O O O
+ 9 y Fyy + 29xF Fxx zeypxryy 2FxFnyy + 26yFnyy
- ZéxFF + 26.sz + 9.2sz = 0 (2 42)
x xy y xy xx ’ °

With the dimensionless coordinates <§ and 77 defined in

(2.21) and with

u‘ = g, (2.43)
F
T = 75- (2044)

a

Equation (2.42) becomes

3p¢2772(T C + T7777 ) + Tég (T77 + g)2 + T7777 (T6 - 77)2

5
- 2T;77 (Tg - 7’])('I‘77 + g) — 3u‘277T - Bu‘zén = 0.
5 77 (2.45)
In this notation the boundary condition (2.9) is

expressed as

(ch - 77)c177 - (T77 + gmf = o. (2.46)

27

If p = u', then Equations (2.45) and (2.24) are
identical. This means-ii =-£:, or {§-= éé- which implies
ln K = an + constant, and consequently g- = constant
during the test. This type of loading is "proportional
straining" in terms of the generalized variables fi< and 6).
Under the condition of proportional straining, there is no
difference between the Hencky and Levy-Mises rigid-

perfectly plastic interaction curves for combined bending

and torsion. This was first noted by Piechnik [24].

2.3.2 Stress-Function Formulation
By eliminating the warping-rate function from the
kinematical relations (2.38), the following compatibility

relation is obtained:

 

8 8
35%.. %:z=9. (2.47)

With use of the stress-strain relations (2.35), and
Equation (2.14) defining ¢;, this compatibility equation

becomes

“STD—85%)] +_§_y[>\__g%] +.Q.=o. (2.48)

Paralleling the procedure used to determine T, the

following relation for A,is derived from the yield

condition:

28

 

0- 15;. é’y
I]: o (2049)
2 I72
[1 - (bx - \[Jy]
Substitution for )\ from (2.49) produces (2.50), a single
second-order partial differential equation for the

unknown #1:

 

 

 

%;[X if] 4, gy[ _%L.§] + 173—€730, (2.50)

where

 

Y
X = 2 2 1/2‘
[1 -“px ““l’lyl
The boundary condition for (2.50) is also

(p = o. (2.51)

Again, if-ég =-€%, Equation (2.50) is identical to
(2.31) and the two theories agree.

Equation (2.50) has been solved numerically by both
Steele [30] and Imegwu [15] and [16]. In his second solu-
tion, Imegwu examines Hill's version of (2.50), which
represents a more general loading condition than considered

here.

29

2.4 Work-Hardening Material Governed
Byfa CeneraIIzequ DeformatiOn Law

Ramberg and Osgood [27] introduced the following

three-parameter stress-strain curves for uniaxial tension:

0'
V5):

 

1 2n
8 = E [l + ]U, (2952)

 

 

where 8 is strain and O stress. The law fits a variety of
metals and includes work-hardening effects. These stress—
strain curves for n = l, n = 9, and n -- 03 are shown in
Figure 3. For n-«- 00 the curve approaches the elastic-

perfectly plastic law.

 

 

 

 

n
105 [— l
O 9
1.0 ~ 00
\/3 k
0.5 b
O l l l
1.0 200 300

ES

 

Figure 3. Ramberg-Osgood stress-strain curves

The uniaxial stress-strain law (2.52) has been

extended by Nadai [6] to the following generalized

30

Ramberg-Osgood J2 deformation theory:

n
2685.1 = [1 + J2]O‘ij, (2.53)

where J2 is the normalized second invariant of deviatoric

stress,
_ 1

The development is again analogous to Section 2.2.
However, the yield condition is no longer required since

the proportionality factor is now specified. Letting

q; = 1 + J2, (2.55)
Equation (2.53) becomes

ZGSij = (bah. (2.56)

The kinematic equations (2.3) and (2.4) remain valid for
the incompressible material. According to the stress-
strain relation (2.56) and Equations (2.3) and (2.4), the

normalized non-zero stresses are

31

 

 

0.21/3
Vic C13Ow
O z _ C) (at
+_$[ay+x] (2.57)
O‘xz__®[at ]
‘17-'55 “gs-Y)
where
(:>='§é29 H ='i:
9
and

2.4.1 Warping-Function Formulation
The shear stresses of (2.57) are introduced into

the equilibrium Equation (2.6a) to obtain

“53%- )1 +gY[_<§15(

The normalized second invariant of stress J2 is

 

] = 00 (2058)

 

040)
Mr?
+
x

 

_E2_[_1
8" d:

 

32

If the stresses from (2.57) are employed in (2.59), then

2 2

_®2 2 2
J2 -«——5 [3p y + (ty + x) ]. (2.60)

CD + (tx — y)

Thus, Equation (2.55) becomes
cp-1+ @211 [322+ (t - )2+ (t +302]n (2 61)
- 13?; H Y x Y y . .
Differentiation and simplification of (2.58) lead

to

CI>(txx + tyy) - cbxvcx - y) - cpyuzy + x) = o, (2.62)

where Cb is obtained from the (2n+1)th degree polynomial

representation of (2.61). Thus,

(find - c132“ - 5n = o, (2.63)
where

)2

S = ®2[3u2y2 + (tx - y + (’(:y + x)2].

The boundary condition on the function t according

to Equation (2.9) is

(tx - y)dy - (ty + x)dx = 0. (2.64)

Equations (2.62) and (2.63) are solved simulta-

neously for the unknown functions t and db.

33

2.4.2 Stress-Function Formulation
In view of the stress-strain relation (2.56) and
Equations (2.14) defining the stress function, the shear

strains are expressed as
k at!) k ; 841
8X2. = Egg—5;, Eyz = - de—é—i-o (2.65)

Since Equations (2.3) and (2.4) remain valid, the
compatibility Equation (2.15) must be satisfied. Substi-

tuting from (2.65) into (2.15) yields

a aw a eat/1 26
—é-’-‘[¢-—a-§]+_-a—i[ W]+—@-=O. (2.66)

Equations (2.14), the first Equation (2.57), and Equation

(2.59) then bring Equation (2.55) into the form

2
Cb: l + [3.2 “21,2 + \IJ: + $3] n. (2.67)

These equations can be brought into a form super-
ficially similar to those for the warping-function formu-

lation. Equation (2.66) takes the form
@(KPXX + ‘11”) + 6px (bx + cpy L/Jy + 2@ = o. ._ (2.68)

Equation (2.67) is expressed as a (2n+1)th degree poly-

nomial,

rt-

34

c133“ - c132“ — sIn = o, (2.69)

where

S = 3 ®2u2y2 + ¢2(1,b)2( + $5). (2.70)

It should be noted that Equation (2.69) is not really of
the same form as Equation (2.63), since S now contains the
unknown<i>. But in the numerical algorithms both equa-
tions are treated in the same manner.

The boundary condition on the function MD is #1: O
on the surface contour.

Equations (2.68) and (2.69) are solved simulta—
neously for the unknown functions 1p and d).

2.5 work—Hardening Material Governed
By a Eeneralized quFIow Law

 

 

According to Prager-Laning type flow theory [6],

[26], the generalized Ramberg-Osgood law takes the form,

° ' 2n+l n-l '
i — O 1
where the dot represents differentiation with respect to
an appropriate time-like variable. Since 69 is assumed to
increase monotonically with time, we may take 9 as the
time-like variable. In the following, the dot refers

specifically to differentiation with respect to 9 .

35

Letting
permits (2.71) to be expressed as

268'

ij = dij + Aoij. (2.73)

Differentiation of the displacement Equations (2.3),

with respect to 9, yields the velocity components,

- x2 + y2) + xy (2.74)

yz + f(x,y,(9).

.3;

From the kinematical relations (2.36), the corres-

ponding strain rates are

m.

 

10
xx = _'2’<y 8xy = 0
° 1 ' ° 1 a?
eyy = --7I<y eyz =-§ [ ay+ x] (2.75)
O . O 1 ap
szz=Ky 8xz=2[ax-y]

where i%=-§€% and F is the derivative of the warping
function with respect to 6 .

36

2.5.1 Warping-Rate Function Formulation
Equations (2.73) and (2.75) are combined to obtain

the following non-zero stress relations:

22 + x022 = BGKy
yz + >\°yz = G(Fy + x) (2.76)
XZ + >\ze = G(Fx -’ Y).

The shear stresses must satisfy the equilibrium

equation

 

50x2 5012
+ = o, (2.77)
a X a Y
and, on the contour surface, the boundary condition,

In order to determine an equation representing
(2.77) in terms of the warping-rate function F, the first
order differential equations (2.76) would have to be
solved explicitly for the stresses as functions of the
strain rates.

Since Equations (2.76) represent a coupled set of
non—linear equations, they will not be solved explicitly.

Rather, in Chapter III, by means of a backwgrds difference

37
scheme for Q , a discrete system of equations will be
cieveloped. This system, in the discrete sense, will be

shown to satisfy Equations (2.76) and (2.77).

2 . 5. 2 Stress-Function Formulation
By virtue of the constitutive relation (2.73) and

— C 2.14) defining (p, the shear strain rates are

. k .
(2. 79)

éyz =5; (Ll/x + )(pr).

Eliminating the warping—rate function from the
Shear strain rate expressions in (2.75) produces the

f o 1 lowing compatibility equation:

88 z agxz
_.‘L_ax -Ty .-. 1, (2.80)

After substitution of the shear strain rates from (2.79)
aInd further simplification, the compatibility equation

( 2.80) becomes
. . . 2
LZ/ch + Lpyy + >\(L/Jxx + Lpyy) + >\x 4}}: + >\yklqu' 7% = 0’

where, according to (2.72) ,

n-l

’_ 2n+l
X-TJZ J20

n.‘

«U
‘L.

MU
RH a

38

From (2.59), the normalized second invariant J2 is

 

 

 

 

 

 

 

— + o
2 V5):
or
J ,_. _Z_Z. 4,2 (pi, (2.83)

 

2 V_k2

In order to determine a system of equations
comparable to the deformation-theory stress-function
Equations (2.68) and (2.69) , J2 must be expressed expli-
Citly as a function of Ry, 1.1/x, and \lJy. In Chapter III,
this will be accomplished, in the discrete sense, through
the use of a backward difference in the 9 variable.

The governing deformation—theory equations for both
the warping function and the stress function are repre-
sented as a system of two equations in terms of two
unknown functions. This formulation for a combined bend-
ing and torsion of a Nadai work-hardening material is
believed to be new. For a numerical analysis approach, it
will be shown in Chapter III that this formulation leads
i=0 a reasonable computational algorithm.

The equations for the flow-theory formulation
Qennot be expressed as two explicit continuous equations
in a form analogous to those of the deformation theory.

But, through the use of a backward difference in 9 , it

r»

le

1
tun
on

39
will be shown that both the warping-function and the
stress-function formulations may be represented, at each
discrete level of 9 , as a system analogous to that of
the deformation theory. In addition, the computational
algorithm for these equations will be shown to be almost
equivalent to the deformation—theory algorithm.

It is generally believed that a flow-theory formu—
lation leads to a much more difficult problem than its
counterpart deformation-theory formulation [2], [6], [18].
The results of the present formulation indicate that both
theories, in the case of combined bending and torsion,
lead to problems in which the solutions are obtained

through comparable numerical computations.

III. NUMERICAL SOLUTIONS FOR

A UNIT SQUARE CYLINDER

3. 1 Symmetry Properties of the Solution
The unit square cross section is geometrically

:5)(mmetric with respect to both the x and y axes (Figure 4).

YA

 

 

 

 

 

 

Figure 4. Coordinate system for
unit square cross section

Since the bending moment acts about the x—axis, the
IThermal stresses are symmetrical with reSpect to the y-axis.
Because the material is isotropic, the normal stresses are
Einti-symmetric with respect to the x-axis. Moreover, only
‘iilne square of the normal stress enters into the calculation
(DEE the stress function and warping function, since the

40

-o u
.D

a, u.
an

LIQv

s: .

'0

at.

e

.

E

‘n
.‘I

41

Inormal stress is considered only in the second invariant

<>f deviatoric stress J2. This implies that the normal

stresses produce only symmetric effects on the governing

equations. Hence, the solutions to the differential equa-

tions have the same symmetry properties as in the case of

pure torsion .
When a square bar is subjected to pure torsion, the

stress function is symmetrical with respect to the x and y

axes, and the warping function is anti-symmetric with

respect to these axes [34].

3. 2 Digcretization and
Q nite Difference Operators

A square mesh of equally spaced horizontal and

Vertical lines is superimposed on the positive quadrant of

the cross section (Figure 5).

Y? —-h-—;

 

3‘

 

 

 

 

 

 

 

 

 

 

 

>
X

Figure 5. Lattice spacings for square mesh

42
The intersections of these lines are called mesh
points or nodal points. The values of x and y at the

nodal points are given by

(i-l)h 1: is N

>4
II

(3.1)

y (j-l)h lEjEN

where i and j are integers, h the lattice Spacing, and N
the number of nodal points in either the x or y direction.
A function g(x,y) defined at the point (x,y) has the dis-
crete analogue g(i,j) defined at the point (i,j).

For the cross section, continuous physical para-
meters are replaced by discrete functions defined only at
the mesh points. Discretization of the problem is accom-
plished through replacing the governing differential equa-
tion by its discrete analogue. Thus, the problem is
reduced to that of solving a system of algebraic equations
for the values of the discrete function.

The spatial derivatives are represented by the
following finite-difference operators:

a) first derivatives at an interior point,

gx [g(i+l,j) - g(i-l,j)]/2h

(3.2)
[g(i,j+l) - g(i,j-1)]/2h,

LO
II

(1

‘1

43

b) second derivatives at an interior point,

gxx = [g(i+l,j) - 2g(i,j) + g(i-1,j)]/h2
(3.3)
gyy = [g(i,j+l) — 2g(i,j) + g(i,j-l)]/h2,
c) mixed second derivative at an interior point,
gxy = [g(i+l,j+l) - g(i-l,j+l) + g(i-l,j—1)
- g(i+1,j-l)]/4h2, (3.4)
d) first derivatives at a boundary point,
9x = [3g(N,j) - 4g(N-l,j) + g(N-2,j)]/2h
(3.5)
gy = [3g(i,N) - 4g(i,N-l) + g(i,N-2)]/2h.

The error for each of the above spatial derivative
operators is of order h2 [29].

In the flow-theory considerations, a function
g(x,y, 6) is defined at the point (x,y,6) and has the
discrete analogue g(i,j,,[) defined at the point (i,j,j).
The derivative with respect to 9 , at the point (i,j,/(),

is represented by a backward difference quotient. Thus,
96 = [g(iaj:/() " g(iajaj'l)]/A9, (3.6)

and this approximation has an error of order A69 [29].

44
Quadrature is performed by the repeated use of the
two-dimensional form of Simpson's-% quadrature formula.
The integral of the function over the area bounded by
(i—l)h:E x E (i+l)h and (j-l)h:E y E (j+l)h is approxi-

mated by

h2
ffg(x,y)dxdy = -—g
A

+ g&+l,j-l) + 4[g(i+l,j) + g(i,j+l) + g(i-l,j) + g(i,j-l)]

g(i+l,j+l) + g(i—l,j+l) + g(i—1,j-l)

 

+ l6g(i,j) . (3-7)

 

The error of this quadrature formula is of order h4 [6].

3.3 Piechnik Equation for
Rigid—Perfectly Plastic Material

 

As was observed in Section 2.3.1, the warping—
function equations for flow and deformation theory of
rigid—perfectly plastic material are identical. Conse-
quently, a solution is provided only for the case of
deformation theory.

The Piechnik Equation (2.24) is quasi-linear in the

unknown function t(x,y) and may be expressed as

Atxx + Zthy + Ctyy + Dty + E = 0 (3.8)

where the non—linear terms are contained in the coeffi-

cients

Cu

Fh.\

45

A = 3u2y2 + (ty + x)2
B = -(t)( - y)(tx — y)2
C = 3u2y2 + (tX - y)2
D = -342y

E = -3u2xy

In view of the boundary conditions (2.9) and the
anti-symmetry property of the warping function, t must

satisfy the additional constraints,

t = -x if y = 0.5 (3.9)

t=0 if x=0 or y=0,

on the first quadrant of the unit cross section.
Following the Gauss—Seidel over-relaxation proce-

dure outlined in Appendix I, the values of tm+1(i,j) for

)th

the (m+1 sweep through the mesh are calculated by

t +1(i,j) = tm(i,j) +w[t

m (i,j) - tm(i,j)] (3.10)

m+l

where QJis the over-relaxation factor and the t

(i,j)

m+1

46

values are furnished by

tm+l(i,j) =

+ (E+fi)tm(i,j+1) + (E-fi)tm

- tm+l(i-l,j+l) - tm(i+l,j-l) + tm+l(i-l,j-l)] + E

where

>l

(III

0|

E

and the discrete values of x and y are determined by

Equations (3.1).

2(A+C)

l

+

[tm(i+l,j) - tm+l(i-l,j)]/2h
[tm(i,j+l) - tm+l(i,j-1)]/2h

2 2

2
3p y + (ty + X)

-[(tx - y)(ty + X)]/2

2 2

3p2y + (tX - Y)

-l.5u2yh

3u2xyh2,

X[tm(i+1,j) + tm+1(i-l,j)]

l(i,j—1) + §[tm(i+1,j+1)

(3.11)

At the boundaries, x = 0.5 or y = 0.5,

the respective derivatives are known from (3.9) and these

are incorporated into the above equations.

For points on

47
the center lines, x = O or y = O, the values of t(i,j)
are known.

For each sweep through the mesh, the cyclic order
is prescribed so that the rows of the unknown matrix
t(i,j) are successively displaced. Thus, for each i, the
index j runs through the range 1 s j;s N, before i is
increased through its range 155i 5 N.

The iterations are continued until

Max (i,j) - tm(i,j) <= 8, (3.12)

tm+l

for a preassigned convergence parameter 8, where the
maximum is taken over all mesh points. More specific
information on the choice of the over—relaxation factor
and the convergence parameter 8 for all algorithms is
provided in Chapter V, Section 5.2.

The stresses are calculated by the finite-
difference analogue of the stress-strain relations (2.13).
The parameter (D is obtained through the analogue of
(2.18). The difference representations of these quanti-

ties are

9 2 2

‘Pk(i,j) = 1UP (1.)) =6 1)) y + éUtx — y)2 + (ty + :02] 1/2

48
ozz(i,j)

-VQ
W-zfiffiv

Ox (i,j) 9[ t(i+1,j) - t(i—l,_L) __ y]

x g 2h , (3.13)
2¢k(i,j) '

 

 

t(i j+l) - t(i j—l)
0 (i,j) 6[ ’ ’ + x]
yz 2h
k = 2(pk(i,j) ’

where x and y are determined by (3.1). The values com-
puted according to (3.13) are normalized stresses.
The normalized stress resultants, bending moment

mn, and torque tn, are determined by the following

integrals:

mn ='fi% Aydzsz

(3.14)
t =-—l j (x0 — yo )dA
n To A yz xz

where MO is the moment for pure bending and T0 is the
torque for pure torsion. The integrals of Equations
(3.14) are approximated by repeated use of Simpson's-%
quadrature formula, Equation (3.7).

The numerical calculations for the solution to the
Piechnik equation are performed by a FORTRAN program,

WARPI. A flow chart showing the sequence of operations

49
for program WARPI is provided in Appendix III, Figure 17.
The listing of FORTRAN statements is given in Appendix IV.

3.4 Warping-Function Equations
for Work—Hardening Generalized

‘gz Deformation eory

Equation (2.62) may be expressed as a non-linear

 

Poisson equation,

txx + tyy - E(tx,ty,x,y) = 0, (3.15)

where

= ¢x(tx - y) + d>y<ty + x)
(13

E

 

and cp is obtained from the polynomial,

d32n+l - $2“ - 6n = o, (3.16)
where

2

2
s = C) [3u2y2 + (tx - y) + (ty + x)2].

Equations (3.15) and (3.16) are treated as two simulta-
neous equations for the unknown functions t(x,y) and
Cb(x,y). Equation (3.15) is a second-order partial
differential equation and the Equation (3.16) is a

(2n+l)th degree polynomial in d).

50
The solution t(x,y) to the system must satisfy the
conditions of anti—symmetry and the boundary condition

(2.9). Thus, the constraints,

t = .f = 005
x y i x
t = 0 if x = 0 or y = 0,

must be satisfied.

Each equation is represented by its corresponding
finite-difference analogue. Using the Gauss-Seidel over-
relaxation procedure (Appendix I), Equation (3.15) is
solved as a Poisson equation.

Since this technique assumes values of the unknown
function t(i,j) at each grid point, the difference esti-
mate of E(i,j) is made at each point using the current
estimates of t(i,j). Thus, in the approximate discrete
sense, E is considered as a function of the spatial
variables. This means, that to some degree, the conver-
gence of the system to the solution t(x,y) depends on the
manner in which E(i,j) converges to the function
.E(tx,ty,x,y).

At each point, after the new estimate for the
function t(i,j) has been determined, a new estimate to

db(i,j) is made by means of the Newton-Raphson method

51

(Appendix II) applied to Equation (3.16). Thus, the two
iteration processes are performed during the same sweep
through the mesh. It appears that the calculation of
these corrected values for the two functions in this
successive manner tends to accelerate the convergence of
the system.

To begin the iteration, initial solutions for
t(i,j) and Cb(i,j) are assumed. These values are cor—
rected by succeeding sweeps through the net. According

to the Gauss-Seidel procedure, the values of tm+l(i,j) for

the (m+1)th sweep are calculated by

tm+l(i,j) = tm(i,j) + Cd[tm+l(i,j) - tm(i,j)] (3.18)
where

Em+1(1,j) = [tm(1+1,j) + tm(i,j+l) + tm+l(i-1,j)

+ tm+l(i,j-1) - Em+1]/4, (3.19)
and

13m+1 = [¢m(1+1,j) - ¢m+l(i-l,j)][tm(i+l,j)

- tm+l(i—1,j) - 2hy] + [Cbm(1,j+1) - d>m+l(i,j-1)]

[tm(i,j+l) - tm+l(i,j-l) + 2hx] /4<bm(1,j). (3.20)

52
The discrete values of x and y are obtained from (3.1).
Modifications in Equations (3.18) through (3.20) are made
at the boundaries, x = 0.5 and y = 0.5, so that the con-
straints of (3.17) are satisfied. Furthermore, on the
center-lines, x = O and y = 0, the values of t(i,j) are
known.

During the same sweep, the values of Cb (i,j) are

m+1
calculated by the Newton-Raphson iteration for a (2n+1)th

degree polynomial

¢§n(i,j)[d>m(1,j) - 1] - 3;“
" _ 3
Chi“ I<i.j)[(2n+1)c1>m(i,j) - 2n]

 

Cbm+1(i,j) = cpm(i,j)

(3.21)
where Sm+1 is the difference representation of
(.92 22 2 2
S =,fi 3p y + (tx - y) + (ty + x) .
The iterations are continued until
Max tm+1(i’j) - tm(i,j)' <= 8, (3.22)

 

with a preassigned convergence parameter 8, where the
maximum is taken over all mesh points.

In general, if the necessary and sufficient condi-
tions stated in Appendix I are satisfied, then the Gauss-
Seidel method will converge for an arbitrary initial

solution to(i,j). However, as shown in Appendix II, the

53
Newton-Raphson method requires an initial estimate of
Cfio(i,j) near the exact solution. Therefore, to ensure
convergence, the process is initiated near C) equal zero.
At this starting point, with C) equal to zero, the exact
values for t(i,j) and Cb(i,j) are known. Consequently,
if C) is increased by a small amount AC), and this known
solution for (D = 0 is used as an initial estimate, the
process may converge. In fact, if AC) is small enough,
one would expect the Newton-Raphson technique to converge
rapidly.

Thus, the mathematical algorithm assumes the solu-
tion to be known at G) = 1 A69. This solution is used as
the initial estimate for the solution at (9 = (I +1)A@.
Iteration is continued until the convergence criterion
(3.22) is satisfied, giving the solution at (1 +1)A®.
Experience indicates that if AC) is kept small enough,
this procedure will provide convergence throughout the
desired range of() .

After a solution is obtained for a given(3 , the
normalized stresses are calculated by the following

relations:

 

HE! {P .4:

. ._ .fl. 1:)"

54

ozz(i,j)

_\/§®
was???

____1?__. = d§(i,j)' (3.23)

a 751.1) ®[t(i,j+l) - t(i,j-1) + x]

“x"?— = cram °

The normalized bending moment and torque for this

 

(D are calculated according to Equation (3.14) using the
two-dimensional Simpson's-% rule.

The calculations for this system of equations are
performed by a FORTRAN program, WARPO. A flow chart
depicting the sequence of operations is shown in Appen—
dix III, Figure 18. The FORTRAN listing of this program
is provided in Appendix IV.

3.5 Stress-Function Equations

for‘WErk—Hardening GeneraIizéd
'2, Deformation Theory

 

 

 

The solution to the stress-function equation for
deformation theory parallels the solution presented in
Section 3.4 for the warping-function equation. Equation
(2.68) is expressed as the non-linear Poisson equation,

we (1

yy + 3(1111)‘, wy,X,Y) = O, (3.24)

where

 

 

 

th

and q: is obtained from the (2n+1) degree polynomial,

®2n+l - c112” - 5n = o (3.25)

where
2 2 2 2 2 2
s=3®uy+CI>(\/JX+LIJY).

By virtue of Equation (2.34) and the symmetry
properties of Vb, indicated in Section 3.1, the following

conditions must be satisfied:
NU = 0 if x = 0.5 or y = 0.5

we Xe
II ll
0 0

if y=00

As in Section 3.4, Equations (3.24) and (3.25) are treated
as two simultaneous equations for the unknown functions
1p(x,y) and Cb(x,y). Equation (3.24) is a second-order
partial differential equation and (3.25) is a (2n+1)th
degree polynomial ian . Each equation is represented by
its finite-difference analogue. Consequently, the Gauss-
Seidel method (Appendix I) is employed on the analogue of

(3.24) and the Newton-Raphson technique (Appendix II) is

 

 

56
used on (3.25). This algorithm is completely analogous
to that applied in Section 3.4 to the warping-function
equations.
The iteration begins with assumed initial solutions
for Vb(i,j) and Cb(i,j). These are subsequently corrected
by successive sweeps through the mesh. The values of

me+1(i93) for the (m+1)th sweep are given by

kpmflid) = (Hung) + wt $m+1<id> - L/Jm(i,j)], (3.27)
where

flung) = [l/Jm(i+1,j) + Lilm(i,j+1) + L/Jm+l(i—1,j)

“ wmu‘id'l’ + Banal/4) (3.28)
and

Em+l = ([cpm(i+1,j) - Cbm+l(i-l,j)][ (11m(1+1,j)

LPm+1(1-1,j)] + [q3m(1,j+1) — cpm+l(i,j-1)][L[Jm(i,j+1)
- ¢m+l(i’j-l)] + 8®h2 /4c:>m(1,j). (3.29)

During the same sweep, the values of CD (i,j) are

m+l
calculated according to the Newton-Raphson iteration,

2n n
cpm (i,j)[cpm(1,j) - 1] - s“1+1

Cbrnfn-l(i,j)[(2n+l)d)m(i,j) - 2n],
(3.30)

 

<Dm+l(i’j) = Cbh(i’3) -

 

57

where Sm is the difference representation of

+1
2 2 2 2 2 2
S = 3(3 p y + (b [1px + (Pyj .
The iterations are continued until

Max 1ph+l(i,j) - 1pm(i,j)| < e, (3.31)

 

with a preassigned convergence parameter s, where the
maximum is taken over all mesh points.

The algorithm for the deformation-theory stress-
function formulation proceeds in a manner identical to
that for the warping function in Section 3.4. Thus, the
computations are initiated at C) equal to AC) and C) is
incremented until the desired stresses and stress result—
ants are obtained for a given range of C).

The normalized stresses are calculated by the
difference analogues of the first equation of (2.57) and
Equations (2.14). These difference representations are,

respectively,

qzz(i,j)

03k

366—39

Oxz(i’j) _ l(J(i,j+l) - IMini-l)
_ 2h

—_T<—_ (3.32)

o z(i,j)

((1 1, ) - Viki-1, )
..i_k_ . .. + 1 2h 1.

 

58
For the stress—function formulation, the torque is
calculated directly from the stress-function 1p(x,y).

From the last equation of (3.14), the normalized torque is

1

Introducing '¢; and ‘¢§ from (2.14) yields

tn = Ti [-ffy\/Jydxdy - ffxwxdxdy] . (3.34)

The integrals of (3.34) are integrated by parts, and since

Vb = 0 on the boundary, the expression for tn becomes

2

The normalized bending moment is calculated by

_ 1
mn ‘fi'; fydzzdxdy. (3.36)

Numerical integration of (3.35) and (3.36) is per-
formed through Simpson's-% rule according to (3.7).

The computations for the stress-function system of
equations are furnished by the FORTRAN program STRO. The
sequential flow chart for the program is given in Appen-

dix III, Figure 19, and the FORTRAN listing in Appendix IV.

 

 

59

3.6 Warping:Function Equations
for Work-Hardening Generalized
_3 ., Flow Theory

 

The warping—rate function F(x,y, 6) depends not
only on the spatial coordinates, x and y, but also on the
time-like variable 6). Thus, the difference analogue of
the function must be described by three integers i, j, 1,

where

(i—l)h

>4
ll

9: 189.

Using the backward difference quotient (3.6) for Ozz’ o

Oyz’ but not for %() and denoting

xz’
Ozz xz z

I _.__ __ o = _ v =

ozz '\/3-k, Oxz k ’ Oyz E

and C9 = «E6, the stresses at the level C») = [AC-D,

according to (2.76), are

0' (1)“ V3A®F<y + Géz(j-1)
ZZ _ [\(jr)

 

A®(Fx-y) + 6122(1 -1)

0322(1) = N1) (3.38)

 

o' (I) = A®(Fy+x) + o§z(’(-1)
yz N1)

 

60

where
m1) =1+A )(1).

The quantity )f’[) is obtained from the backwards dif-

ference representation of (2.72). Hence,

NI) = 1 + 223$ Jg‘1(1)[.12(1)- 32(1-1n, (3.39)

where J2(,€) is the normalized second invariant of stress
at the [ACE level.
For convenience, the notation,

T = A®F, (3.40)

is introduced. The shear stresses of (3.38) are substi-

tuted into the equilibrium Equation (2.77) to give

a g; — A®y + 0&2‘1‘1’
“'53? [\(1)

 

+

 

8T
+ AC)x + O' ( -l)
a [51 12 I J = 0, (3.41)

By IQCZ Y

which after differentiation and simplification becomes

2 2
A(/)[_8_§._§.;g] - AM) [33.18, + (wad-1)]

ax

 

 

aT ' aOJ'(2(/ -1)
_Ay(j)[ay+8®x+oyz(j-1)] + [\(j)[ ax

 

5092(j-l)] z 0

BY (3.42)

61
The unknown function.1N(/() is obtained from Equation
(3.39).

The boundary condition (2.9) takes the form

 

Lg; — A®y + 0322(1 -1)] dy
5T . ..
+ [51’ + A®x + O'yz(1-l)] dX — O. (3.43)

For a unit square cross section, Equation (3.43) and the

anti-symmetry properties of Section 3.1 imply

TX = A®y if X = 005
T = 0 if x = 0 or y = 0.

In Equation (3.44), the conditions

II
0
II
0
0
U]

O£z(,(-l) if x

and

N
C
II
0
0
U1

(5)./2(1-1) if y
have been used.
The three equations, (3.42), (3.39), and (3.43),

represent a discrete system in terms of the (3 variable

and a continuous system in terms of the spatial variables

62
x and y. In addition, the entire system represents an
initial value problem, requiring the additional condition
that the solution is known for (3 equal zero.

These equations now represent a system which is
discrete in the time—like variable C) and continuous in
the spatial variables x and y. At each discrete level of
C), the system may be considered as a continuous system
very similar to the continuous system given for the
deformation theory (Section 3.4). The solution to this
system at a given C) level will be somewhat equivalent to
the solution of the deformation system. Since, in
general, the solution is required for a range of C), the
two systems may be solved by similar algorithms, which
will depend on incrementing AC). This formulation,
through the backwards difference on the (3 variable,
yields a system of flow-theory equations which are essen-
tially like the counterpart deformation-theory equations,
and are solvable through numerical calculations of com-
parable difficulty.

The technique for solving this system of equations
is very similar to that used in the case of deformation
theory. However, now not only the iteration but also the
solution to the system at the ,(AC) level depends on the
solution at the previous level.

If the solution is known at the previous level,

then the continuous equations at the present level take

9%!

'(‘J

an:

U)
’\

8:9

fun:

63

the form,

Txx + Tyy + E(Tx,Ty,x,y) = 0, (3.45)

where

[\(1) 3T [\(1) air
E=-—Kx-r7-T[—5-§-A®y+0;(z(1-l)]--%T7-y[-5§

--6v (1-1) 0" ( -1)
+A®x+d§z(1-l)] + Oxz ayzj .

 

 

ax + By
The variable AJ,() is obtained from the (2n+1)th degree
polynomial
A2184”, ) - A371) - m1) = o, (3.46)
where
8(1) =39; sn‘1(j)[5(j) - [(2% )32(1—1>J, (3.47)

and the parameter 5(j ) is given by

2
5(1) = [WA®Ry+oéz(j-l)] + [—%}T-E-A®y

2

2
+ G'XZ(1—1)] + [3; + A®x + O§z(1-l)] .

 

At a given C) level, Equations (3.45) and (3.46)
are treated as two simultaneous equations in the unknown

functions T(x,y,/() and ,A(x,y,/€). Each is represented

64
by its corresponding difference analogue, and simulta-
neously (3.45) is solved as a non-linear Poisson equation
(Appendix I) and (3.46) as a (2n+1)th degree polynomial
(Appendix II).

To begin the iteration, initial solutions for
T(i,j,/€) and ,A(i,j,/() are assumed to be the final
solutions at the previous level. These values are cor-
rected by successive sweeps through the mesh. The values

Tm+1(i,j,/[) for the (m+1)th sweep are furnished by

Tm+l(i’j’/C) Tm(i’j’1) + wETm+l(i’j’/C) - Tm(i’j’1)]’

(3.48)
where
fm+l(1,j,j) = [Tm(i+1,j,1) + Tm(1,j+1,,()
+ Tm+l(1-1,j,1) + Tm+l(1,j—1,1) - Em+l]/4, (3.49)

and
Em+l = [([Am(i+1’j91) - Am+l(i-1,j,1)][Tm(i+1,j,1)
- Tm+1(i-1,j,1) - 2hA®y + 2ho)'(z(,{7—l)] + [Am(i,j+1,1)

— Am+1(1,j-1,/g)][1m(1,j+1,1) - Tm+l(i,j-l,1 )

+

211869;: + 2ho§z(1 -1)] /4Am(1,j,1)] - 2h[q;cz(i+l,j,/(-D

— c§z(1-1,j,1-1) + G§z(i,j+l,1-l) - o§z(i,j-1,1-1)l.

65
At the boundaries, the constraints of Equations
(3.44) are incorporated into the above equations.
During the same sweep, the values (\m+1(i’j”1)

are calculated by the Newton-Raphson iteration. Thus,

Am+l<i,j.1) = Am(i.j,1)

A:n(i,j,j)[/\m(i,j,j) - 1] - ((m+1
Aifi-I<i.j.1)[(2n+1)Am(1,j,j) - 2n]

 

(3.50)

where Km+ is the difference analogue of (3.47).

1
Experience indicates that a one—to-one correspon-
dence between the number of iterations for the Gauss-
Seidel and Newton-Raphson methods is not sufficient to
give general convergence. Thus, at each point an addi-
tional convergence requirement is placed on the Newton-
Raphson iteration of j\m+1(i,j, I).
Temporarily, the notation [\ (i,j,/[) is

introduced. At each point in the mesh the Newton-Raphson

m+1,n

iteration is continued on /\ (i,j,/() until

m+1,n
Am+l,n(i’j’1) " Am+l,n—l(i’j’1) < 51, (3.51)

(for a preassigned convergence parameter 81, at that
particular point. '
It is worth noting that the primary difference

between the flow-theory calculation and deformation-theory

66
calculation (Section 3.4) is this additional iteration.
From a computer standpoint, these additional calculations
represented a very small difference in the two solutions
for the combined bending and torsion problem.
After inequality (3.51) is satisfied, the (m+1)th
sweep is continued. The number of sweeps are continued

until

Max Tm+1(1,j,j) — Tm(i,j,1) <2 8, (3.52)

for a preassigned convergence parameter 8, where the
maximum is taken over all mesh points.

As in the deformation-theory solution, the compu—
tations are begun for the undeformed bar. However, in
flow theory a particular path for the forward integration
is prescribed. The iteration at any level takes the pre-
vious solution as an initial estimate and iterations are
continued until condition (3.52) is satisfied.

To better illustrate the method of solution, the
curvature I< and unit angle of twist 69 space is con-
sidered (Figure 6).

Any point (/<,(9) in this space may be reached by
any number of paths pl, p2, p3, etc. The only restriction
is that 69 be a monotonically increasing function. Speci-
fication of dK/de as a function of e , with the initial

condition, K'= 0 when 6): 0, defines the load path.

67

X

(K.9)

 

Figure 6. Arbitrary integration paths
for a given deformation (K , 9)
Experience indicates that if AC) is kept small

enough the method will converge at each level of C), in a
manner analogous to that of deformation theory. However,
in flow theory an additional reason for keeping AC) small
is that according to the difference representation, the
error in the forward integration is of order AC).

At each level, the normalized stresses are calcu—

\

lated by

68

. ‘\/‘3A®ky + Oéz(i,j,1-l)
022(1’3J) = Add.) I

 

(3.53)

T(i+l,j,,( )2;T(i-l,j,1) _ A®y + 0322(1’3’1'1)
Ari-'13)] ,)

 

O
x o
N

p.

C.
Xe

I

T(iij+lgjy)2;T(lfij-121) + A®x + O§Z(i,j,1-l)
116.131) °

 

Q
P.
U.
K
V
l

The normalized bending moments and torque are deter-
mined by substitution of (3.53) into (3.14) and carrying
out the integration by Simpson's-% rule.

The computations for this system of equations are
furnished by a FORTRAN program WARPLV. The flow chart and
FORTRAN listing of WARPLV are shown in Appendix III,

Figure 20, and Appendix IV, respectively.
3.7 Stress-Function Equations

for WEEk-HardeningGéneraIizéd
J, EIow Theory

In flow theory, the stress-function W(x,y,‘6) is a
function of the three variables x, y, and 9 . The dif-
ference analogue is represented by the integers i, j, and
’[as defined by (3.37).

Through application of the backwards difference on

6, in the same manner as in Section 3.6, Equation (2.81)

becomes

69
A(1)[prx(1)+ WWW” + [\x(j)IIJx(/)) + Aydnpyd)
+ 2463 - Bond-1) - (AMI-l) = 0. (3.54)

where .A(/() is again defined by (3.39). But the norma-

lized second stress invariant J2(/() is now given by

V3A®ky + 022(1-1) 2
J2(»() = [\(1)

Equation (3.54) is expressed as the non-linear

 

4 (113251) + (113(1).

(3.55)

Poisson equation,
(bud) + WWW) + E( \JJX,(py,x.y) = 0. (3.56)

where

a -.- [Ax<1> (w) + Ay<1>wy<1> + 249- (4.4-v

- L(Inge-1.)] //\(j).

The function 1\(/() is obtained from the (2n+1)th degree

polynomial,

A2“*1(,(7)- firm) — («1): o, (3.57)

x<1>= 39-33“- sn‘1(j)[8(j) - A2(j)32(1-1)], (3.58)

70

and the parameter 5(1 ) is provided by

s<1>=£1f33®ky + Oéz(1-l)]2

+ A2(1)[Kl/fc(1)+ 1123(1)]. (3.59)

Since the stress function (p is still defined by
(2.14), the boundary conditions on.1p given by (3.26) are
valid for any level of C). The condition of the
unstressed bar, #lidentically zero throughout the cross
section, is assumed as the initial condition at C) equal
zero.

Equations (3.56) and (3.57) are represented by
their discrete analogues and solved simultaneously by an
algorithm identical to that employed in Section 3.6.

At a given level 1A6), the values ¢m+l(i,j,j)
for the (m+1)th sweep are furnished by the Gauss—Seidel

over-relaxation (Appendix I),

\P+1(1,j,1) \lJm(1.j,1) + wt¢m+1(i.1.1> - \Pmuawl

m
(3.60)

where

¢m+1(isj)1) = [wm(i+1,j,1) + me(i,j+l,1)

+ libm+l(i-l’j’/C) + lemma-1.1) + amp/4. (3.61)

and

71

m+

E 1 = ([Am(i+l,j,1) - Am+1(i-l,j,j)][ %(i+l,j,1)
¢m+l(i-1ajs,€)] + [Am(i:j+191) " Am+1(iaj'l:1)]
[tlJm(1,j+1,1) - ¢m+l(i,j-1,1)J + 8A®h2 + 4[\/J(1+1,j,1-1)

+ ‘lJ(i,j+1,/(-l) + l/J(i-l,j,/(-l) + (Pug-1,14)

— 41,1/(1,j,/( -1) )/4Am(i.j.1)-

At the boundaries, the constraints of Equations
(3.26) are incorporated into the above equations.

During the same sweep, the values A_ (i,j,/() are

m+l

calculated by the Newton-Raphson iteration (Appendix II),

Am+l(i’j’1) = Am(iaj91)

A§n(1,j,j )[Am(1,j,1) - 11- Km+1

Am (1.j,j)[(2n+1)/\m(i,j,j) - 2n]

 

where Km is the difference analogue of (3.58).

+1

Again, an inner iteration is performed on
j\m+l(i,j,,() at each spatial point, until condition

(3.51) is satisfied. The (m+1)th sweep is then continued.

The sweeps are continued until

Max
m

 

(11‘4de) .. \Dm(i,j,j)| e 6, (3.63)

72
for a preassigned convergence parameter s, where the
maximum is taken over all mesh points.

As in Section 3.6, a load path is specified in the
curvature-twist space. This specification indicates the
dependency of the stresses and stress resultants on the
load history.

At each level, the normalized stresses are deter-

mined

‘\/—3A®I.<y + O‘éz(i,j,/€-l)

 

 

 

“éz‘iJJ’ = AflsL/{H—
“322(11311) ___ W(1,3+1,£)2% ¢(i,j-l,,() (3.64)
O§z(i,j,/€) =_‘/J(1+1,j,1)2; W(1-1,j,1).

The normalized torque and bending moment are
obtained by substitution of the above results into Equa-
tions (3.35) and (3.36). These integrals are evaluated by
Simpson's-% rule, (3.7).

The computations for this system of equations are
performed by a FORTRAN program STPLDT. The flow chart and
FORTRAN listing of STPLDT are presented in Appendix III,
Figure 21, and Appendix IV, respectively.

The numerical results for the unit square bar calcu-
lations are given in Chapter V. In order to demonstrate

the feasibility of applying this technique to other

73
coordinate systems, a circular bar of unit radius is

considered in Chapter IV, using polar coordinates.

IV. NUMERICAL SOLUTION FOR A

CIRCULAR BAR OF UNIT RADIUS

4.1 Transformation of Equa-
tions to Polar Coordinates

 

A point P, defined by the polar-coordinates (r,a),
may be transformed into rectangular coordinate (x,y)

through equations
x = r cos a and y = r sin a. (4.1)
The inverse form of relations (4.1) is

r = (x2 + y2)1/2

and a = arctan«¥. (4.2)

If u is a function of r and a, then according to
the chain rule and the inverse relations (4.2), the rec-
tangular coordinate derivatives are determined by the

following formulas:

an au sinoc au
"56 C°S “‘JF’T‘EEZ

u u cos a u
-€;§ sin a'%%5 +-——E—--€;I.

Similar expressions may be developed for the second deri-

(4.3)

vatives azu/ 8x2, 8211/ 3Y2. and 3211/ Ex ay.

74

75
Substituting relations (4.3) into the Piechnik
Equation (2.24) yields

2 3

3p r 2)2

2 2
sin a(r trr + rtr + tad) + rtrr(ta + r

2 2
- 2tr(rtra - ta)(r + ta) + rtr(rtr + tau)

3 5

- Buzr sin g(rtrsin a + tacos a) - 3u2r sin a cos a = 0.

(4.4)
From equation (2.18), the parameter (D given in

terms of r and a is

t

e2
r

2 1/2
+ r + Zta:] .

(4.5)

w = k ¢ = -g 3u2rzsin2a + t: +

 

 

The non-zero stresses (2.13) are transformed to

 

Ozz V3.6
k

= pr sin a
V316 H

 

 

0‘ 6
z _ . COS a
.—%— "2WF; trSln a + ta-—fF—— + r cos a (4.6)
0x2: 9 tcosa-tM—rsina
E 2<Dk r a r '

 

 

Equations (4.3) permit the deformation—theory
Equations (2.62) and (2.63) to be transformed to their
respective polar-coordinate representations. These equa-

tions become

76

t
1 1 Choc a __ _
Cb trr+?tr+:2t0(0( - CI)rtr- (pa 0’
(4.7)
and
C1D2n+l _ ¢2n _ Sn = O, (4.8)
where
t 2 2
S = ()2 3u2rzsin2a + t3 + -—% + r2 + Zta ] . (4.9)

 

 

Equations (4.7) and (4.8) may be treated as two simulta-

neous equations for the unknown functions t(r,a) and

@(rmn).

The stresses are determined by the polar represen-

tation of Equations (2.57). These normalized stresses

become ,

0'

 

22 = lggpr sin 0.

g.

 

 

O'
2 ‘Cl . cos a
—%—»= CD tr51n a + ta __?_— + r cos a (4.10)
O'
xz -.£2 _ sin a _ .
._E_ _ CD trcos a ta‘_—F_- r Sin a .

 

 

The normalized bending moment and torque relations

of (3.14) are transformed into their polar-coordinate

77
representation. That is,

.1

n M J[r Sln a ozsz
0 A

(4.11)
-1
t. — TP'J[(r sin a Oxz - r cos a Gyz)dA
o A
where dA is now given by rdadr.
The boundary conditions for both the Piechnik
Equation (4.4) and the work-hardening deformation Equa-
tion (4.7) are given by the transformation of the general

boundary condition (2.64). This equation may be expressed

d dx
(tx - y) 3% - (ty + x) E = 0. (4.12)

 

 

 

 

 

 

Figure 7. Relationship between unit normal
direction cosines and coordinate angle

78

From Figure 7 the direction cosines are given by

«%§ = cos a and -%§ = -sin a,

and when x and y are on the surface of the unit circular

cross section,
x = cos a and y = sin a.
Thus, (4.12) becomes

t cos a + t sin a = 0. (4.13)
x Y

However, Equation (4.13) is recognized as the partial
derivative of t with respect to r. Consequently, the
general boundary condition (2.64) reduces to the following

condition for the case of a unit circular cross section:

This boundary condition (4.14) is appropriate for
both the rigid-perfectly plastic and the deformation—
theory work-hardening formulations.

4.2 ‘Polar—Coordinate Discretization
and—Finite-Difference Representation

 

A mesh of equally spaced concentric circles and
radial lines is superimposed on the circular cross section
(Figure 8). Again the symmetry arguments of Section 3.1

apply to the solution. Thus, in determining the solution,

79
only one—quarter of the cross section needs to be

considered.

y!)

 

 

 

 

Figure 8. Lattice spacing for
polar-coordinate mesh

Since the spacing between concentric circles is h
and the angle between radial lines is Aa, the values of r

and a at each nodal point are given by

(i-l)h 1 _<. i _<_ N

I"!
II

(4.15)
(j-1)Aa l f j.§ N.

9
n

The spatial derivatives of a function g(r,a) are repre-
sented by the following finite-difference operators:

a) first derivatives,

80

= g(i+l,j) — g(i-l,j)

 

 

 

 

gr 2h
(4.16)
- g(i.j+1) - g(i,j-1)
96 ’ 2(56) ’
b) second derivatives,
_ g(i+l,j) - 29(i,1) + g(i-l,j)
grr ‘ 2
h
(4.17)
- 3
ad (Aa)2

c) second mixed derivative,

_ g(i+l 141) - g(i-l,j+l) + g(i-l j-l)-g(iHHj€D
gra ‘ ’ 4h(Aa) ’ ’
(4.18)

The integral of the function g(r,a) over the region
(i-l)h E r fi(i+l)h and (j-l)rAa fire 5 (j+1)rAa is

approximated by Simpson's-% rule for polar coordinates,

 

‘1]g(r,a)rdadr = (Aa)3r(i) [g(i+l,j+l) + g(i—l,j+l)

+ g(i+l,j—l) + g(i—l,j-l) + 4[g(i+l,j) + g(i,j+l)
+ g(i—l,j) + g(i,j-l)] + l6g(i,j)] . ' (4.19)

4.3 Piechnik Equation for
Rigid-Perfectlyplastic Material
The polar-coordinate representation of the Piechnik

Equation (4.4) is expressed as

81

At + Bt + Ct + Dt + Et + F = 0, (4.20)
rr ra aa r a

where the coefficients are given by

A = BuzrssinZa + r(t0L +r2)2

2
B = —2rtr(r * ta)

C = 3u2r35in2a + rt;

2 3

2 .
E — 22tr(r + ta) — 3p r Sin a cos a

F = ~3u2rssin a cos a.

By virtue of the anti—symmetry property of t(r,a)
and the boundary condition (4.14), the following condi-
tions must be satisfied on the first quadrant of the unit

circle:

t = 0 if a = O or a =

min

(4.21)

Equation (4.20) may be replaced by its respective

finite-difference analogue. The Gauss-Seidel over-

82
relaxation procedure, outlined in Appendix I, applied to
this system of algebraic equations, gives the numerical

solution. Thus, the tm l(i,j) for the (m+1)th sweep is

+

given by

tm+l(i,j) = tm(i,j) + (uEEm+l(1,j) - tm(i,j)], (4.22)
where

Em+1(1,j) = [(i+5)tm(1+1,j) + (A-D)tm+l(i-1,j)

+ (C+E)tm(i,j+l) + (C-E)tm+l(i,j-l) + EEtm(1+1,j+1)
-tm(i+l,j-l) + tm+l(i-l,j-l) - tm+l(i-l,j+l)]

+ E] /[2(K+E)], (4.23)

with the coefficients having been previously calculated

according to

 

 

 

 

- 2
A — (AOL)2 [3p2rssin2a+-r 2h + r2 ] ,
- 2 2h

 

 

 

r2 + tm(i,j+l) - tm+l(i,j-1)
2(Aa) ’

83

 

 

 

 

 

 

 

 

= H r Sln a + r 2h 3
5 _ hm0021.2 tm(i+l,j) - tm+l(i-1.i) 2

' 2 2h ’
E h2(Aa)[2 (tm(i+l’j) - tm+l(i-l’j))

' '—“§" 26

t (i j+l) - t (i j-l)
2 m ’ m+l ’ 2 3 .
r + 2(Aa) + 3p r Sin a cos a] ,

F =-h2(Aa)2(3u2rssin a cos a).

The iterations are continued until the condition,

Max tm+l(i,j) - tm(i,j) 5 6, (4.24)

is satisfied for a preassigned s, where the maximum is
taken over all nodal points.

The normalized stresses are calculated from the
discrete analogues of (4.6). The normalized bending
moment may then be obtained by successively applying
Simpson's-% rule (4.18) to Equations (4.11). The compu-
tations for this system of equations are performed by the
FORTRAN program POLARR. The flow chart and FORTRAN list-
ing for POLARR are shown in Appendix III, Figure 22, and

Appendix IV, respectively.

84

4.4 Warping-Function Equation
fBr Work-Hardening Generalized
En Deformation Theory
Equation (4.7) is treated as a non-linear Poisson

equation expressed in polar coordinates. That is,

 

trr +-% tr +-—% tad + E(tr,ta,r,a) = O, (4.25)
r
where
cbt -®°‘t°‘-c1>
r r r2 a
E: ,
CID

and Q) is obtained from the (2n+1)th degree polynomial,

 

2n+1 2n 2 2 2 . 2 2
q; — c1) - [C9 (3p. r Sit) on + tr
t 2 2 n
a 2 _
+ .5 + r . 26“) J - 0. (4.26)

 

Equations (4.25) and (4.26) may be considered as two
simultaneous equations for the unknown functions t(r,a)
and Cb(r,a). In addition, the function t(r,a) must

satisfy the boundary conditions,

t=0 if a=0 0r C1=

min

(4.27)

85
Each of the simultaneous equations may be repre-
sented by its respective finite-difference analogue.
Applying the Gauss—Seidel and the Newton—Raphson methods
to these difference analogues yields the following
sequence for correcting the initially assumed values for
t(i,j) and Qb(i,j). During the (m+1)th sweep, the values

tm+l(i’j) are calculated by

tm+l(i,j) = tm(i,j) + (0[Em+1(1,j) - tm(i,j)], (4.28)
where tm+l(i,j) is given by

tm+l(i,j) = (K+B)tm(i+1,j) + (Kimmie-1m

+ C[tm(i,j+l) + tm+1(i,j-l)] - Em+1 -/[2(X+E)]. (4.29)

The coefficients are determined according to

 

 

A = (AOL)2r2
E = h2
2

5 = h(Aa) r

.___§___

(4.30)

E _ h2(Aa)2r2 cbm(1+l,j) - cbm+l(i-l,3))
m+l ‘ cpm(i,j) 21.

+

 

 

(tm(1+1,j) - tm+1(i-1,j)

Cbm(i,j+l) - dbm+l(i,j-l))
2h

2(Aa)

 

 

 

tm(i,j+l) - tm+l(i,j-l) _ 1)
2(A0L)r2

 

86

Furthermore, the corrected values of d) (i,j) are

m+l

calculated during this sweep according to

 

2n . n
(1 j)[<1> (1 j) - 11- s
<13... (i,j)[(2n+l)CI>m(i,j) - Zn]

(4.31)

where the parameter S is given by the (m+1)th discrete

m+l
representation of S in Equation (4.9).

The sweeps continue until

Max tm+l(i’j) - tm(i,j) «=8, (4.32)

for a preassigned convergence parameter 8, where the
maximum is taken over all nodal points.

The normalized stresses may now be obtained from
the difference analogues of (4.10). As a result of
Equations (4.11), these stresses determine the normalized
stress resultants.

The calculations for this system of equations are
performed by the FORTRAN program POLARO. The flow chart
and FORTRAN listing in POLARO are shown in Appendix III,
Figure 23, and Appendix IV, respectively.

The polar-coordinate solution was developed only
for these two cases. However, due to the similarity
between the algorithms for the deformation and flow
theories in the case of a square bar, it is believed that

a similar polar-coordinate algorithm could be developed

87
for flow theory. The purpose here is to demonstrate the
feasibility of the technique applied to circular geometry

rather than to compare the two theories.

V. RESULTS AND DISCUSSION

5.1 Effect of Mesh Refinement

 

Steele [30] observed, in his numerical solution of
the Hill-Handelman equation, that a 9 x 9 mesh on the
entire unit square cross section provides torque and
moment values which vary by less than 0.5 per cent from
those for a 13 x 13 mesh. Thus, assuming the finite-
difference solution converges to the continuous solution
as the lattice spacing approaches zero, the normalized
torque and moment provided by the 9 x 9 mesh appear to be
valid to at least two decimal places.

The present investigation utilizes a 12 x 12 cross-
sectional mesh. As a consequence of the symmetry pro-
perties, the solution must be obtained for only a 7 x 7
set of lattice points or 49 discrete points. In order to
estimate the effect of grid size, the number of grid spac-
ings was doubled. This latter system, a 24 x 24 mesh,
requires the determination of the solution for a 13 x 13
set of lattice points or 169 discrete points. Two solu-
tions were obtained for the Piechnik rigid-perfectly
plastic equation with N = 13. A comparison between the
normalized torques and bending moments for the two mesh
sizes is shown in Table l.

88

89

Table 1. Effect of mesh refinement on
normalized torque and moment
(calculations at N x N points)

Per cent difference
2

N u tn mn Atn Amn
7 0.333 0.8294 0.6388
13 0.333 0.8350 0.6370 +0.67 —0.28
7 1.333 0.6417 0.8230
13 1.333 0.6452 0.8188 +0.54 -0.51

The results of Table 1 indicate that the actual
error inthe torque and moment would tend to cause dis-
agreement in the third decimal place. Thus, with the 49
discrete point solution, at least two place accuracy is
assumed in the normalized torques and bending moments.

The variations in stresses for the above calcula-
tion on the Piechnik equation are presented in Table 5
of Appendix V. The stress results tend to show a dif-
ference in the second decimal place.

Greenberg, Dorn, and Wetherell [6], in their
numerical solution to the plastic torsion problem, used
an 8 x 8 mesh on one quadrant of the cross section. With
this mesh size, a solution is required for 81 discrete
points, and accuracy of the calculated torques was assumed
to be valid to three decimal places.

Using N =3B, u = 1.00, and n = 12, computer runs
were performed for both the deformation-theory and the

flow-theory warping-function equations. With each run,

90
the difference between the normalized torques and bending
moments for N = 13 and those for N = 7 was of the same
order as that shown for the Piechnik equation in Table 1.

In view of the large number of calculations
required for the combined bending and torsion problem,
and since the major emphasis of this study is confined to
the torque and moment relationships, the remaining calcu-
lations use a discrete system of 49 points on the first
quadrant of the cross section. The comparisons between
calculated stresses and stress resultants according to
different theories assume the numerical accuracy of the
order presented in the above discussion.

5.2 Effect of Over-Relaxation Factor
'EHH_CFSIEE_3f—50nvergenc6*?arameter

In the case of linear second—order equations,
formulas are available for selecting an optimum relaxa-
tion factor [32]. Generally, these formulas are dependent
on the mesh size h and the geometric shape of the cross
section. Apparently, formulas of this nature are not
available for the non-linear cases.

Since most of the computer runs were for a system
of 49 discrete points, various relaxation factors were
considered in an experimental analysis. Data, showing the
effect of the over-relaxation factor a), was determined
for the Piechnik equation, the warping—function equation,
and the warping-rate function equation. This information

is presented in Figures 9a, 9b, and 9c.

91

140—

120—

100—

I 80—

60-

40—

20~

 

1 1 1 1 I L 1 l 1 1
1.0 1.1 1.2 1.3 1.4 1.5 L6 1.7 L8 1.9 2.0

a)

 

Figure 9a. Effect of relaxation factor (0 on
the number of iterations I for the Piechnik
equation with N = 7 and u = 1.00

l40~
120
100
I 80
60

40

 

 

0 L l J l L I 1 1 l .4
1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0

0)

Figure 9b. Effect of relaxation factor Q) on
the number of iterations I for the warping-
function deformation-theory equations with
N = 7, u = 1.00, and n = 2

92

140 .

120

100

I 80

60

40

 

20

 

0 1 1 1 1 1 1 1 1 1 1
1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0

U)

 

Figure 9c. Effect of relaxation factor cu on
the number of iterations I for the warping-
function flow-theory equations with
N 3 7, p. = 1000, and n = 2

In the case of the Piechnik equation, Figure 9a,
the system failed to converge after 250 iterations, for
a): 1.9. With the other two equations, the system con-
verged for the entire range of a). In all cases, the
results indicate that the rate of convergence is quite
dependent on the values of the relaxation factor. Values
of (0 from 1.5 to 1.6 appear to give reasonable conver-
gence rates for the values of u, n, and C) tested.

The value of 1.55 for (U was used in all calcula-

tions except for the experimentation just described. This

93
implicitly assumes that the optimum over—relaxation
factor is independent of u, n, and C). As Figures 9b and
9c indicate, this probably is not exactly the case, since
the optimum (0 changes when n and u are held constant and
C)is allowed to vary. However, a run for n = six and
u = zero for the case of deformation theory gave the
optimum a) as approximately 1.45. In addition, with the
optimum over-relaxation factor, these deformation and flow
equations tend to converge faster than the Piechnik
equation.

For the case where N = 13, the number of discrete
points is 169, and the optimum over-relaxation factor is
1.80. Using the optimum over—relaxation factors, 1.55 for
N = 7 and 1.80 for N = 13, the following comparisons were
observed with the warping-function deformation—theory
equations for u = 1.00 and n = 12 with 0 5 C9 5 4.00:

1) With the increase in grid lines, the average number

of iterations increased from 23 to 38.

2) The execution time increased from 0.82 minutes to

4.27 minutes for the respectives runs.

Because of the large number of computer runs, it
was felt that the most effective use of the computer time
would be realized if the convergence parameter required
the computed variable to have only four significant non—
zero digits after the decimal place. Since only first

derivatives of this function are needed, the accuracy

94
should be sufficient to allow the making of reasonable
comparisons. The convergence parameter is given in each
program listing (Appendix IV), under the variable name El.
In each program E1 is set equal to the desired constant.
5.3 Results of the Piechnik Equation
for Rigid—Perfectly PIastib Material

The Piechnik equation was solved numerically for
the unit square and unit circular cross sections. The most
important aspect of each solution is the interaction curve
between the normalized torque and bending moment.

The interaction curve for the unit square cross
section is shown in Figure 10. These results are compared
with those obtained by Steele [30] and Imegwu [15] through
a numerical solution of the Hill-Handelman equation. The
numerical results are essentially equal, but in the pre-
sent study the Piechnik equation has been solved through-
out a greater range of p. This is possible since the
Piechnik equation does not have the line singularity pre-
sent in the Hill-Handelman equation. In addition,

Figure 10 contains the upper and lower bounds for a unit
square section, obtained by Steele [30] by means of energy
principles developed by Hill and Siebel [13].

Because of the line singularity in the Hill-
Handelman equation, nodal points cannot be placed on the
x-axis. Steele [30] and Imegwu [15], [16] begin the grid

lines at a distance of-% from both the x and y axes. If

95

    
 

 

 

upper bound
lower bound
t
P— {}——D—Piechnik equation
—- -<}—-C6Hill-Handelman equation 2‘
‘3.
1“ 1
—- \
| 1
1
F.
1, l 1 l J l I I l -
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8. 0.9 1.0
m ,

n

Figure 10. Normalized torque vs. bending moment
interaction curve for a unit square bar

96
the Piechnik equation is solved with N = 10, then the
stresses at selected spatial points may be compared with
those obtained by Steele at these points. Steele computed

the stresses for p =-—l and p =-—3. The Piechnik equation
3

was solved with N = 10 for thesevgame values of u. A com-
parison between the stresses given by Steele and those of
the Piechnik solution for p =-—£ are given in Table 6 of
Appendix V. Because of the fingr grid, the results of the
Piechnik equation should be more accurate. However,
general agreement between the stresses for the two solu-
tions is observed.

A comparison between the interaction curve for the
unit circular bar and unit square bar is illustrated in
Figure 11. These results seem to indicate that the inter-
action curve is somewhat insensitive to the geometric
cross section. Prager and Onat [20] obtained a similar
result for the case of combined plastic bending and an
axial load. In his more limited analysis, Imegwu [15]
came to the conclusion that such an insensitivity appears
to exist also for the case of combined bending and torsion.

Piechnik [24] solved his equation by assuming the
bending to act as a perturbation on the pure plastic
torsion. Due to mathematical difficulties, this solution
was confined to a unit circular cross section. A com-
parison between the numerical and perturbation solutions

is presented in Table 2. Piechnik states that his results

97

 

 

 

0.4 - ‘D—dD—square bar

0.3 - ~{F—0-circular bar

0.2 t

0.1 r

o 1 l 1 I l) L, l l L
0.1 0.2 (L3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
m
n

Figure 11. Comparison of the normalized torque vs.
bending moment interaction curves for the
square and circular bar

98

should be reasonably accurate for p less than 0.25.

The

values shown in Table 2 indicate a favorable comparison

within this range of u.

Table 2.

u

0.10
0.15
0.20
0.25
0.30

radial cross section

Numerical

tn mn
0.9929 0.1329
0.9845 0.1945
0.9735 0.2518
0.9603 0.3050
0.9454 0.3542

Comparison between numerical
and perturbation solutions for unit

Perturbation
tn mn
0.9927 0.1328
0.9846 0.1951
0.9747 0.2529
0.9647 0.3240
0.9562 0.4107

The results for both of these numerical solutions

of the Piechnik equation agree with those determined

through energy, perturbation, and numerical methods, by

other researchers.

Thus, even though no mathematical

justification is available to ensure convergence of the

system to the correct result, the above comparisons

strongly suggest the method indeed converges to a satis-

factory solution.

5.4 Generalized J

 

Deformation—Theory,Solutions

The combined bending and torsion problem for the

generalized J

deformation theory was solved numerically

by both a warping-function and a stress-function approach.

For the case of the unit square bar, the two formulations

provide essentially the same result.

99

This generalized J deformation law was also used

2
in the solution of the unit circular bar in polar
coordinates. The polar-coordinate solution shows that the
mathematical algorithm appears to be a tractable means to
solve the problem when the coordinate system corresponds
to the cross-sectional geometry.

5.4.1 Warping-Function

Solution for Unit Square Bar

The warping-function deformation-theory equations
were solved for the values of u and n as shown in Tables 7
of Appendix V. The resulting normalized torques and bend-
ing moments are also given in these tables. For each
particular combination of H and n a given range for C) is
specified. The range of C) and the given AC) may be
obtained from inspection of the data in each table.

To illustrate the general effect of the work-
hardening on the torques and moments, the normalized
torques and bending moments vs. the unit angle of twist,
for selected values of the Ramberg-Osgood exponent n, are
shown in Figures 24 of Appendix VI.

Due to the large number of calculations, it is not
convenient to present the entire stress output. However,
Table 8 of Appendix V shows the typical stress output
obtained by the warping-function formulation for some

values of u, n, and® .

100
5.4.2 Stress-Function
Solution for Unit Square Bar

The solution to the stress-function equations
serves as a check on the warping-function solution and was
developed so that a comparison could be made with the
above results. The equations were solved for the combina-
tions of u and n as presented in Tables 9 of Appendix V,
and the normalized stress resultants obtained are also
given in these tables.

Calculations of the stresses corresponding to the
same values of u, n, and C) as given in Table 8 (warping-
function results of Section 5.4.1) are presented in
Table 10 of Appendix V. A comparison between the two
stress calculations may be obtained by examining these two
stress tables.

Since the results are only assumed to be accurate
to two decimal places, the stress resultants for the
warping and stress-function formulation are essentially
equal. The computational algorithms are also essentially
the same. Thus, it appears neither method offers an
advantage, although, from the standpoint of error, the
stress-function formulation would be expected to be
better, since the values of the unknown function are pres-
cribed on the boundary, while in the case of the warping
function, the normal derivatives are prescribed at the

boundary.

101

The warping-function approach has the advantage of
giving directly the axial displacements, which might be of
physical interest. For this reason, and because it repre-
sents a more novel approach in plasticity, the major
effort in this study was made with the warping-function
approach.
5.4.3 Warping-Function Solution
for a Unit Circular Bar

The polar—coordinate deformation-theory warping-
function equations for the circular bar were solved for
the values of u and n as shown in Tables 11 of Appendix V.
The normalized stress resultants are presented graphically
in Figures 25 of Appendix VI.
5.4.4 Approximation of the Elastic-
Plastic Boundary for an Elastic-
Perfectly Plastic Material

The generalized J work-hardening deformation-

2
theory approach provides a method of approximating the
elastic-plastic boundary for an elastic-perfectly plastic
material. As was pointed out in Section 2.4, as n-a-oo
the uniaxial Ramberg-Osgood stress-strain curve approaches
the elastic-perfectly plastic curve represented by the
curve consisting of two straight segments in Figure 3.
When the tensor version, Equation (2.53), of the Ramberg-

Osgood law is assumed to apply over the whole cross sec-

tion in a deformation-theory analysis with no unloading,

102
it is reasonable to expect that the solution for n large
will approximate the elastic-perfectly plastic solution.
The region where J2 approaches the constant unity approxi-
mates the perfectly-plastic region,axu1the boundary of
this region is taken as the approximation to the elastic-
plastic boundary.

On the basis of letting n become large, the
elastic-plastic boundary was approximated for the square
section for u = 0 and for p = l, and for the circular
section with u = 1. These boundaries for various values
of G) are shown in Figures 12.

By the sand hill-soap film analogy, the boundary is
known for the case of pure torsion. Figure 12a has the
appearance of the elastic-plastic boundary predicted by
this familiar analogy. The boundary for the combined
bending and torsion problem is not known. However, an
indication that the solution is approximating a perfectly—
plastic material is given since the second stress deviator
invariant, J2, approaches unity. Since J2 approaches its
greatest values along the y-axis, J2 is plotted as a
function of distance along the y-axis in Figures 13.
5.4.5. Comparisons of the
Deformation—Theory Solutions

A general comparison between the stresses and

stress resultants calculated by the warping function

(Tables 7 and 8) and by the stress function (Tables 9

103

 

C9 = 6.00
@ = 4.00
+ _@ = 3.00

 

 

//\\

Figure 12a. Approximation to elastic-plastic
boundary with u = 0.00 and n = 12

 

 

 

 

 

1\\ C3 = 2.50
Mh® = 3.00
—@ = 4.00

_.|_

 

 

\\l_,/"”'* 7““‘-.__.//
// \

\

 

 

 

 

Figure 12b. Approximation to elastic-plastic
boundary with u = 1.00 and n = 12

104

ll
l-‘
C
U1
C

II
N
O
O
O

 

 

 

 

 

 

II
p
o
O
O

 

 

Figure 12c. Approximation to elastic-plastic boundary
for circular section with p - 1.00 and n - 15

and 10) indicate that they agree within the limits of the
assumed numerical accuracy. In addition, for the case of
pure torsion, each method confirms the results presented
by Greenberg, Dorn, and Wetherell [6]. As a final check,
if n becomes large, then the moments and torques approach
those given by the solution to the Piechnik equation.
Typical comparisons for both the unit square and circular
bars are shown in Table 3. These arguments offer physical
assurance that the results are indeed converging to

reasonable values.

N

mHXMIm so a .m> h unmeum>ca vacuum uwueameuoz .mH wusoeh

0.6 .o.m .o.m .mé u© 6.6 .o.m .m.~ .o.m n ® 6.6 .06 .06 6.... n @
me n 6 .oo.H u 1 NH 1 6 .oo.H u 1 «H n 6 .oo.o n 1
coewuwm Hmasuufiu flue soauuom ouwsqm Abe cowuumm mumsqm Ame
:5 m h a
nu
1. oo.a om.o

 

 

 

 

 

 

 

106

Table 3. Comparison between deformation-theory
solution approximation to perfectly-plastic
material and the Piechnik solution
Deformation theory Piechnik equation

Square Bar

n = 12, p. = 1000, Q = 4.00 p = 1000
tn mn tn mn
0.6834 0.8093 0.6874 0.7885
Circular Bar
n a 15, p. = 1000, G) =3 4000 H = 1000
tn mn tn mn

0.7119 0.7855 0.6936 0.7488

In the case of the unit circular bar, expressed in
polar coordinates, the warping function offers an advantage
over the stress-function formulation. That is, the govern-
ing equations in both cases are singular at the origin.
However, with the warping function, the singularity is
removable since the warping is known to be zero at the
origin. Thus, the warping-function approach provides a
numerically tractable way to approach the problem through
the use of polar coordinates.

The flow-theory equations are solved in Section 5.5
for the same parameters as the deformation-theory equa-
tions. As will be shown, the two theories appear to give
the same results under the condition of proportional

straining in terms of the generalized variables K and 9 .

107
5.5 Generalized J
mgns
The generalized J2 flow-theory analysis was also
formulated for both the warping and stress functions.
Both methods give essentially the same results for the

stresses and the normalized torques and bending moments.

5.5.1 Warping-Function Solutions

The warping-function system of equations was solved
for various values of n and u; the resulting normalized
stress resultants are given in Tables 12 of Appendix V.
Typical values for calculated stresses are presented in

Table 13 of Appendix V.

5.5.2 Stress-Function Solutions

The stress-function system of equations was solved
for the various values of u and n; the resulting normalized
torques and bending moments are given in Tables 14 of
Appendix V. Table 15 of Appendix V gives calculated
values for the stresses which compare with those for the

warping-rate function calculations given in Table 13.

5.5.3 Comparisons for Flow Theory

A comparison between the torques and moments for
the warping function (Tables 12) and the stress function
(Tables 14) shows that the two formulations provide

results which agree within the assumed numerical accuracy.

108
An examination of Tables 13 and 15 shows a similar
agreement between the stresses. Thus, both methods pro—
vide essentially the same results.
From a computational point of view, neither method
offers an advantage, since the mathematical algorithms
are practically equivalent.

5.6 Comparison Between Flow
and DeformatiOn Theory

 

 

5.6.1 Conditions Where
the Theories Agree

Piechnik [24] showed that for a rigid-perfectly
plastic material, the condition-é; =-€% produces the same
governing equation for flow theory as for deformation
theory. The solution to the Piechnik equation is given by
the interaction curve between torque and moment. In the
case of deformation theory, this curve results from 0 f-é;
5300. Similarly, in the case of flow theory, it results
from 0:5 Jg-EEOO. Thus, the interaction curves for both
theories are identical.

In the work-hardening materials, the solution to
the system of equations is no longer an interaction curve
but rather the relationships of the moments and torques to
the deformation parameters K7 and 9 . These relationships
depend on the Ramberg-Osgood constant n. A comparison

between the results presented for deformation theory

(Section 5.4) and flow theory (Section 5.5) indicates that

109
the two theories appear to agree under the condition

K' I(

9 It should be noted that the solutions to the flow
theory represent a forward integration in the C) variable.
This is accomplished through a backwards difference
quotient in AC). The error in the system is of order AC)
and it may tend to accumulate with increasing C). Thus,
the two theories would not be expected to agree as well
for C) large as for C) small.

In all cases tested, with N = 7, the calculations
for deformation theory and for flow theory agree to two
decimal places when C) is small. As C) assumes its maxi-
mum value in the specified range, there is a little less
agreement. This difference is believed to be due to an
accumulation in the error of the backwards difference
quotient.

To investigate these numerical differences,
computer runs were made with N = 13, u = 1.00, and n = 12
for both the deformation and the flow-theory warping-
function equations. In addition, with the flow theory the
increment in C) was AC) = 0.0625. The results of these
calculations with the two theories agree to three places
initially and gradually decrease to two places as C) takes
on its maximum value. This indicates the theories produce
the same results if«é§ = constant and that any difference

in results can be attributed to numerical error rather

110
than to actual differences in the predictions of the
theories.

It is generally believed that the solutions for a
deformation theory and its counterpart flow theory will
agree only under the conditions of radial loading in the
stress space. Greenberg, Dorn, and Wetherell [6] showed
the two theories provided the same numerical result for
pure torsion even though the stresses were non-radial.
The present study shows that under the loading condition
«f? = constant, where f< and (9 are generalized strain
parameters, the deformation and flow theories provide
essentially the same numerical results. In this more
general situation, the loading is such that the stresses
follow a non—radial load path.

As was remarked in Section 1.1, a radial stress
= Co’

13 ii
is a function which increases monotonically with the time-

load path means a where dij are constants and C

like variable. Thus, for the combined bending and tor-

sion, a radial stress load path implies

0' O O
22 XZ Z
22 xz yz

O O 0
where 022’ Oxz’ Oyz may be taken as the respective stresses

when C1) equal AG).2

 

2This choice of the value of (:>at which the
stresses Oij are taken is arbitrary.

111
To better illustrate the non—radial nature of the

combined bending and torsion problem, consider

37": c2 (5.2)

where the CB are functions of (9. A radial loading is
indicated if C1 = C2 = C3. If C1’ C2, and C3 are plotted
as functions of the time-like variable 6 , then the dif-
ference between the curves will tend to indicate the degree
of non-radial stress loading.

Figures 14 show the functions CB plotted as func-
tions of (9 for u = 0 and for u = l and for two choices of
the Ramberg-Osgood exponent, n = 2 and n = 12. The
stresses are taken at the spatial point where i = 4 and
j = 5. The results for u = 0 confirm the results reported
by Greenberg, Dorn, and Wetherell.

An analysis of these graphs shows that the stress
paths tend to be more non-radial when u is smaller and n
larger. However, in the torque, moment, and stress

results, it was found that the two theories agree under

the condition-é; = constant, and there appears to be no

 

O l l l

1 (3?
l 2 3 4

Figure 14a. Functions CB vs.-§éL showing non-radial
stresses for u = 1.00 and n = 2 at the

 

discrete point i = 4, j = 5

 

 

Figure 14b. Functions CB vs.-§éL showing non-radial
stresses for u = 0.00 and n = 2 at the
discrete point i = 4, j = 5

113

 

 

 

Figure 14c. Functions CB vs. 9&1 showing non-radial
stresses for u = 1.00 and n = 12 at the
discrete point i = 4, j = 5

 

 

6

Figure 14d. Functions CB VS"§E— showing non-radial
stresses for p = 0.00 and n = 12 at the
discrete point i = 4, j = 5

114
dependence on the values of n and p, insofar as agreement
between the two theories is concerned. The agreement
between the numerical results from the two theories seems
to suggest that, although the condition of a radial stress
path is a sufficient condition for integrability of the
flow-theory equations to those of deformation theory, it
may not be a necessary condition. Otherwise, with non—
radial stress paths, one would expect more disagreement
between the two theories than that observed with this
study.

Budiansky [1] has proposed a condition under which
the deformation theory is an "acceptable competitor" to
the flow theory even under the condition of non-radial
stress paths. For the generalized J2 deformation theory,

this condition is [6]

a1< tan-lj/Zn-l (5.3)

where

1 013011
1 1 I72 .1 '. I/Z'
[Uklokl] [opqopq]

(5.4)

 

a = cos-

Budiansky and Mangasarian [2], [18], for the case
of an infinite plate with a circular hole, showed that

generalized J deformation theory and generalized J2 flow

2
theory give comparable results even if the stress paths
are non-radial. However, in all cases the Budiansky

criterion was satisfied.

115

Calculations with the combined bending and torsion
problem indicate that under the condition ۤ.= constant,
the Budiansky criterion is satisfied. Thus, on the basis
of Budiansky's physical arguments, perhaps one could
expect the results of the two theories to be comparable.
Nevertheless, this condition does not answer the question
of integrability posed earlier in this section.

In Section 5.6.2 examples will be given where the
theories do not agree. This disagreement between the
theories will be shown to depend on how much the load

variable-$51 is allowed to vary during the loading.

5.6.2 Examples of Load Paths Where
the Theories Do Not Agree

To further compare the results of the deformation—
theory and the flow-theory calculations, the parameter-$5;
was allowed to vary during the loading. As was shown in
Chapter III, Section 3.6, a point in the generalized K—Q
space may be reached by any number of paths. For the
paths A, B, C, and D shown in Figure 15, the values calcu-
lated at the point (3.00,3.00) for the normalized torque
and moment are shown in Table 4. These calculations were
made with n = 2 and with n = 12. This table shows the per
cent difference between the normalized torques and moments
and those for load path C. In this case, path C repre-

sents the loading where é; is constant. The values for

the deformation-theory calculations are also included in

3.00—

116

(3.00,3.00)

 

 

 

2.00-
GK?
._E_
1.00
0.00 ’ ' ad
1.00 2.00 3.00
Figure 15. Paths chosen for comparing the deformation

and flow-theory calculations at the point (3.00,3.00)

Table 4. Comparison of torques and moments for
different load paths (see Figure 15)

Per Cent Dif—
ference from C

Path tn mn torque moment
Ramberg-Osgood exponent n = 2
B 00711 0.743 +11009 " 6077
C 0.640 0.797
C (deformation) 0.642 0.807 + 0.3I + I.38
D 0.593 0.836 - 7.34 + 4.52
Ramberg-Osgood exponent n = 12
C 0.659 0.781
C (deformation) 0.656 0.789 - 0.45 + 1.11
D 0.625 0.807 - 5.16 + 3.33

117
this table; deformation theory predicts the same results
for all paths.

For the paths tested, the differences are not
extremely large; however, they are of such magnitude that
they cannot be attributed to numerical error. Thus, under
the condition of non-radial loading in terms of the
generalized strain parameters l< and 6 , a difference in
the predictions of the flow and deformation theories is
observed. This difference becomes greater as the path in
the F<- 6 space becomes more non-radial.

The Budiansky criterion, Equations (5.3), was
checked with deformation-theory calculations for paths A
and B with the Ramberg-Osgood exponent n = 2. For n = 2,
the maximum permissible value for the angle a is 65.9“.
With path A the maximum calculated a was 76.2“ while with
path B the maximum calculated a was 63.9“. Thus path A
violates the Budiansky criterion while path B does not
violate this condition. These deformation-theory calcu-
1ations deviated by as much as 14.8 per cent of the flow-
theory values with path A while the maximum deviation was
about 10.0 per cent of the flow-theory values with path 8.3
Although no definite conclusion may be drawn from this

particular analysis, it would seem that if the Budiansky

 

3When calculated as a percentage of the deformation—

theory values, the percentages for paths A and B were as
high as 17.7 per cent and 11.1 per cent, respectively.

118
criterion is violated with a deformation theory then the
results certainly should be questioned.

The development of the Budiansky criterion is
predicated on the assumption that one wishes to set a
limit on when the deformation theory may be considered as
an "acceptable competitor" to a counterpart flow theory
[1]. This implicitly assumes that one would prefer the
flow-theory solution, but that generally it is easier to
obtain the deformation-theory solution [2], [18]. The
results of the present study indicate that the flow-theory
analysis is comparable to the numerical calculations
required for the deformation theory. If the flow theory
is the preferable theory, then it may be more reasonable

to proceed with a flow-theory analysis.

5.7 Conclusions

The following general conclusions relating to the
combined bending and torsion problem are made as a result
of this investigation:

1) The numerical solution of the Piechnik equation
confirms the results of other investigators. But,
from a computational point of view, the Piechnik
equation has better convergence properties, and,
consequently, all points on the interaction curve
may now be calculated.

2) In both the flow-theory and the deformation-theory

generalized J2 work-hardening material calculations,

3)

119
the stress-function equations and the warping—
function equations give comparable results with
approximately the same amount of numerical compu-
tation. Thus, neither method appears to offer any
computational advantage. However, with more com-
plicated geometric boundaries, the stress-function
approach would be preferable.
The two work-hardening theories, deformation and
flow, appear to agree under loading conditions
where a constant ratio is maintained between the
generalized strain parameters K and 9 even though
the stress load path is non-radial. This strongly
suggests that a radial stress load path may be only
a sufficient condition for integrating the flow-
theory equationstx>those of deformation theory, but

not a necessary one.

4) As the strain path is made non—radial in terms of

5)

the generalized variables K and 9 , differences
are observed in the results of the two theories.
It is customarily believed that the flow theory
should be the preferable theory [12], [17]. Thus,
for this type of non-radial straining, the flow-
theory results should be preferred to those of the
deformation theory.

The numerical analysis calculations for the flow-

theory solution are comparable to those for the

120

deformation theory. In all solutions, both flow

and deformation, the times required for convergence

at any given (9 level are comparable. The flow—
theory calculations take slightly longer due to the
additional iteration mentioned in Section 3.6, but
this iteration is not significant in terms of the
total time. The more important consideration
between the theories is the size of the increment

AC). In flow-theory computations, this must be

kept small because of the error, and, consequently,

more calculations are required for a given range

of C).

It is generally believed that even with numerical
analysis the flow-theory solutions are much more difficult
to obtain than the deformation-theory solutions. With the
iteration procedures of the present study, this was not
the case for combined bending and torsion, and it is
possible that similar procedures could be used in other
problems to make flow-theory competitive with deformation

theory from the standpoint of computational difficulty.

VI. POSSIBLE EXTENSIONS OF THE RESEARCH

It may be possible to extend the research along
two primary avenues: (l) to apply the technique to other
related problems in mechanics, and (2) to establish a
mathematical criterion which would ensure the convergence
of the iterative schemes.

6.1 Related Mechanics Problems
6.1.1 More General
Loading Conditions

The equations should be extended to a more general
loading situation. Hill [11], [12] extended Handelman's
work [7] to the general case of combined torsion, bending
about two orthogonal axes of symmetry and uniaxial exten-
sion. In a similar manner, the work-hardening formula-
tions could be generalized to include these loadings.

In this generalization, the displacement Equations

(2.3) would be replaced by the following displacements:

u = -%K1xy + %K2(y2 - x2 - 222) -%6x - eyz
v = %K1(x2 - y2 - 222) - 3211(ny - )2” by + 6x2 (6.1)

w = Klyz + szz + 62 + f(x,y, 9),

121

122
where I<1 and I<2 are curvatures, 69 unit angle of twist
and 6 unit extension. These more general displacements
represent the sum of the displacements for each loading
considered separately.

Again, the elastic analogy assumption on the
stresses would be made and an incompressible material
assumed. The development should parallel that of Chap-
ter II, but, in the solution, the manner in which the
curvatures and the extension depend on the unit angle of
twist would have to be specified. Numerical solutions
could be obtained with algorithms similar to those used in
this study. Since under the more general loading the
effect of the unit angle of twist on the normal stress is
greater, extreme care should be maintained in incrementing

the unit angle of twist variable.

6.1.2 General Boundaries

From a more practical viewpoint, it may be reason-
able to consider other geometric cross sections. The
present analysis treated only the square and the circular
bars. In each of these cases, the algorithms exhibit
reasonable convergence rates. For these examples, the
normal derivatives on the boundary are parallel to the
grid lines, and neither the warping function nor the
stress function offer a computational advantage. In the

case of an irregular boundary, however, the stress-

123
function approach might give an advantage, since the value
of the stress function is prescribed on the boundary while
the normal derivative is prescribed for the warping
function.

In the extension to a more general boundary, the
major problem would be computer programming. The program
must consider the different possibilities at each boundary
point. In addition, an irregular geometric cross section
may tend to cause some difficulty in the convergence of
the algorithms. The convergence properties of the Gauss-
Seidel method applied to the linear Poisson equation are

known to depend on the boundary geometry [32].

6.1.3 Other Theories

A number of attempts have been made to determine
mathematical stress-strain relations. Many of these are
summarized by Osgood [21]. Some of these theories may be
preferable to the ones used in this study. Since the pre-
sent formulation considers the constitutive relation as a
separate equation, the algorithm could be adapted to these
other theories, and the results of different theories
could be compared.

Perhaps one of the more interesting theories would
be one due to Prager [12], assuming the tensile stress-

strain curve as

o = Y tanh [23) . (6.2)

124
Prager [26] and Hill [12] give the counterpart generalized
flow theory for this stress-strain curve. From a physical
standpoint, one of the objections to this relation is the
difficulty in defining the parameters [21]. However, in
the case where the material has a well-defined yield point
and the actual stress-strain curve approaches the elastic-
perfectly plastic situation, this formulation provides a
good fit to the uniaxial physical data [12]. Thus, using
Prager's relation for the whole body provides a numerical
method of determining displacements when the material
behavior approaches that of an ideal elastic-perfectly
plastic medium.

Hill [12] gives an equation for determining the
warping in the case of pure torsion of a prismatic bar
with an elastic-perfectly plastic material. Hodge [14]
has applied this equation to some special cross sections.
In Hodge's formulation, the deformations must be such that
the elastic strains may be neglected. This restriction
would not be necessary in the numerical calculations with
a smooth stress-strain curve such as Prager's hyperbolic
tangent law.

The Ramberg-Osgood curve also approaches the
elastic-perfectly plastic case when n approaches infinity.
But, in the plastic region, the perfectly-plastic stress
is approached from above, whereas with Prager's method it

is approached from below. Consequently, in this special

125
case, it would seem that Prager's relation offers some
advantage over that of Ramberg and Osgood for approxi-
mating an elastic—perfectly plastic curve by a smooth
curve.
6.1.4 Work-Hardening Equations Not
Related to Combined Bending and Torsion

The iteration procedure developed in this study
permits the material non-linearity to be confined to a
single equation. This non-linearity, resulting from the
constitutive relation, is treated by the Newton-Raphson
iterative procedure. The governing equation, i.e., the
equilibrium equation with the warping function or the com-
patibility equation with the stress function, is still
treated numerically in a manner analogous to the corres-
ponding linear elasticity problem.

It may be possible to develop iterative procedures
similar to those used in this study for other work-
hardening problems not related to combined bending and
torsion. For example, the governing work—hardening equa-
tions for an Airy's stress function in plane stress or
plane strain problems may be developed. This development
would assume the entire body to be governed by work-
hardening plasticity equations so that no undetermined
elastic-plastic boundary is involved. The governing equa-
tion would be a fourth-order partial differential equation

which would be treated in a manner similar to the linear

126
biharmonic equation. The non—linear material properties
would be contained in an equation similar to the polyno-
mial equation considered in this study. The computational
algorithm should somewhat parallel that employed with the
stress-function approach of this investigation.

As in elasticity, the Airy stress-function approach
would be suitable only for problems with traction boundary
conditions. For displacement boundary conditions or mixed
problems, it might be possible to develop a similar proce-
dure for equations analogous to the Navier equations of
elasticity. Another possible class of problems are those
of axial symmetry including, but not limited to, torsion
of bars of non-uniform cross section. Also, no considera-
tion has yet been given to combined torsion and flexure
with transverse loads, introducing dependence on the

z—coordinate.

6.1.5 Experimental Verification

The validity of any engineering analysis must be
measured by the degree to which it actually predicts the
behavior of physical phenomena. To determine if this
method produces reasonable results, one would have to com-
pare these analytical results to experimental measurements.
Since the deformation and flow theories produce different
results under certain loading conditions, there is some

question as to the physical validity of the methods. It

127
should be possible to develop experiments which would
provide greater insight into the application of the
theories.

It would be quite difficult to subject a prismatic
bar simultaneously to arbitrary pure bending and torsion
where the rates of curvature and rates of twist could be
independently controlled. Hill and Siebel [13] have
reported on an experiment with a steel circular bar where
the ratio of moment to torque was held constant.

The problem of combined torsion and tension should
exhibit similar differences between the theories. It
should be feasible to design an experiment where the rate
of extension and rate of twist are independently controlled
to determine the relationship between tension and torsion.
This type of experiment may be a possible way of checking
the physical validity of the numerical procedures.

Ramberg and Osgood [21] suggest that their formu-
lation will provide a reasonable fit to a variety of
materials. For the uniaxial stress conditions, the three
parameters are easily obtained [21].

In order to obtain a valid mathematical stress—
strain relation for a given material, it may be necessary
to consider other theories such as those suggested in
6.2.3. If a reasonable fit to a given material is made,
then the above experiment should provide additional infor-
mation concerning the range of application for both the

deformation and the flow theories.

128
6.2 Mathematical Investigations

6.2.1 Bound on the Theta Increment

Experience indicates that the algebraic system for
both the flow and deformation theories converges at
reasonable rates if the increment in the angle (9 is
small. Thus, at any given level of G) , the convergence
is possible only if the initial solution is near the
actual solution. It appears that the governing factor on
the "nearness" of the initial solution is the Newton-
Raphson iteration on the (2n+1)th degree polynomial. That
is, if the initially estimated roots provide a good
approximation, then the system tends to converge.

The object of an additional investigation might be
to determine a bound on the increment AC). Previous
results show that this bound depends on the Ramberg-Osgood
exponent n and on the curvature-twist (or curvature rate-
twist rate) ratio p.

By experimenting with different increments, a bound
depending on n and u could be established experimentally.

Mathematically, the term Sm (for example, see Equation

+1
3.21) must be bounded for a given increment. With this
bound, a limit may be established on the degree of change

in the root<i>(i,j) at an arbitrary point (i,j).

129
6.2.2 Effect of Newton-
Raphson Iteration
The effect of the Newton-Raphson iteration should

be further studied. It appears that an equation of the

type.

A(u)uxx + B(u)uxy + C(u)uyy + D(u)ux + E(u)uY + F(u) = 0,
(6.3)

solved by iteration will tend to converge if the estimates
of the non-linear coefficients A(um), B(um), C(um), D(um),
E(um), and F(um) approach their respective limits at a
faster rate than um approaches u. This seems to indicate
that, if values of um approximate u in a relatively coarse
manner, some technique must be utilized by which the
coefficients are predicted quite accurately. If such
techniques can be developed, it seems reasonable to expect
the system to converge.

In a certain sense, the Newton-Raphson iteration
takes rather poor values of tm(i,j) and predicts quite
accurate values for Em(i,j). Thus, this iteration scheme
acts as a predictor in determining more accurate esti-
mates. It may be possible that other predicting tech-

niques may provide better convergence rates.

6.2.3 General Convergence
The question of general convergence must be con-

sidered from two standpoints: (1) will the discrete

130
system converge to the continuous system as the lattice
parameter approaches zero, and (2) will the iterative
method for the algebraic system converge. In non-linear
problems, very little theory is available in regard to
either of these two problems.

Recently, Parter and Greenspan [22], [23] con-
sidered these convergence problems with a mildly non-
linear second-order partial differential equation. They
considered the Dirichlet problem, on the union of R and S,
with

Vzu = f(u) on R

(6.4)
u = g(x,y) on S,

where R is the interior of the closed path S. In addition,
the function f(u) must be such that both its first and
second derivatives exist and are always greater than or
equal to zero.

Parter [22] considered the convergence of the dis-
crete system to the continuous system. Both authors [23]
treat the convergence of a number of different algorithms
for the solution of the discrete system. In each analy-
sis, the most important restriction appears to be the
bounds on the function f(u) and its first and second
derivatives.

It may be that the equation considered in this

study,

131

0216 + E(tx,ty) = o, (6.5)

also possesses similar non-linear properties. The problem
is to establish certain bounds on the function E(tx’ty)‘
If this is accomplished, then perhaps a procedure analo-
gous to that of Parter and Greenspan could aid in the
establishment of a general convergence requirement of the
system. In addition, with the warping-function formula-
tions, the analysis requires the consideration of the

Neumann boundary conditions.

APPENDIX I

SOLUTION TO FINITE-DIFFERENCE EQUATIONS BY

THE GAUSS-SEIDEL OVER-RELAXATION METHOD

This method, also known as Liebmann's method or the
method of successive displacements, is considered in most
books treating iterative methods for solving linear alge-
braic equations. Additional details may be obtained by
examining references such as Fox [5], Forsythe and Wasow
[4], Varga [33], etc.

To illustrate the procedure, the method will be
applied to the solution of the discrete representation of
the Poisson equation, where the equation is defined on a
square region R. However, the method may be applied to
other partial differential equations.

The Poisson equation is

uxx + uyy = F(x,y) (I-l)

where either u or the normal derivative is specified on
the boundary.

An N x N set of grid lines with lattice dimension h
is superimposed on the region R. Equation (I-l) is repre—
sented by its discrete finite-difference analogue. Thus,
at all mesh points

132

133
u(i+l,j) + u(i,j+l) + u(i-l,j)
+ u(i,j-l) — 4u(i,j) = h2F(i,j), (1-2)

where on the boundary either u(i,j) or the normal deriva—
tives are known.

If u(i,j) is known at each boundary point, then the
number of points where u(i,j) is unknown is (N — l)2.
However, if the normal derivatives are specified, then
expressing (I-2) along with the derivative approximation
permits the determination of u(i,j) on the boundary as a
function of the interior points. In this latter case, the
system can be solved within an arbitrary constant. In the
more general case, where u(i,j) is specified at some
boundary points and the normal derivative is specified at
the other points, the rank of the system will be equal to

the number of unknowns.

To begin the Gauss-Seidel over-relaxation proce-
dure, an initial solution uo(i,j) is assumed. These
values are recalculated during successive sweeps through
the mesh. A particular cyclic order for replacing u(i,j)
in any one sweep is specified.

If um(i,j) represents the values for the mth sweep,

then the values um+l(i,j) for the (m+1)th sweep are calcu-

lated by the formula,

um+1(i,j) = um(i,j) + QJ[fim l(i,j) - um(i,j)], (1-3)

+

134
where (U is a relaxation factor and fim+l(i,j) is the
unextrapolated value given by Equation (I-2). Thus,
fim+l(i,j) depends on its nearest neighbors and on the
cyclic order. If the cyclic order is such that for each
i, all values of j are computed before i is increased,

then Gm+l(i,j) is given by

fim+l(i,j) = [um(i+l,j) + um(i,j+l) + um+l(i-l,j)
+ um+l(i,j-1) - h2F(i,j)]/4. (1-4)

The relaxation factor (D may take on values

If (0 <= 1 the system is under-relaxed and if a) =6 l the
system is over-relaxed. For linear systems, various
theorems are available for estimating the optimum a) [4],
[32].

The iterations, or sweeps, are continued until some

convergence criterion is satisfied. A common criterion is

Max um+l(i’j) - um(i,j) <: 8, (I-6)

where the maximum is taken over all mesh points.

From the above, it appears that convergence depends
on the initial uo(i,j), the relaxation factor a), and the
specified cyclic order. However, the general question of

convergence depends on the properties of the system of

 

135
equations specified by (I-2). The system represented in
matrix form is

Au = B. (I-7)
The matrix A may be expressed as the sum

A=T-T—T (I-8)

D L U’

where TD is the diagonal of A such that ‘TD‘-l # 0, TL and

TU are respectively strictly lower and upper triangular
matrices, whose entries are the negatives of the entries
of A below and above the main diagonal. A necessary and

sufficient condition for convergence is that all the roots

of the polynomial,
P( A) = [( >( + (11- 1)1‘D - an )er + TU) = 0, (I-9)
have modulus less than one [33].4

The method may also be applied to the general

quasi—linear second-order equation,

Au + Bu + Cu + DuX + Euy + F 0, (I-10)

XX XY yy

where A, B, C, D, E, and F are functions of u, ux, uy, x,

and y. In this latter example, the (m+1)th sweep consi-

ders the equation

Amuxx + Bmuxy + Cmuyy + Dmux + Emuy + Fm = 0, (I-ll)

 

4Condition follows from Varga [33], Section 3.4,
page 75. Presented in this form by C. S. Duris, class
notes, Math 852, Winter Term, 1965.

 

136
where the subscript m denotes that these coefficients are
evaluated at each iteration as functions of um, x, and y.
The above method is applied with this additional consi-
deration that as m —>-OO , the solution to (I-11) con-
verges to the solution to (I—lO). At the present time,
general theorems giving the conditions for convergence are

not available.

APPENDIX II

NEWTON-RAPHSON ITERATIVE METHOD FOR

SOLVING A NON-LINEAR EQUATION

The general problem is to find the solution or

zeros of the equation

In general f(x) is a non-linear function, and the problem
reduces to that of finding the zeros of (II-l) through
some iterative sequence. The Newton-Raphson iterative
method is a procedure for determining a single root (5 if
an initial estimate of the root is near enough to g'.

The technique may be applied to either real or
complex roots. However, in the following the interest is
only in determining a single real root. In this case, the
method corresponds to using the tangent to the curve f(x),
at the last determined point[xi,f(xi)] to evaluate a new
approximation x.+

i l'

the new approximation to xi+

From Figure 16, it is observed that
1’ according to the tangent
approximation, is

x1+1 = x1 " TITS??? (II'Z)

138

f(x)“

 

Figure 16. Estimation of new
approximation by tangent

A more detailed discussion of the properties of
(II-2) and the iterative method may be obtained from Todd
[32] and Hildebrand [10]. However, the convergence is
dependent on the error of the initial approximation.
Figure 16 illustrates that if this approximation is near
enough to the root €:, then the method will always
converge.

If the method converges to fir, and k is taken

sufficiently large, then the error at the k + 1 iteration

is [10]

~ 1.5.45) 2

(5 - xk+l) ~ - .2 (f- xk) . (II-3)
Thus, the error in xk+1 is proportional to the square of
the error in xk. Iterative methods having this property

are called second-order processes and should provide

reasonable convergence rates.

APPENDIX III
BASIC FLOW CHART PATTERNS

The basic flow pattern for each program is shown
in the following flow charts. The details of the calcu-
lations may be obtained by examining the equations listed
on each flow chart. Additional information on the
specific equation of each program is contained in the

program listings of Appendix IV.

139

140

Hmm<3 Emumoum mom cumuuma uumnu 30am uemmm .sH wuzmam

‘

 

 

 

R... 2...;
2 :L

as .65 “em mwc om mmmmmmem
s » BZHma
aszm mesmzoo

 

 

 

 

 

 

 

 

 

 

 

 

 

 

\) ANH.mv .Om Anafivp AMH.mV .0m
2 momHmmHeem BZHmm mmmmmeem
// mozmomm>zou meomzoo

 

 

 

 

 

 

 

e

AHH.mV .Om moemx

 

 

 

 

// An.wvu mmadzHemm .2 .«d .1
mesmsou q<H9H2H sum o<mm

 

 

 

 

 

 

omm<3 Emumoud Mom cuwuema uumnu 30am UHmmm .ma mufimam

 

®6+@u®

 

141

9mm

 

 

 

Z Z

 

 

 

 

ae®n® /HV\18«moemx/ s
ex/ \

66 .6 16H.mc om mummmmem
u , c .6
asza 6 6 asza
6956206

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ANN.mv .0m 6 a . AMN.mv .CQ
.6 \\ aamHmmHe<m in vamemfi as» mmmmmmem
/rmozmomm>zoo mesmzou.

 

 

 

 

 

 

 

 

e

moamx

 

 

 

AHN.mv

 

 

.Amfi.mc .mom mme<zaemm

 

 

6966200 .45 .1 o<mm

 

 

 

 

E
ahoHvAU “Afiaflvu. WHBHZH Bum ad “2 n ® o®q

 

 

 

omem Smumoum wow Gumwumm uumnu 30am vammm .ma unseen

 

 

142

 

 

 

 

 

z
. /
a mu ".AVUVIIWWAHNO n dogma
» 1\.»
6e .6» 1mm me: on 66666696
eszm , e . u BZHmm
nebmzoo .

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Aam.mc .om . . . 1mm.mc .66
z eanmmHe<m an «cmemwn eye. , mmmmmmem
/rmozmomm>zoo mesmzoo
WL .AOm.mc moemx
Am~.mv .mom mme<zHemm . . .e .
“n.63nv .hfi.acs7 aeaeHzH 966 6..u .amuaemwe
mesmsoo ‘

 

 

 

 

 

 

 

>qmm¢3 Emumoud mom sumuumm 30am uwmmm .om wuzmam

 

®<+®u®

 

143

Bum

 

 

 

Z Z

 

 

mefiv u mo

 

 

no u moamx /

K .

h?
66 .6» new mcc om mmmmumem
e . u azamm
BZHmm mammzoo

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ANm.mv .OH :2 a a . Amm.mv .0m
2 \ SmHmmHefi Q. . e<eeMHWmn fie mmmmmmem
/ 66266665.on mesmsou
mw .Aom.mc . \\\( 6096s
Amq on omom m. a mmaéHm—wmm a a a5 a
22262 I??? 6533 666 ”.32. ® @<
webmzoo .. . .X a<mm

 

 

 

 

 

 

 

Baamem Emumoum now cumu#ma unmzu 30am uwmmm .Hm musmam

 

qu + AU nAuv

 

144

9mm

 

 

 

 

2 z

 

 

 

 

me® "@Ifl KS n 85 / w
x A \

:2 .cu Amm.mvc.cm mmmmmmem
BZHmm _ a . u , BZHmm
mapmzou

 

 

 

 

 

 

 

 

 

 

 

 

 

 

\1 Amm.mv .Om Awm.mv .Om
«omHmmH94m Awpm.fla<;fiwpn.av mmmmmmam
z/ruuzmwmm>zou aszm mabmzou

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

    
 

moamzd

E
om «Zao ® a®<
1V x a<mm

Amm.mv
.Aam.mv .mom mma<zaamm
Q43; .A 3.23% .5353 9mm
meamzou

 

 

 

 

 

 

 

 

145

mm<qom Emuvoum mom cumuuma unmcu 30am oammm .mm musmwm

 

 

 

 

 

 

 

    

a c “HH.¢V .ou ‘
e . u :5 .cu mmmmmmam
BZHmm mapmzou BZHmm

 

 

 

 

 

 

 

 

 

 

 

 

z «omHmmHa<m w BZHmm mmmmmmam
/Kmuzmwmm>zoo mepmzou

 

 

 

 

 

 

 

 

 

 

 

 

Ammfiwwavmm mme<zHamm .zmmwwu.1
mammzou q<HaH2H 9mm
a<mm

 

 

 

om<qom Emumoum Mom aumpuma unmnu 30am uammm .mm musmwh

 

mud + 0V "AU

 

146

8mm

 

 

 

Z

 

 

 

 

 

z
~e® n® / Raoumoamx
K»)

L\ .w
as .c AHH ¢V om mmmmmmam
u as .cu aszm
aszm mapmzou

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Amm.¢v .Om . . . a u Aoa.vv .Om
2\\~amHmmHa<m An avmemwn av mmmmmmam
/rmuzmwmm>zou masmzou

 

 

 

 

 

 

 

 

 

 

 

MBDmZOU

.Aam.¢v fl mosmx
Amm.¢v .mom mme<zHemm .: . .5 .
a a a 2. mu 6
Ah avnv Ah Hvu A<H9H2H sum .qu .1 n<ww

 

 

 

 

 

 

APPENDIX IV
FORTRAN PROGRAMS

The listings for all FORTRAN programs are provided
in this appendix. The purpose of each program and the
page listings are as follows. WARPI (pages 148-50) solves
the Piechnik equation for a unit square section. WARPO
(pages 151-54) solves the warping-function deformation-
theory equations for a unit square cross section. STRO
(pages 155-58) solves the stress-function deformation-
theory equations for a unit square cross section. WAPPLV
(pages 159-63) solves the warping-function deformation-
theory equations for a unit square cross section. STPLDT
(pages 164-67) solves the stress-function deformation-
theory equations for a unit square cross section. POLARR
(pages 168—70) solves the Piechnik equation for a unit
radial cross section. POLARO (pages 171-72) solves the
warping-function deformation-theory equation for a unit
radial cross section.

All programs follow the coordinate representation
of Piechnik [24] while the theory follows Hill [12].
Program variable names generally agree with those of the
theory. Major exceptions are T for ND, LAMDA for<i>, and
TAU for the shear stresses.

147

148

DPOGRAM WARP!
TYPE REAL LAMDA
DIMENSION T‘25025)¢Tl(25¢25)9T2‘25925)OSIGMAZ‘ZSOZS)CTAUYZ‘25625,0
5 TAUXZ‘ZSOZS)
1 READ 39PHIQBETA0NQITERQKSTOP
3 F0RMAT‘2F10059315)
H = 005 / (N'lo) 5 IT 3 l 5 E1 = loE-6 5 PH! = pHI**2
DO 5 I'ION 5 DO 5 J=19N
SIGMAZKIQJ) 3 000 5 TAUYZ‘IOJ) = 000 5 TAUXZ(IQJ) 3 000
T1(IQJ) = 000 5 T2‘19J) 3 000
5 T(toJ) = 000
6 DO 99 I =19N 5 DO 99 J=10N
CHI 9 (1-1)*H 5 ETA = (J-1)*H
IECI*1)17999017
l7 IF(I-N)21019021
l9 IF(J-N)25923025
CASE I=N J=N
23 T(IQJ) 3 (T(I~10J)+T(IQJ-l))/Zo 5 GO TO 99
25 [F(J‘l)35099035
CASE I=N J NOT 31 OR N
35 5X 3 (3*T(IOJ)-4*T(I‘1OJ)+T(I‘20J))/(2*H)
SY = (TCIQJ+1)-T(IQJ‘1))/(2*H)

A = 3*pHI*ETA**2

B = 000

C = 3*pHI*ETA**2 + (SY-CHI)**2

D = -3*PHI*ETA*H/2 5 E = 3*PHI*ETA*CHI*H**2

T(IQJ)' (2*C*(T(I‘loJ)‘H*ETA)+(A+D)*T(IoJ+1)+(A-D)*T(IqJ-l)‘B+E)/
5 (2*(A+C))
GO TO 99
21 IF(J-l)29099929
29 IF(J-N)31033931
CASE I NOT = 1 OR N J=N
33 5X 8 (T(I+10J)+T(I-10J))/(2*H)
SY = (3*T‘IOJ"4*T(ICJ-1’+T(IQJ‘2)’/(2*H)
A 8 3*pHI*ETA**2 + (SX+ETA)**2 5 B = 000
C = 3*PHI*ETA**2
D 3 ‘3*PH1*ETA*H/2 5 E 3 3*9HI*ETA*CHI*H**2
TCIQJ) = (C*(T(I+l0J)+T(I-10J))+(A+D)*(2*H*CHI+T(IQJ‘I))
5 +(A’D)*T(IoJ-I)+B+E)/(2*(A+C))
GO TO 99
CASE I NOT a 1 09 N J NOT a 1 OR N
31 SX = (T(I+19J)‘T(I-10J))/(2*H)
SY 8 (TCIOJ+1)-T(I9J”1))/(2*H)
A 3 3*PHI*ETA**2 + (SX+ETA)**2
B 3 -(SX+ETA)*(SY-CHI)/Zo
C = 3*PHI*ETA*§2 + (SY—CH!)**2
D = -3*PHI*ETA*H/2 5 E = 3*PHI*ETA*CHI*H**2
T‘IOJ’ 3 (C*(T(I+l0J)+T(I‘19J))+(A+D)*T(IQJ+1)+(A-D)*T(IQJ-1)+
5 B*(T(I+19J+1)‘T(I+19J‘1)’T(1-10J+1)+T(I-19J-1))+E)/(2*(A+C))
99 CONTINUE 5 X = 0.0
DC 301 I=10N 5 DO 301 J=10N
A] a ABSE‘ABSF(T(IOJ)’-ABSF(T2(IQJ)3)
IF(AI‘X)30103059305

()

149

305 X=A1
301 CONT1NUE
IFCX‘E11200051‘51
51 DO 53 1310N 5 DO 53 JEION
TZCIQJ) 3 T110J1
T1119J) 3 T1119J1 + BETA*(T(19J)-T1(10J1)
53 T(19J) = T1(10J)
1T 3 1T + 1 5 1F(1T-1TER)694OOO4OO
200 PRINT 201.1T98ETA9PH10‘1TC19J101310N19J'10N)
201 FORMAT (1H19* TC19J) VALUES FOR ITERATION * 9 1100//* RELAXATION
5FACTOR *9F10050* PHI VALUE *9F1005o//(7E1708 )1
GO TO 112
400 PRINT 40101TOBETAOPHI
401 FORMAT11H19* FAILED TO CONVERGE AFTER * 01100*1TERAT10NS*

S ///* BETA EQUAL*0FIO.5o* PHI EQUAL*¢F10¢51
GO TO 105
112 DC 199 1=loN $ 00 199 J=10N

CH1 3 (I*1)*H 5 ETA 3 (J-1)*H
IFCI~11117v1079117
107 1FKJ*1)10991990109
CASE I=1 J81
110 SX 8 T(1+10J1/H
SY : T(1oJ+11/H 5 GO TO 151
109 1F¢J¢N111191139111
CASE 131 J NOT 8 1 OR N
111 SK 3 T11+1cJ1/H
SY 8 (T(19J+1)-T(10J-1))/(2*H1 5 GO TO 151
CASE 1:1 J=N
113 5X 8 T11+10J)/H
SY 8 CH1
GO TO 151
117 1F(1-N)12101190121
119 1F(J~N)12591230125
CASE 1 8N J=N

123 SX = - ETA
SY 2 CHI
60 T0 151

125 1FCJ-1113591370135
CASE 1=N J81
137 5X 8 ~ ETA
SY 3 T119J+11/H
GO TO 151
CASE 1=N J NOT 8 1 OR N
135 5X 8 - ETA
SY 8 ( T110J+1)‘ T(19J-1))/(2*H)
GO TO 151
121 1F(J-1)12901279129
CASE 1 NOT 3 1 OR N J31
127 5X 3 1 T11+10J)“ T(1-10J11/(2*H)
SY 3 T(19J+11/H
GO TO 151
129 1E(J-N)13101330131
CASE 1 NOT 3 1 09 N JaN

133

131

151

199
411
413

415
900

911

901

903 FORMAT11HOQ* TORQUE EGUALS *0F10080* BEN01NG
5 0F1OOB)

150

SX = ( T(1+19J)~ T(1*1.J))/(2*H)

SY = CHI
60 TO 151
CASE NOT 8 1 0R N J NOT = 1 OR N

1
5X = ( T11+10J1“ T11’19J))/(2*H)
SY 8 ( T(19J+1)- T(19J-1))/(2*H1

LAMDA 8 SORTF( 2.25*PH1*ETA**2 + 0075*($Y**2-2*CHJ*SY+CHI**2 +
s SX**2 + 2*ETA*SX + ETA**2))

SIGMAZC10J) = 105*SQPTF(PH11*ETA/LAMDA
TAUXZ‘IOJ) I31005*(SX-fi’ETA1/LAMDA)‘I'SGPTF1301
TAUYZC19J) 8(005*(SY*CH1)/LAMDA)*SORTF(3C)
CONTINUE

PRINT 4110((51GMAZ(10J)01=1QN)QJ=10N)
FORMAT(1HO¢* SIGMA Z STPESSES§//(7E1708))
PRINT 4139((TAUXZC10J)91=19N)0J310N)
FORMAT¢1HO¢* TAU XZ STPESSES* // (7E17081)
PRINT 4159((TAUYZ‘10J)01:10N19J810N)
FORMAT11H09* TAU Y2 STPESSES* //(7E1708))
DO 911 1:19N

DO 911 JrloN

CH1 3 (1-1)*H 5 ETA 3 (J-1)*H

T1(19J) = ETA*S1GMAZ(10J1

TCIOJ) 8 ETA*TAUXZ(10J) - CH1*TAUYZ(10J)
TORO 3 0.0 5 BEMO I 000

DO 901 1=ZQN12

DO 901 J=20N02

T090 = TOPO + T(1-19J‘1)+T(1+19J-1)+T(1-1vJ+1) + 40*(T110J-1)+
S T(1+19J)+T(1.J+1)+T(1*1.J)1 + 16¢*T(I¢J)
BEMO = BEMO + T1(I‘1oJ-1)+T1(1+19J-1)+T1(1+1oJ+1)+T1(1-10J+1)
5 + 4.*(T1(19J‘1)+T1(1+1cJ)+T1(IoJ+1)+T1(1~1.J)) + 16o*T1(19J)

HSQ = (H**2)/9o
TOPO = 12*TOPG*HSQ
BEMO = 16.*HSQ*BEMO
PP1NT 903vTORQoBEMO

105 IFCKSTOP12150 1 9215
215 STOP 5 END

+T(1+19J+1)

MOMENT EOUALS

*

151

PROGRAM WARPO
TYPE REAL LAMDAOJZOLAMDAXOLAMDAY
DIMENSION T125925)9LAMDA(25925)0T1(2592511J21250251
s .sxGMAZ(2s.25).TAux2(25.25).TAUYZ(25.25)
5 QBEN(25925)0TOR(25025)
1 READ 39PH1¢ BETAQDTHETAvTHETAMcNoNEoITERcKSTOP
3 FORMAT (4F100504151
PRINT AONEQNQBETAQPHI
4 FORMAT11H19* DEFORMATION THEORY OUTPUT FOR NE EOUAL*0150* N EOUA
5L*9159* BETA EOUAL*0F6.39* MU EOUAL*9F6.3¢/ 0
5* IT THETA TORO MOMENT SIGMA 1 TAUXZ 1 TAUYZ 1 SIGMA
5 2 TAUXZ 2 TAUYZ 2 SIGMA 3 TAUXZ 3 TAUYZ 3*)
M = N-l 5 H = 0.5/M 5 IT = 1 5 KP = 5 5 E1 = IoE-T
NEZ 2*NE 5 NE2N1= NEZ - 1 5 THETA 8 DTHETA
PHI PHI**2
DO 5 I=loN
DO 5 J=19N

STNIIQJ) = 000 5 STX(19J) = 0.0 5 STYIIOJI 3 0.0
T(19J) = 0.0 5 TIIIQJ) = 000
SIGMAZIIQJ) 3 0.0 5 TAUXZ(19J) = 0.0 5 TAUYZIIQJ) 3 0.0
JZIIQJ) = 0.0
5 LAMDACIQJ) = 1.0

112 DO 199 1=IQN
DO 199 J=19N
X = (I-I)*H 5 Y = (J-1)*H
IECI-1111791070117
107 IEIJ-1110991990109
109 IEIJ-N111191130111
CASE I=1 J NOT 3 1 OR N
111 SX = TII+19J)/H
SY = IT‘IOJ+1)-T(IOJ‘1)’/(2*H)
J2(19J1 = (THETA**2)*(30*PHI*Y**2 + (SX+Y)**2 + (SY-X)**2)
LAMDA(19J1 = LAMDACIOJ) - (LAMDA(IQJ)**NE2*(LAMDA(I1J1‘10) ‘
5 JZIIQJ)**NE) / (LAMDA(IQJ)**NE2N1 *((2*NE+1)*LAMDA(IQJ) “ 2*NE11
GO TO 199
CASE I=1 J=N
113 SX = T(I+loJ)/H
SY = X
J2(10J1 = (THETA**2)*(30*PHI*Y**2 + (SX+Y)**2 + (SY-X)**2)
LAMDAIIQJ) = LAMDA(10J) - (LAMDA(IQJ)**NE2*(LAMDA(I9J)-Io) ‘
5 J2(10J)**NE) / (LAMDA(IQJ)**NE2N1 *((2*NE+1)*LAMDA(IQJI - 2*NE11
GO TO 199
117 IF(I“N)12191190121
119 IF(J-N)12591230125
CASE I=N J=N
123 TBAR = (TII’IQJ)+TIIoJ-1))/20
TIIQJI = T(19J) + BETA*(TBAR - T(IOJ)1
J2(IOJ) = (THETA**2)*(3.*PHI*Y**21
LAMDAIIQJ) = LAMDA(19J) — (LAMDAIIOJ)**NE2*(LAMDA(IQJ)“101 ‘
5 J2(19J1**NE) / (LAMDAIIQJ)**NE2N1 *((2*NE+11*LAMDA110J) - 2*NE11
GO TO 199
125 IFIJ-1113501379135
CASE I=N J=1
137 SX = «Y

-3

152

SY 8 TIIQJ+11/H
J2(IOJ) = (THETA**2)*(30*PHI*Y**2 + (SX+Y)**2 + (SY‘X)**2)
LAMDAIIOJ) 8 LAMDAIIOJ) - (LAMDA(10J)**NE2*(LAMDA(IOJ1‘101 “
5 J2(1¢J)**NE) / (LAMDACIQJ)**NE2N1 *((2*NE+1)*LAMDA(19J) ~ 2*NE11
GO TO 199
CASE I=N J NOT = 1 OR N
135 SY=(T(I¢J+1)‘T110J*1))/(2*H1
LAMDAY = (LAMDA(10J+1) ’ LAMDA(IoJ-1))/(2*H1
E = LAMDAY*(SY-X)*(H**21
TBAR 8 (LAMDA(IQJ)*(2*(T(I-19J)-H*Y)+ T(IoJ+1)+ T(19J‘11) -E)
5 /(4*LAMDA(IOJ))
T(IoJ) 8 T(IoJ) + BETA*(TBAR - T119J1)
J2(I¢J) 3 (THETA**2)*(3¢*PHI*Y**2 + (SY‘X)**21
LAMDA(10J) 8 LAMDAIIQJI - (LAMDA(19J)**NE2*(LAMDA(IQJ)“101 ‘
5 J2(19J)**NE) / (LAMDAIIQJ)**NE2N1 *((2*NE+1)*LAMDA(IQJ) ‘ 2*NE11
GO TO 199
121 IF(J“1)12901279129
CASE 1 NOT = 1 OR N J=1
127 SX‘= ( TCI+1¢J)- T(I-1oJ))/(2*H)
SY = T119J+1)/H
J2(10J) = (THETA**2)*(30*PHI*Y**2 + (SX+Y)**2 + (SY-X)**2)
LAMDA110J) a LAMDA(1¢J) - (LAMDA(IOJ)**NEZ*(LAMDA(IQJ)-101 ”
5 J2(1¢J)**NE) / (LAMDA(19J)**NE2N1 *((2*NE+1)*LAMDA(IOJI “ 2*NE11
GO TO 199
129 IF1J~N113191330131
CASE 1 NOT ‘-'-'-' 1 OR N J=N
133 SX = (T(I+19J)-T(I-10J))/(2*H)
LAMDAX = (LAMDA(I+19J) - LAMDA(I-10J))/(2*H)
E = LAMDAX*(SX+Y)*(H**2)
TBAR = (LAMDACIQJ)*(2*(T(IoJ-1)+H*X)+TII+10J)+T(I-10J)I-E)
5 /(4*LAMDA(I¢J))
TIIQJ) ‘ TIIQJ) + BETA*(TBAR - TIIQJ)1
J2(19J) = (THETA**2)*(30*PH1*Y**2 + (SX+Y)**2)
LAMDA(19J) = LAMDA(19J) - (LAMDA(1QJ)**NE2*(LAMDA(IvJ)-101 “
5 J2(10J)**NE) / (LAMDA(IoJ)**NE2N1 *((2*NE+1)*LAMDA(19J) - 2*NE11
GO TO 199
CASE 1 NOT = 1 OR N J NOT = 1 OR N
131 SX = (T(I+19J1-T(I-10J)1/(2*H)
SY 8 (T119J+1)-T(19J‘1))/(2*H)
LAMDAX = (LAMDA(I+19J) - LAMDA(1-10J))/(2*H1
LAMDAY = (LAMDA(19J+1) ‘ LAMDA(IQJ~1))/(2*H)
E =(LAMDAX*(SX+Y) + LAMOAY*(SY-X)1*(H**2)
TBAR = (LAMDA(10J)*(T(I+1oJ)+ TCIQJ+11+ T(I-19J)+ T(IQJ‘1)) ‘E)
5 /(4*LAMDA(I¢J))
T(IoJ) 8 T(IvJ) + BETA*(TBAR - TIIQJII
JEIIOJ) = (THETA**2)*(30*PHI*Y**2 + (SX+Y)**2 + (SY-X)**2)
LAMDA(10J) = LAMDA(10J) - (LAMDAIIOJ)**N52*(LAMDA(IOJI‘10) ‘
5 J2(19J)**NE) / (LAMDA(19J)**NE2N1 *((2*NE+1)*LAMDA(IQJ1 ~ 2*NE11
199 CONTINUE
IF(1T-ITER)179179411
17 X = OOO
DO 61 I‘lQN
DO 61 J=19N

I

63
61

16
51
400

401

404
405

317

407

153

A1 = ABSF(ABSF( T(19J))‘ABSF(T1(IOJ)1)
IF(A1‘X)61963063

59F10o5o/(7E17o811

PRINT 4059 ((LAMDA(IQJ)01=10N19J=19N1
FORMAT(1HO¢* LAMDA VALUES* /(7E17oB))
DO 317 I=10N

DO 317 J=loN

J2(1¢J) = J2(19J) / LAMDA(19J)**2
PRINT 4079 IIJ2110J101=19N19J=10NI
FORMAT!1HO¢* J2 VALUES * /(7E17.B))

COMPUTE NORMALIZED STRESSES

200

207

209

211

213

217

219

DO 299 I=loN

DO 299 J=IoN

X = (I-I)*H 5 Y = (J-1)*H
IF(1-1)217o207o217
IE(J*1)2O992999209
IF(J-N121102139211

CASE I=1 J NOT = 1 OR N

SX = T(I+1oJ)/H

SY : (T(IoJ+l)-T(IcJ-1))/(2*H)
GO TO 251

CASE 181 J=N

5X 8 T(1+19J)/H

SY = X

GO TO 251

IF(I-N)22102190221
IF(J-N)22592239225

CASE I=N J=N

223 SX = -Y

SY 2 X
60 TO 251

225 IF(J-1)23592379235

. CASE I=N J=1
237 SX

= -Y
SY = T119J+11/H
GO TO 251

CASE I=N J NOT = 1 OR N

235 SX = ‘Y

221

227

SY = (T(IqJ+1)-T(IoJ-1)1/(2*H)
GO TO 251
IF(J-1)229¢2279229

CASE 1 NOT 8 1 OR N J:1

sx s'(T(I+1.J)-Tc!-1.J))/(2*H)

NE EQUAL *
BETA EQUAL

X = A1

CONTINUE

IF(X-E1)400.16016

DO 51 I=loN

DO 51 J=19N

T1(19J) = T(IoJ)

IT = IT + 1 5 GO TO 112

PRINT 401QITQNEITHETAOPHIoBETAO‘IT‘19J191=10NIOJ=19NI
FORMAT(1HOQ* T(I.J) VALUES FOR ITERATION*91100*
5oIlOo/* THETA EQUAL *9F10059* PHI EQUAL *9F10.59*

*

154

SY a T(IoJ+1)/H
GO TO 251
229 IEIJ‘N123102330231
CASE I NOT = 1 OR N JBN
233 SX = (T(1+19J)‘T(1-10J))/(2*H)
SY a X
GO TO 251
CASE I NOT = 1 OR N J NOT a 1 OR N
231 SX = (T(I+1oJ)-T(I-19J))/(2*H)
SY = (TIIoJ+l)-T(IoJ~1))/(2*H)

251 SIGMAZ(I.J) = SORTF(3.*PHI)*THETA*Y/LAMDA(IoJ)
TAUXZ(I¢J) = THETA*(SX+Y)/LAMDA(19J)
TAUYZ(I9J) = THETA*(SY-X)/LAMDA(I¢J)

299 CONTINUE

4001 PRINT 4IIo((SIGMAZ(IcJ)oI=IoN)oJ819N)

411 FORMAT(1HO.* SIGMA Z STRESSES * /(7E17.8))
PRINT 4130((TAUX2119J191=19N19J819N1

413 FORMATIlHOo* TAU XZ STRESSES * /(7E1708)1
PRINT 4159((TAUYZ(19J101=10NI9J819NI
415 FORMAT(1HOQ* TAU YZ STRESSES */(7E1708))

4002 DC 911 IzloN
DO 911 J=10N
X = (1-1)*H S Y = (J-1)*H
BEN(19J, SIGMAZIIOJ1*Y
911 TORIIOJ) : Y*TAUXZ(19J) - X*TAUYZ(I.J)
T090 z 0.0 5 BEMO : 0.0
900 DC 901 1:29N92
DO 901 J=29N92
TORO = TORG + TORII-loJ‘l)+TOR(1+19J-1)+TOR(1-10J+1)+
S TORC1+19J+11 + 40*(TOR(1+19J1+TORI10J+11+TOR(1-10J)
$ +T09110J‘1)’ + 160*TOR(10J1
901 BEND = BEMO + BEN(I-10J-11+BEN(1+10J-11+BEN11-10J+1)
5 +BEN(1+19J+1) + 40*(BEN(1+10J1+BEN(IOJ+11+BEN(I‘19J)
5 + BENIIQJ‘I)’ + 160*BEN(ICJ)
H50 3 H**2 /9o
T090 3 120*TORQ*HSQ
BEMO = 16.*H50*BEMO
PRINT 9039TORQQBEMO
903 FORMAT11H09* TORQUE EQUALS *9F10050* SENDING MO‘ENT EQUALS*9
$ F1008)
500 THETA = THETA + DTHETA
IF‘THETA‘THETAMISOIOSOI9411
501 IT 2 1 5 GO TO 112
411 IFIKSTOPI413¢1Q413
413 STOP 5 END

V

a
U

(1

1

155

PROGRAM STRO
TYPE REAL LAMDAOJZOLAMDAXQLAMDAY
DIMENSION T(25925)sLAMDA<25o2510T1¢2592510J2(250251

S vSIGMAZIZSoZS)oTAUXZ(25925)oTAUYZ¢250251

READ 30PHI. BETAoDTHETAqTHETAMoNgNEoITERoKSTOPoKODE

3 FORMAT¢4F10¢50515 1

5

S J211¢J)**NE) / (LAMDA(IqJ)**NE2N1 *((2*NE+1)*LAMDA(19J1

$ J2(1.J1**NEI / (LAMDAIIoJ)**NE2N1 *((2*NE+1)*LAMDA(IoJ1

M a N-l s H = 0.5/M 3 IT = 1 s KP : 5 s E1 = 1.Eu6
NEZ = 2*NE s NEZNIB NE2 — 1 s PHI = PHI**2

THETA x DTHETA s 00 5 1:1.N s 00 5 JlloN

TIIoJ) = 0.0 s TIIIoJ) = 0.0

SIGMAZIIoJ) = 0.0 s TAUXZIIoJ) : 0.0 s TAUYZ(I.J) a 0.0
J2111J1 9 0.6

LAMDAIIoJ) = 1.0

112 DO 199 I=1oN S DO 199 J=1oN

X = (1-1)*H S Y = (J-1)*H
IFCI-1111701070117

107 IFIJ-11109011149109
CASE I=1 J21
1114 5x = 0.0 5 SY = 000

E = (LAMDAX*SX + LAMDAY*SY + 2*THETA)*(H**21/LAMDA(IoJ)
TBAR 8 (TII+1.J) + T(1¢J+1) + T11+1¢J1 + T119J+11 + E)/4o
TIIoJ) = TIIoJ) + BETA*(TBAR - TIIQJII $ GO TO 199

109 IF(J-N)11191139111
CASE I=1 J NOT = 1 OR N
111 SX = 0.0

SY 8 (T(IvJ+1)-T(10J-1))/(2*H)

LAMDAX 3 0.0

LAMDAY = (LAMDACI9J+1)-LAMDA(IoJ-11)/(2*H1

E = (LAMDAX*SX + LAMDAY*SY + 2*THETA)*(H**2)/LAMDA(IoJ)
TBAR 8 (T(1+10J) + TIIoJ+11 + T(I+10J) + TIIoJ-11 + E1/4o
T‘IOJ’ 8 T(IoJ) + BETA*(TBAR * TCIOJ1)

JZIIQJ) = 3*PHI*(Y**2)*(THETA**2) + (LAMDA(I¢J)**2)*(SX**2 +SY**2)

LAMDAIIQJ) a LAMDAIIQJ) - (LAMDAIIoJ)**NE2*(LAMDA(IQJ)-101 ‘

s JZIIoJ)**NE1 / (LAMDA(19J)**NE2N1 *(¢2*NE+1)*LAMDA(1qJ) - 2*NE11

GO TO 199

CASE 181 J=N
113 SX 8 0.0

SY 8 (3*TIIvJ)-4*T(IoJ-1)+T(1oJ-211/(2*H)

JZIIQJI = 3*pHI*(Y**21*(THETA**2) + (LAMDA(I¢J)**2)*CSX**2 +SY**2)

LAMDAIIoJ) = LAMDA(I¢J) - (LAMDA(19J1**NE2*(LAMDA(19J1-101 ”

GO TO 199

117 IFtI-N)!21.119.121
119 1F(u-N)125.123.125
CASE I=N J=N
123 5x = 0.0 3 SY no.0
J2(I.J) = 3*PHI*(Y**2)*ITHETA**2) + (LAMDA(I.J)**2)*(SX**2'+$Y**21

LAMDAIIQJ) I LAMDA(I¢J) “ (LAMDACIQJ)**NE2*(LAMDA(10J1-101 ‘

GO TO 199

125 IF(J-111359 1379135

CASE I‘N J=1

137 SX 8 (3*T119J1-4*T(1-1.J)+T(I-20J))/(2*H1

- 2*NE11

- 2*NE))

156

SY 8 OoO
J2(IoJ) = 3*PH1*(Y**2)*(THETA**2) + (LAMDA111J1**21*ISX**2 +SY**2)
LAMDA(10J) 8 LAMDAIIOJ) - (LAMDAII9J)**NE2*(LAMDA(I¢J)-1o1 “

$ J2(IoJ)**NE) / (LAMDA110J)**NE2N1 *1(2*NE+1)*LAMDA(IOJ) - 2*NE11

GO TO 199
C CASE I=N J NOT = 1 OR N
135 SX = (3*TIIoJ1-4*T(I-1oJ)+T(I-29J))/(2*H) 5 SY 8 000

JEIIQJ) 3 3*pHI*(Y**2)*(THETA**2) + (LAMDA(19J)**2)*(SX**2 +SY**2)
LAMDAIIQJ) = LAMDAIIQJ) - (LAMDAII0J1**NE2*(LAMDA(IQJ1*101 ‘
3 J2CIOJ)**NE) / (LAMDA119J)**NE2N1 *((2*NE+11*LAMDA(19J) ‘ 2*NE11
GO TO 199
121 IFIJ-1112901279129
C CASE 1 NOT = 1 OR N J31
127 SX = I T11+19J)’ T11-10J11/(2*H) 5 SY = 0.0
LAMDAX 3 (LAMDA(1+19J) “ LAMDAII’19J11/(2*H1
LAMDAY 3 0.0
E 8 (LAMDAX*SX + LAMDAY*SY + 2*THETA)*(H**2)/LAMDACI9J1
TBAR 8 (T(1+19J) + T(10J+1) + T(I*19J) + T(10J+1) + EI/4o
TIIQJ) 3 TIIOJ) + BETA*(TBAR ~ T‘IOJ’)
J2119J1‘2 3*pH1*(Y**2)*(THETA**2) + (LAMDA110J)**2)*(SX**2 +SY**2)
LAMDAIIQJ) = LAMDAIIQJI “ (LAMDAIIQJ)**NE2*(LAMDA(I9J1-1o) ‘
$ J2(IQJ1**NE) / (LAMDA(IQJ)**NE2N1 *((2*NE+1,*LAMDA(IOJ) ’ 2*NE11

GO TO 199
129 IFIJ~N)13101330131
C CASE I NOT = 1 OR N JBN

133 5X 3 (T(I+19J)-T(1-19J1’/(2*H1
SY 8 (3*T(19J1 -4*T(IOJ‘1) + TIIOJ“2))/(2*H)
J21IOJ) = 3*pH1*(Y**2)*‘THETA**2) + (LAMDAIIQJ)**2)*(SX**2 +$Y**2)
LAMDAIIQJ) 3 LAMDACIIJ) - (LAMDA(10J)**NE2*(LAMDA(IOJ1‘101 ‘
$ J2(10J1**NE1 / (LAMDAIIQJ)**NE2N1 *((2*NE+11*LAMDA(IQJ) - 2*NE11
GO TO 199
C CASE 1 NOT 8 1 OR N J NOT = 1 OR N
131 5X 8 (T(1+19J1-T(1-10J11/(2*H)
SY = IT‘IOJ+11-T(10J‘111/(2*H1
LAMDAX 3 (LAMDA(I+10J) - LAMDAII‘10J11/12*H1
LAMDAY = (LAMDA(19J+1) ‘ LAMDA(IoJ-1))/(2*H1
E 3 (LAMDAX*SX + LAMDAY*SY + 2*THETA)*(H**2)/LAMDA(IOJI
TBAR 8 (T(1+10J) + T(10J+1) + TII‘IOJ) + T‘IOJ’I) + E1/4o
T‘IOJ) 3 TIIOJ) + BETA*(TBAR - T‘IOJ)’
JZIIOJ) 8 3*9H1*(Y**21*(THETA**2) + ILAMDA119J)**2)*ISX**2 +SY**2)
LAMOA‘IOJ) = LAMDA(10J) ‘ (LAMDA(10J1**NE2*(LAMDA(ICJI‘IO) ‘
5 J2119J1**NE) / (LAMDA‘IQJ)**NEZN1 *((2*NE+11*LAMDA(10J1 “ 2*NE11
199 CONTINUE
IFIIT‘ITERII79170411
17 X 3 00°
0° 61 1=ION 5 DO 61 JBIQN
A1 8 ABSFCABSFC TCIQJ’I‘ABSF(TI(IOJ111
IFIAI‘X161063963
63 X = A1
61 CONTINUE
1F‘X‘E114OOQ16016
16 DO 51 I=10N
DO 51 J=10N

D

157

51 TICIOJ) 3 TIIOJ1
1T = IT + 1 5 GO TO 112
400 pRINT ‘01QITONEQTHETAOPHIOBETA!((T‘IOJ’QI‘1QN10J319N’
401 FORMATI1H00* T‘IOJ’ VALUES FOR ITERATION*01100* NE EQUAE *
$01100/* THETA EQUAL *9F10059* PHI EQUAL *0F10050 * BETA EQUALS*
$9F10050/17E17081)
IFIKODE120004049200
404 pRINT 4059 ((LAMDA(IQJ)OI=1ON10J=19N1
405 FORMAT11H09* LAMDA VALUES* /(7E1708)1
DO 317 1:10N $ 00317 J=19N
317 J2‘10J’ = J2119J) / LAMOA‘1OJ)**2
pRINT 4070 ((J2(10J191=19N)CJ=10N)
407 FORMATCIH00* J2 VALUES * /(7E1708))
COMPUTE NORMALIZED STRESSES
200 DO 299 1:10N
DO 299 J=10N
X = (1'1)*H S Y = (J-1)*H
IFII-1121702070217
207 IF1J‘1120992999209
209 IFIJ-N121192130211
CASE 1:1 J NOT = 1 OR N
211 SX = 0.0
SY 8 (TIIQJ+1)-T(19J‘1)’/(2*H1
GO TO 251
CASE 1‘1 J=N
213 SK 3 000
SY 8 (3*T(10J) “4*T119J“1) +T(IOJ-2))/(2*H)
GO TO 251
217 IFII-N122102190221
219 IE1J‘N122592230225
CASE IzN J=N
223 SX = 0.0 5 SY 3 000
GO TO 251
225 IF(J*1)23592370235
CASE 1=N J31
237 SX 8 (3*T119J) ‘4*T(I‘10J) + T11’2QJ11/(2*H1
SY 3 000
GO TO 251
CASE 1=N J NOT = 1 OR N
235 SX = (3*T110J) -4*T(1‘19J)+T(1*20J))/(2*H)
SY = 000
GO TO 251
221 IFIJ-1122992270229
CASE 1 NOT = 1 OR N J=1
227 SX 2 (T(1+10J1“T(I*10J)1/(2*H)
SY 9 000
GO TO 251
229 1F¢J‘N)23102330231
CASE 1 NOT = 1 OR N J=N
233 SX = 000
SY = (3*T110J) “4*T110J“1)+T(1¢J‘21)/(2*H)
GO TO 251

158

CASE 1 NOT 8 1 OR N J NOT = 1 OR N
231 SX (T(I+19J)-T(I~19J))/(2*H1
SY = (T(IoJ+1)-T(I¢J‘1)1/12*H1
251 SIGMAZIIQJ) = SORTF¢3¢*PHI)*THETA*Y/LAMDA(IOJ1
TAUXZ‘IQJ) = ~SY
TAUYZIIOJ) = SX
299 CONTINUE
IFIKODEISOOOvSIIOoSOOO
5110 ERINT 5110((SIGMAZIIQJ)QI=1oN19J=loN1
511 FORMATI1H09* SIGMA Z STRESSES * /(7E17o8))
PRINT 41301(TAUXZ(I¢J)oI=19N19J=19N)

413 F09MAT<1H0.* TAU xz srnessss * /(7E17.B))
PRINT 415.((TAUYZ(1oJ)oI=1oN).J-loN)
415 FORMATI1HO.* TAU Y2 57955525 ”/17617.811

5000 DO 911 1:10N
DO 911 J=19N
X = (1‘1)*H S Y = (J-1)*H
SIGMAZ‘IQJ) = Y*SIGMAZ(I¢J)
911 TAsz11QJ) = Y*TAUXZ(I¢J) ‘ X*TAUYZ(19J)
TORQ = 000 $ BEMO 3 0.0
TORQZ = O o O
900 D0 901 1:29NO2
DO 901 J=20N02
TORQZ 3 TORQE + T11‘10J‘11+T(1-10J+1)+T(I+19J“1)+T(I+19J+11
5 + 4*(T(1+1¢J)+T(10J+1)+T(1~19J)+T(IoJ-1I) + 16*T119J1
TORO = TORQ + TAUXZCI-I9J-1)+TAUXZ(1+10J~1)+TAUXZ(I‘19J+1)+
S TAUXZII+10J+11 + 40*(TAUXZ(I+10J1+TAUXZIIQJ+1)+TAUX2(I‘19J1
S +TAUXZIIoJ-1II + 160*TAUXZ(IQJ1
901 BEMO = BEMO + SIGMAZII‘IQJ‘II+SIGMAZ(I+19J*11+SIGMAZII~19J+II
S +SIGMAZII+19J+11 + 40*(SIGMAZ(1+1OJ)+SIGMAZ(10J+11+SIGMAZ(1‘19J)
5 + SIGMAZIIoJ-111 + 160*SIGMAZIIQJ)
HSQ = H**2 /9o
TORQ2 = 24.*TOR02*HSQ
TORQ = 120*TORQ*HSQ
BEMO = 16.*HSQ*BEMO
PRINT 9039TOR098EM09TORQ2
903 FORMAT(1HOO* TORQUE EQUALS *0F10080* SENDING MOMENT EQUALS*9
S F10.89* TORQZ EQUALS * oFlOoB)
500 THETA 8 THETA + DTHETA
IF(THETA-THETAM150195010411
501 IT = 1 $ GO TO 112
411 IFIKSTOP1413010413
413 STOP 5 END

159

PROGRAM WARPLV
TYPE REAL LAMDAQJZQLAMDAXQLAMDAYQJelOLAMDAI
DIMENSION T(25925)9LAMDA(25925)9T1(2592519J212562519LAMDA1(25925)
$ QSIGMA2125925)oTAUX2125025)9TAUYZ(25925)QTOR125925)
$ 9J21(2592519TAUXZI(2592519TAUYZI(25.25)OSIGMA21125925)OBEN125925)
1 READ 39PH19 BETAQDTHETAOTHETAMoNoNEOITEROKPSQKSTOP
3 FORMAT(4F10059515 1
READ 49 PHIloTHETAloDTHETAIQKPSI
4 FORMAT (3F1005915)
M = N‘I 5 H = 0.5/M 5 IT
NE2 = 2*NE $ NE2N1= NE2 -
E2 = IoE-4 S KODE = 0
THETA = DTHETA $ PHI 3 PHI**2 $ PHII = PH11**2
DO 5 1:19N
DO 5 J=10N
TIIQJ) = 000 $ T1(19J) = 000 $ TAUY21(19J) = 000
SIGMA21(IQJ) = 000 S LAMDA1(IOJ) 3 000
SIGMAZ(19J) = 0.0 $ TAUXZ(19J) 8 000 $ TAUYZI+OJ1 ‘ 0.0
J2(IQJ) 3 0.0 5 J21(IOJ) = 0.0 $ TAUXZI(I¢J) = 000
5 LAMDA‘IQJ) = 100
PRINT 800NEOPHIODTHETA
80 FORMAT11H19* WARPING FUNCTION FLOW THEORY CALCULATIONS FOR NE EO
$UAL *01100 * PHI EQUAL*9F10¢59* DTHETA EQUAL *9F1005)
112 DC 199 1:19N
DO 199 J=19N
X = (I~1)*H $ Y = (J-1)*H
IFI1“1)11701079117
107 IF(J-1)10901990109
109 IF(J-N)11191139111
CASE I=1 J NOT = 1 OR N
111 SX = T(1+1oJ)/H
SY = (T(19J+11-T(10J*1)I/(2*H)
GO TO 169
: CASE I=1 J=N
113 SX = TII+19J)/H
SY X*DTHETA
GO TO 169
117 IF(I~N)12101199121
119 IF(J*N)12591239125
CASE I=N J=N
123 SX = ‘DTHETA*Y
SY = DTHETA*X
DXTXZ = (30*TAUX21110J) '4.*TAUXZI(I‘19J) + TAUX21(I-29J))/(2*H)
DYTYZ = (30*TAUYZI(IOJ) - 40*TAUY21110J’1) + TAUYZIIIOJ‘211/12*HI
LAMDAX = (30*LAMDA(IcJI-40*LAMDA(1‘19J1+LAMDA(I*29J))/(2*H)
LAMDAY = (30*LAMDA(19J1“40*LAMDA(IoJ-I)+LAMDA(IoJ-2))/(2*H)
E = (DXTXZ + DYTYZ ‘ (LAMDAX*((SX + DTHETA*Y) + TAUXZICIQJII
S + LAMDAY*((SY - DTHETA*X1 + TAUYZIIIQJ))1/LAMDA(10J)1*(H**2)
TBAR = (20*(T(1-10J1-H*DTHETA*Y)+20*(T(IQJ‘I)+H*DTHETA*X1+E)/4o
TIIQJ) = TIIQJI + BETA*(TBAR * T(IQJ))
GO TO 169
125 IFIJH1113591379135
C CASE I=N J=1

= 1 S KP = 1 S E1 8 1.E~e
1 $ NEN1 2 NE - 1

I 1

1|

(1

‘

160

137 SX 3 ‘Y*DTHETA
SY 8 TIIQJ+11/H
GO TO 169
CASE I=N J NOT = 1 OR N
135 SX = ‘DTHETA*Y
SY = (T‘IOJ+1) ‘ TIIoJ-111/12*H)
DXTXZ 8 (30*TAUXZI(IQJ) -4o*TAUX21(I’19J) + TAUXZI(I-20J)1/(2*HI
DYTYZ 3 (TAUY21(IOJ+11 ‘ TAUYZIII.J‘1))/(2*H)
LAMDAX = (30*LAMDA119J)-40*LAMDA(1‘19J1+LAMDA(1’29J11/(2*HI
LAMDAY = (LAMDA(19J+11 ‘ LAMDAIIoJ-III/(2*H)
E = (DXTXZ + DYTYZ * (LAMDAX*((SX + DTHETA*Y) + TAUXZlIIOJI)
S + LAMDAY*((SY - DTHETA*X) + TAUYZ1(I¢JI))/LAMDA(19J)1*(H**2)
TBAR 8 (20*(T(I-I.J)-H*DTHETA*Y)+T(I9J+1)+T(IoJ-1)+E)/4o
T(19J) 3 T(1¢J) + BETA*(TBAR - T(10J))
GO TO 169
121 IF(J*1)12901279129
CASE I NOT = 1 OR N J=1

127 SX = ( T(1+1¢J)- T(I-1oJ))/(2*H)
SY = T(19J+1)/H
GO TO 169
129 1F(J*N)131¢1339131
CASE I NOT = 1 OR N J=N
133 SX = (T(I+19J)-T(I-19J1)/(2*H1
SY = DTHETA*X

DXTXZ = (TAUX21(I+19J) ‘ TAUXZIII-l.J))/(2*H)

DYTYZ = (3o*TAUYZI(IoJ) - 4o*TAUYZI(IoJ‘1) + TAUYZIIIOJ‘ZII/(2*H)
LAMDAX = (LAMDA(I+19J) - LAMDA(I-10J))/(2*H)

LAMDAY = (30*LAMDA(IOJ)-40*LAMDA(IQJ‘II+LAMDAIIOJ‘2))/(2*HI

E = (DXTXZ + DYTYZ - (LAMDAX*((SX + DTHETA*Y) + TAUXZIIIQJI)

$ + LAMDAY*((SY - DTHETA*X) + TAUYZI(I¢J)))/LAMDA(I¢J)1*(H**21
TBAR = (T(I+1oJ)+T(1*1.J)+2-*(T(IoJ-1)+H*DTHETA*X)+E)/4o
TIIqJ) = T(IoJ) + BETA*(TBAR - T(IQJ))

GO TO 169

CASE 1 NOT = 1 OR N J NOT = 1 OR N
131 SX = (T(I+1.J)~T(1-10J))/(2*H)

SY = (T(19J+1) ‘ TIIoJ-111/(2*H1

DXTXZ = (TAUX21(I+19J) - TAUXZIII-leJI/(2*H)

DYTYZ = (TAUYZI(I¢J+1) - TAUYZI(1vJ‘1))/(2*H)

LAMDAX = (LAMDA(I+10J) - LAMDAII-19J11/(2*H)
LAMDAY = (LAMDA(10J+1) - LAMDAII.J-1))/(2*H)
E = (DXTXZ + DYTYZ - (LAMDAX*((SX + DTHETA*Y) + TAUXZ1(IOJ11

$ + LAMDAY*((SY - DTHETA*X) + TAUYZIIIQJ111/LAMDA(19J)1*(H**2)
TBAR (T(1+19J) + T(19J+11 + T(1‘19J) + T119J“1) + EI/4o
T(19J) T‘IQJ) + BETA*(TBAR “ T(19J))
169 J2110J) = (DTHETA*SQRTF(30*PHI)*Y + SIGMA21119J))**2 + (SX+DTHETA*
5 Y + TAUXZIIIOJI)**2 + (SY‘DTHETA*X + TAUYZI‘IOJ))**2
170 IEINEN1117391719173
171 LAMDA(19J) = LAMDAIIQJ) ‘ (LAMDA(IOJ1**NE2*(LAMDA(10J) ‘ 10) -
$ (INE2+Io)/201*(J2(IQJ) - (LAMDA(IOJ)**2)*J21(IQJIII
S / (LAMDAIIQJ)**NE2N1*((20*NE+101*LAMDA(IQJ) - 20*NE) + (20*NE+10)
S *J21(IQJ)*LAMDA(IOJ)1
GO TO 177
173 LAMDA(IQJ) = LAMDAIIQJ) - (LAMDAIIoJ)**NE2*(LAMDA(IoJ) ’ 1.) -

177
175

199

17

63
61

16

51

400
502

504

402
401

405 FORMAT¢1H09*

317 J2(I9J) =

407 FORMAT11H09*

S

5

((NE2+19)/29)*J2(19J)**NEN1*(J2(I9J)
$ / (LAMDA(19J)**NE2N1*((29*NE+19)*LAMDA(IOJ)
S *(J2II9J)**NEN1)*J21(19J)*LAMDA(19J11

- LAMDA1II9J11

IFIABSFILAMDA(I9J)

LAMDA1(I9J) = LAMDA(19J)

GO TO 170

CONTINUE
IF(IT‘ITER)1791794111

X = OoO

DO 61 1319N

DO 61 J=19N

A] = ABSFIAB$F( T(I9J))‘ABSF(T1(I9J)))
IF(A1-X)61963963

X = A1

CONTINUE

IF(X-E114OO916916

DO 51 I=19N

DO 51 J=19N

T1(19J) = T(I9J)

IT = IT + 1 $ GO TO 112
IF(KP-KPS)50294029502

DO 504 I=19N

DO 504 J=19N

J2(I9J) = J2(I9J)/LAMDA(I9J)**2
GO TO 200

161

(LAMDA(I9J)**2)*J21(19J))1
- 29*NE1 + (29*NE+191

* E2119991759175

PRINT 4019IT9NE9THETA9PHI9C(T119J191819NI9J319N1

FORMATI1H09*
91109/*
PRINT 4059

DO 317 1819N
DO 317 J=19N
J2(I9J)
PRINT 4079
J2

COMPUTE NORMALIZED STRESSES
200 D0 299 I=ION

00 299 J=19N
x = (1-1)*H s Y a
IF(I-l)21792079217

(J—1)*H

207 IF(J-1)20992999209
209 IFIJ-N121192139211

CASE I=1 J NOT = 1 OR N
211 SX = T(I+19J1/H
SY = (TII9J+1)-T(I9J-111/(2*H)
GO TO 251
CASE I=1 J=N
213’SX = T(I+19J1/H
SY = X*DTHETA
GO TO 251

217 IF(I-N)22192199221
219 IFIJ-N122592239225

CASE 1=N

J=N

/ LAMDAII9J)**2
((J2119J191=19N)9J=19N)
VALUES * /17E1708)1

TII9J) VALUES FOR ITERATION*91109* NE EQUAL *

THETA EQUAL *9F10959* PHI EQUAL *9F10959/(7E1798))
((LAMDA(I9J)9I=19N19J=19N)

LAMBDA VALUES* /¢7E17981)

162

223 SX = *Y*DTHETA
SY 8 X*DTHETA
GO TO 251
225 IFCJ~1123592379235
CASE I=N J=1

237 SX = -Y*DTHETA
SY = T119J+11/H
GO TO 251

CASE ISN J NOT = 1 OR N
235 SX = -Y*DTHETA
SY = (T(19J+1)-T(19J‘1)1/(2*H)
GO TO 251
221 IEIJ-1122992279229
CASE 1 NOT = 1 OR N J=1
227 SX = (T(1+19J)-T(I‘19J))/(2*H)
SY 8 T119J+11/H
GO TO 251
229 IE(J~N)23192339231
CASE I NOT = 1 OR N J=N
233 SX = (T(1+19JI‘T(I‘19J))/(2*H1
SY = X*DTHETA
GO TO 251
CASE 1 NOT = 1 OR N J NOT = 1 OR N
231 SX = (T(I+19J)“T(I*19J))/(2*H)
SY = (TII9J+11‘T(I9J-1)1/(2*H1

251 SIGMAZ(I9J) = ISQRTF(39*PH1)*DTHETA*Y + SIGMAZI(19J))/LAMDA(I9J)

TAUXZII9J) g (SX+DTHETA*Y + TAUXZI(I9J))/LAMDA(I9J)
TAUYZII9J) = (SY-DTHETA*X + TAUYZI(I9J))/LAMDA(I9J)
299 CONTINUE
IF(KP-KPS)50095059500
505 PRINT 4119((SIGMAZ(I9J)91=19N19J819N1
411 FORMAT(1H09* SIGMA Z STRESSES * /(7E1798))
PRINT 4139((TAUXZ(I9J)9I=19N)9J=19N)

413 FORMATI1H09* TAU XZ STRESSES * /(7E1798)1
PRINT 4159((TAUYZ(I9J)9I=19N)9J319N)
415 FORMAT(1H09* TAU YZ STRESSES */(7E1798)1

DO 911 1=19N
DO 911 J=19N
X = (I-l)*H 5 Y = (J-1)*H
BEN(I9J1 = SIGMAZ(I9J)*Y
911 TOR(19J) = Y*TAUXZ(I9J) - X*TAUYZ(I9J)
TORQ = 0.0 $ BEMO = 090
900 DO 901 1:29N92
DO 901 J=29N92
TORQ = TORQ + TORII‘19J*1)+TOR(1+19J~1)+TOR(I-19J+1)+
$ TOR(I+19J+1) + 49*(TOR11+19J1+TOR(I9J+11+TOR(I*19J)
S +TOR(I9J-1)) + 169*TOR(I9J)
901 BEMO = BEMO + BEN(1'19J‘1)+BEN(I+19J-1)+BEN(I-19J+1I
$ +BEN(I+19J+1) + 4.*(BEN(1+19J)+BEN(I9J+1)+BEN(I-19J)
s + BEN(I9J-1)) + 169*BEN(I9J)
HSQ z H**2 /99
TORQ = 129*T0RQ*HSQ
BEMO = 169*HSQiBEMO
PRINT 9O39TORQ9BEMO

163

903 FORMAT11H09* TORQUE EQUALS *9F10989* SENDING MOMENT EQUALS*9
S F10981

KP=O
500 THETA = THETA + DTHETA
IF(THETA ‘ THETA11501950195001
5001 1F(KODE)50039500595003
5005 KODE = I 5 THETA z THETAI + DTHETAI
KPS = KPSI $ DTHETA = DTHETAI $ PHI 3 PHII $ KP 8 O
5003 IEITHETA-THETAM150195019411
501 IT 2 1 5 KP = KP +1
DO 601 1=19N
DO 601 J=19N
TAUX21119J1 = TAUXZ‘IOJI
TAUYZIIIoJ) = TAUYZ(IOJ)
SIGMAZI‘IOJ) = SIGMAZ‘IOJ)
601 J21(19J) = J2(I9J1
GO TO 112
4111 PRINT 411291T
4112 FORMAT11HO9* FAILED TO CONVERGE AFTER *91109* ITERATIONS*)
411 IFIKSTOP)4I3919413
413 STOP 5 END

Cl

C.

C

164

PROGRAM STPLDT
TYPE REAL LAMDA9J29LAMDAX9LAMDAY9J21OLAMDAI
DIMENSION T(25925)9LAMDA(25925)9T1(2592519J2(25625)
$ 9SIGMAZI25925)9TAUXZI25925)9TAUY2125925) 9LAMDA1125925)
S 9J21(2592519SIGMA21(25925)9T2(2592519BEN(259251
1 READ 39PHI9 BETA9DTHETA9THETAM9N9NE9ITER9KPS9KSTOP 9K0DE
3 FORMATI4F10959615 I

M = N-1 5 H = 0.5/M S IT = 1 $ KP = 1 5 E1 = 19E“6
NE2 = 2*NE S NE2N1= NE2 - I S NEN1 = NE - 1

PHI 3 PHI**2 $ THETA = DTHETA

DO 5 1=19N 5 DO 5 J=19N

SIGMAZI(19J) = 090 5 T2119J) 3 000 S J21‘19J) = 090
T(I9J) = 090 5 T1(I9J) = 090 S J2(I9J) 3 000

SIGMAZ(I9J) = 090 S TAUXZ(I9J1 3 090 S TAUYZ(I9J) 3 090
LAMDA1(I9J) = 190
5 LAMDA¢I9J1 = 190
PRINT 809NE9PHI9DTHETA
80 FORMATIIH19* FLOW THEORY CALCULATIONS FOR NE EQUAL *91109 *
SPHI EQUAL *9F10959 * DTHETA EQUAL *9F1098)
IFIKODE14911294

4 PRINT 6
6 FORMAT (1H09* ITERATION THETA TORQUE MOMENT *)
1.12 DO 199 1=19N $ DO 199 J=19N
X = (1-1)*H S Y = (J-1)*H

IF(I-1>11791079117
107 IFIJ—11109911149109
CASE 1:1 J21
11145x=0.0 s sv=0.0
STAR 2 T2(I+1.J) + T2¢I.J+11 + T211+19J1 + T2(I.J+1) - 4.*T2(I.J1
E = ((LAMDAX*SX + LAMDAY*SY + 2.*DTHETA)*(H**2) - STAR1/LAMDAII9J)
TBAR a (T(I+19J) + T(I9J+1) + TII+I9J1 + TII9J+11 + E1/4.
T(I9J) = T(I9J) + BETA*(TBAR - T(I9J)) 3 GO TO 169
109 IFIJ-N111191139111
CASE I=1 J NOT I 1 OR N
1.11 SX = 0.0
SY = (TII9J+1)-T(I9J-1)1/(2*H1
LAMDAX 3 OoO
LAMDAY = (LAMDA(I9J+1)-LAMDA(I9J-1))/(2*H)
STAR = T2<I+19J1 + T2(19J+1) + T2II+19J) + T2(I9J-I) - 4.*T2(19J)
E = ((LAMDAX*SX + LAMDAY*SY + 2.*DTHETA)*(H**2) - STAR)/LAMDA(19J1
TBAR = (T(I+1.J) + T(I.J+1) + T(1+19J) + T(I9J-11 + E1/49
T(I9J) = T(I.J) + BETA*(TBAR - T(1.J)1 s GO TO 169
CASE I=1 J=N
113 SX = 0.0
SY = (3*T1I9J)-4*T(I9J-1)+T(I.J—2))/(2*H) 3 GO TO 169
117 IF(1-N)12191199121
119 IF(J-N)12591239125
CASE I=N J=N
123 SX = 0.0 5 SY = 0.0 $ GO TO 169
125 IF(J—1)13591379135

(I

165

CASE I=N J21
137 SX = (3*TII9J1-4*T(I‘19J)+T(I-29J)1/(2*H)
SY 2 090 5 GO TO 169
CASE I=N J NOT = 1 OR N
135 SX = (3*T1I9J1-4*T(I-19J)+T(1-29J))/(2*H)
SY = 0.0 $ GO TO 169
121 IF(J-1112991279129
CASE I NOT = 1 OR N J=1
127 SX = ( T(I+19J)- T(I-19J))/(2*H) $ SY = 090
LAMDAX 8 (LAMDA(I+19J) - LAMDA(I~19J))/(2*H)
LAMDAY 3 090
STAR = T2(I+19J) + T2(I9J+1) + T2(1“19J) + T2(I9J+1) - 49*T2119J)
E 8 ((LAMDAX*SX + LAMDAY*SY + 29*DTHETA)*(H**2) - STARI/LAMDAII9J)
TBAR = (T(I+19J) + T(I9J+1) + T(I*19J) + T(I9J+1) + E1/49

T(I9J) 3 T(I9J) + BETA*(TBAR - TII9J11 5 GO TO 169
129 IF(J-N)I3191339131
CASE I NOT = 1 OR N J=N

133 5x 2 (T(I+19J)-T(I-19J))/(2*H)
3v 2 (3*Ttt.d) 34*T1I.0-11 + TII.0-211/(2*H) s 00 T0 169
CASE I NOT = 1 on N J NOT = 1 on N
131 sx = (T(I+19J)-T(1-19J))/(2*H)
SY = (TII9J+1)‘T(19J“1)1/(2*H1
LAMDAX = (LAMDA(I+19J) - LAMDA(1-19J))/(2*H1
LAMDAY = (LAMDA(19J+I) - LAMDA<1.J~1)1/¢2*H1
STAR z T2I1+19JI + T2(I9J+1) + T2<1~1.J1 + T2(19J-l) - 4.*T2(I9J1
E = ((LAMDAX*SX + LAMDAY*SY + 29*DTHETA)*(H**2) - STARI/LAMDA(19J)
TBAR = (T(I+I9J) + T(1.J+1) + T(I-19J) + T(I9J~1) + E1/4.
T(I9J) = T(19J) + BETA*(TBAR - T(I9J))
159 J2(I9J) = (DTHETA*SORTF(3.*PHI)*Y + SIGMAZI(19J)I**2 +
s (LAMDAII9J)**2)*(SX**2 + sv**2)
170 IFINEN1)1739I719173
171 LAMDA(I9J) = LAMDA(I9J) - (LAMDA(I9J)**NE2*(LAMDA119J) - 1.) -

$ ((NE2+19)/29)* (J2(19J) - (LAMDAII9J)**2)*J21(I9J1))
$ / (LAMDAII9J)**N52N1*((29*NE+I.)*LAMDA(I9J) ‘ 29*NE) +129*NE+191
$ *J21(I9J)*LAMDA(I9J)) 5 GO TO 177

173 LAMDAI19J) = LAMDA(19J) - (LAMDAII9J)**NE2*(LAMDA(I9J1 ’ 19) -

$ ((NE2+Io)/20)*J2(19J)**NEN1*(J2(I9J) ‘ (LAMDA(I9J)**2)*J21(I9J1I)
S / (LAMDA(I9J)**NE2N1*((29*NE+19)*LAMDA(I9J1 ‘ 29*NE) +120*NE+10)
S *(JZIIoJ)**NEN1)*J21(I9J)*LAMDA(I9J)1
177 IF(ABSF(LAMDA(19J) ‘ LAMDAIII9J11“ E1119991759175
175 LAMDA1(19J) = LAMDAII9J) 5 GO TO 170
199 CONTINUE
1E(IT-ITER11791794111
17 X = 000
D0 61 1=19N 5 DO 61 J=19N
A1 = ABSFIABSFI T(I9J))‘ABSF(T1(19J)1)
IFIAI-X161963963
63 X = A1
61 CONTINUE
IFIX¢E11400916916
16 DO 51 1=19N 5 DO 51 J=19N
51 TIIIOJI = T110J)

166

IT = IT + 1 5 GO TO 112
400 IFIKP-KP5150294029502
502 D0 504 1=19N $ DO 504 J=19N S Y 8 (J-1)*H
J2(I9J) = J2II9J1/LAMDAII9J)**2
504 SIGMAZ(19J) = (SQRTF(3.*PHI)*DTHETA*Y + SIGMAZI(I9J))/LAMDA(I9J)
GO TO 500
402 KP = 0
IF(KODE)4O6940009406
4000 PRINT 4019IT9NE9THETA9PHI9((T(I9J)9I=19N)9J=19N)
401 FORMAT(1H09* T(I9J) VALUES FOR ITERAT10N*91109* NE EQUAL *
$91109/* THETA EQUAL *9F10959* PHI EQUAL *9F10.59/(7EI798))
PRINT 4059 ((LAMDAII9J)9I=19N19J=19N1
405 FORMATI1H09* LAMBDA VALUES* /(7E1798)1

406 D0 317 1=19N $ DO 317 J=19N
317 J2(I9J) = J2<I9J1 / LAMDAII9J)**2
IF(KODE)20094089200

408 PRINT 4079 ((J2(I9J191=19N)9J319N)
407 FORMAT(1H09* J2 VALUES * /(7E1708)1
COMPUTE NORMALIZED STRESSES
200 D0 299 1=19N 5 DO 299 J=19N
X = (1-1)*H 5 Y = (J-1)*H
IE(I-1)21792079217
207 IE(J-1)20992999209
209 IEIJ-N)21192139211
CASE 1:1 J NOT = 1 OR N
211 SY = (TII9J+1)~T(I9J‘1)1/(2*H) $ SX = 000
GO TO 251
CASE 181 J=N
213 SY = (3*TI19J) -4*T(I9J-1) +T(19J-2))/(2*H) $ SX
GO TO 251
217 IFII~N122192199221
219 IF(J“N)22592239225
CASE I=N J=N
223 SX 3 000 5 SY = 000 $ GO TO 251
225 IEIJ“1)23592379235
CASE 13N J31
237 SX = 13*T119J1 ‘4*TII“19J1 + TII-2OJII/(2*H1 5 SY 3 000
00 T0 251
CASE IBN J NOT = 1 OR N
235 SX = (3*T1I9J) -4*T(1-19J)+T(I~29J))/(2*H) 5 SY = 000
GO TO 251
221 IF(J‘1)22992279229
CASE I NOT = 1 OR N J=1
227 SX = (T(I+19J)-T(I-19J))/(2*H) 5 SY 8 000
GO TO 251
229 IF(J*N)23192339231
CASE 1 NOT = 1 OR N J=N
233 SY = (3*T119J) -4*T119J‘1)+T(I9J-2))/(2*H) S SX = 000
GO TO 251
CASE I NOT = 1 OR N J NOT = 1 OR N
231 SX ITII+19J1‘TII‘19J11/(2*H1
SY (T(19J+1)~T(I9J-1)1/(2*H1

11
O
o
O

1111

167

251 SHGMAZIIOJ) = (SORTF139*PHI)*DTHETA*Y + SIGMAZIII9J11/LAMDAII9J)
TAUXZII9J1 3 allSY
TAUYZII9J) = SX
299 CONTINUE
IFIKODE13OO99O993OO
909 PRINT 5119((SIGMAZIIOJ)9I=I9N19J=19N1
511 FORMAT11H09* SIGMA Z STRESSES * /(7E1708)1
PRINT 4139((TAUXZ(I9J)9I=19N)9J‘19N)

413 FORMAT11H09* TAU XZ STRESSES * /¢7E1798))
PRINT 4159((TAUYZII9J19I=19N19J819N1
415 FORMAT(1H09* TAU YZ STRESSES */(7E1798))
300 DO 911 1=19N $ DO 911 J=19N
X = (1-1)*H 3 Y = (J-1)*H

BENII9J) 3 Y*SIGMAZ(I9J1
911 TAUXZ‘I9J) = Y*TAUXZ(I9J) ‘ X*TAUYZIIQJ)
TORQ = 090 $ BEMO : 090 S TORO2 = 090
900 D0 901 1:29N92
DO 901 J=29N92
TORO2 = TORQ2 + T(I-19J-1)+T(1-19J+1I+TII+19J*1)+T(1+19J+1I
5 + 4*(TI1+19JI+TI19J+11+T(I-19J)+T(I9J-1I) + 16*TIIOJ)
TORQ = TORQ + TAUXZII‘I9J-1)+TAUXZ(1+19J‘11+TAUXZ(I‘19J+11+
$ TAUXZ(I+19J+1) + 49*(TAUXZII+19JI+TAUXZII9J+1)+TAUXZ(I’10J)
5 +TAUXZII9J‘111 + 160*TAUXZIIOJ1
901 BEMO = BEMO + BENII’IOJ’I) + BENII+19J-1) + BENII-19J+1)
$ + BENII+I9J+11 + 40*(BEN(1+19J1 + BENII9J+1) + BENII‘IOJI
S + BENII9J-11) + 160*BENIIOJ1
HSO = H**2 /9o
TORQ2 = 240*TORO2*HSO
TORQ = 129*TORQ*HSO
BEMO = 169*HSQ*BEMO
IF(KODE19O4990299O4
902 PRINT 9039TORO9BEMO9TOR02

903 FORMAT(1H09* TORQUE EQUALS *9F10989* SENDING MOMENT EQUALS*9
$ Flo-89* TORQZ EQUALS * 9F1098)
GO TO 500

904 PRINT 9069IT9THETA9TOR029BEMO 9TORQ
906 FORMAT(11294F1296)
500 THETA = THETA + DTHETA
IFITHETA-THETAM150195019411
501 IT = 1 $ KP = KP +1
DO 601 1=19N S DO 601 J=19N
SIGMAZIII9J) = SIGMAZ(I9J)
T2(I9J) = T(I9J)

J21(19J) = J2(I9J)
601 CONTINUE 5 GO TO 112
4111 PRINT 411291T
4112 FORMATIIHO9* FAILED TO CONVERGE AFTER *91109 * ITERATIONS*1

411 IF(KSTOP)413919413
413 STOP 5 END

1

168

PROGRAM POLARR
TYPE REAL LAMDA
DIMENSION T(25925)9T1(25925) 9 SIGMAZI25925)9TAUY212592510
S TAUXZI25925) 9TOR(25925) 9 BENI25925)QR(251 9 THA125)
1 READ 39PSI9BETA9N9ITER9KSTOP
3 FORMAT(2F109593151
PRINT 59PSI9BETA
5 FORMAT(*1*9* TIIoJ) VALUES FOR CIRCULAR BAR OF UNIT RADIUS WHERE
SPSI EQUAL*9 F10959* BETA EOUAL*9F1005)
PSI = PSI**2 S M = N-l S MI 3 N+1
H = Io/M S DTHA = 195708 / M S 1T8 1 S E1 = 19E*6
DO 7 1=19N
THAII) = (1‘1)*DTHA
7 RII) 3 (I-I)*H
DO 9 I=I9M1 S DO 9 J=19N
SIGMAZ(I9J1 = 000 S TAUYZ(I9J) = 0.0 S TAUXZCI9J) a 0.0
T1(10J) = 000
9 T(I9J) = 000
10 DO 11 1:29N S DO 11 J=29M
A = (39*PSI*(R(I)**5)*SINF(THA(J))**2 + RIII*(((T(19J+1)‘T(IOJ‘III
S/(29*DTHA)1 - R(11**21**2) * (DTHA**2)
B = (20*RII) * (IT(I+19J)' TII‘l9J))/(29*H))*(R(I)**2 * (ITIIOJ+1)
S - TIIQJ-III/(29*DTHA)))1*IO¢25*DTHA*H1
C 8 (30*PSI* (RII)**3) * (SINFITHAIJ))**21 + R11) * ((TII+19J1
S - TII‘19JI) /(29*H11**2)*1H**21
D =(IRII)**2) * ((T(I+19J) ‘ TII‘IOJ))/(20*H))**2)*(005*H*DTHA**2)
E = (“29*((T(I+19J1 ‘ TII-19J)1/(20*H1)*(R(11**2 - ((T(10J+1) ‘
S TII9J~1I)/(2o*DTHA1))-3o*PSI*(R(I)**3)*SINF(THA(J11 *
S COSFITHAIJIII * (005*DTHA*H**21
F = (39*PSI* (R(1)**5) * SINFITHA(J)) * COSFITHA(J111
S *(H**2*DTHA**2)
TBAR 3 I(A+D)*T(I+19J) + (A-D)*T(1-19J) + (C+E)*T(IOJ+1)
S + (C‘E1*T(I9J‘1) + B*(T(I+19J+1) “ T(I~19J+1) + TCI‘19J‘1)
S - TII+19J-1)I + F) / (29*(A+C))
11 T(19J) = TII9J) + BETA*(TBAR * T(19J1)
DO 13 J=19N
13 T(M19J) = TIMoJ) S X = 000
DO 301 1=19N S DO 301 J=19N
A1 = ABSF(ABSFIT(I9J11-ABSF(T1(I9J))1
IFIA1‘X130193059305
305 X=A1
301 CONTINUE
IFIX‘E11200951951
51 DO 53 1:29N S DO 53 J=29N
53 TIIIQJ) = TIIQJ)
IT 8 IT + 1 S IFIIT‘ITER110940094OO
200 PRINT 2019IT9BETA9PSI9(1T119J191319N19J319N)
201 FORMAT (1H09* T(19J) VALUES FOR ITERATION * 9 1109//* RELAXATION
SFACTOR *9F10059* PSI VALUE *9F10959//(7E1708 )1
GO TO 112
400 PRINT 40191T9BETA9PSI
401 FORMATI1H19* FAILED TO CONVERGE AFTER * 91109*ITERATIONS*

I I

169

S ///* BETA EOUAL*9F10.59* PSI EQUAL*9F10951
GO TO 1105
112 D0 111 1:29N S DO 111 J=19N
IF(J-1110391019103
CASE J=1
101 STHA = TII9J+1)/DTHA
SR 8 (T(I+19J) - T(I‘19J)1/(2o*H) S GO TO 109
103 IF(J-N)10591079105
CASE J=N
107 STHA = ‘TII9J-11/DTHA
SR 3 (T(I+19J) - T11‘19J11/(29*H) S GO TO 109
CASE J NOT = 1 OR N
105 STHA = (T(I9J+1) ~ TII9J-1))/(29*DTHA)
SR = (T(I+19J) - TII-19J11/(29*H)
109 LAMDA = SORTFI39*PSI *(R(I)**2)*(SINF(THA(J)1**2) + SR**2 +
S (STHA/RII))**2 + R(II**2 - 29*STHA)
SIGMAZI19J) = (SQRTFI3.*PSI)*R(I)*SINF(THA(J))I/LAMDA
TAUX2119J) = (SR*COSF(THA(J)) “ STHA*SINF(THA(J))/R(9) + R(I)
S *SINF(THA(J)))/LAMDA
TAUYZII9J) = (SR*SINF(THA(J)) + STHA*COSF(THA(J))/R(I) - R‘I)
S * COSFITHAIJIII/LAMDA
111 CONTINUE
PRINT 4119((SIGMAZII9J)91319NI9J=19N1
411 FORMAT(1H09* SIGMA Z STRESSES*//(7E179811
PRINT 4139((TAUXZII9JI9I=19N)9J=19N1
413 FORMATIIHO9* TAU XZ STRESSES* // (7E179811
PRINT 4159((TAUYZII9J)91=19N)9J=I9N)
415 FORMATI1H09* TAU YZ STRESSES* //(7E1798))
DO 113 1=19N
DO 113 J=19N
BENII9J) = SIGMAZ(I9J)*(R(I)**2)*SINF(THA(J)1
113 TORIIOJ) (RII)**2)*(SINF(THA(J)1*TAUXZII9J) - COSFITHAIJII
S *TAUYZ(I9J))
TORQ = 0.0 5 DEMO 2 000
D0 115 I=29N92
DO 115 J=29N92
TORQ = TORQ + TORII+19J+11 + TOR(I~19J+1) + TOR(I-19J“1) +
s TOR(I+1.J-1) + 4.*(TOR(1+1.J) + TORII.J+1) + TOR(1-19J) +
S TORII9J-11) + 169*TOR(I9J)
115 BEMO = BEMO + BENII+19J+II + BENII‘19J+I) + BENII-IoJ‘I) +
S BENII+19J-1) + 49*(BEN(I+19J) + BENII9J+11 + BENII‘19J) +
S BENII9J-111 + 169*BENII9J)
H50 3 (DTHA*H)/9o
TORQ = 4o*TORQ*HSQ/(29094)
BEMO = 30*BEMO*HSQ
PRINT 9039TORO9BEMO

903 FORMAT11H09* TORQUE EQUALS *9F10989* BENDING MOMENT EQUALS

S QFIOOBI
1105 IFIKSTOP12159 1 9215
215 STOP S END

*

(I

(1

170

PROGRAM POLARO
TYPE REAL LAMDA 0 J2
DIMENSION T125925)9T112592519J2125925)9SIGMAZIZ5925)OTAUY2125925)9
S TAUXZ‘25925) 9TORI25925) 9 BENI25925)9R(25) 9 THA125)
S 9 LAMDA125925)
1 READ 39PSI9BETA9DTHETA9THETAM9N9NE9ITER9KSTOP9KODE
3 FORMATI4F109595151
IFIKODE129492
PRINT 69PSI9BETA9NE S GO TO B
FORMATIIHI9§ DEFORMATION THEORY - CIRCULAR BAR FOR PSI EQUAL*O
SF10959* BETA EQUAL*9F10959* NE EQUAL*9159//* ITERATION THETA
S TORQUE MOMENT*)
4 PRINT 59PSI9BETA9NE
5 FORMAT(*1*9* T(I9J) VALUES FOR CIRCULAR BAR OF UNIT RADIUS WHERE
SPSI EQUAL*9 F10.5.* BETA EQUAL*9F1095 9* NE EQUALS*9151
8 NE2 8 2*NE S NE2N1 = 2*NE - 1
P51 8 PSI**2 S M 8 N‘l S M1 8 N+1
H = Io/M S DTHA 3 1057OB / M

010

IT 8 1 S E1 8 10E‘6 S THETA 8 DTHETA
DO 7 I819M1
THAIII = (1-1)*DTHA

7 R11) = (1811*H
DO 9 I819M1 S 00 9 J819N
SIGMAZII9J) = 000 S TAUYZII9J1 = 000 S TAUXZIIQJ) 3 OoO
T1(19J1 3 090 S LAMDA(19J1 = 100

9 TIIoJ’ = 090
10 DO 11 1:29N S DO 11 J=19N

IF(I-N)6O396OI96O3
601 IFIJ-1160796059607
CASE I8N J81
4505 SR 8 OoO
STHA 8 T119J+1)/DTHA S GO TO 625
607 IFIJ-N161196O99611
CASE I8N J=N
609 SR 8 OoO
STHA 8 *TII9J-11/DTHA S GO TO 625
CASE I8N J NOT 8 1 OR N
611 SR 8 090
STHA 8 (T(I9J+1) - T(I9J-l))/(29*DTHA1
A 8 (RII)**2)*(DTHA**2) S C 8 H**2
E (((LAMDAII9J+1) - LAMDA(I9J-111/(29*DTHA)1*(STHA/(RIII**2) -
S IoOI/LAMDAII9J)1*(RII)**2)*(H**21*(DTHA**21
TBAR 8 (29*A*T(1*19J) + C*(T(I9J+1) + T(I9J-1)) - E)/(29*(A+C))
TII9J) 8 TII9J) + BETA*(TBAR - T(19J))
GO TO 625
603 IFIJ-1161596139615
CASE I NOT 8 N J81
613 SR 8 (T(I+19J) - TII‘19J11/(29*H)
STHA 8 T(I9J+1)/DTHA S GO TO 625
615 IFIJ-N161996179619
CASE 1 NOT 8 1 OR N J=N
617 SR 8 (T(I+19J) - TII‘I9J))/(2.*H)
STHA 8 ‘T(I9J-1)/DTHA S GO TO 625

171

CASE 1 NOT 8 N J NOT 8 1 OR N

619 A 8 (R(11**2)*(DTHA**2)

D a 0.5*R(11*H*(DTHA**2)
C 8 H**2
E 8 ((((LAMDA(I+19J) - LAMDA(1819J))/(2.*H))*((T(1+19J) 8

S T(1819J))/(29* H ))+((LAMDA(19J+1) 8 LAMDA(I9J81))/(2o*DTHA1)*
$ ((T(I9J+1) 8 T(19J81))/(29*DTHA))/(R(I)**2)8(LAMDA(I9J+11 8
SLAMDA(I9J81))/(29*DTHA)1/LAMDA(19J))*((RCI)**2)*(H**2)*(DTHA**211

TBAR 8 ((A+D)*T(1+19J) + (A-D)*T(I819J) +C*(T(19J+1) + T(I9J81)1 -

S E)/(29*(A+C),

T119J1 8 T(19J) + BETA*(TBAR - T‘I9JI)
SR 8 (T(I+19J) 8 T(I819J))/(29*H1
STHA 8 (T(I9J+1) 8 T(I9J81))/(29*DTHA)

625 J2(19J) 8 (THETA**2)*(39*PSI*R(1)**2*SINF(THA(J))**2 + (SR’COSF

11

305
301

51
53
200

202
201

402
4005
406

4006
404

400
401

112

101

s (THA(J)) 8 STHA*SINF(THA(J))/R(I) + R(1)*SINF(THA(J)))**2 + (SR*
5 SINF¢THA(J)) + STHA*COSF(THA(J))/R(I) 8 R11)*COSF(THA(J)))**2)

LAMDA(I9J) 8 LAMDA(19J) 8 (LAMDA(I9J)**NE2*(LAMDA(I9J1819) 8

5 J2(IqJ)**NE) / (LAMDA(19J)**NE2N1 *((2*NE+11*LAMDA(19J) 8 2*NE))

X 8 00°

88 476 IM6P9

A1 8 ABSECABSF‘TKI9J118ABSFIT1(19J)1)
1E(A18X)30193059305

X=A1

CONTINUE

IF(X8E112OO951951

DO 53 I‘Z9N

DO 53 J829M

T1(I9J) a TII9J)

IT 8 IT + 1 S IF(IT‘ITER)IO9AOO9400
IECKODE14005920294005

PRINT 2019IToBETA9PSI9THETA9((T(I9J)9I=19N)9J=19N)
FORMAT (1HO9* T(I9J) VALUES FOR ITERATION * 9 1109//* RELAXATION

SFACTOR *9F10959* PSI VALUE *9F10959* THETA EQUALS *9F10959/I
S (7E1798 ))

PRINT 4029((LAMDACI9J)9I819N)9J819N)

FORMATI1H09* LAMBDA VALUES *//(7E1798))

DO 406 I819N S DO 406 J819N

J2(19J) 8 J2(19J)/LAMDA(I9J)**2

IE(KODE)112940069112

PRINT 4049((J2(I9J)91819N)9J819N1

FORMAT(*O*9* J2 VALUES * //(7E1798))

GO TO 112

PRINT 4Ol9IT9BETA9PSI

FORMAT¢1H19* FAILED TO CONVERGE AFTER * 91109*ITERATIONS*

S ///* BETA EOUAL*9F10.59* PSI EQUAL*9F10.5)

GO TO 1105

DO 111 I829N

DO 111 J819N
IF(J-1)10391019103

CASE J81

STHA 8 T(I9J+1)/DTHA
SR 8 (T(I+19J) 8 T(1-19J))/(2.*H)
GO TO 109

(I

172

103 IE(J8N110591079105

CASE J=N

107 STHA 8 8T(I9J-1)/DTHA
SR 8 (T(I+19J) 8 T(I819J))/(2.*H)
GO TO 109

CASE J NOT 8 1 OR N
105 STHA 8 (T(I9J+1) 8 TII9J8111/(29*DTHA)
SR 8 (T(I+19J1 8 T(I819J))/(29*H1
109 SIGMAZII9J) =((SORTFC39*PSI)*Q(I1*SINF(THA(J)1’/LAMOA(I9J1)*THETA
TAUXZ(I9J) 8((SR*COSE(THA(J)) 8 STHA*SINF(THA(J))/R(I) + R11)
S *SINF(THA(J)1)/LAMDA(I9J))*THETA
TAUYZII9J) 8((SR*SINF(THA(J)) + STHA*COSE(THA(J))/R(I) 8 R‘I)
S * COSE(THA(J)))/LAMDA(19J))*THETA
111 CONTINUE
IE(KODE)4OIO9401194010
‘4011 ERINT 4119‘(SIGMAZ(I9J)91819N19J819N)
411 FORMAT(1HO9* SIGMA Z STRESSES*//(7E17¢8})
PRINT 4139((TAUXZ(I9J)9I819N)9J819N)
1413 FORMAT11H09* TAU X2 STRESSES? // (7E179811
PRINT 4159((TAUYZII9JI9I819N)9J819N1
415 FORMATIIHO9* TAU YZ STRESSES* //(7E1708))
‘4010 DO 113 1819N
DO 113 J819N
BEN(I9J) 8 SIGMAZCI9J)*(R(I)**2)*SINE(THA(JII
113 TOR(I9J) 8 (R(I)**2)*(SINF(THA(J)1*TAUXZII9J) - COSE(THA(4))
S *TAUYZ(I9J))
TORQ 8 090 S BEMO 8 090
DO 115 I829N92
DO 115 J829N92
TORQ 8 TORQ + TOR(I+19J+1) + TOR(I819J+1) + TOR(1819J81) +
S TOR(I+19J81) + 49*(TOR(I+19J) + TORII9J+1) + TOR(I819J) +
S TOR(19J81)) + 169*TORIIqJ)
115 BEMO 8 BEMO + BEN(I+19J+1) + BEN(I819J+1) + BEN(I819J-1) +
S BEN(I+19J81) + 49*(BEN(I+19J) + BEN(I9J+1) + BEN(I819J) +
S BENII9J81)) + 16.*BEN(I9J)
H50 8 (DTHA*H)/9o
TORQ 8 49*TORQ*HSQ/(29094)
BEMO 8 39*BEMO*HSO
IF(KODE)9O4990299O4
904 PRINT 9089IT9THETA9TORQ9BEMO
908 FORMATI2X9I1092X9F100692X9E109692X9F1096)
GO TO 500
902 PRINT 9039TORQ9BEMO -
903 FORMAT(1HO9* TORQUE EQUALS *9F10989* BENDING MOMENT EQUALS
S 9E1098)
500 THETA 8 THETA + DTHETA
IFITHETA 8 THETAM)501950191105
501 IT 8 l S GO TO 10
1105 IE(KSTOP)2159 1 9215
215 STOP S END

{-

APPENDIX v
TORQUE, MOMENT, AND STRESS OUTPUTS

Table 5. Effect of grid dimension on Piechnik e uation

warping function and stresses for u = 1/ 3
(lower values for finer mesh)
Warping function: t x 10-2
i l 2 3 4 5 6 7
j 1 3 5 7 9 11 13
1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
2 0.0000 0.5028 0.7515 0.6005 0.2114 -0.3191 -0.9613
3 0.0000 0.5036 0.7243 0.5591 0.1647 -0.3681 -1.0111
3 0.0000 0.9032 1.4392 1.3673 0.7540 -0.2326 -1.4970
5 0.0000 0.8853 1.3786 1.2682 0.6365 -0.3565 -1.6220
4 0.0000 1.3055 2.1900 2.3449 1.7040 0.3830 -1.4691
7 0.0000 1.2702 2.0980 2.1946 1.5140 0.1777 -1.6758
5 0.0000 1.7619 3.0814 3.5996 3.1341 1.6625 -0.7165
9 0.0000 1.7139 2.9621 3.4021 2.8740 1.3722 -1.0060
6 0.0000 2.3023 4.1675 5.2016 5.1393 3.7911 0.9913
11 0.0000 2.2485 4.0263 4.9563 4.8026 3.4076 0.6383
7 0.0000 2.9435 5.4874 7.2137 7.8358 7.0955 4.0434
13 0.0000 2.8951 5.3334 6.9167 7.3898 6.5419 3.7365
Normal stress: ozz/ngik
i 1 2 3 4 5 6 7
j 1 3 5 7 9 11 13
1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
2 0.5016 0.5353 0.5702 0.4287 0.2765 0.1898 0.1398
3 0.4914 0.5242 0.5656 0.4167 0.2687 0.1856 0.1374
3 0.5180 0.5462 0.6206 0.6552 0.5487 0.4052 0.2996
5 0.5126 0.5430 0.6238 0.6515 0.5362 0.3954 0.2937
4 0.5235 0.5467 0.6158 0.7077 0.7329 0.6291 0.4828
7 0.5209 0.5460 0.6198 0.7114 0.7256 0.6150 0.4726
5 0.5217 0.5404 0.5983 0.6941 0.7980 0.8165 0.6870
9 0.5207 0.5409 0.6024 0.6998 0.7981 0.8037 0.6732
6 0.5151 0.5297 0.5760 0.6584 0.7783 0.9092 0.8890
11 0.5147 0.5307 0.5805 0.6550 0.7819 0.8999 0.8792
7 0.5054 0.5162 0.5513 0.6148 0.7112 0.8746 0.9943
13 0.5052 0.5171 0.5565 0.6235 0.7202 0.8527 0.9988

173

174

Table 5 (continued)

Shear stress: lo /k|
xz

 

i 1 2 3 4 5 6 7
j 1 3 5 7 9 11 13
1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
2 0.8651 0.8242 0.6105 0.2621 0.0935 0.0296 0.0000
3 0.8709 0.8256 0.5771 0.2381 0.0871 0.0280 0.0000
3 0.8553 0.8296 0.7247 0.4939 0.2328 0.0769 0.0000
5 0.8586 0.8303 0.7126 0.4676 0.2163 0.0721 -0.0000
4 0.8520 0.8344 0.7698 0.6254 0.3880 0.1501 0.0000
7 0.8536 0.8344 0.7638 0.6093 0.3677 0.1403 0.0000
5 0.8531 0.8405 0.7966 0.7010 0.5200 0.2507 0.0000
9 0.8538 0.8401 0.7927 0.6915 0.5046 0.2382 0.0000
6 0.8571 0.8480 0.8169 0.7509 0.6205 0.3663 0.0000
11 0.8574 0.8474 0.8135 0.7441 0.6118 0.3641 0.0000
7 0.8629 0.8565 0.8343 0.7884 0.7050 0.4768 0.0000
13 0.8630 0.8559 0.8308 0.7817 0.6932 0.5199 0.0000
Shear stress: Oyz/k|
i l 2 3 4 5 6 7
j 1 3 5 7 9 11 13
1 0.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1 0.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
2 0.0000 0.1870 0.5498 0.8646 0.9564 0.9814 0.9902
3 0.0000 0.2088 0.5891 0.8773 0.9593 0.9822 0.9905
3 0.0000 0.1153 0.2994 0.5716 0.8030 0.9110 0.9541
5 0.0000 0.1254 0.3210 0.5974 0.8159 0.9157 0.9559
4 0.0000 0.0696 0.1679 0.3287 0.5588 0.7627 0.8757
7 0.0000 0.0748 0.1801 0.3501 0.5816 0.7759 0.8813
5 0.0000 0.0382 0.0862 0.1637 0.3047 0.5200 0.7266
9 0.0000 0.0406 0.0935 0.1791 0.3291 0.5453 0.7395
6 0.0000 0.0158 0.0308 0.0524 0.0957 0.1980 0.4577
11 0.0000 0.0162 0.0348 0.0638 0.1195 0.2399 0.4765
7 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
13 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

175

Table 6. Comparison between stresses computed from
Piechnik e uation and HilluHandelman equation for
p = 1‘V5 (lower values due to Steele)

Normal stress: ozz/‘V3 k

i 2 4 6 8 10
j 1 2 3 4 5
2 0.497 0.496 0.246 0.137 0.091
1 0.348 0.452 0.277 0.139 0.077
4 0.526 0.622 0.629 0.444 0.296
2 0.507 0.626 0.612 0.455 0.628
6 0.531 0.612 0.738 0.729 0.543
3 0.554 0.633 0.727 0.712 0.487
8 0.524 0.586 0.712 0.865 0.819
4 0.540 0.601 0.735 0.844 0.815
10 0.510 0.554 0.648 0.802 1.000
5 0.530 0.571 0.668 0.837 1.000
Shear stress: '0 /k‘
xz
i 2 4 6 8 10
j 1 2 3 4 5
2 0.856 0.492 0.110 0.027 0.000
1 0.907 0.553 0.097 0.027 0.000
4 0.847 0.718 0.385 0.112 0.000
2 0.849 0.664 0.339 0.116 0.000
6 0.846 0.778 0.583 0.252 0.000
3 0.829 0.746 0.532 0.229 0.000
8 0.852 0.809 0.693 0.424 0.000
4 0.841 0.788 0.651 0.381 0.000
10 0.860 0.832 0.761 0.598 0.000
5 0.848 0.821 0.744 0.548 0.000
Shear stress: 10 /kl
yz
i 2 4 6 8 10
j 1 2 3 4 5
2 0.137 0.715 0.963 0.990 0.996
1 0.238 0.700 0.956 0.990 0.997
4 0.076 0.311 0.675 0.889 0.955
2 0.103 0.370 0.687 0.881 0.963
6 0.039 0.142 0.339 0.636 0.840
3 0.053 0.190 0.407 0.651 0.873
8 0.015 0.050 0.113 0.269 0.574
4 0.020 0.073 0.170 0.349 0.580
10 0.000 0.000 0.000 0.000 0.000
5 0.000 0.000 0.000 0.000 0.000

176

0N00.H
5000.H
5800.H
0000.H
H080.H
8H80.H
00mo.H
8Hm0.a
00N0.H
NONO.H
H8H0.H
8500.H
~000.H
mm00.0
0m00.0
0850.0
mm00.0
NH00.0
m5m0.0
00N0.0
0000.0
0050.0
0080.0
0000.0
0005.0
0000.0
50H0.0
00H0.0
m0H8.0
mmam.0
NO0N.0
H80H.0
0000.0

:
u

NH n C

NH50.H
H500.H
5N00.H
H000.H
N000.H
H080.H
0N80.H
0000.H
5000.H
N8N0.H
m5H0.H
0000.H
0H00.H
0000.0
0000.0
5N50.0
0000.0
0580.0
0Nm0.0
00H0.0
0800.0
0000.0
05m0.0
0000.0
0005.0
NN00.0
0mH0.0
00H0.0
m0H8.0
mmam.0
N0om.0
H80H.0
0000.0

085H.H
0HOH.H
55NH.H
0HOH.H
0m50.H
Nm80.H
8000.H
0H50.0
H0m0.0
H000.0
8NN0.0
0005.0
0000.0
0H00.0
0000.0
050N.0

0000.0

HH5N.H
H8MN.H
~00H.H
N80H.H
50HH.H
0800.H
NOH0.H
0m00.0
N800.0
0080.0
5055.0
0H00.0
0000.0
mm08.0
5m0m.0
000H.0

0000.0

00.0 n 1 spa: coauomm mumsvm m mo mfimhamcm Snows»
ICOHmeuommU coawucsm10CH00m3 Mom mmsvuou UmNHHmEuoz

00.0
05.5
00.5
00.5
00.5
05.0
00.0
0N.0
00.0
05.0
00.0
00.0
00.0
05.8
00.8
0N.8
00.8
05.m
00.0
00.0
00.m
05.0
00.0
00.0
00.N
05.H
00.H
0N.H
00.H
05.0
00.0
0N.0
00.0

00

.M5 GHQMB

177

00H0.0 0000.0 0000.0 0000.0 0005.0 00H0.H 0000.0
05H0.0 0000.0 00H0.0
80H0.0 0000.0 0H00.0 0000.0 0000.0 0000.0 0000.0
00H0.0 0000.0 0508.0
HHH0.0 H550.0 50H0.0 0800.0 0000.0 0000.0 0000.0
0000.0 .8050.0 0000.0
8000.0 0000.0 0000.0 H850.0 0080.0 5000.0 0050.8
0000.0 0000.0 0500.8
HH00.0 0000.0 0H00.0 0000.0 00H0.0 0500.0 0000.8
H000.0 H000.0 00H0.8
0800.0 0000.0 0000.0 0000.0 5000.0 8080.0 000H.8
0000.0 0080.0 0500.0
0500.0 H000.0 0000.0 0000.0 5H00.0 0000.0 0005.0
0000.0 H000.0 0000.0
8050.0 0800.0 H050.0 00H0.0 0H00.0 0005.0 0050.0
0050.0 HOH0.0 050H.0
0000.0 5800.0 0500.0 8805.0 5008.0 00H5.0 0000.0
0000.0 0005.0 00H0.0
8080.0 5055.0 0500.0 0005.0 5008.0 0000.0 0000.0
0000.0 H005.0 0508.0
0000.0 0H85.0 50H0.0 0005.0 0008.0 5000.0 0000.0
0000.0 50H5.0 0000.0
5H08.0 5000.0 H058.0 0000.0 0050.0 0880.0 0050.0
0008.0 0080.0 0500.H
HOH8.0 8500.0 0008.0 0000.0 0000.0 0058.0 0000.0
H050.0 H500.0 0000.H
0000.0 0008.0 5000.0 0008.0 H000.0 H000.0 000H.H
0050.0 0000.0 0500.0
00H0.0 00Hm.0 00H0.0 00H0.0 800H.0 0050.0 0005.0
800H.0 0800.0 0000.0
000H.0 HO0H.0 000H.0 HO0H.0 080H.0 500H.0 0050.0
H800.0 0050.0 050H.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
GE cu :6 c0 as aw .L$w
0 n a 0 u c H n c Q

00.0 n 1 sud: coavumm mumsvm m 00 mammamcm huomzuICOHpmeuommU
coauUQDMImcwaum3 Mom mucmeoe mafivcmn 0cm m030000 wwNHHmeuoz .95 magma

178

0000.0
0500.0
0000.0
0000.0
8000.0
0505.0
0005.0
8005.0
8005.0
0005.0
0005.0
0055.0
0055.0
0055.0
0805.0
8005.0
0005.0
0085.0
8005.0
0005.0
0005.0
8050.0
0000.0
0000.0
0800.0
0000.0
0008.0
0000.0
5000.0
0000.0
0880.0
0050.0
0000.0

C
E

00

"

8000.0 0000.0
0000.0 5000.0
0050.0 0000.0
0050.0 0500.0
0050.0 0800.0
8000.0 0000.0
0800.0 0505.0
5000.0 0805.0
8000.0 0005.0
0000.0 0005.0
0080.0 5005.0
0080.0 0055.0
0000.0 0055.0
0000.0 0505.0
0000.0 0005.0
0000.0 0005.0
0000.0 0085.0
8000.0 0005.0
0050.0 0805.0
5500.0 8005.0
0000.0 0000.0
5000.0 0000.0
5008.0 0080.0
0088.0 8000.0
0008.0 0000.0
0000.0 0000.0
0000.0 0008.0
0000.0 0000.0
0000.0 5000.0
0000.0 0000.0
0800.0 0880.0
0000.0 0050.0
0000.0 0000.0

cu GE

: 0

84 u

0000.0
0000.0
0050.0
5050.0
0050.0
0000.0
0800.0
5000.0
0800.0
0080.0
0880.0
0000.0
8000.0
0800.0
0000.0
0000.0
0000.0
0000.0
0000.0
5000.0
0000.0
0000.0
0058.0
0888.0
0008.0
0000.0
0000.0
0000.0
0000.0
0000.0
0800.0
0000.0
0000.0

0000.0 0000.0
0000.0 0000.0
5000.0 0050.0
0000.0 0000.0
5805.0 0000.0
0005.0 0000.0
0005.0 0800.0
0805.0 0500.0
0005.0 0000.0
0005.0 0500.0
8050.0 0000.0
0000.0 0008.0
0580.0 0008.0
0008.0 0000.0
0000.0 0000.0
0880.0 0800.0
0000.0 0000.0
:5 a»
0 n c

5000.0 8005.0
5000.0 8800.0
0580.0 0050.0
0000.0 0000.0
0000.0 0080.0
0005.0 0000.0
8505.0 0000.0
0005.0 0050.0
0000.0 0080.0
0500.0 0000.0
8000.0 5008.0
0500.0 8008.0
0008.0 0800.0
0000.0 0000.0
0000.0 8800.0
0880.0 0000.0
0000.0 0000.0
:5 a»
0 n C

5080.0 5005.0
0000.0 0005.0
0000.0 0500.0
5000.0 0000.0
0000.0 0000.0
0005.0 0500.0
0085.0 5050.0
0505.0 0880.0
0800.0 0000.0
0000.0 8058.0
0000.0 0508.0
0000.0 0000.0
8008.0 0000.0
0000.0 0000.0
0000.0 8000.0
8000.0 0000.0
0000.0 0000.0
as £0
0 n c

1 £u03 c00vumm mumsUm m we m0m00mcm anomnuI000vmeuommU

:00uUGSMI0c00um3 uom mucmeos UCHUQmQ 0cm mmzvuon 00N00meuoz .05 m0nma

000.8
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0

179

0800.0

5008.0

8000.0 8008.0 0000.0 8008.0 0000.0
0000.0 0008.0 0500.0
0500.0 0008.0 0000.0 8008.0 8550.0 0008.0 0050.0
5000.0 0008.0 0000.0
0000.0 0008.0 0000.0 0008.0 0880.0 0500.0 0005.0
0000.0 0008.0 0500.0
5000.0 0000.0 0000.0 0800.0 5000.0 5050.0 0000.0
0000.0 8000.0 0000.0
0000.0 0500.0 0000.0 0800.0 0850.0 0000.0 0000.0
0800.0 0000.0 0508.0
0050.0 0050.0 0050.0 8050.0 0000.0 8000.0 0050.0
5050.0 0050.0 0000.0
0000.0 8800.0 8000.0 0000.0 0805.0 0000.0 0000.0
0000.0 8500.0 0500.0
0080.0 5080.0 0000.0 0080.0 8005.0 0000.0 0000.0
0800.0 0080.0 0000.0
0000.0 0000.0 0500.0 0000.0 0005.0 8050.0 0000.0
0500.0 5000.0 0500.0
0005.0 0000.0 0055.0 0000.0 5080.0 0000.0 0050.0
0005.0 0500.0 0000.0
0085.0 0000.0 0805.0 0050.0 0000.0 8000.0 0005.0
5005.0 0000.0 0500.0
5050.0 0080.0 0500.0 0080.0 0000.0 0000.0 0000.0
0500.0 0000.0 0000.0
0050.0 0000.0 0000.0 0800.0 0008.0 0050.0 0000.0
0800.0 0000.0 0508.0
0008.0 0000.0 0008.0 0000.0 8000.0 8500.0 0050.0
0000.0 0000.0 0000.0
5000.0 0800.0 5000.0 0800.0 0000.0 0500.0 0000.0
0000.0 0050.0 0500.0
0880.0 0000.0 0880.0 0000.0 0080.0 0000.0 0000.0
0050.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0

CE Cu. CE Cu. 0.05 Cu, Ilmwm
0 nl... 0H 0 N C H nu c

00.0 n 1 £003 c00900m mumsvm m 00 m0mh0mcm xuomspIQOHHmEuommv
COHvocnmI0C0mum3 mom mucmEOE 0:0Ucmn 0cm mmsvuou U¢N00meuoz .05 m0nt

Second invariant of stress:

\lO‘U‘ioDUONI-‘L-b 000.5me0.

800.5me0.

Table 8.
deformation—theory with n
and ® = 3.00

i

i

i

i

1 2 3
0.000 0.030 0.132
0.259 0.274 0.327
0.618 0.617 0.621
0.868 0.863 0.850
1.064 1.057 1.037
1.231 1.224 1.202
1.383 1.376 1.355

Normal stress: Ozz/‘VE k

1 2 3
0.000 0.000 0.000
0.406 0.403 0.391
0.627 0.627 0.625
0.741 0.744 0.754
0.812 0.818 0.834
0.860 0.867 0.885
0.892 0.898 0.916

Shear stress. onz/k|

1 2 3
0.000 0.000 0.000
0.307 0.293 0.251
0.474 0.457 0.405
0.566 0.549 0.498
0.636 0.621 0.574
0.701 0.687 0.645
0.766 0.754 0.718

Shear stress: Io /k|
yz

1 2 3
0.000 0.172 0.363
0.000 0.160 0.333
0.000 0.122 0.258
0.000 0.086 0.184
0.000 0.054 0.116
0.000 0.025 0.054
0.000 0.000 0.000

800.5me0-

J

180

2
4

0.314
0.433
0.646
0.840
1.011
1.168
1.320

0.000
0.365
00611
0.761
0.856
0.915
0.947

0.000
0.185
0.317
0.410
0.489
0.568
0.651

0.560
0.515
0.414
0.304
0.195
0.092
0.000

5

0.528
0.585
0.712
0.851
0.989
1.127
1.273

0.000
0.322
0.575
0.753
0.876
0.953
0.991

0.000
0.111
0.205
0.285
0.361
0.443
0.540

0.726
0.685
0.583
0.450
0.302
0.148
0.000

1.00,

6

0.733
0.760
0.824
0.905
0.993
1.093
1.218

0.000
0.275
0.516
0.714
0.872
0.986
1.046

0.000
0.047
0.093
0.139
0.191
0.256
0.352

0.856
0.826
0.741
0.612
0.443
0.236
0.000

Typical stress output for warping-function
2: P =

7

0.923
0.936
0.968
1.008
1.052
1.101
1.180

0.000
0.231
0.447
0.644
0.822
0.978
1.086

0.000
0.000
0.000
0.000
0.000
0.000
0.000

0.961
0.940
0.877
0.770
0.612
0.379
0.000

181

mmv0.H
ma¢0.H
¢0m0.H
bvm0.H
00m0.H
b0N0.H
MNNO.H
0bH0.H
0NHO.H
N500.H
ma00.H
0000.0
H000.0
0000.0
mmbm.0
0m0m.0
0Nmm.0
00vm.0
mem.0
HHHm.0
0Nmm.0
5000.0
bmm0.0
mm00.0
vab.0
0Hm0.0
0vH0.0
vmmm.0
Hmav.0
mva.0
000N.0
0¢0H.0
0000.0

mth.H
NO0H.H
¢0NH.H
000H.H
NMb0.H
0m¢0.H
>000.H
mmbm.0
vomm.0
0N00.0
00N0.0
0005.0
vm00.0
N0mm.0
N¢0¢.0
NOHN.0
0000.0

a»

mu:

bth.H
VNVN.H
mm0N.H
HNOH.H
00HH.H
HN50.H
¢NNO.H
0000.0
0HH0.0
0500.0
0055.0
vvm0.0
0000.0
mmm¢.0
000m.0
mHON.0
0000.0

00.0
05.5
00.5
mm.b
00.5
05.0
00.0
0N.0
00.0
05.0
00.0
0N.m
00.0
mh.v
00.v
mm.v
00.v
mb.m
00.m
mm.m
00.m
mb.N
00.N
mm.m
00.N
mb.H
00.H
mN.H
00.H
mb.0
00.0
mm.0
00.0

lu..l
mo

00.0 n 1 spa: coauumm mumaUm m mo mamhamcm whoonu

Icoauwsuommv coapucsmlmmmuum Mom mmsvuog Umuwameuoz

.mm magma

182

>0N0.0
0va0.0
5000.0
0005.0
0005.0
0055.0
0005.0
0H0>.0
0mmn.0
0005.0
N050.0
0¢N0.0
Hbv0.0
Homv.0
0000.0
m¢¢H.0

0000.0

C
E

00.H

0000.0
N000.0
NH00.0
HHb0.0
0000.0
v0¢0.0
0000.0
vNH0.0
0000.0
0000.0
0mN0.0
NN>¢.0
bmov.0
0NHm.0
000N.0
b¢0H.0
0000.0

G
u

0050.0
0H00.0
0mv0.0
0v~0.0
0m00.0
0005.0
0005.0
m0mb.0
bm00.0
0000.0
NHH0.0
0000.0
050v.0
0b0m.0
0HON.0
vaa.0

0000.0

C
E

v0ab.0
0m0>.0
Nb00.0
0000.0
b0v0.0
m0N0.0
b¢00.0
m0>0.0
vmvm.0
H¢H0.0
0m5v.0
00Nv.0
0500.0
bv0N.0
>00N.0
0¢0H.0
0000.0

a
u

000.¢
050.0
00>.m
0N0.m
000.0
mbm.m
00N.m
0NH.m
000.0
0N0.N
00>.N
0N0.N
000.N
mum.m
00N.N
0NH.N
000.N
050.H
00>.H
0N0.H
000.H
0bm.H
00N.H
0NH.H
000.H
050.0
005.0
0N0.0
000.0
05m.0
00N.0
0NH.0
000.0

0.0

1 nud3 coauumm GHMSUm m mo mamhamcm h000£¢l¢0dumauomoo

coauUGDMImm000m Mom mucmeos 0GHUnon Ucm mmsvuou vaHHmEuoz .30 magma

Second invariant of stress:

\lO‘U‘hb-wNI—‘t—h

Table 10.

i

1

0.000
0.253
0.622
0.875
1.070
1.237
1.387

183

Typical stress output for stress-function
deformation-theory with n = 2, p =
and ® = 3900

2

0.030
0.269
0.622
0.869
1.063
1.229
1.380

3

0.131
0.325
0.625
0.856
1.043
1.207
1.359

Normal stress: 022/ V3 k

\lmtflPWNI-‘L—I-

Shear stress:

\lO‘iflbWNl-‘L—h

i

i

1

0.000
0.407
0.624
0.736
0.808
0.856
0.888

1

0.000
0.295
0.483
0.578
0.646
O. 710
0.773

2

0.000
0.404
0.624
0.740
0.813
0.862
0.894

 

2

0.000
0.283
0.465
0.561
0.631
0.696
0.762

«In/kl

3

o. 000
0.392
0.622
0.750
0.830
0.881
0.913

3

0.000
0.246
0.410
0.508
0.583
0.654
0. 725

Shear stress: loyz/k|

\lO‘U‘l-bUJNI-‘Lb

i

1

0.000
0.000
0.000
0.000
0.000
0.000
0.000

2

0.174
0.161
0.125
0.089
0.056
0.026
0.000

0.362
0.334
O. 264
0.190
0.120
0.056
0.000

J2

4

0.307
0.436
0.650
0.845
1.015
1.173
1.323

0.000
0.365
0.609
0.758
0.853
0.912
0.944

0.000
0.184
0.321
0.417
0.498
0.577
0.657

0.554
0.514
0.420
0.311
0.201
0.096
0.000

5

0.523
0.584
0.715
0.855
0.993
1.131
1.275

0.000
0.323
0.573
0.750
0.872
0.950
0.990

0.000
0.112
0.208
0.290
0.368
0.453
0.543

0.723
0.684
0.586
0.457
0.310
0.154
0.000

1.00,

6

0.729
0.758
0.825
0.907
0.997
1.097
1.218

0.000
0.275
0.515
0.712
0.868
0.982
1.047

0.000
0.048
0.094
0.141
0.195
0.266
0.350

0.854
0.824
0.742
0.617
0.453
0.249
0.000

0.933
0.949
0.982
1.016
1.052
1.098
1.180

0.000
0.228
0.441
0.639
0.822
0.981
1.086

0.000
0.000
0.000
0.000
0.000
0.000
0.000

0.966
0.947
0.887
0.779
0.612
0.367
0.000

184

vm00.H
HH00.H
h0b0.H
N050.H
bm>0.H
0H50.H
H000.H
H000.H
0N00.H
0000.H
0000.H
HH00.H
0bv0.H
vmv0.H
mbmo.H
NNmo.H
d0m0.a
NON0.H
mma0.a
>000.H
Hb00.0
H500.0
0vb0.0
0000.0
0000.0
0¢H0.0
H0h0.0
m0H0.0
0NN>.0
0H00.0
H0bm.0
050H.0
0000.0

0mH0.H
0H50.H
00N0.H
H00¢.H
>0mv.a
H0hm.a
0HNm.H
0HON.H
500H.H
0mNH.H
mmv0.H
¢N00.0
0000.0
00Hb.0
0500.0
0000.0
0000.0

c
#

H u G

oo.m
ma.»
cm.»
34.
cos
2.6
8.6
mud
cod
mad
8....
mm.m
cod
2.4
84
mm...
84
2..“
Sim
mm.m
oo.m
ms.m
om.~
mmxm
oo.m
3.4
84
mm;
84
3.0
36
36
8.0

LI
mm

00.0 n 1 spa: cowuumm “waauuau m mo mwmhamcm huomnu

Icoaumeuommn coaaucaMI0CHQHm3 Mom mwsvuou nouaameuoz

.MHH OHQMB

185

0005.0
0005.0
0005.0
0005.0
0H05.0
5055.0
0055.0
0555.0
0055.0
0055.0
0055.0
0055.0
0005.0
0005.0
0005.0
0005.0
0505.0
0005.0
0005.0
0005.0
0H05.0
0005.0
0505.0
55H5.0
0005.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
050H.0
0000.0

00.H n

0HH5.0
00H5.0
0005.0
0505.0
0005.0
0005.0
0005.0
HH05.0
H000.0
0500.0
5000.0
0000.0
0000.0
0000.0
0000.0
0050.0
0050.0
0000.0
H000.0
0000.0
0000.0
0000.0
0000.0
H000.0
0H00.0
0000.0
0500.0
0000.0
0050.0
0H00.0
050H.0
0000.0
0000.0

HOH0.0
0000.0
0000.0
0005.0
0005.0
0005.0
0005.0
0055.0
0005.0
5005.0
0005.0
0005.0
0000.0
0000.0
0000.0
H000.0

0000.0

C
E

0

0005.0
0005.0
0005.0
50H5.0
00H5.0
0005.0
H000.0
H500.0
0050.0
5H00.0
0000.0
00H0.0
H050.0
5000.0
0050.0
050H.0
0000.0

c
p

C

0000.0
0000.0
0000.0
0000.0
00H0.0
5000.0
0005.0
5005.0
0055.0
0005.0
0005.0
50H5.0
0000.0
5000.0
0H00.0
H000.0

0000.0

C
E

0

0505.0
0H05.0
0005.0
5005.0
0H05.0
00H5.0
5005.0
HH00.0
0550.0
0000.0
5500.0
0500.0
0000.0
0000.0
0000.0
050H.0
0000.0

:
u

C

000H.H
000H.H
000H.H
0H50.H
0000.H
0000.0
0000.0
00H0.0
0000.0
0000.0
5005.0
5000.0
0000.0
0HH0.0
0000.0
0000.0

0000.0

:8

H

HH00.0
0000.0
0000.0
H550.0
0000.0
00H0.0
0555.0
0005.0
5500.0
0000.0
0000.0
0000.0
0050.0
0000.0
0000.0
0050.0
0000.0

:
p

C

1 ngflz cowuumm amasuuau m mo.mamhamcm muomzulcoaumeuommv
coauucsmI0CH00m3 mom mucmeoe mcaocmn 0cm mwSUHOMfiUmNHHmEnoz

000.0
000.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0
050.0
005.0
000.0
000.0
050.0
000.0
00H.0
000.0
050.H
005.H
000.H
000.H
050.H
000.H
00H.H
000.H
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0

INF
0w

.QHH mHQMB

186

mmmo.H mmmo.a mbmH.H mmmm.H oo.m
osmo.a memo.a mo¢H.H m¢vm.a ms.u
«cmo.a mmmo.H mmmH.H mmmm.H om.s
mHNo.H mm~o.a o¢NH.H boom.a mm.s
bmao.a ammo.H ONHH.H oumH.H oo.s
mmao.fl mmmo.a mmmo.H mmmH.H ms.m
NNHO.H Hmao.a mmmo.H HovH.H om.m
mmoo.a mmao.a Hmbo.a s¢NH.H mm.m
swoo.a mmoo.H ommo.a bmoa.a oo.m
mooo.a fimoo.H N¢¢o.fl oomo.H ms.m
mmmm.o mmmm.o mmmo.H memo.H om.m
soma.o Hmmm.o mmao.a mmmo.H mm.m
m0mm.o 0mmm.o «mam.o Nuoo.a oo.m
¢msm.o mubm.o wssm.o Namm.o ms.v
oaum.o mamm.o Nmmm.o afimm.o om.v
bumm.o Homm.o mumm.o mmmm.o mm.¢
mmmm.o mmvm.o Hmam.o ¢omm.o 00.0
wmvm.o 05mm.o bmmm.o mmmm.o m>.m
bmmm.o vam.o mumm.o mmmm.o om.m
waam.o mmom.o Novm.o Hamb.o m~.m
«mom.o mmmm.o moam.o Hmms.o oo.m
Hflbm.o Hmmm.o maps.o mwmh.o mu.m
Ha¢m.o ¢Hmm.o NH«>.o vam.o om.~
mmom.o mmmb.o moo>.o mo¢m.o m~.m
mmmb.o mmvs.o m¢mm.o mmmm.o oo.~
o¢mo.o momm.o mmom.o mm¢m.o mb.H
cmam.o ooam.o Hm¢m.o mqu.o om.a
mowm.o mmflm.o mvsv.o nmmv.o m~.H
Nuav.o H>H¢.o mmmm.o mumm.o oo.H
mmam.o mNHm.o mmom.o omm~.o mb.o
vmo~.o «mom.o abom.o mmmH.o om.o
avoa.o avoa.o o¢oa.o «moa.o m~.o
0000.0 0000.0 oooo.o oooo.o oo.o
a» an cu cu. .%m
NH u c m u a m u a H u a

00.0 n 1 spa: cowuumm mumswm m mo mfimhamcm huomnu
Isoam coapoc30I0caaum3 How mmswuou omuaameuoz .mNa magma

187

50H0.0 0000.0 0000.0 0000.0 0H05.0 00H0.H 0000.0
0000.0 5H00.0 00H0.o 0H00.0 0H00.0 0000.0 00H0.0
0HH0.0 0050.0 00H0.o 0500.0 0000.0. 0000.0 0000.0
0000.0 0050.0 00H0.o 0000.0 0000.0 H000.0 0500.0
0500.0 0H50.0 5OH0.0 0050.0 0000.0 0H00.o 0000.0
0000.0 0000.0 0500.0 0050.0 0500.0 0000.0 0000.0
0H00.0 0000.0 0000.0 0000.0 0000.0 00H0.o 0050.0
H000.0 H000.0 H000.o 5000.0 0000.0 0000.0 0500.0
0000.0 0000.0 0000.0 0000.0 00H0.o 0000.0 0000.0
0000.0 5000.0 0H00.0 0000.0 0500.0 0000.0 00H0.0
0000.0 0000.0 0500.0 0000.0 0000.0 0000.0 000H.0
5000.0 0000.0 0000.0 0000.0 0050.0 0000.0 0500.0
0H00.0 0000.0 H550.o 0000.0 0000.0 0000.0 0005.0
0050.0 0000.0 0H50.o 00H0.0 5H00.0 0005.0 0000.0
0H50.0 00H0.0 5000.0 00H0.0 0000.0 0005.0 0050.0
0000.0 0000.0 0500.0 0005.0 00H0.0 H505.0 050H.0
5000.0 0005.0 0000.0 0505.0 0000.0 00H5.o 0000.0
5000.0 0505.0 0000.0 0055.0 0050.0 0500.0 00H0.0
0H00.0 0055.0 0000.0 0005.0 0500.0 0H00.0 0000.0
H000.0 0005.0 05H0.0 H005.0 0000.0 0000.0 0500.0
00H0.0 H005.0 0000.0 00H5.0 05H0.0 0000.0 0000.0
0000.0 HHH5.0 0000.0 0000.0 0000.0 0050.0 0000.0
0050.0 0000.0 H000.0 H000.o 0050.0 5000.0 0050.H
0500.0 0H00.0 0000.0 H000.o 0000.0 0000.0 0500.H
0000.0 5000.0 0H00.o OH00.0 H000.0 0000.0 0000.H
0050.0 0000.0 0000.0 0000.0 0000.0 5000.0 0000.H
5000.0 0000.0 0000.0 H000.o 0000.0 5050.0 000H.H
0050.0 0000.0 0000.0 5000.0 0000.0 0000.0 0500.0
00H0.0 0000.0 00H0.0 00H0.o 000H.0 0050.0 0005.0
000H.0 0000.0 000H.o 0000.0 000H.o 00H0.o 0000.0
000H.0 000H.0 000H.0 000H.O 000H.o 000H.0 0050.0
H000.0 H050.o 0000.0 H05o.o 0000.0 0050.0 050H.o
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
CE GD. GE CU. CE GU. IIWII

m u C m N c H nu AH m 0

00.0 n 1 £003 coauumm mumswm a mo mamhamcm huoonulonm
coauUGSMImcaaumz mom mvcmeoe madncmn vam mmsvuou voNHHwEHoz .nma magma

 

188

0000.0
0005.0
5005.0
0005.0
0005.0
0005.0
5005.0
0005.0
0005.0
0555.0
0055.0
5055.0
0005.0
0005.0
0005.0
0005.0
0005.0
0005.0
0005.0
5005.0
0000.0
5050.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
5000.0
0000.0
0000.0
0050.0
0000.0

C
E

00

0500.0
0000.0
0000.0
0050.0
5050.0
0050.0
0000.0
0000.0
0000.0
5000.0
0000.0
0000.0
5000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0050.0
0000.0
0000.0
0000.0
0000.0
0000.0
0500.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0

0000.0
0000.0
0000.0
0005.0
0005.0
0005.0
0005.0
0005.0
0005.0
0555.0
5055.0
0005.0
0005.0
0005.0
0005.0
0005.0
0005.0
0005.0
0005.0
0005.0
0000.0
0000.0
0000.0
0500.0
0000.0
0500.0
0000.0
5000.0
5000.0
0000.0
0000.0
0050.0
0000.0

C
E

0

00.0 n 1 £u03 c00uumm mumavm m mo m0m00mcm hu00£UI300M
GOHuUGDMImc0aum3 Mom mucmeoa 0:00:09 0cm mmavuou ©0N00meuoz

0000.0
0000.0
0000.0
0050.0
0050.0
0050.0
0500.0
0000.0
0500.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0050.0
0000.0
0000.0
0000.0
0000.0
5000.0
0000.0
0000.0
0000.0
0000.0
0000.0

c
p

C

5500.0
0000.0
0000.0
0000.0
0000.0
0005.0
5005.0
0005.0
0005.0
0055.0
0055.0
5005.0
0005.0
5005.0
0005.0
5005.0
0005.0
0005.0
5005.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0050.0
0000.0

C
E

0

0000.0
0500.0
0000.0
0050.0
0050.0
0000.0
5000.0
0000.0
0000.0
0500.0
0000.0
5000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0500.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0

a
u

C

0050.0 00O5.o
0000.0 0005.0
0000.0 5000.0
0500.0 5000.0
0500.0 0550.0
0000.0 0000.0
0000.0 0000.0
5500.0 0000.0
0005.0 0000.0
0005.0 0000.0
0055.0 0000.0
0505.0 0500.o
0505.0 0000.0
0005.0 0000.0
0005.0 0000.0
0005.0 0000.0
0000.0 0000.0
0500.0 0000.0
0500.0 0000.0
0000.0 0000.0
0000.0 0000.0
0050.0 0000.0
0500.0 0000.0
0000.0 5000.0
5050.0 0000.0
5000.0 0500.0
0000.0 0000.0
0000.0 0000.0
0050.0 0000.0
5000.0 0000.0
0000.0 0000.0
0050.0 0000.0
0000.0 0000.0
GE up
0 n C

0000.0 0005.0
0000.0 0005.0
0000.0 0005.0
0000.0 0000.0
0050.0 0000.0
0000.0 0050.0
0000.0 0000.0
0000.0 5500.0
0000.0 0000.0
0005.0 0000.0
mubb.o mmom.o
5005.0 0000.0
0005.0 0050.0
0005.0 0000.0
5005.0 0000.0
0050.0 0000.0
0500.0 5000.0
0000.0 0500.0
0000.0 0000.0
0000.0 0500.0
0000.0 0000.0
0000.0 0000.0
0500.0 0550.0
0000.0 0000.0
0000.0 0000.0
0000.0 0000.0
0000.0 0000.0
0000.0 0000.0
0000.0 0500.0
0000.0 0000.0
0000.0 0000.0
0000.0 0000.0
0000.0 0000.0
GE a»
0 n C
.000 m0QMB

000.0
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0
050.0
005.0
000.0
000.0
050.0
000.0
000.0
000.0

«mm

189

0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.H
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 H0H0.0 00H0.0 0NH0.0 0000.0 0000.0 0NHO.H
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.H
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.H
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.H
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.H
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 00NH.H
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 Hmbm.0 0000.0
0000.0 0000.0 0000.0 00H0.0 H000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0HON.0 H0H0.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 00HN.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 000H.0 0000.0
0000.0 0000.0 H000.0 0000.0 0000.0 000H.0 0000.0
0000.0 000H.0 0000.0 000H.0 0000.0 000H.0 0000.0
0000.0 0000.0 0000.0 000H.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 000H.0 0000.0 NbHH.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 000H.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0
CE cu GE G» as cu lwl.

0 n a 0 u a H n a 0m0

00.0 n 1 £003 c0000mm mumdvm 0 mo m0m>0mcm hu0m50I300m
c00uoc50I0C0QH03 how mucmEOE 0c00cmn 0cm mmsvho0 0000006uoz .000 w0QmB

1

Table 13.
flow-theory with n

Second invariant of stress:

\lmmbWMl-‘L-b

i

1

0.000
0.257
0.608
0.852
1.043
1.206
1.354

2

0.030
0.272
0.608
0.847
1.036
1.199
1.347

3

0.130
0.324
0.613
0.836
1.018
1.178
1.327

Normal stress: ozz/N/B k

\lmU'lobmei-‘K—h

Shear stress: loxz/kl

i

1

0.000
0.404
0.617
0.730
0.800
0.847
0.878

2

0.000
0.401
0.619
0.733
0.805
0.853
0.884

3

0.000
0.390
0.618
0.743
0.822
0.872
0.903

3

0.000
0.254
0.406
0.499
0.574
0.644
0.716

3

0.361
0.329
0.256
0.186
0.118
0.055

1 l 2

J'

1 0.000 0.000
2 0.306 0.293
3 0.475 0.458
4 0.565 0.549
5 0.635 0.620
6 0.699 0.685
7 0.764 0.752
Shear stress: ‘0 /k|

yz
i l 2

J'

1 0.000 0.172
2 0.000 0.158
3 0.000 0.120
4 0.000 0.087
5 0.000 0.055
6 0.000 0.026
7 0.000 0.000

0.000

90

= 29 P
and @= 3000

J2

4

0.307
0.427
0.638
0.828
0.994
1.146
1.294

0.000
0.364
0.606
0.752
0.845
0.903
0.934

0.000
0.189
0.318
0.411
0.490
0.567
0.649

0.554
0.508
0.411
0.306
0.199
0.095
0.000

5

0.516
0.575
0.702
0.839
0.973
1.107
1.248

0.000
0.324
0.571
0.745
0.865
0.942
0.981

0.000
0.113
0.207
0.287
0.363
0.444
0.535

0.718
0.677
0.577
0.449
0.305
0.150
0.000

6

0.717
0.745
0.811
0.891
0.977
1.076
1.196

0.000
0.276
0.514
0.708
0.863
0.974
1.034

0.000
0.048
0.093
0.138
0.190
0.261
0.355

0.847
0.816
0.733
0.608
0.444
0.243
0.000

Typical stress output for warping-function
= 1000,

7

0.904
0.919
0.953
0.993
1.032
1.080
1.158

0.000
0.232
0.446
0.638
0.816
0.970
1.076

0.000
0.000
0.000
0.000
0.000
0.000
0.000

0.951
0.930
0.869
0.765
0.605
0.374
0.000

191

0N00.H
0000.H
0000.H
0N00.H
0bN0.H
HON0.H
00H0.H
00H0.H
m000.H
0m00.H
0500.0
0000.0
mm00.0
0050.0
5000.0
0500.0
0000.0
0000.0
00H0.0
HN00.0
HN00.0
0500.0
NON0.0
NHO>.0
bm0b.0
Hm00.0
0500.0
HOHm.0
00H0.0
0mam.0
000N.0
000H.0
0000.0

00.0 I 1 and: coauuom mumavm m mo mflmhamcm knows“
Izoam coauUGSMImmmuum mom mmavuou nmuaameuoz

b00H.H
H0mH.H
0NOH.H
mHmH.H
00HH.H
000H.H
0m00.H
0050.H
0000.H
0000.H
0mm0.H
00H0.H
0000.0
0000.0
bH00.0
mH00.0
m0H0.0
0000.0
00b0.0
0m00.0
0NH0.0
0055.0
0N0h.0
0H0>.0
m0m0.0
0000.0
b00m.0
mm>0.0
000m.0
N00m.0
whom.0
m00H.0
0000.0

cu

0 n c

0Hbm.H
NNmm.H
HmmN.H
0MHN.H
mm0H.H
bNbH.H
mamH.H
HomH.H
NO0H.H
0000.H
NN00.H
00m0.H
0NHO.H
0000.0
5000.0
0Hm0.0
0H00.0
N000.0
mbm0.0
0N00.0
m00h.0
0bmb.0
0000.0
0m00.0
0000.0
0000.0
m000.0
00N0.0
0000.0
0m0m.0
0b0H.0
0NOH.0
0000.0

a»

an:

00.0
mu.b
00.5
mw.b
00.5
05.0
00.0
mm.0
00.0
0b.m
00.0
mm.m
00.0
05.0
00.0
mm.0
00.0
mb.m
om.m
mm.m
00.m
mb.m
00.N
mm.m
00.N
0>.H
00.H
0N.H
00.H
05.0
00.0
0N.0
00.0

IUF
mo

.m0H magma

192

0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
0000.0 0000.0 0000.0 0000.0 000.0
CE G00 CE Cu. u
m n C N n C Q 0

00.0 n 1 £003 co0uumm mumsvm m mo m0m00mcm huomsul300m
COHuUGSMImmmuum uom mgcmeoe 0:00:09 cam mmSUHOp nmn00msuoz .n00 m0nma

Second invariant of stress:

000.0me0.

Normal stress:

00501000100.

Shear stress:

~403thtum)Hu0

Table 15.

i

i

i

1

0.000
0.251
0.612
0.857
0.105
1.211
1.358

1

0.000
0.405
0.616
0.725
0.795
0.842
0.874

1

0.000
0.294
0.482
0.576
0.645
0.708
0.771

193

Typical stress output for stress-function

 

flow-theory with n
and (9 = 3.00

Shear stress: loyz/kl

0ChUthuNHHh»

i

1

0.000
0.000
0.000
0.000
0.000
0.000
0.000

2 3
0.031 0.132
0.268 0.324
0.611 0.616
0.852 0.840
1.041 1.022
1.204 1.183
1.351 1.331

arm/V? k

2 3
0.000 0.000
0.402 0.390
0.617 0.616
0.729 0.740
0.801 0.818
0.849 0.868
0.880 0.899

ze/kl

2 3
0.000 0.000
0.282 0.245
0.464 0.409
0.559 0.506
0.630 0.582
0.695 0.652
0.759 0.723

2 3
0.175 0.363
0.162 0.335
0.126 0.265
0.090 0.190
0.056 0.120
0.026 0.056
0.000 0.000

= 2: H

J:2

4

0.307
0.430
0.642
0.831
0.996
1.149
1.297

0.000
0.363
0.603
0.749
0.842
0.900
0.932

0.000
0.183
0.319
0.415
0.496
0.574
0.654

0.554
0.514
0.420
0.311
0.202
0.096
0.000

= 1000,

5

0.520
0.581
0.707
0.842
0.975
1.110
1.250

0.000
0.322
0.568
0.743
0.863
0.940
0.979

0.000
0.111
0.206
0.288
0.366
0.450
0.540

0.721
0.682
0.585
0.456
0.310
0.154
0.000

0.721
0.750
0.816
0.894
0.980
1.077
1.194

0.000
0.275
0.512
0.706
0.860
0.973
1.036

0.000
0.047
0.093
0.140
0.194
0.265
0.348

0.849
0.820
0.738
0.614
0.450
0.248
0.000

0.919
0.935

0.967

0.999
1.033
1.078
1.158

0.000
0.228
0.438
0.633
0.814
0.972
1.076

0.000
0.000
0.000
0.000
0.000
0.000
0.000

0.959
0.940
0.881
0.774
0.608
0.365
0.000

APPENDIX VI

GRAPHICAL REPRESENTATION 0F TORQUES

AND MOMENTS

Since all work-hardening formulations provide

essentially the same numerical results when-43 = constant,

the results are shown graphically for only té: deformation-
theory warping-function formulation. These figures
include both the unit square and unit radial solutions.

In all graphs, values of n were chosen so that the
general effect of the work—hardening on the torques and
moments is illustrated. More specific information on the

exact values of the torques and moments may be obtained

from Appendix V.

194

195

12

 

Figure 24a.“ Torque vs. unit afgle of twist for
warping-function deformation-theory solution
of a square section with u a 0.00

196

om.o a 1 nuaz coauuwm mHMSUm m mo coausHOm knowsuIGOHnnEHOMMU
coauucaMImcaaumz mom umazp mo wamcm vac: .m> pcmEOEI05UMOB

nu!
mo

 

-—4

 

m.o

m.o

o.H

.nvm ouaoam

 

N.o

m.o

v.0

m.o

0.0

5.0

m.o

m.o

 

197

oo.H u 1 nfiaz coauumm mumavm m mo :OHvsHOm whoosulcoaumeuommv

coauucsmlmadaums Mom umflzu mo mamas uans .m> unmeoaloavuoa

 

.uvm wusmah

 

.1H.0

lm.o

1m.o

.lm.o

 

m.o

¢.o

m.o

0.0

5.0

m.o

 

o.H

198

oo.m u 1 Suez coauumm mumsvm m mo coau3H0m huowzylcoaumauommo
COHHUGSMIacflaumz mom umfizu mo mamcm “fins .m> “c0605I090u08 .pvm whamfim

a. a.

 

 

 

0.N m.H 0.H m.0 0.N m.H 0.H m.0
a l _ _ 0 _ _ _ _ 0
l H.0 |.H.0
1 N.0 I.m.0
I m.o Ilm.0
H-
ivoo lfioo
1 80 m.-. 480
as a
l 0.0 u.m.0
l 5.0 1.5.0
1 m.o no.0
m 1m o :m o
H l 0.H |.0.H

 

199

 

 

 

Figure 25a. Torque vs. unit angle of twist for
warping-function deformation-theory solution
of a circular section with p a 0.00

200

oo.H u 1 nuas coauumm Hoaauuau m mo cowv3H0m mucosulcodumenommp
cowpuc3MIGCHQH03 Mom uwakp mo mamcm was: .m> usuaoenmzvuoa .Qmm unawam

gum—l.
mo

 

H.0
l N.0

L 0.0

 

ma.‘

L 0.0

l 0..n

l H...”

 

lml
mo

 

~
—1
—I
I—

 

H.0

«.0

m.0

1.0

m.0

0.0

0.0

0.0

0...”

H...”

9.

10.

BIBLIOGRAPHY

Budiansky, B., "A Reassessment of Deformation Theories
of Plasticity," J. Appl. Mech., Vol. 29, Trans ASME,

Budiansky, B., and Mangasarian, O. L., "Plastic Stress
Concentration at a Circular Hole in an Infinite Sheet
Subjected to Equal Biaxial Tension," J. Appl. Mech.,
Vol. 27, Trans ASME, 1960, pp. 59-64.

Drucker, D. C., "A More Fundamental Approach to
Plastic Stress-Strain Relations," Proceedings of the
First U.S. National Congress of AppIIed MeChanics,
m: 1951, pp. 414-727.

 

Forsythe, G. B., and Wasow, W. R., Finite Difference
Methods for Partial Differential Eguations, Jofin WIIey
and Eons, Inc., New York, N.Y., 19 .

Fox, L., Numerical Solution of Ordinary_§nd Partial
DifferentIaI Equations, AddiSonJWesley Publishing 5o.,
Tnc., Reading, Mass., 1963.

 

Greenberg, H. J., Dorn, W. 5., and Wetherell, E. H.,
"A Comparison of Flow and Deformation Theories in
Plastic Torsion of a Square Cylinder," PLASTICITY,
Proceedings of the Second 5 m osium on NavaI Struc-
tural Mechanics, Editédfiby E. H. Lee and P. 3.

Eymonds, Pergamon Press, Inc., 1960, pp. 279-96.

 

 

Handelman, G. H., "A Variational Principle for a
State of Combined Plastic Stress," Quart. Appl. Math.,
V01. 4, 1944, pp. 351-530

Handelman, G. H., Lin, C. C., and Prager, W., "On the
Mechanical Behavior of Metals in the Strain-Hardening
Range," Quart. Appl. Math., Vol. 4, 1947, pp. 397-407.

Hencky, H., "Zur Theorie plastischer Deformationen und
der hierdurch im Material hervorgerufenen Nach-
spannungen," Z. angew. Math. Mechanik, Vol. 4, 1924,
pp. 323-34.

Hildebrand, F. B., Introduction to Numerical Analysis,
McGraw-Hill Book Company, Inc., New.York, N.Y., I956.

201

11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

202

Hill, R., "A Variational Principle of Maximum Plastic
Work in Classical Plasticity," Quart. J. Mech. Appl.
Hath., V01. 1, 1948’ pp. 18-28

Hill, R., The Mathematical Theory of Plasticity,
Clarendon Press, Exford: 1950.

Hill, R., and Siebel, M. P. L., "On the Plastic Dis-
tortion of Solid Bars by Combined Bending and Twist-
ing," J. Mech. Phys. Solid, Vol. 1, 1953, pp. 207-14.

Hodge, P. G., "On Torsion of Plastic Bars," J. Appl.
Mech., Vol. 16, Trans ASME, 1949, pp. 399-405.

Imegwu, E. 0., "Plastic Flexure and Torsion," J. Mech.
Phys. Solids, Vol. 8, 1960, pp. 141-46.

Imegwu, E. 0., "Combined Plastic Bending and Torsion,"
J. Mech. Phys. Solids, Vol. 10, 1962, pp. 277-82.

Lubahn, J. D., and Felgar, R. P., Plasticity and Creep
of Metals, John Wiley and Sons, Inc., New York,’NIY.,
I961.

 

Mangasarian, O. L., "Stresses in the Plastic Range
Around a Normally Loaded Circular Hole in an Infinite
Sheet," J. Appl. Mech., Vol. 27, Trans ASME, 1960,
pp. 65-73 .

von Mises, R., "Mechanik der festen KSrper im
plastisch deformablen Zustand," GBttinger Nachrichten,
math.-phys. Klasse, 1913, pp. 582-92.

Onat, E. T., and Prager, W., "Limit Analysis of
Arches," J. Mech. Phys. Solids, Vol. 1, 1953,
pp. 77-89.

Osgood, W. R., "Stress-Strain Formulas," J. Aero.
Science, Vol. 13, 1946, pp. 45—48.

Parter, S. V., "Mildly Nonlinear Elliptic Partial
Differential Equations and their Numerical Solution.
I," Numerische Math., Vol. 7, 1965, pp. 113-28.

Parter, S. V., and Greenspan, D., "Mildly Nonlinear
Partial Differential Equations and their Numerical
Solutions. II," Numerische Math., Vol. 7, 1965,

pp. 129-46.

24.

25.

26.

27.

28.

29.

30.

31.

32.

33.

34.

203

Piechnik, S., "The Influence of Bending on the Limit
State of a Circular Bar Subjected to Torsion," Arch.
Mech. Stos., Vol. 13, 1961, pp. 77-106.

Piechnik, S., and Zyczkowski, M., "On the Plastic
Interaction-Curve for Bending and Torsion of a Circu-
lar Bar," Arch. Mech. Stos., Vol. 13, 1961, pp. 669-92.

Prager, W., "An Introduction to the Mathematical
Theory of Plasticity," J. Appl. Phys., Vol. 18, 1947,
pp. 375-83.

Ramberg, W., and Osgood, W., "Description of Stress-
Strain Curves by Three Parameters, NACA TN No. 902,
July, 1943.

Salvadori, M. G., and Baron, M. L., Numerical Methods

in Engineering, Prentice—Hall, Inc.,-Englewood CIITFE,
N.J., O

Southwell, R. V., Relaxation Methods in Theoretical
Ph sics, Clarendon Press, Oxford, 1946.

 

Steele, M. C., "The Plastic Bending and Twisting of
Square Section Members," J. Mech. Phys. Solids,
V01. 3’ 1954, pp. 156-66.

Timoshenko, 3., and Goodier, J. N., Theory of Elasti-
city, McGraw-Hill Book Company, Inc., ew YOrk, N.Y.,

Todd, J. (Editor), Survey of Numerical Analysis,
McGraw—Hill Book Company, Inc., New York, N.Y., 1962.

Varga, R. 8., Matrix Iterative Anal sis, Prentice-
Hall, Inc., EnglewoodCIIffs, N.J., 1962.

Wang, C. T., Applied Elasticity, McGraw-Hill Book
Company, Inc., NewIYork, N.Y., 1953.