A NONLINEAR, SEMI-ANALYTICAL FINITE ELEMENT

ANALYSIS FOR NEARLY AXISYMMETRIC SOLIDS

MANSOURSEDAGHAT
Facultad de lngenieria Civil, Universidad de Conception, Conception-Chile

and

LEONARDR. HERRMANN
Department of Civil Engineering, University of California, Davis, CA 95616,U.S.A.

(Received 14 June 1982:receioed for publication 3 August 1982)

Ah&act-The presented research is concerned with the development of the theory and accompanying computer
program for a semi-analytical finite element analysis of non axisymmetrically loaded, nearly axisymmetric solids.
The theoretical basis of the method together with numerical procedures for handling boundary conditions, rigid
body motion and the iterative solution process are described. A finite element program for evaluating the analysis is
discussed, evaluated for effectivenessand appliedto severalexamples.
The range of applicationof the analysis for inhomogeneous, orthotropic, nonlinear and nearly axisymmetric
bodies is demonstrated by a series of examples. The substantial savings in computer time and memory as compared
to a conventional finite element analysis is discussed.

1.~TRODU~ION range of application of the semi-analytical linite element

Standard finite element methods are capable, in principle, method to a wider class of problems. The first restriction
of handling one, two and three-dimensional, linear and is partially relaxed so that “nearly axisymmetric” solids
nonlinear problems. The solution cost however, in- can be analyzed. The second restriction is lifted com-
creases greatly with each added dimension. The exces- pletely, so that the material properties are no longer
sive cost and computer memory requirements generally required to remain constant in the circumferential direc-
make three-dimensional finite element stress analyses tion, thus permitting the analysis of general non-
unfeasible except on the largest machines. However, two homogeneous, anisotropic and nonlinear materials. The
particular classes of three-dimensional bodies, i.e. solids relaxation of the third restriction will be the subject of a
of revolution and infinite prismatic solids can be future paper.
analyzed at a fraction of the cost and with more con- The theoretical basics for the analysis is presented in
venience if semi-analytics finite element methods are Sections 2 and 3, five examples demons~ting the
used. capability of the method are presented in Section 4 and
Semi-analytical finite element methods use a Fourier finally conclusions are drawn in Section 5.
series representation in one coordinate direction and a
standard finite element approximation in the other two. 2.PROBLEMSTATEMENTFOR ~~OMOGENEOUS,ANISOTROPIC
Substantial savings in computational cost are achieved MATERIALS
because the resulting sets of simultaneous equations for This section is concerned with the development of a
the node point unknowns for each of the several terms semi-analytical finite element analysis when the material
in the Fourier series are uncoupled. The cost of solving is anisotropic and/or inhomogeneous. For a linear elastic
the several sets of uncoupled equations is far less than material write the stress-strain relation in matrix form,
solving the one very large set obtained from a conventional i.e.
three-dimensional finite element analysis,
Early applications of the method involved the analysis
of plates and boxes subject to flexure[l], prismatic bars where the transpose of a column vector (tr)= is written,
and concrete box bridges[2,3], axisymmetric solids of (II),
revolution and shetls of revolution with non-symmet~c
loads[4,5]. U~o~unately the range of application of the ifl) = (an a;g,o;? or**CM,a&)
semi-analytical method is Iimited. For solids of rev- (ef = (a %, t;9 Yrz*Y&Y
%A (2)
olution it is required, of course, that the geometry of the
body be axisymmetric; in addition, it is required that the and Dii = Dii(r, 8, z). It is to be noted that even though a
material properties be constant in the circumferential material is homogeneous, if it is orthotropic with prin-
direction, and the boundary condition type (i.e. specified cipal axes of x, y and z, that its properties expressed in
force or specified displacement) should not vary in the cylindrical coordinates are theta dependent [7,8].
circumferential direction. (The magnitudes of the The potential energy II of an elastic body with initial
specified boundary conditions of course can vary in the stresses may be written as
circumferential direction).
The second of these restrictive assumptions was par-
tially lifted in f7,8] to allow for the analysis of solids of
revolution with rectilinear orthotropic properties.
The purpose of the present study is to extend the

385
390 M. SEDACHAT
and L. R. HERRMANN

where (E) are the strains, [II] the matrix of elastic

coefficients, {uO}the initial stresses, {F} body forces,
(8) displacement vector, {P} the applied surface tractions
and d V, dA and A, respectively denote the differential
elements of volume and surface area, and the surface
area over which {P} is applied.
The displacement components in cylindrical coor- E, = % T,,, = eqmlT,,,
dinates are u,, ue, and u,; their theta dependence is
expanded in Fourier series: Y’_=(~+~)T,,,=e,,,,T,,,
MI MI

& = c U”,, cm(m,0) + ml_(l

WI,=0
z u,,sin(mzO) YrH = (-zuml+~-~
r J
T rtt=e,naT ,,,

ue= 2 v,, sin(m,~)t~~~v,~cos(m,~) (4) T ,,, = ehn,T ,,>. (10)

“,=&I

u, = j$w,,, cos (m,f?)

mj=0
t 2 w,? sin (m20).
m2=-0 The body forces are also expanded in Fourier series

To simplify the resulting equations the function T, is

introduced to represent both the cosine and sine func- F, = F,,,,T,,,
tions, i.e. define: F, = F”,,,T ,)I

cos (me) when m 2 0 F: = FJT,,, (II)

T,(8) = sin (me) (5)
when m < 0’ where
It follows that:
dT,_--mT
do ,,,.
Lw
FH, = f F,,T n, d0, etc. (13)
For the same purpose the following definition will be mI 0
needed later.
Because of the theta dependence of the 4, coefficients,
substituting eqns. (8), (10) and (11) into the expression
for potential energy does not yield the usual uncoupled
Fourier expansion. To overcome this difficulty the matrix
Introducing the summation convention that the [D] is decomposed into two parts:
repeated index m requires summation from -Mz to MI
eqn (4) takes on the simplified form
[D(r, 0, z)l = [D(r,z)l +[D(r, 0, z)l. (14)
u, = u,,Tm
The components of [o] are called the base properties
ue = v,T r,,
and in most cases will be obtained by averaging the
uz = w,,T,,,. (8) properties in the theta direction. However, to prevent
coupling between the stiffn_ess matrices for the m and
In general u, is a function of (r, 0, z) while u, is a -m terms in the series, [D] must be selected sp that
function of (r, z), etc. For convenience the two zero fi,, = &=O, i= l-+4. The components of [D] are
order terms of eqn (4) have been written as a single term called the deviation properties which, as their name
in eqn (8), even though in the conventional semi-analy- indicates, represent the deviations from the base values.
tical method they represent the uncoupled phenomena of With eqn (14), an iterative solution to the problem is
torsionless and torsional axisymmetry respectively. used; for iteration N, the stress-strain relation is written
Depending on the symmetry of the body and the applied in the form:
loads not all terms from -Mz to M, will actually be
present for a given problem; denote the number of
[# = [d(r, z)]{6}N +[D(r. 0, z)lN ‘{EIN ‘. (15)
non-zero terms as NT. For simplicity in the development
of the finite element program eqn (7) is written in the
form The possibility of [d] being dependent on the iteration
number is introduced for later use in nonlinear problems.
NT
For each iteration the first part of eqn (15), where the
u, = x ulTtncI,,etc.
,=I
(9) base properties are independent of theta, leads to the
usual series of uncoupled two-dimensional problems (and
where m(l) is the particular Fourier number correspond- associated uncoupled stiffness matrices). The theta
ing to term I in the series; in most of the remaining dependent deviation properties and strain estimates
paper, however, reference is given to eqn (8). {E}~-’ contribute to the load matrices in a form similar
The strain-displacement relations become (the terms to the contributions of initial stresses.
ei,,,are introduced for later use and are the elements of a In view of eqn (IS), eqn (3) for the potential energy
vector {e},). takes on the form:
 Analysis for nearly axisymmetric solids 391

The element of area dA is given by:

Fl=
I( ;(e)N[D]{e)N
” t (eymN-'IelN-'
dA= Jdtdq. (24)
-(F){S}“]dV-I (IWIN dA (16)
AS Expressions for fi, Gi and the Jacobian J may be found
in a standard text on Finite Element Theory (e.g. see [9]).
This expression can be decomposed into two parts; for a The partial derivatives of eqn (22) are of the form:
finite element analysis the first term
+ = $,lJi”F,, etC
I
n, = $~)~[b]{e}~ dV (17)

3.2 Element stiffness marrix

leads to the stiffness matrix, and the second term The element stiffness matrix is derived by considering
the contribution of the element in question to eqn (17).
llL=-
II -(E)N[B]N-‘{e}N-’
” Substituting eqn (25) into eqn (10) and the results into
eqn (17), integrating with respect to 0 and making use of
the orthogonality properties of the Fourier functions, and
+(F){6}N]dVt/ PXSIN dA (18) differentiating with respect to the node point unknowns
A
gives the coefficients of the element stiffness matrix.
Using standard finite element matrix notation[9] the 3 x 3
to the load matrix.
matrix [Km]ij which gives the terms in the element
3. FINITEELEMENT
ANALYSIS stiffness matrix [&I (four Fourier term m) resulting
from considering nodes i and j of the element is:
3.1 Element type
The Fourier series representation developed in Section
2 will be used in conjunction with a four node iso- [Knlij = L l-1 I’, [&JiT[~ll[MJ~dZdv (26)
parametric element approximation for the displacement
coefficients in the r-z plane to develop a semi-analytical
finite element analysis. Details concerning the program to The components of [fil are the elastic base properties
evaluate the analysis and its use may be obtained by previously defined and the matrix [B,,,], has the following
writing to the Authors. The shape functions are: form:

Ni(tv a
rl) = f1 + t&5(1+ 7%) (19)
F,
N
-! mN
0 0
0
where 0’ -7
(27)
&=(l,-l,-1,l) t&Jr =
qi = (191, -l,- 1). (20) ‘: N, F,-+ 0

The coordinate transformation and the approximations 0 G, -TN,

of the Fourier coefficients of the displacements of eqn (8)
are written in the form
The spatial integration in eqn (26) over the domain of
the element is carried out as usual by four point Gaus-
r = i riNi
i=, sian quadratures.
The system stiffness matrices (one for each term of the
Z = f: ZiNi (21) Fourier series) are obtained by a direct combination of
i=, the element matrices. They are functions of the base
properties and the Fourier number m, but most important
and are independent of the iteration number N. Therefore at
iteration one, all the system stiffness matrices for several
urn= $,
ui”Ni the terms included in the Fourier expansion are for-
mulated, reduced and stored on disk. At any subsequent
iteration the new load matrices are evaluated, the reduced
urn= $,
Vi”‘Ni (22) stiffness matrices are recalled from disk and backsub-
stitutions are performed to obtain the new estimate of
w,,, = i WimNi the displacement coefficients.
i=l Formation and reduction of the system stiffness
matrices are a large portion of the total computational
where U,“, Vim and Wi” are respectively the Fourier effort required for the analysis, therefore once the first
coefficients of the node displacements in the r, ~3and z iteration is completed, subsequent iterations require only
directions. The partial derivatives of the shape functions a fraction of the cost of the first.
with respect to r and z are denoted by E and Gi:
3.3 Element load matrix
The element load matrix is derived from the part of the
(23) potential energy given by eqn (18) which is the sum of
392 M. SEDAGHATand L. R. HERRMANN

three terms. The contributions from the last two terms theta dependence of the fi,, terms are approximated in it
are identical to those from a conventional semi-analytical piecewise constant manner and analytical integration is
finite element analysis; the reader is referred to Refs. [4, used. For nonlinear (strain dependent) material proper-
8, 91 for details. The contribution from the first term is ties numerical integration is used. The number of non-
new and requires detailed consideration. trivial, independent terms in eqn (29) equals (NP)(NE)
Substituting eqn (10) into the first term &., of eqn (18) (NT t l)(N7’)/2. The quantity NP is the number of
and letting the index m correspond to summations in- elements in the r, z grid that have unique non-trivial H
volving terms from iteration N and the index n to variations of their elastic properties and NE is the num-
iteration N- 1 (e.g. u,,,~T,,, = a,T,,, and aZ-‘T, = ber of elastic coefficients in [D] that have unique non-
a.T,), dropping the superscript N - 1 from the dij terms, trivial theta dependence. The number of terms in the
and integrating over the element, the following expres- Fourier series is denoted by NT. The quantity NE varys
sion is obtained: from 1 for an isotropic material with constant ~8and
variable E to 21 for a completely general, in-
homogeneous anisotropic material. One of the keys to
the development of an efficient (in computer time and
space) program for the evaluation of the analysis is to
avoid calculating and storing zero and dependent terms
in the [C,,,]” matrices.
+ &esnT-.) + e2,T,di[i12e,,Tn+ . .I
rJded[dq. (28) 3.4 Rectilinear orthotropic properties
Geometric solids of revolution made of orthotropic
The only theta dependent terms in eqn (28) are the materials are becoming increasingly important in
Fourier functions and the deviation material properties engineering applications. Unless the material is cur-
dye’. Matrices [Cm]. of integrals of the above quan- vilinear orthotropic with one principal direction being the
tities are defined by the components circumferential direction, the conventional semi-analy-
tical finite element analysis is not applicable and an
expensive general three-dimensional finite element
(29)
analyses is required. The range of application of the
semi-analytical finite element analysis is expanded in this
where k= m for i= l-4 and k= -m for i=5,6; and section to include rectilinear orthotropic materials (see
I = n for j = 1+4 and I= - n for j = 5,6. Now recall (see also [7,8]).
eqns 10 and 25) that the quantities e,, are functions of The transformation of stresses from rectangular to
the nodal unknowns U,‘“, etc. The partial derivatives of cylindrical coordinates can be expressed as follows
II,, with respect to Uj”, Vim and Wj”’ gives the con- (where the two systems have a common z axis and theta
tribution to the load matrix; is measured from the x axis)
-
I-

gr co?0 sin’ 0 0 0 2 sin e cos e 0 G

0% sin’ % co? % 0 0 -2sin ecos 0 0 ov
c* = 0 0 I 0 0 0 o, (31)
rrZ 0 0 0 cos 0 0 sin e ls 7x2
r@ - sin e cos 0 cos e sin e 0 0 cos’ 6 - sin2 e 0 G”
1%J L 0 0 0 - sin 0 0 cos 0 J j T,,,
-

or in condensed form

{a]= [ Tl{d (32)

[B,liT[Cmln14nrJdZdv. (30)
but

{g’}= [DI]{E’} (33)

The repeated “n” subscript requires summation while where [D’] is the matrix of rectilinear orthotropic
“m” does not. For linear materials the coefficients of material properties; thus:
[C,,,],, are independent of the iteration process hence
they can be computed once and stored for use in the {o1= ]rl]D’Ke’]. (34)
computation of the load matrices at every iteration. For
materials that are nonlinear due to strain dependent It is easy to demonstrate that the strain transformation
properties, the coefficients of [D] and hence of [D] has the form:
change from iteration to iteration, therefore it is neces-
sary to evaluate [C,]” at every iteration. The terms in 1e’J= Vl’le]. (35)
{E}”are evaluated from the results of iteration N - 1.
In the program written to perform the analysis, two Substituting eqn (35) into eqn (34) gives:
subroutines are used for the evaluation of the integrals of
eqn (29). For linear but inhomogeneous materials the {a] = [71[D’1[7-lT{eI. (36)
 Analysisfor nearly axisymmetricsolids 393

Comparing eqns (36) and (1) shows: u, = (8, t ~4,) cos 0 t (8, - z&) sin 0
us = r& - (6, t sin 8 (6, - cos 0
[Dl = mmm (37)
= S, r+? cos t r& 8. (41)
Expanding out this transformation gives:
prevent any the rigid motions from occuring
D,, = 0; , cod 6 + D& sin4 0 t . . . , etc. (38) arbitrary and, thus, assuring uniqueness of
solution requires (refer to 4 and the
Since the Dii coefficients are now functions of theta, following be specified:
the iterative semi-analytical finite element analysis is (a) m= w0 must specified least at
required. While the integrals of eqn (29) could be evalu- point and must be specified least at point for
ated analytically it is found to be more convenient to rf 0.
evaluate them by numerical integration. They are, For m 1: u, u, must specified least at
however, still independent of the iteration process and point and a second either u,, or w,
are computed only once and stored on disk for use in the be specified.
evaluation of the load matrices at every iteration. (c) For =- L, or must be specified least
at point and a second either u-, Do, or
3.5 Nonlinear properties must be specified.
An important factor present in many analyses is the The of points each of cases “b” “c”
strain dependence of the material properties. For such may selected arbitrarily fact they be the
cases the deviation properties defined by eqn (14) are point) as as the following condition satisfied. If
now functions of the iteration number N.-The integrals the the specified displacements
defined in eqn (29), are functions of the Dif coefficients and the body motion are written matrix
and have to be reevaluated every iteration, or possibly form, determinant of coefficient matrix be
every few iterations. In fact two iteration processes are non-zero. For example that condition “b” is
occuring simultaneously, i.e. the iteration to account for satisfied specifying u, at points (ra, and (rp,
the material nonlinearities and the iteration to account The relationship between displace-
for the coupling of the Fourier terms. The program ments u,, ulp and rigid body terms S,
written to evaluate the analysis updates these integrals #Y is:
(i.e. accounts for the material nonlinearity) every
ITERUP iteration. Thus, for ITERUP specified as 1 it
updates the integrals after each iteration. To date no
attempt has been made to determine an optimum range
of values for ITERUP; for the nonlinear example
presented in Section 4 the analysis was done with
ITERUP = 1. 1 za =z,-z,fOorz,fz,.

3.6 Continuity conditions for points on z axis

I1 I ZD
(43)

The cylindrical coordinate displacements u, ue and u, As a second example assume that “b” is satisfied by
of a point on the z axis (r = 0) can be written in terms of specifying u, and w, at two points which yield:
rectangular components u, u and w as:

u,=ucost?tusin0
us=-usinfVtucos0 (39)
Thus r, f 0. Of course if the corresponding term is not
u, = w. included in the expansion of eqn (4) no specification is
required.
Comparing these expressions to eqns (4) and (8) shows
that for nodes with r = 0 the following conditions must be 3.8 Block iterative solution of simultaneous equations
specified: The general form of the simultaneous equations to be
solved for Fourier term m are
for m = 0: u0 = v. = 0
for[m[=l:u,=-mu,andw,,,=O [K,]{S,N} = {FmN}- M2 5 m I M, (45)
forIml>l: u,=o,=w,=O. (4) where [K,,,] is the system stiffness matrix, {SmN}the
vector of Fourier coefficients of the nodal displacements
The program written to evaluate the analysis automa- for iteration N and finally {FmN}is the load vector. The
tically specifies these conditions. evaluation of {FmN} requires a knowledge of all NT
solution vectors from iteration N - 1, {Sz-‘}, and the
3.7 Suppression of rigid body motion matrices of integrals [&In, eqn (29).
The prevention of possible rigid body motion (in order Since memory economy was one of the objectives of
to assure a unique solution for the displacements) is the present development Jacobi iteration is used to ac-
another area where special care is needed. Consider count for the coupling between the Fourier terms. Jacobi
arbitrary rigid body translations S,, S,, S, parallel to the iteration was selected for simplicity; future work will
x, Y, z axes and rotations I&, &, dL about the same axes. investigate possible computational advantages of Gauss-
The corresponding displacements in cylindrical coor- Seidel iteration[lO] and for large computers, the simul-
dinates for a point located at (r, 0, z) are: taneous direct solution of all equations. In the block

CAS Vol. 17, No. 3-F

394 M. SEDAGHAT and L. R. HERRMANN

Jacobi iteration process, at a given iteration N, eqn (45) and the hoop stress around the circular hole ih:
is solved (using only back substitution for N > 1) for
each Fourier term and a total solution vector {A”} = ‘fl,, = P( 1+ 2 tan’ i#~).
({SE,), . ‘ * {S&j) is formed and stored on disk. In the
next iteration the computation of each load matrix It is easily seen that:
{F,N+‘j requires the use of the total sohrtion vector {A”}
from the previous iteration. Once all the load matrices
{FE-“} are updated a new system of equations is solved (49)
for each Fourier term m and a new estimate of the total
solution vector is obtained.
This process is continued until convergence is If d is much larger than R, then the edge of the plate
achieved as measured by a convergence criterion in- has little or no effect on the stress distribution around the
corporated in the program. The first order norm, LI{AN}, hole. On the other hand if d is almost equal to R, the
of a vector is the sum of the absolute values of the hole is very close to the edge of the plate and the body is
components. The convergence criterion of not “nearly axisymmetric”, thus, the finite element
method developed in the previous section would not be
appropriate. Choosing d/R = 4 the ratio of eqn (49) is
approx. f indicating that enough of the edge effect
where c is a specified small quantity, is used. A value of remains to make the problem meaningful.
E = 2 x lo- 3 was found to be appropriate for most of the The semi-infinite plate is approximated as a nearly
examples studied. axisymmetric solid as shown in Fig. 2 where the radius p
is chosen so that the stress o, at p is small. For p = 6R
4. ~UME~CAL EXAMPLES
the value of cr, is 0.03P which is considered to be small.
4.1 ~x~~pie I. Levi-in~njte plate containing a circ~iur Thus a geomet~c solid of revolution is obtained for
hole which the typicai element shown in Fig. Z(b) has the
A semi-infinite thin plate containing a circular hole following variation of properties in the theta direction:
subjected to a uniform internal pressure P is considered,
(see Fig. 1). An analytical solution for this problem can ForO~tYS6,and8,t8,SBc:2~:E=E,,. (SO)
be found in [It]. Along the edge of the plane the stress
cr, is For 8, < t?< & f 0,: E = 0.

The elastic coefficients Dii are written as

(47)
Dij = Oi,
t bi (51)

with & representing the base properties and Dij the

deviation properties. The base modulus is taken to be the

(a) Finite Element Grid With I4 Eisnnts

(a) Semi Infinite Plate With Hole

(b) RsgioRs A and B

(b) Details of the Hole
Fig. 2. “Nearly axisymmetric” idealization of Example 1 and
Fig. 1. Example 1. finite element grid.
 Analysisfor nearlyaxisymmetric
solids 395

average rn~~us, i.e. Analytical Valur ( 1.U)

A
B=&*B,+E,(27r-&)I. (52)
X x (1.12)
1.1 -
*
The deviation modulus is then simply computed using X .
eqn (51); for the fictitious material fi = - Eo. The value X
l
of Poisson’s ratio Y for the fictitious material is unim- %nax
P .
portant and taken to be 0.25.
A typical finite element grid with 10elements in region
A and 4 in region B is shown in Fig. 2(a).RegionsA and 14 Element Grid
l

B are definedin Fig. 2(b).The z axis is located at point C x 24 Element Grid

of Fig. 1 and is perpendicularto the plane of the paper, it 1.0 I I , I I I
is the axis of rotation for the body. For the choice of the 12 3 4 5 6 7
axis B= 0 (Fig. 2), there is symmetry of defo~ation
about the line B= r//2,hence, the Fourier terms that need Wumbsr of Fourier Terms
to be included in the expansion are the even terms for Fig. 3. Comparison of maximum hoop stress from the finite
the first series of eqn (4) and the odd terms for the elementanalysis to the exact value--Example 1.
second. Had, instead, 6 been measuredfrom the vertical
then ah terms from the first series and none from the
second would be required.
The loading consists of upon pressure P inside the while those in region B do (either Eo or 0). From a
hole. Thus, the only non-zero node point loads (applied consideration of the anti-symmetry of the deformation
in the r direction at nodes 1 and 2) correspond to the about the lines 0 = 0, 1ri4* +.2n it is easily seen that the
zero order term of the first series of eqn (4); their values only non-zero terms in eqn (4) are ml = 0,4,8.. .
are: The boundary conditionsfor St. Venant’storsion con-
sist of specifiedzero displacementsin the r and 19direc-
tions for nodes 1,3,5.. .21, specified8 ~splacements (of
Fr = FQ= (2~~)(~)(~) = 7&P,for m = 0. (53) 4 r) for nodes 2,4,6.. .22 and a zero z~ispla~ementat
node 1 to prevent rigid body motion.
The specified boundary conditions consist of the The finite element analysisgives values of 78Zat values
specified radial loads at nodes 1 and 2, and to prevent of ~9specifiedby the user. Numericalintegrationover the
rigid body motion, rollers along the bottom of the body cross-sectionalarea of the shaft gives the total torque T
and ue = 0 at nodes 29 and 30. required to produce the specific twist, i.e.,
Fire 3 shows the predicted vahtes of the ma~um
hoop stress for two finite element grids and for varying T= tZT&r, 8) dr de. (54)
numbers of Fourier terms; also shown is the analytical flA
value. The variation of u@along the edge of the hole for
two different iinite element grids is shown in Fig. 4. For The torsional constant rrlG# is expressed as k,(2ar and
the 24 element grid the maximumerror is less than 2% is evaluated for the four grids shown in Fig. 6 and is
while for the coarser grid it is less than 4%. tabulated in Table 1; the exact value is taken from [12].
For all four of the analyses two Fourier terms were used
4.2 ~xarnp~~ 2. Torsion of a square shaff in the expansionfor displacements.The predicted values
A square shaft, Fig. 5(a), subjected to torsion is con- are in good agreementwith the exact value with the error
sidered to further ilhtstrate the capabilities of the for the last two grids beingless than 2%.The variationof
method. The shaft is idealized as a nearly axisymmetric the predicted shear stress at 0” and at 45” compared to
body as shown in Figs. 5&c). Elements in region A the theoretical values are plotted in Figs. 7 and 8 for all
have no variation of properties in the theta direction four analyses.

- Theoretical Result
x 14 Elements Grid, 5 Fourier Terms
l 24 Elements Grid, 5 Fourier Terms

0 10 20 30 40 50 60 70 60 90 -8

Fii. 4. Comparisonof the exact and finite element predictionsfor u8-Exmple 1.

396 M. SEDAGHA~
and L. R. H~;KKMJ~N

Table 1. Predicted values of the torsional constant (exact value 0.1406) and number of itrratlan\ requjred-
Example ?

Finite element Predicted Number of iterations

grid number value of k, to converge

I .I098 5

2 .134x 8

3 .I430 8

4 .13x7 X

A Region A Reglon B

rlc
I
T

l--a----l
I
I

Grid I

(a) Cross Section of Actual

Shaft Subjected to s
Torque T Grid ::2

Grid LJ3

(b) Cross-Section of Idealized ‘Nearly’

Axisymmetric Solid
Grid 1: 4
Fig. 6. Four finite element representations of the square shaft of
Z
Repion A Retion B Example 2.

.
aI 22
_L2, 4 vergence problems were encountered when the number
1.0
I I llllll 21
_r of Fourier terms was increased to three. It is thought that
the presence of the artificial discontinuities introduced
T1 3 into the problem, by the approximation of the properties
(c) Finite Element Grid in region B as piece-wise constant (see Fig. 9), is the
Fig. 5. Problem configuration and Finite element idealization-
cause of the poor convergence. The artificial discon-
Example 2. tinuities appear to magnify the importance of the higher
order Fourier terms. Thus, in this case the inclusion of
additional terms in the analysis could actually decrease
It should be noted that example 2 is not the type of the accuracy of the analysis, i.e. the solution would begin
problem that is best suited to the method since the area to approximate the behavior of the fictitious shaft (Fig. 9)
of interest (region of maximum stress) is at the edge of rather than the real one (Fig. 5). What is clearly indicated
the body which is the region most affected by the near is that an increase in number of Fourier terms should be
axisymmetric representation. Thus the accurate predic- accompanied by a corresponding increase in number of
tion of shear stress is quite encouraging. However, con- elements to better approximate the actual geometry.
 Analysisfor nearly axisymmetricsolids 397

t
- Tbcmtical Values

o Grid 1
5.0
o Grid 2
l Grid 3
4.0 /* 1 Grid 4

tTa37gr
3.0
--I T 0”
/O

2.0 *6

1.0

/RX0 , I I I
0 )r la
0 0.5 1.0
Fig. 7. Variationof rer at 0”vs r-Example 2.

- Theoretical Values

0 Grid 1

r/a

Fig. 8. Variation of r8: at 45”vs r-Example 2.

4.3 Example 3. Stress distribution in an orthotropic plate For analysis purposes a circular region (radius R) of
with a circular hole the plate sufficiently large so that the disturbance of the
An example of an infinite orthotropic plate with a small hole is not felt at the outer boundary is considered.
circular opening subjected to a uniform stress uX= P at A ratio of R to a of 10 was found to be appropriate (see
infinity is shown in Fig. 10. Lekhnitskii et a1[13] give a Fig. 11).The outer edge of this zone is subject to stresses
complete analytical solution for this problem.

7 I

Fig. 9. Approximation to the square shaft as defined by piece-

wise constant properties.
r[
Fig. IO. Orthotropicplate with circular hole-Example 3.
17
398 hf. SEDAGHATand L. R. HERRMANN

2 3 4 5 6
- I a

QFinit@ ElementResults
.25-

.M l

.
l---R--4
Fig. 11. Finite element idealizationof Example3.
Fig. 12. Stress concentrationnear the circular hole-Example 3.

4+os(2H) material however is assumed to be isotropic, nonlinear

elastic. For illustrative purposes, the following simple
strain dependence of the elastic isotropic parameters is
TM= - f sin (28) (55) assumed

which correspond to a uniform tensile stress of P in the

(57)
x direction.
The following orthotropic properties for pfywood are
taken from the Wood ~andbook[4], units are in
lo-” inZ/lb (see Section 3.4).

0.83 - 0.06 -0.06

l.b67 -0.06 ii 0 0
[a]-’ = 1.5 0 0 0
(56)
(symm) 20.0 0 0
14.25 0
15.0

It is assumed that the plywood is oriented with principal where

axes of x, y, and z as shown in Fig. 10.
Six terms are incIuded in the Fourier series expansion
and eleven elements are used to represent the r-z cross-
section. Convergence occurred in 7 iterations. Predicted &, = e, - h. (58)
deformations (multiplied by the factor lO%aPD;:)) at
the hole are compared to theoretical values in Table 2. The maximum and minimum principal strains are
The stress concentration at the hole is ~ustrated in denoted as B, and e3. For computations purposes the
Fig. 12. To estimate the value of the hoop stress (divided following numerical values are used:
by P) at the edge of the hole and 6 = O”,linear extrapoli-
tion is used and gives -0.684 which compares favorably E&a;, = 750, v<,= 0.3. (59)
with the theoretical value of -0.706. At 90” the finite
element value is 5.82 as compared to the theoretical
value of 5.57. The tensile stress P is taken to be (r,/2 and a value 11a
was selected for R to define the solid of revolution.
4.4 Example 4. Analysis of a nonlinear efastic plate The boundary stresses and Fourier expansion used in
The plate shown in Fig. 10 is again considered. The the previous example are still valid. The base properties

Table 2. Compa~sonof the defo~ation at points A and B-Example 3

Deformation Theoretical Finite element

of point PZSUltS results

A 5.472 5.456

B -1.420 -1.423
 Analysisfor nearly axisymmetricsolids 399

are computed using En and ZJ,,.At each iteration, after 4.5 example 5. Three-~imensionu~ a~af~~is of a can-
the first, strains are computed, principal strains are fonnd ti~@ve~ beam
for all elements, new E and v values are evaluated, AU of the previous examples are two-dimensional in
deviation properties are computed, the integrals of eqn nature and thus, the semi-analytical method only
(29) are computed, new load matrices are formed and a required one-dimensional finite element representations
new solution estimate is obtained. (i.e., a single row of elements). To fully illustrate the
No analytical solution is available for this problem, usefulness of the mod&ad semi-anaiytical approach, a
hence a two~im~~sional, nonlinear, plane stress kite thr~dimension~ problem is analyzed. The problem
element program is used to analyze the problem in order considered is selected not for its inherent value nor
to provide results for comparison. The program uses because it is particularly appropriate for the modified
successive substitution to handle the nonlinearities. The semi-analytical approach but because a simple ap
grid used for the two-dimensional plane stress analysis is proximate analytical solution (beam analysis) exists for
shown in Fig. 13, 132 elements are considered to be comparison. The cantilever beam of square cross-section
adequate. By contrast 14 elements and four Fourier of Fig. 14 subject to an imposed tip deflection S is
terms were found to be su~cient for the semi-an~~i~al analysed. This nearly axisymmet~c body is idealized as
finite element analysis. shown in Figs. S(b) and 15. For clarity the r-t cross-
The convergence is slow in both cases, each analysis section of the idealized body, shown in Fig. 15, is not
requiring 14 iterations for convergence. A comparison of drawn to scale. The grid shown in Fig. 15 includes 70
the hoop stress, normalized with respect to P, at three elements, the 20 elements in the right-hand two columns
points around the hole is shown in Table 3; overall the have material properties that are theta dependent similar
comparisons are quite satisfactory. to those in region B of example 2.
Recognizi~ the axes of deformation symmetry and
antisymmetry shows that only the odd terms of the first
series of eqn (4) are present; 3 terms were included in the
154 155 154 153 152 151 150 analysis.
The convergence of the solution is very fast, three
iterations being sufficient to meet the convergence cri-
terion. Figure 16 shows the transverse shear stress dis-
t~butions (at z = 5a and y =O) predicted in the first
iteration and after convergence has occurred; the later is
compared to that given by beam theory.
The bending stress distributions for the same location
are shown in Fig. 17. The cantilever beam is analyzed a
second time for the case of a load P applied to the tip of
the beam. The predicted detIection at the centerline vs
that given by beam theory is shown in Fig. 18. The
results are particularly encouraging when it is considered
that only three iterations and 3 Fourier terms were
needed for the analysis. This example was run on a
Burroughs 7800 computer (the other examples were all
run on a LSI-I l/23 microcomputer) in order to compare
1 14 27 40 53 66 l9 92 I05 116 l3l 144 the cost with the anticipated cost for a regular three-
Fig. 13. Grid for the two-dimensional finite element analysis of dimensional finite element analysis. The semi-analytical
Example 4. analysis using the 70 element grid shown in Fig, 17, 3

Table 3. Comparison of q)P at three points around the hole-Example 4

Location of Tw~imensioMl Sf2mianalytiC~l

point 8%3lpiS analYSiS

O'L -1.06 -1.11

450 .95 .91

900 1.41 1.49

y /------F-lord
Fig. 14. Cantilever beam subject to tip deflectio=Example 5.
400 M. SEDAGHATand L. R. HERRMANN

- Beam Theory Predtctlon

x/a X Finite Element Predictron After
Convergence Has Occured
! 0 Finite Element Predrction
After First lteratron
81 82 83 64 85 86 87 88 1

Fig. 16. Transverse shear stress distribution at midspan-

Example 5.

1 2 3 4 5 6 7 8 -r
3EI
i 0.5a 0.207a p13”x
fl I- I y=x=o

Fig. 15. Finite element grid for cantilever beam-Example 5.

I
Q 0.5 1.0

terms in the Fourier series and 3 iterations to con-

-0.5 -
vergence required 1.57 minutes CPU time, the total
charge including memory and printing was less than
$4.00. Although no actual figures are available for a -1.0 -
conventional three-dimensional finite element analysis of 0 Finite Element Prediction
such a structure, experience with large two-dimensional - Beam Theory Prediction
problems would indicate that it would be at least an Fig. 18. Prediction deflection at centerline-Example 5(b).
order of magnitude more.

5. CONCLUSIONS plication of the conventional semi-analytical finite ele-

By partially relaxing a number of the assumptions ment method has been extended to include a much larger
required to formulate the method, the range of ap- class of problems.

xia
l Finite Element Prediction
- Beam Theory Prediction

2P2
3E8a”z I z=E/Z,y=O

Fig. 17. Bending stress distribution at midspan-Example 5.

 Analysisfor nearly axisymmetricsolids 401

The requirement that the material properties be con- 4. E. L. Wilson, Structural analysis of axisymmetric solids.
stant in the circumferential direction, is completely AZAA 3.3,2219-2274 (1%5).
removed, thus the range of application is extended to 5. J. A. Stricklin and.J. C. De Andrade, Linear and nonlinear
include completely non-homogeneous and/or, orthotropic analysis of shells of revolution. Pfoc. 2nd Conf. Mahix
Mehods in Strut. Me&., Air Force Inst. Tech., Wright-
(rectilinear principal axes) and/or nonlinear materials. By
Patterson A. F. Base, Ohio (1%8).
partially relaxing the requirement that be geometry be 6. L. R. Herrmann, Three-dimensionalelasticity analysis of
axisymmetric the range of application is further extended non-axisymmetricallyloaded solids of revolution. Dtpt. of
to include nearly axisymmetric solids. Civil Engin., Rep. 68-12-1, Univ. of Calii. Davis Campus
Because the size of any one of the system of equations (1%8).
that must be solved is relatively small the analysis may 7. C. G. Pardoen, Improved structural analysis technique for
be performed on small computers including microcom- orthogonalweavecarbon-carbonmaterials.AIAA J. 756-761
puters. The resulting analysis is computationally very (1975).
efficient and permits the practical analysis of a large class 8. C. G. Pardoe.n,A. D. Falco and J. G. Crose, Axysymmetric
of two- and three-dimensional linear and nonlinear prob- stress analysis of axisymmetric solids with rectangularly
orthotropicproperties.AUA J. 14,1419-1426(1976).
lem with only modest computer facilities. 9. 0. C. Zienkiewicz,ThePi& E/ementMethod in Engineering
Science. McGraw-Hi, New York (1971).
10. B. Camahan, H. A. Luther and J. 0. Wilkes, Applied
Numerical Methods, Oxford, Wiley, New York (1%9).
RWERENCES 11. G. N. Savin, Stress Concentration around Holes. Per8amon
1. Y.K. Cheung, The 6nite strip method in the analysis of Press, Oxford(l%l).
elastic plateswith two oppositesimplysupportedends. Proc. 12. S. Timoshenkoand J. Goodier, Theoryof Elasticity, 3rd Edn.
Inst. Cio. Engng 40, l-7 (1968). McGraw-Hi, New York (1970).
2. 0. C. Zienkiewicz and E. Hmton, A note on a simple thick 13. G. S. Lekhnitskii,(Translatedby S. W. Tsai and T. Cheron),
finite strip. Inr. J. Num. M&k Engng 11,905-909 (1972). Anisotropic Plates. Gordon and Breach Science Publishers
3. 0. C. Zienkiewicz,The finite prism in analysis of simply (1%8).
supportedbridge boxes. Proc. Inst. Ciu. Engng 53, 147-172 14. Wood Handbook, U.S. Department of Agriculture Hand-
(1972). book, No. 72 (1955).