.

LAMS-2745

CIC-14 REPORT COLLEGti~i ‘

REPRODUCTION

COPY

LOS ALAMOS SCIENTIFIC LABORATORY

OF THE UNIVERSITYOF CALIFORNIAo LOSALAMOS NEW MEXICO

... . . -. .—. _=.,—-

.

.

FOR

..-—

.— ,.. . .=. .

.-

.—

.-

DIFFERENCE EQUATIONS P“LANE THERMAL ELASTICITY

—. .—. —.- ..-. — .— . . . ..- .

.

. . .

.———a— ,. .

—

...

— ..—

.- ,.-
t ... . ...–. -. .-

LEGAL NOTICE

Thisreportwas preparedas an accountofGovernment sponsoredwork. NeithertheUnitedStates,northe

Commission,noranypersonactingon behalfoftheCommission:

A. Makesanywarrantyorrepresentation, expressed

or implied,withrespectotheaccuracy,completeness,or

usefulnessoftheinformationcontainedinthisreport,or

thattheuseofanyinformation, apparatus,method,or process disclosedinthisreportmay notinfringeprivately

ownedrights;or

B. Assumes any liabilities withrespecttotheuse

of,or fordamagesresultingfrom theuseofanyinformation,apparatus,method,or processdisclosedinthisreport.

As usedintheabove,“personactingonbehalfof the

Commission”includesany employeeor contractorofthe

Commission,ofemployeeofsuchcontractor,totheextent

thatsuch employeeor contractoroftheCommission,or

employeeof such contractorprepares,disseminates,or

providesaccessto,any informationpursuantto hisemploymentorcontractwiththeCommission,orhisemployment withsuchcontractor.

PrintedinUSA. Price$1.75. Availablefrom the

OfficeofTechnicalServices

U. S.DepartmentofCommerce

Washington25,D. C.
-1 —.

I

LAMs-2745

MATHEMATICS AND COMPUTERS

TID-4500(17thEd.)

LOS ALAMOS SCIENTIFIC LABORATORY

OF THE UNIVERSITYOF CALIFORNIA LOS ALAMOS NEW MEXICO

REPORT WRITTEN August1962

REPORT DISTRIBUTED: October16,1962

DIFFERENCE EQUATIONS FOR PLANE THERMAL ELASTICITY

by

GeorgeN. White,Jr.

ContractW-7405-ENG. 36 withtheU. S.AtomicEnergyCommission

All LAMS reports are informal documeots, usually prepared for a speokl purpose and prhnarIly prepared for use within the Laboratory rather than for

genernl dfetrilaition. Tbls report has not been edited, reviewed, or verified

for accuracy. AU LAMS reports expre.ee the vfews of the autbore as of the

time they were written and do not necessarily reflect the opinions of the LoB

Alamos Soientiflo Laboratory or the final ophion of the autbore on tbe subject.

1

1

ABSTRACT

Difference equations for the displacement components of a heated,

inhomogeneous, orthotropic material are obtained from a variational

principle. Including boundary terms, the matrix of the resulting

linear system of equations is positive semi-definite, a property which

assures solution by many iterative techniques. The equations are

solved for an infinite hollow cylinder of isotropic material in which

the temperature varies inversely with the radius, whose inner surface

is stress free, and whose outer surface is fixed.

3

INTRODUCTION

For many problems in plane thermal elasticity which cannot be

treated by classical analysis, a numerical method of general applicability

consists of replacing the linear differential equations and boundary

conditions by linear difference equations. The matrix of the system of

linear equations obtained by this direct differencing has usually such

a complex structure that the choice of a successful method for solving

the system is extremely difficult. In this paper it is shown that the

uncertainties of direct differencing can be avoided by deriving the

difference equations Nom a variational principle. For a wide class of

thermal elasticity problems, the variational

system with a positive semi-definite matrix,

methods to be used with confidence.

approach yields a linear

which permits iterative

The use to which a code based on the variational difference equations

was to be put required that the basic equations describe anisotropic,

inhomogeneous mterials in multiply connected regions. A formulation of

the therul elasticity problem in terms of the displacement vector met

5
these requirements as well as allowing boundary conditions of all types

to be imposed in a simple manner. The Ufferential eqwtions and

boundary conditions in terms of displacement components are derived by

the variation method in the first section of this report. The second

part contains a discussion of the stress-strain relation for orthotropic

materials, a type of anisotropic substance of sufficient generality for

the purposes of the code. The third part is a derivation of difference

equations from the variational principle. The last section contains a

comparison of an example calculation with the exact solztion for the

stresses and displacement in an infinite hollow cylinder of isotropic

mterial in which the temperature varies inversely with radius.

6
I. THE BASIC DIFFERENTI&L EQUATIONS

For therml stress problems a variational principle based on the

Helmholtz free energy is a convenient means of obtaining both the

differential equations

ment components.

The free ener~ F

and the difference equations for the displacers assumed to be a function of the temperature

T and the strain tensor e. .(i

lJ ,j s 1,2,3). From the combined first and

second laws of thermodynamics for isothermal processes, the rate of

change of the free energy eqwls the ne~tive of the rate at which

work

From

that

W is being done by the stresses cr.:

lj

dF dW —=- —=

deiiox

dt dt “ij dt

(I 1)

the assumption of the dependence of F on T and e.., it is evident

lJ

*A1l vectors and tensors w-jll be referred to a rectan@ar Cartesian

coordinate system. The suumation convention for repeated indices will

be employed; e.g.,

a. .b lJ ik = aljblk + a2jb2k+ a5jb5k.

7
a

ij ‘%ij)T” (I 2)

Consider the following variational principle for the case of no body

forces:

‘Jmv -$; ‘“im = 0)

(I 3)

v

St

where the variations of the displacement components Ui are carried uut

isotherxmlly andT~

of the surface S’.

yields

o

are stress vector components

Integration by parts and the

prescribed on a ~rt

use of Greenis formula

= “0* (I 4)

s -s1 v

The surface integral is split as follows:

J

Ti 6Ui ~ =

J

Ti 6U1 ~ +

J

Ti &l~ ~ +

J

T; ~U~ ~. (1 5)

s d Sy St

The first integral is over that part of surface where the Usplacements

are specified and vanishes with the bui appropriate there. The second

integral, over a symmetry boundary on which no work is done, also

vanishes. The third part is cancelled by the corresponding term in

(I k). Thus the Equation (I k) reduces to

* IY.

2

1 ?s

% >s Xj
J (I 6) aij,j

5ui dV = O.

v

For suitable functions a

ij and arbitrary bu., we obtain by standard

1

arguments the equilibrium equations,

aij,j

= o. (1 7)

The boundary conditions are the natural conditions for surface

tractions,

Ti = T;(%P),

(a,@, are surface parameters); specified displacements,

u. =

1

U;(%P);

and symmetry conditions,

Tf &2~ = O.

(I 8)

This last condition canbe simplifiedby observing that for any symmetry

boundary the norml component of the displacement must vanish. Since

the tangential component of the displacement is still arbitrary, we can

conclude that the tangential component of the stress vector must also

9
vanish if no work

symmtry boundary

natural.

The e~ation

is to be done on a symmetry surface. Thus the

conditions are tied - - partly specified and partly

for displacements is next found from (1 7) by =w=S-

‘ng‘he‘tresses ‘ij

in terms of the displacements through a stressstrain relation. For therrralstress problems in which the temperature

changes and strains are not large, the free energy can be approxinmted

by an expansion(1)

The

the

The

F(T,eiJ) = F(TO) + Cij eiJ(T-To) + l/2cijkj eij ekj. (I 9)

coefficients Cij are anisotropic coefficients of thernal e~nsion;

cijk~ are anisotropic elastic constants. The stress from (I 2) is

(T-TO) + cijkj~~” (I 10) aij = Cij

equilibrium equations (I 7) my now be written in terms of the displacements,

“ijk~ %@) ,j i- [Cij(T-To)l,j= O.

The boundary condition for a specified stress vector can also be

(I U.)

written in terms of the dii.splacements.For an outward nornal with

components n , the stress vector components are

i

10
Ti=a c n.(T-To) + c

ijnj=ij J ijk~ W,4 ‘j “

(I 12)

For symmetry boundaries the normal component of the displacement must

vanish,

u =uin i= o,

n

(I 13)

and the tangential component of the stress vector must also vanish,

‘t

= a.. n. t = O, lJJi (I 14)

where ti is any vector lying in the tangent plane. If the material is

isotropic, it can be shown that (I 13) and (I 14) together in@y that

aut

&T=o” (I 15)

If the =terial is anisotropic, only (I 14) can be used.

Differential Equations (I XL) plus the boundary conditions (I 8),

(I 12), (I 13), and (I 14) constitute a displacement formulation of

the thermal elastic problem for anisotropic ~terials. In the next

section the equations will be specialized for a mterial with orthotropic symmetry and for the case of plane strain.

11
II. THE STRESS-STRAIN RELATIONS

Several as_pectsof the stress-strain relations (I 10) need to be

discussed. First the form of the relation specialized for materials

with orthotropic symmetry is to be determined. Also the positive

definiteness of that part of the free energy quadratically dependent

on the strains must be established. Finally the inequalities relating

the elastic constants are to be obtained as conditions for positive

definiteness.

A material is said to have orthotropic symmetry if under reflection

through any of three orthogonal planes the mterial tensor components,

do not change their values. cijk~’ In a single reflection a component

of a tensor has a change in sign if the index of the effected coordinate

appears an odd number of times. Since all coordinates m.y be reflected

in an orthotropic material, all components with an odd nuniberof indices

must vanish if the tensor is to be unaffected by reflection. Taking

this into account and usin~ coordinate lines which are the intersection

of the symmetry planes, one has for the free energy

12
F=

+

+

+

+

+

+

2

*[cuu eu

+ C1122 ’11 ’22 + c2211 ’22 ‘XL + cl.133’11 ’33

c3311 ’33 ‘u + c1212

2 2

’12 + c2112 ’21 ’12 + c2121 e a

c1221 ’12 ‘a + C2222

2

C2323 ’23 + C3223 ’32

2 2

C3333 ’33 + C1313 ’13

2

’22 + C2233 ’22’33 + C3322

’23 + C3232 ‘;2 + C2332 ’23

’33’22

’32

2

+ C1331 ‘13e31 + ‘3131 ’31

c3113 ’31’13 I+C ~ eU(T-To) + C22 e22(T-To)

C33 ’33(

T-TO) +F(TO). (II I)

The symmetry of the strains and elastic constants reduces (II 1) to

F .+[C 2 2

111.1‘1.l+ C2222 ’22 + C3333 ‘;3 + 2c1122 ’11 ’22

+ 2C1133 ‘u. ’33 + 2C2233 ’22 ’33 + 4c1212 ‘:2+ 4=1313 ‘:3

2]+C

+ 4C2323 ’23 ~ =ll(T-TO) + C22 e22(T~o) + c33 =33(T-TO)

+ F(TO).

(II 2)

The stresses obtainedby (I 2) from (II 1) are

13
%3.= cm-l ‘I-1.+ C1122 ’22 + c1133 ’33 + CU(T-TO)

’22 = C1122 ’11 + C2222 ’22 + C2233 ’33 i-C22(T-TO)

a33 = cl..133’11 + C2233 ’22 + C3333 ’33

-i-C33(T-TO)

= 2C

a12 1212 ’12

=2C

a13 1313 ’13

.2C

a23 2323 ’23”

(II 3)

These equations take a more familiar form when solved for the strains

in the following manner:

‘1 ‘u = all - ’12 a22 - ’13 a33

+ El (xl(T-To)

‘2 ’22 = -V21 all + a22 - ’23 a33 ‘E2

‘3 ’33 ‘ ‘V31 an -’32 ’22+ a33 ‘E3

a12

’12 = —2V12

a2(T-To)

a3(T-To)

.2.2,

’23 %?3

14

(II 4)
The Ei and v

ij are Young’s moduli and Poisson’s ratios, the V. . are

lJ

shear moduli, and the (Liare coefficients of therml expansion. By

regarding the c

ijk~ as coefficients obtained from solving (II 4) for

the stresses, one finds

Cllll =

C1122=

C1133 =

C2222 =

C2233 =

C3333 =

C1212 =

C2323 =

C1313 =

‘~(1-v~3v32)/A

‘2(v12+vl#52 )/A

313 12v23)/A

E (V +V

‘2(1-v~3v3~)/A

‘3(V23+V21V13)/A

‘3(1-v~2v2~)/A

’12

’23

h3

= E1(V21+V31V23 )/A

= E1(V51+V21V52)/A

= E2(V32+V12V31)/A

-1

cll = ‘f~ ~(1-v23v32) + ‘2 a2(v12+v52v15)

c22 = +1 al(V~+V31V23) + E2 ~2(1-V13V31)

+ ‘3 ‘3(V13+V23V12))/A

+E cZ(V )}/ 3 3 23+V21V13 A

-1

C33 = -t ‘l ”l(v31+v21v32) ‘E2a2(v32+v31v12) ‘E3 a3(1-v12v21)j/A

A=l-v

12’21 - ‘13’31 - ‘23V32 - ‘12’23’31 “VZLV32V13° (II 5)

15
Symmetry of the coefficients in the equations for the stresses, (II 3),

implies symmetry of the coefficients in (II 4) as follows:

’12 ‘2 =

V~E1> V13E3

‘12V23V31 = ‘21V32V13”

Since the properties of

= ’31 ‘1’ ’23 ‘3 = ’32 ‘2’

(II 6)

the differential and difference equations

for the displacement components, (I 11), are governed in part by the

magnitudes of the elastic coefficients, it is important to find

additional information about these coefficients. The requirement of

stable thermodynamic equilibrium in terms of the free energy provides

the needed facts. For the free energy given by (II 2) the requirement

of stability is met if the quadratic part of F is positive definite.

(2)

Necessary and sufficient conditions

quadratic form are the positiveness

of the 6 x 6 mtrix associated with

cijk~ are

for positive definiteness of the

of the principal

the form. These

(a) cl~~ c~22y C3333‘ C1212’ C1313‘ C2323 >

2

(b) Cllll CW22 - cll~ >0>

2

CKU1 C3333 - c~33>o’

2

C2222C3333 - c2233>0’

minor determinants

conditions for the

o,

16
(c) ‘~l ‘2222 C3333 +

The three inequalities in (b)

and (c) is independent of (a)

2

2c1122 C2233 cl.133“ Cllll C2255

2

C2232 cu33 - C3333 CL22 > 0“

=e symmetric forms equivalent

and (b). The inverted form of

(II 7)

tocll>O,

the stressstrain

yields

relations (11 4) must also have a positive definite matrix and

the inequalities

1

- ’12 ’21’ 1 - ’23 ’32’ 1 - ’13 ‘31>0’

A>O. (II 8)

A necessary condition of positive definiteness is that the sum of the

coefficients of the xmtrix be positive. Application of this condition

to the first three equations of (II 4) with T = To implies that the

bulk Mdulus, k, defined by

be

‘1

1 ‘kk 1 1 1 .=— =. 4LL L__..’21 ~-k-k>.

k -P ‘1 ‘~2+fi El El E2

3 ‘3 ‘3

(II 9)

positive and therefore

*2(1-v31-v32) + El *3(1-v=-v23) +E2E3(1-v12-v13) >0. (II 10)

17
Applied to the second order principal minors, the same arguments give

‘i + ‘j ~ ‘ij ‘j ~ ‘ji Ei >0 (i, j not summed). (II IQ

Inequalities (II U) in conjunction with the symmetry conditions (II 6)

give inequalities for the individual Poisson’s ratios in terms of the

Young’s moduli,

~> ]Vijl.

i

(11 12)

From (II 6) it should be noted that v and v. always have the same sign.

ij Ji

If the orthotropic elastic constants of a particular material in

terms of the c.ljkklor E ) V. J v. i

satis~ the inequalities (II 7) or lj lj

(II 8), it can be concluded that the material tensor is positive definite.

If the constants do not satisfy the inequalities, the constants are

thermodynamically inconsistent and may not have been determined correctly

by experiment, or the model which assumes the free energy to have the

form (II 2) may

mterial cannot

be found.

As is well

not be physically adequate. In the latter case the

be orthotropic, and a new form of the free energy must

known, positive definiteness

sufficient for the uniqueness within a rigid

value problems of elasticity.(3)

of the free ener~ is also

body mtion to the boundary

The differential equations (I Xl), for the displacement components

18
in an orthotropic material can be written explicitly

(Cllll’-h,l),l +

+ (c ) 1122 ‘2,2 ,1

+ (C1133 ‘3,3),1

(c~212ul,2) ~+

9

+ (c ) 1212 ‘2,1 ,1

+ (c ) 2233 ‘3,3 ,2

(c )+ 1313 %,3 ,1

+ (c ) 2323 ‘2,3 ,2

+ (C1313 ‘3,1),1

(C1ZL2 %,2),2+ (C1313 %,3),3

+ (c ) 1212 U2,1 ,2

+ (c ) -I-(C1l[T-TOI),l = 0, 1313 ‘3,1 ,3

(c ) 1122 U1,1 ,2

+ (c2222 U2 a) a + (C2323 U2,3),3

>)

+ (c ) + (C22[T-TOI),2 = 0) 2323 ‘3,2 ,3

(c ) 1133 ‘1,1 ,3

+ (c ) 2233 ‘2,2 ,3

+ (c ) +(C 2323 ‘3,2 ,2 )

3333?>3 ,3

+ (C33[T-TON3 = 0= (II 13)

If Equations (II 13) are written in the general form (I 11), the

terms containing the highest derivatives are

Cijkl?h,~j” (II 14)

19
Since the quadratic form associated with the coefficient mtrices is

the free energy, which is assumed positive definite with respect to

the strains, and if ri#d motions are excluded, the form is also positive definite with respect

The differential eqwtions

elliptic.(4)

to gradients of the displacement components.

associated with this form are therefore

In the case of plane strain, which will be investi~ted in detail,

the

The

The

strains satisfy

‘i3 = o.

stress-strain

= Cmlul,l +-CD( T-TO),

’11 + CXL22U2,2

+ CW(T-TO) ,

’22 = C1122 %,1 + C2222U2,2

(II 15)

relations for an orthotropic nm.terialreduce to

’12 = c1212%,2+u2,1) “

equilibrium equations in terms of displacements

(II 16)

are

m
(cull %,1),1+ (c1212 %,2),2

+ (c ) + [C1l(T-TO)],l=O, 1212 U2,1),2+ (C1122U2,2 ,1

) + @l122 %,1 ,2 (’1212 %,2 ,1 )

+ (C1Z2 ‘2,1))1 + @2222 ‘2,2),2 + [C22(T-TO)1,2 = 0. (II 17)

The Helmholtz free energy is for plane strain

F= +(cllll.e~+ c2222e~2+ 2cll=eU e22+4c1a2e~2)

+ ClleU(T-To) +Cme22(T-To) +F(TO) (II 18)

In the next section the free energy (II 18) will be used to obtain

difference equations which approximate the differential equations (II 17)

and the boundary conditions (I 8), (I 12), (I 13), and (I 14), appropriately modified for the plane case.

a
III. THE DIFFERENCE EQUATIONS lKIRPIANE PROBIEMS

The difference equations for the plane problems will be derived

froma variational principle equivalent to (I 3). In addition to the

advantage of giving a positive semi-definite matrix, this approach

offers a uniform

points. For the

method for forming the equations for the boundary

current IBM 7090 code the boundary of the plane region

is approximated by straight segm?nts consisting of ~id lines as shown

in Figure 1. Future codes will approximate the boundaries by straight

lines which are not necessarily grid lines. Difference equations for

both boundary treatments till be presented although the current scheme

will be examined in more detail. As before the coordinate surfaces are

taken to be parallel to the planes of symmetry

mterial.

Using Greek indices to signify integers 1

variational principle as

in the orthotropic

or 2, we can write the

JT[ t51

Z%@7b ‘a,@”7,5 + ca@ua,p IJ (T-TO) dA - b TauadS = O,

R s

(III 1)

22
where the first integral

is over that part of the

is over the entire plane region and the second

boundary curve on which the stress vector is

specified. The coefficients of the material appearing in (III 1) are

those in (II 16) and (II 18). The variation, as before, is mde on

the displacements at constant temperature.

To obtain difference equations from (III 1) it suffices to approximate the integrals by sums and to approxi-te the derivatives in the

integrand by differences. Application of the necessary conditions for

a stationary value of the sum with

@elds the difference equations.

For the first scheme in which

respect to the displacement components

only square mesh cells appear, the

free energy of a single

(i+l,j), (i+l,j+l), and

properties evaluated at

cell with corner points denoted by (i,j),

(i,j+l) is a basic quantity. With material

the point (i)j) and with the derivatives at

(i,j) obtained from fitting Ua bya plane throu@ (i,j), (i+l,j), and

(i,j+l), the energy in the

by

cell with sides of length h is approximated

2C

{ ‘Y J

Uli+l.-uli. U,i.+l-u,i.

l122iJ h

‘) h

/. /

‘2i+lj-u2ij.2

,h, )]

(

Uli+l.-u

lij’’)(T-T)

[ )

‘2ij+l-u2i~

+ Clli h o ij + c22i h

(T-To)i~k2.

/

(III 2)

23
Because the mterial tensor is positive definite with respect to the

strains, (III 2) is semi-definite with respect to the displacement

components. The energy in this cell could as well have been approximated

by quadratic forms evaluated at the other corners (i+l,j), (i+l,j+l),

and (i,j+l). A symmetric and still positive semi-definite form is

obtained by averaging these four contributions to give the energy

centered at (i+$,j+$)

1

[

)2+ ac~aaij(%i+lj-%ij %-+ j~ 1 = 2 Cllllij(uli+lj-%ij )%ij+l-%ij)

+

+

+

-1-

-1-

i-

-t

(u c2222ij 2ij+l-u2ij)2+C -i+lj (uli+lJ-ulij)2

2c1122ii-lj(%i+lj‘“lij)(u2i+lj+l-u2i+lj)

c2222i+lj(u=+lj+l-u2i+1j )2 + Clmi+lj+l(%i+lj+l-%i+l)z

( 2c1122i+lj+l %i+lj+l%ij+l )(u2i+lj+1-ua+lj )

(u c2222i+lj+l 2ii-lj+l-u2i+lj)2+ c ( X1.31ij+l%.i+lj+l-%ij+l )2

2c1122ij+l(%i+lji-1-%.ij+l)‘“2ij+l-u2ij)

(u c2222ij+l 2ij+l-u2ij‘2 + cl~2ij(~ij+l-~ijw2i+lj-u2ij ‘2

24
-1-

+

+

(“li+lj+l-%i+lJ )2 c1212i+lj “+u2i+lj-u2ij

(u c121*i+lj-1-lli+lj+l-%.i+lj+u2i+lj+l-u2ij+l

)2

\*1

(u cl~2ij+l lij-i-l-5ijw2i+lj+l‘“2ij+l~ J

. .) +, + ‘[c~ij(%i+lj-’&J

+

+

+

h[cui+lj (%i+lJ ‘-~~j)

(u .)](T-To)ij c22ij 2ij+l-u2iJ

.)l(T-To)i+lj

+ c22i+lj(“2i+lj+l-u2i+lJ

h[c ( lli+lj+l %i+lj+l-%ij+l )

h[c ( llij+l %i+lj+l-ulij-t-l

)i-

(U .)l(T-To)i+lj+l

+ c22i+lj-i-l2i+lj+l-u2i+lJ

(u

c22ij+l 2iji-1-u2ij)l(T-To)ij+l” (111 3)

The difference equations apart from the prescribed boun~ry terms are in

every case obtained by minimizing a sum of cell ener@es like ~~ j+”

It is convenient, therefore, to obtain the partial derivatives of

%+$ j+ with respect to ~ij andu2ij* A.ninterior Point (ij) ~s

neighboring cells with energi-es~~ j~) ~~ j~) ~~ j~~ ad

%+ j-$” Taking partial derivatives with respect to Uaij of these

four energies, we have eight basic quantities from which difference

equations can be formed in a systematic way. The eight partial

derivatives are

25
) - ‘11 ?Pij ‘“2ij+l-u2ij

.

*

-“

= ( “l.lJlij %i+lj ‘Ulij

)

ij

.( “llUi+lJ %i+lj-ulij )-c ~22i+lj(u2i+lj+l-ua+lj )

“l~zij (%_ij+l-%_ij+u2i+lj-u2ij)

“121_2ij+l(%ij+l-ulij+u2i+i ji-l-u2ij+l)

‘hcllij(T-To)ij -hclli+lj(T-To) i+lj~

. .

*

.

=

2ij

“1122ij(%ii-lj-’’i)i)- c2222ij(u2ij+l-u2ij)

( “1122ij-i-l%i+lj+l-%-ij+l

)-c 2222ij+l(u2ij+1-u2ij)

“l~2ij (%ij+l-%ij+u2i-t-lj-u2ij)

)

“1212i+lj (“li+lj+l-%i+lJ“+u2i+lj-u2ij

-hc22ij(T-T~)~j ‘hc22ij+l(T-T~) ij+l)

W

-- =

clll.lij(%ij-%i-lj ) + c~i-lj (%ij-%i-lj )

ij

) + C~22i-lj(u2i-lj+l-u2i-lJ “1122ij(u2ij+l-u2ij

.)

( “1212ij %ij+l-%ij+u2ij-u2i-lj )

( “1212ij+l %ij+l-ulijW2ij-l-l-u2i-lj+l)

‘hcllij(T-To)ij +hclli-lj(T-To) i-lj)

26
*

“-”

(u

= “2222ij 2ij+l-u2ij)-c 2222ij+1(u2ij+l-.u2ij)

2ij

( “1122ij %ij-%.i+lj )- c~122ij+#11ij-uli-lj+l)

“12121J..(Ulij-i-l-%ij+u2ij-u2i-lJ.)

(u “1212i-lj li-lj+l-uli-lj+u2ij-u2i-lj)

-hc22ij(T-To)ij -hc22ij+1(T-To) ij+l,

a%_l.. ~

( ) + cl~i-lj(%ij-uli-lJ =

Cllllij %ij-%i-lj .)

ij

“1122ij(u2ij ‘“2ij-l)+C ~22i-lj(ua-lj -u2i-lj-1)

( “1212ij “llij-ulij-l+u2ij-u2i-lj)

(u “l~2ij-l Iij-ulij..l+uaij-l-uai-lj-l )

‘hcllij(T-To)ij + cHi-lj(T-To)i-lj,

1

*

-“-

= (u c2222ij 2ij-u2ij-l)+C 2222ij-l(uaj-u2ij-1)

2ij

( “1122ij %ij-uli-lj )*C (u l122ij-1 lij-l-%i-lj-l )

“l~aij (~ij-ulij-l+uaij-uai-lj )

(

“l~ai-lj ~i-lj-uli-lj-l+uaij-u~-lj )

+hc22ij(T-To)ij + hc22ij-l(T-To)ij-1,

2’7
aq-@L .

-=& (is

“lllli j hi-l j‘hi j) - clllli+lj%i+lj%i j)

“1122i j

(u2ij-u2ij-1) - c112a+1j (Lla+lj-u2i+lj-1)

( “1212ij “lij-%ij-l+u2i+lj-u2ij )

“1212ij-l@li j-”lij-l+%+lj -1-”2ij-l)

-hcUij(T-To)ij - hclli+lj(T-To)i+lj ‘

1

%+

--

=

c2222ij(u2ij-u2ij-1) +

2ij

(U “1122ij li+lj-ulij

)+

(u c2222ij-1 2ij-”2ij-l)

c1122ij-l(%i+lj-l-%i j-l)

(III 4)

(u “1212ij lij-%.ij-l+u2i+lj-u2ij)

(u “1212i+lj li-i-lj-%ii-lj-l+u2i+lj-u2ij )

+hc22ij(T-To)ij + hc22ij-l(T~o)ij-1 ●

If the total free energy, Q, of the system is written as the sum

of all the individual cell energies,

(III 5)

the difference equations for the point (ij) are

28
(III 6)

unless the point is on

that can contribute to

derivatives (III 4) of

a line of symmetry. Since the only terms of Q

the expressions (III 6) are the eight partial

the energies in the four adjacent cells of an

interior point, the general difference equations for an interior point

(ij) are

(III 7)

To obtain natural boundary conditions (specified stress vector) it

is only necessary to note which types of cell ener~ can contribute

and form a sum similar to (III 7). For

boundary with the =terial on the left.

to the difference equations are %$j+

equations are, therefore,

where ~a are averaged components of the

boundary of the two cells.

For a corner point on the boundary

example, consider a vertical

The only energies which contribute

and ~~3~. The difference

(III 8)

specified stress vector on the

the cell energy contributing

might be ~~j ~. The difference equations would be

29
. 1

*

.- -

- T#h = o.

aij

(III 9)

All the other specified stress cases can be treated in a similar way.

The boundary conditions at an internal line of symmetry can be

obtained by minimizing the energy with respect to the tangential displacement vector only. As the most simple

vertical symmetry boundary. For this case

varied. The difference equation is simply

example consider a~in a

%.ij = 0’ andu2ij ‘s ‘0

be

Cl. J

‘hich‘ith ‘hij =

O form the boundary conditions.

(III 10)

For the more

an angle through

‘n=%%+

Variation of the

general situation where the symmetry line runs at

the mesh, the ~ and U2 components are related by

‘2n2 = 0 “

(III 11)

energy in a cell must take this constraint into

account. If% is taken to be the independent variable, only the

variation of the energy with respect to U1 is taken. Consider again

the case where the symmetry line is replaced by a segment of mesh lines

and the cell whose ener~ contributes terms to the difference eq~tion

is centered at i~j-$. The difference equation is then

30
(III 12)

The tangent vector at the point ij which approximates a boundary point

has components

‘1

= -n2, t2 = rj-.

Equation (III 10) is therefore

(III 13)

1

*

“-- ta = o.

ai j

(III 14)

Thus (III 11) and (III 14) are the approximate boundwy conditions for

the case of a symmetry boundary not coincident with a mesh

view of Equation (III 9) it is evident that (111 14) is an

to the natural part of the symmetry boundary conditions (I

line. In

approxhmtion

14)

‘t = aa#3n(3ta= 0”

Comparison by a Taylorts expansion of the interior difference

equations (III 7) with the differential equations (II 17) shows that

the error is of order h2. On

difference equations gives an

error must be ascribed to the

mesh cells are required to be

the boundary lack of symmetry in the

error of order h. Psrt of the boundary

distorted boundary which is used when all

square. In order to reduce this latter

error, a modification of the boundary treatment will now be presented.

31
The treatment proposed is based on an appmximtion of the boundary

by chords between points of intersection of the boundary curve with

mesh cells, Figure 2. In order to keep the number of possible configurations to a XELnimumit will be assumed that the mesh size is so smll

that the boundary curve does not pass through the sides of any cell.mre

than twice. That part of a boundary celJ.which belongs to the interior

is then a triangle, a qymlrilateral, or a pentagon.

The quantity needed for the difference eqmtions of points on or

adjacent to the boundary is the &ee energy associated with such

irregular shaped areas. The energy in a square celll_was calculatedly

averaging an energy associated with each of the four vertices. If each

of the four vertices is assured to account for the energy in a triangle

consisting of the vertex in question and its two neighbors, then the

sum of the energy of all four possible such triangles is twice that of

the whole cell. The sum of the triangular energies is therefore

divided by two. A similar approach is to be used for the irregular

regions.

The derivates of the displacements needed for the energy in a

triangular region can be obtained by passing a plane representingua

through the three vertices of the triangle. The displacements in the

plane can

u=

a

be represented

aij X+ll

ay 7 a

by an eqression of the term

(III 15)

32
where i and j are indices locating what might be called the pr~ry

vertex, and where (i+lj), (ij+l) are the other vertices. The derivative

of Ua at the point (ij) is

and

ij

= a~, ‘a,$

(III 16)

the ener~ associated with this triangle becomes

<j =

+

+

{[

aij aij ..

(ai3)2+ 2c~22ij 11 22

+ Cllll.ij 11 + c2222ij@:)2

4C1212ij‘a~~+a~)2 1

aij(T~o)ij +cllija~(T~o)ij }Aij,

Cllij 11

(111 17)

where A<q is the area of the triangle. Since there are three vertices

-w

to the triangle, a more symmetric form of

. tj + t+lj +4“+1 4+$+ 3

●

the

the

For a quadrilateral the treatment is

energies of the four triangles formed

vertices are averaged,

~Q = ij + i+-lj+tj+l +tl-lj+l

2

.

the enerfy is

(III 18)

the same as that for a square,

bypassing diagonals through

(III 19)

33
The pentagonal boundary region presents a larger number of possible

subdivisions for computing the energy tiom its component triangles.

From the coding viewpoint the simplest expedient would be to readjust

the approximate boundary so that all pentagons are altered into quadrilaterals. Next in order of complexity wouldbe a somewhat arbitrary

division of the pentagon into a quadrilateral and a triangle. A third

possibility, hating less bias, would be to obtain three triangles

associated with each of the five vertices by passing lines from the

given vertex to all the others. Fifteen triangular energies wouldbe

computed in this way whose sum was five times the energy of the pentamn.

If the vertices were labelled as in Figure 3, and the energies associated

with vertex one written q

123’ ‘134’ ‘145’

the total energy of the ce13.

becomes

‘T

qp =

1

T

1/5 2q~23 + q~~k + 2q;45+ $45 + 2q”45 + $15 + 2’$34

To obtain

one forms sums

+ 2q~45 + 423 1●

(III 20)

the difference

of terms like

equations for regions mde of triangles

(III 21)

34
This ecyzationin turn is formed from the derivatives of the energies

in triangles according to (III 18). Thus the basic derintives are

(III 22)

From elementary considerations the quantities a

*

appearing in the

Equations (III 16) and (III 17) are

‘aii-lj(x2i+lj-x2ij )

- ‘lij(x2ij+1-x2i+1j) /2Aij ,

~ij ‘ =Iu (x

22 I_ 2ij+l li+lj-xlij)-u (

2i+lj %ij+l-xlij )

(x

- ‘2ij li+lj-xlij+l)]/a~jJ

ij -

a12 = 1%ij+l(xli+lj-xlij) - ‘li-l-lj(xlij+l-xlij)

)]/aij> + hi j‘Xli+lj‘Xiij+l

1

ij -

aa =

.(x2ij+l-x2ij) - u2ij+l(x2i+lj-.x2ij)

‘2ij+J

- ‘2ij(x2ij+l-x2i+lj)]/aij . (Im 23)

35
l?rom(III 17),

+ cllij(T-To)i3

But

_..

‘Xli+l ) j-xlij+l ~ , =-

ai j

oil.

(x2i.+l-x2i+1.)

=-

2A.,

8

a2”

.LJ

36
Thus

{

[ .(x cllllij %i+lJ 2ij+l-x2ij)- %ij+l(x2i+lj ‘Xpij )

+

+

+

+

(x .)IIX ‘lij 2ij+l-x2i+lJ

]8

2i+lj-x2ij+l al

[ .(x c1122ij %i+lJ 2iji-1-x2ij)-u ~ij+l(x2i+1j-x2ti)

‘%ij‘x2ij+l-x2i+lj)]‘Xiij+l-%i+lj ] 8a2

..[U (x cl.1221J 2ij+l li+lj-xlij) - u2i+1j(x1ij+l-xlij)

( ‘2ij %i+lj-xlij+l )I[x2i+lj-x2ij+l1 bd

..[U (x C22221J 2ij+l li+lj-%ij ) - u2i+1j(x1ij+l-~ij)

( ‘2ij KLi+lj ‘x2ij+l)1[x~ij+~-xli+lj] bti

8C [ (x 1212ij ‘+lij+l li+lj-xlij) - ~i+lj(xli J+l-xlij)

%ij(%i+lj-~i j+l) +

(x ‘2ij+l 2i+lj-x2ij)-

‘[(~ij+l-xli+lJ .)5 ~ +

‘2i-i-lj(x2ij+l-xaj)

(x ‘2ij 2ij+l-x2i+lj)1

)5 ] (x2i+lj-x2ij+l @

+C ~ij (T-To)ij(X2i+1j-X2ij+l)bd

+ Cppij(T-T~)ij (~ij+l-%i+lj )5a J/=.jj)

37
%, ‘{[%i+lJ(x2ij+l-x2ij) - ‘lij+l(x2i-i-lj-x2ij)

- %ij ‘x2i~+l-x2i-t-lj )]“llllij ‘x2i+lj-x2ij+l)5al

(Xlij+l-xli+lj)b~l + ‘“2ij+l(xli+lj-xlij) + c1122ij

(. ) - ‘~ij(~i+lj-xpi~+l

- ‘2i+lj ll.ij+l-xlij )1

x [c~122iJ(x2i+1j-x2ij+l)5d + c2222ij(~ij+l-xli+1j )8@l

+ 8C~a2ij [ulij+l(xli+lJ-xlij)- uli+lj(xlij+l-xliJ)

(x

- ‘lij li+lj-%_ij+l)

(x

- ‘2ij+l 2i+lj-x2ij)

x[(xlij+l-xi+lJ .)bd +

+ Cllij(T-To)ij(X2i+1j-x2ij+1)ad

(x

+ ‘2i+lj 2ij+l-x2ij)

(x

- ‘2ij 2ij+l-x2i+lj)1

(x2i+lj.-x2ij+l)8@]

+ c22ij(T-T~)ij‘hi j+l-ai+lj )8a2}/2Aij. (III 25)

Expressions for a

~+lj/buaij ‘d b~j+l/~ai3 can be ob~fned in a

similar mmner to form the com@ete derivative of the cell ener~

%+j+ /

1 au It is apparent at this time that the price paid for a

aij”

more careful boundary treatment is not only a much greater complexity

in the coefficients of the difference equations but also a proliferation

38
of possible combinations of cell arrangements at the boundary with an

attendant difficulty in the logic of coding any automrtic routines for

setting up the mesh and the difference equations.

39
IV. NUMERICAL CAICUIATIOI?OF THE STRESSES IN AN INTTNITE HOLLOW CYLINDER

In order to test the difference equations of plane elasticity

derived in Section III, the following example problem was considered.

An infinite hollow cylinder of isotropic mterial was assumed to have

a temperature distribution varying inversely with the radius. The

inner surface was stress free; the outer was held fixed.

For the numerical calculation a section of the annulus between 0°

and 60° was chosen for the region of integration. Thus there was one

symmetry boundary coincident with the x-axis, and the other symmetry

boundary was at 60°. Because the basic mesh cells for the difference

equations are square, the displacement specified, 60° symmtry, and

stress free boundaries were zig-zag. In Fig. (hA) the segment of the

annulus and in Fig. (4B) its approximation by zig-zag and straight

boundary sections are shown. The boundary for the numerical solution

is identified by nunibers1 - 9 and letters A - U. The letter “O”

identifies interior points, the letter ‘Y’ exterior points. The numbers

and letters on the boundary tell the autonwtic set-up code which type

40
of boundary condition (stress,

and which of the cell energies

equations for the point. Fig.

displacement, or symmetry)

are going to contribute to

(4) is a reproduction from

a plot with a Benson-Lehner @20.

The analytic solution to the hollow cylinder problem

coordinates,

u

r—=.

a f-) l+V w - 1)

[

1+(1+)(:)*

l-v aTira 1+(1-2V)(:)2 -

—=- %

E

In the above formulae a is the inner radius, b the outer,

perature at a. The

usual significance.

are

~ = Ur Cos e,

=U ‘2 r ‘in “

is in effect

the difference

film taken of

is in polar

(Iv 1)

T1 the temisotropic elastic constants, E, v, and u, have the

The cartesian components correspondingto (IV 1)

41
’12

— = sin (3cos

E (~ -%1

’22 . ar ‘e 2

T

~ sin20 + ~ Cos e.

The components ~ and U2 of (IV 2)

placements calculated from general

(Iv 2)

are to be compared with the dlsdifference equations (III 7) and the

various boundary difference equations (III 8), (III 9), etc. The

stresses of (IV 2) are to be compared with stresses obtained by tifferencing the numerical displacements to form approximate strains and by

use of the stress-strain relations (II 16) with isotropic constants.

The approximate strains, computed from displacements at the corners of

square mesh cells, are centered at the centroid of each cell. The

stresses obtained are therefore not to be compared with the boundary

stresses of the analytic solution but with analytic stresses slightly

renrwed from the boundary.

The numerical problem

parameters:

a = 5.0 cm,

b . 100. cm,

h= 0.25 cm,

Ti = 5000°C,

v = 0.3,

a = 8 x 10-6/OC

was solved for the following values of the

(Iv 3)

42
The method used for solving the linear system of equations was

successive over-relaxation with extrapolation.(5) No attempt was

mde to optimize the extrapolation parameter u in this particular calculation, althou~ the value of u.),equal to 1.6, was chosen on the

basis of an earlier set of problems. The computer program consists of

a mesher, which will set up the boun~ry and interior points and their

difference equations, a relaxation code for the displacements with a

check on residuals, and a stress calculator. All of these programs

were written by L. Stein in Ivy, a compiler-assembler at LASL for the

IEM7040 and IBM 7090. The number of calculated mesh points was 691.

Running time until the maximum residual was reduced to 10-8

was about

25 minutes or 750 cycles on the IBM 709.

The results of the calculation of displacement are compared with

the analytical radial displacement in Figures 5-7. The calculated

stress is compared with the analytic stress on each of the four boundaries and in the interior on a radius of

=gures 8-22. While the displacemmt is

boundary and in the interior, the stress

particularly on the boundary.

Considering the fact that the

to satisfy the same boundary conditions

45°. These are shown in

sufficiently accurate on the

is naturally less accurate,

zig-zag boun~ry

that are correct

was assumed

for the true

boundary, the loss of accuracy at the boundary was not unexpected.

Smoothness of the stress and displacement in the interior suggested

that values could be extrapolated to the zig-zag boundaries which could

43
then be used to obtain an improved solution particularly near the

boundary where the effect would be greatest. Adjustment of the displacement boundary conditions at the outer radius was performed with

the improved results in the stress shown on part of that boundary in

Fi~es 17-19. Adjustmnt of stress free boundary was tried at the

worst point on the bounda~ with the results shown in Figures U-13.

To be really effective an extrapolation scheme of this type shouldbe

mde part of the program for the computing machine. Also an investi@tion

should be made to determine whether such an iteration of bountiry velues

is an always convergent process.

Another procedure for handling difficult areas near the boundary

is made feasible by the existence of L. Stein’s meshing code. The subregion in which difficulty is encountered can be removed from the rest

of the problem and its nuniberof calculation points increased. Boundary

values for the artificial boun~ry created can be obtained from a first

run for the real region of integration. The automatic mesher makes such

a refined calculation straightforward, whereas set-up by hand for a

large number of points is tedious and time-consuming.

It might be thou@t that a differencing scheme for the boundary

points based on interpolation to the correct boundmy location would

be an answer to the difficulties encountered on the boundary. Such an

approach was taken in an early attempt on this problem, but the iteration process failed to converge. The total matrix of the &Lfference

equations was too co~lex to analyze, and that approach was abandoned

44
for the variational derivation. An attempt to obtain interpolation

formulae which would necessarily leave the m.trix of total equations

positive definite should be made. Still to be tried are the more complicated energy relations based on triangular regions described in the

latter part of Section III.

It is hoped that further investigation of the techniques at the

boundary will yield a practicable and simple solution to the problem

of accuracy at the buun&ary which is an inherent difficulty with

elliptic equations. In the meantime the present scheme based on square

cells should provide displacements of sufficient accuracy and, with the

adjustment of boundary values and boundary mesh refinement, stresses

of useful accuracy as well.

45
REmcEs

1. Landau, L. D. and Lifshitz, E. M., ‘?l?heoryof Elasticity”

(Addison-WesleyPublishing Co., Inc., Reading, Massachusetts),

pages 36-41 (1959).

2. Callen, H. B., “Thermodynamics” (John Wiley and Sons, Inc., New

York), pages 226-231 (lg60).

3* Sokolnikoff, I. S., “Mathematical Theory of Elasticity” (McGrawHill Book Company, Inc., New York), Second Edition, pages 86-87

(1956).

4. Courant, R. and Hilbert, D., “Methoden Der Mathematischen Physik”

(Julius Springer, Berlin), Second Volume, pages 141-142 (1937).

5. Forsythe, G. E. and Wasow, W. R., “Finite-differenceMethods for

Partial Differential Eqyations” (John Wiley and Sons, Inc., New

, York), pages 242-266 (lg60).

46
Figure 2

a- -- —— ---- ___

v

I

<

\

I

,\

\

I

‘.

.

1

‘N 1

\

I

Y

I

;

\

I

\

/’ \

\

\

(

/

=-—_—- -

Approximation of boundary by grid lines.

/

\

\

\

Approximation of boundary by chords.

Figure 3

2 1

/“y / ,/ I

, 0’

/ /1

I

/’ //

/

I

/

4 5

Pentagonal boundary cell.

47
o

w

48
. . ...> >. >>.. .>** a.000m9w -- w-w-w W**>*>*

. .

49
o

IA

CD

II

0
%JG

CDEb

a)

51
8

a)

‘lG

52
53
e=~w

“-

0

is

E)

a)

&

0

9+
55

0

(n

i!

-P

i?

0
d

56
%

-’=

5=IW

57
?2

8

al
+.

3

a)

II

is

59
%?

Uj

k!

al

1%

d

60
%3

CD

IIk.m!

-P

al

c-u

i48

a)

l-l

w“ b
62
u.

‘=

b=lw

63

l-l

cd
,%.

?2

2)

a)

s’

x

+

64
-“”= x

~=1-

II.0

9-I i$

L1

%-l 0

LI
?2

!3

a) o

q+

66
x

e=l”’

67