Piecewise hierarchical p-version axisymmetric shell element for non-linear heat conduction in laminated composites

Compw~rs & S~rucrvns Vol. 47. No. I. pp. i-18. Printed in Gnat Britain.

1993 0

004s7949193 56.00 + 0.00 1993 Pergamon Press Ltd

PIECEWISE HIERARCHICAL P-VERSION AXISYMMETRIC SHELL ELEMENT FOR NON-LINEAR HEAT CONDUCTION IN LAMINATED COMPOSITES A.

BOSE

and K. S. SURANA~

Department of Mechanical Engineering, The University of Kansas, 3013 Learned Hall, Lawrence, KS 66045, U.S.A. (Receiued

17 January 1992)

Abstract-This paper presents a three-node axisymmetric shell finite element formulation for non-linear axisymmetric steady-state heat conduction in laminated composites where the temperature approximation for the laminate is piecewise hierarchical and is derived based on the p-version. The temperature approximation for the element is developed first by establishing a hierarchical temperature approximation for each lamina of the laminate and then by imposing inter-lamina continuity conditions of temperature at the interface between the laminas. The approximation functions and the nodal variables for each lamina are derived directly from the Lagrange family of interpolation functions of order p{ and “p,,. This is accomplished by first establishing one-dimensional hierarchical approximation functions and the corresponding nodal variable operators in the e and “q directions for the three- and one-node equivalent configurations that correspond to pt + 1 and “p, + 1 equally spaced nodes in the 5 and kq directions and then taking their products. The nodal variables for the entire laminate are derived from the nodal variables of the laminas and the inter-lamina continuity conditions of temperature. The element formulation ensures Co continuity or smoothness of temperature across inter-element as well as the inter-lamina boundaries. Weak formulation of the Fourier heat conduction equation in the cylindrical coordinate system with temperature-dependent thermal conductivities, internal heat generation, film coeflicients and radiation effects is constructed and the individual lamina matrices and equivalent nodal heat vectors are derived using this weak formulation and the hierarchical temperature approximation for the laminas. Inter-lamina continuity conditions are used to construct the transformation matrices for the laminas. These matrices permit transformation of the lamina degrees of freedom to the laminate degrees of freedom. Using these transformation matrices individual lamina matrices and the equivalent heat vectors are transformed and then summed to obtain the laminate element matrix and equivalent heat vectors. There is no restriction on the number of laminas and their lay-up pattern. Each layer can be generally orthotropic. The material directions and the layer thickness may vary from point to point within each lamina. The geometry of the axisymmetric shell element is defined by the nodes located at the middle surface of the element and the lamina thicknesses. Due to the temperature dependence of thermal conductivities, internal heat generation, film coefficients and radiation parameters; the resulting discretixed equations of equilibrium are non-linear. These equations are solved using the modified Newton-Raphson method. The tangent matrix is derived using the non-linear equations of equilibrium. The modified tangent matrix is symmetric and provides excellent convergence during the equilibrium iteration process. Numerical examples are presented to demonstrate the accuracy, efficiency and overall superiority of the present formulation and the extremely good and fast convergence of the iteration process. Comparisons are made with the available results in the literature.

NOTATION

bilinear functional of a test function u and temperature T derivatives of kth lamina approximation functions incremental length along an element boundary represents a typical element equivalent nodal heat vector for the entire mesh equivalent nodal heat vector for the kth lamina

WV, T)

PC1 dT

if} {YI

t Author addressed.

WI

to

whom

all correspondence

should

be

FHI 4

Rice)

assembled conductivity matrix for the entire mesh conductivity matrix for k th lamina thickness of kth lamina quadratic functional for an element, e Jacobian matrix for the kth lamina matrix of thermal conductivities thermal conductivities in the global r, z directions a typical lamina in laminated composite linear functional of test function v total number of lamina direction cosines of an outward normal to the element boundary parabolic shape functions for geometry

2 PW,

A. BOSEand K. S.

‘rtll

PN,l P&l

T, To aT aT aT aT ---ar ’ a2 ‘ar”af

temperature approximation functions for k th lamina retained elements of kth lamina temperature approximation functions eliminated elements of kth lamina temperature approximation functions polynomial order in the axial (<) direction of the laminas polynomial order in radial (9) direction of the k th lamina internal heat generation per unit mass distributed heat flux per unit area heat flux due to convection per unit area heat flux due to radiation per unit area radial coordinate radial coordinate at node i mean radius temperature temperature at node I’ equilibrium temperature for no radiation equilibrium temperature for no convection specified temperature on boundary of an element derivatives of temperatures with respect to r, z, r’ and z’

kT 3 ” ar ’ dZ”T. I ,.... ar2 apq aem i kr, v Vi

hierarchical kth lamina nodal variables at a node i total thickness of laminate thickness of kth lamina at a node j test function unit normal to the middle surface of the shell at a node i axial coordinate radiation heat transfer coefficient convective heat transfer coefficient emissivity Stefan constant natural coordinates in the axial and radial directions of the laminate natural coordinates in the axial and radial directions of kth lamina eliminated degrees of freedom for k th lamina retained degrees of freedom for kth lamina domain of an element element boundaries on which distributed heat flux, convection or radiation act boundaries on which temperature T is specified vector of generalized nodal variables for an element e coordinate in the thickness direction density INTRODIJCXION

The development of composite material technology has integrated the material design and the component design processes. The flexibility to ‘tailor-make’ a material for a specific application has provided

SURANA

greater design flexibility. Although the advantages are obvious, designing with composites is not without problems. Their inherently complex material structure presents enormous difficulties in simulating the component behavior which is essential during the design process. Numerical simulation of temperature distribution in laminated composites is of much concern, interest and importance. The first application of the finite element method to steady-state heat conduction was presented by Zienkiewicz and Cheung [l]. Since then, Wilson and Nickel1 [2], Zienkiewicz and Parekh [3], Emery and Carson [4], Donea [5], Bruch and Zyvoloski [6], Podovan [7, 81, Rathjen and Jiji [9], Comini et al. [lo], Tham and Cheung [l 11, Banasek [12], Azevedo and Wrobel[13], Fan and Tsai [14], Goldak et al. [15], Kleiber and Sluzalec [16], etc., have presented various aspects of the finite element formulations, solutions procedures, accuracy, material behavior and stability studies related to the heat conduction process. The polynomial approximation concept (papproximation, p-version, p-extension) in the finite element method was first proposed in 1970. A selected number of references are cited here [17-281. A more comprehensive list of references of p-approximation technique can be found in a recent publication edited by Babuska et al. [29]. Although the p-approximation approach has been applied to the two-dimensional field problems there is virtually no published or reported work for axisymmetric shell heat conduction in composite laminates using this approach. In a recent paper, Surana and Kalim [30] (also [31]) presented an axisymmetric shell element formulation for heat conduction where nodal temperatures and nodal temperature gradients (in the element thickness direction) were retained as nodal variables at the element nodes. Many commercial finite element systems such as ANSYS [32], NASTRAN [33] and NISA [34] report two-node axisymmetric shell elements with constant temperature through the shell thickness. In another recent paper Surana and Phillips [35] (also [36]) presented a formulation for the three-dimensional curved shell element where temperatures and temperature gradients in the shell thickness direction were considered as nodal degrees of freedom permitting linear temperature distribution through the shell thickness. In laminated composite shells the assumption of linear temperature distribution through the shell thickness is hardly adequate. The approach where each lamina of the laminated composite is modelled using axisymmetric solid elements is highly impractical and may even lead to numerical difficulties because a composite laminate may contain as many as 20 or more individual lamina each with a thickness of under 1 mm. These modeling practices and analysis methods either lead to extreme inefficiency or begin to breakdown in severe cases where the laminate is either very thin or very thick with many laminas. Accurate and efficient heat conduction studies in such

p-Version axisymmetric shell element applications require simpler geometric models and higher order and variable temperature approximation in both radial and axial directions of the element. In recent papers Surana and Orth [37-391 presented an axisymmetric shell element formulation for isotropic and laminated composite applications where the temperature approximation in the shell thickness (transverse) direction could be of arbitrary polynomial order p. Surana and Abusaleh [40] presented a three-dimensional curved shell element formulation with higher order hierarchical temperature approximation in the thickness direction. These formulations have proven to be extremely useful for thick shell and laminated composite applications while maintaining their effectiveness for very thin shells. These formulations with hierarchical temperature approximation in the thickness direction, although very effective in many applications fail to perform optimally when concentrations and singularities, etc., are present. In such cases a variablep-level control is needed for both longitudinal (r) and transverse (q) directions of the axisymmetric shell element. To remedy these shortcomings Surana and Orth [41,42] presented a completely hierarchical axisymmetric shell element for laminated composites where the temperature approximation for the laminate in the 5 and q directions is described by polynomials of arbitrary orders pt and pII which can be chosen independently to achieve maximum rate of convergence. The formulation is completely hierarchical, i.e., the approximation functions and the nodal variables corresponding to lower p-levels (p[, p,,) are a subset of those corresponding to higher p-levels (pt + 1, ps + 1). The approximation functions and the nodal variables for the element are derived directly by using Lagrange interpolating polynomials of order pc and p,, in the 5 and r~directions. A major drawback .of the formulation presented by Surana and Orth [42] is that the temperature approximation in the transverse direction of the laminate is described by a single higher order polynomial whereas the physics of heat conduction in laminates requires a piecewise behavior, i.e., for each lamina the temperature distribution in the transverse direction may be a higher order curve with inter-lamina continuity of the temperature but discontinuity of the slope of the temperature in the transverse direction at the interlamina boundaries. The formulation of Surana and Orth [42] satisfies continuity of temperature at the inter-lamina boundaries but also enforces continuity of first and higher order derivatives of temperature in the transverse direction at the inter-lamina boundaries which is not in conformity with the physics of heat conduction in laminates. The motivation of the work presented here is to remedy the shortcomings of Surana and Orth [42] and to develop a formulation which is totally in conformity with the physics of heat conduction in laminated composite axisymmetric shells and to provide temperature dependence of thermal conductivities, film

3

coefficients internal heat generation and radiation effects. This paper presents a three-node axisymmetric shell finite element formulation for axisymmetric steady-state heat conduction in laminated composites where the temperature approximation for the laminate is piecewise hierarchical and is derived based on the p-version. The temperature approximation for the element is developed first by establishing a hierarchical temperature approximation for each lamina of the laminate and then by imposing inter-lamina continuity conditions of temperature at the interface between the laminas. The approximation functions and the nodal variables for each lamina are derived directly from the Lagrange family of interpolation functions of order pe and “p,,. This is accomplished by first establishing one-dimensional hierarchical approximation functions and the corresponding nodal variable operators in the l and kq directions for the three- and one-node equivalent configurations that correspond to pc + 1 and “pq+ 1 equally spaced nodes in the 5 and kq directions and then taking their products. The nodal variables for the entire laminate are derived from the nodal variables of the laminas and the inter-lamina continuity conditions of temperature. The element formulation ensures Co continuity or smoothness of temperature across interelement as well as inter-lamina boundaries. Weak formulation of the Fourier heat conduction equation in the cylindrical coordinate system is constructed and the individual lamina matrices and equivalent nodal heat vectors are derived using this weak formulation and the hierarchical temperature approximation for the laminas. Inter-lamina continuity conditions are used to construct the transformation matrices for the laminas. These matrices permit transformation of the lamina degrees of freedom to the laminate degrees of freedom. Using these transformation matrices, individual lamina matrices and the equivalent heat vectors are transformed and then summed to obtain the laminate element matrix and equivalent heat vectors. There is no restriction on the number of laminas and their lay-up pattern. Each layer can be generally orthotropic. The material directions and the layer thickness may vary from point to point within each lamina. The geometry of the axisymmetric shell element is defined by the nodes located at the middle surface of the element and the lamina thicknesses. Due to the temperature dependence of thermal conductivities, internal heat generation, film coefficients and radiation parameters; the resulting discretized equations of equilibrium are non-linear. These equations are solved using the modified Newton-Raphson method. The tangent matrix is derived using the non-linear equations of equilibrium. The modified tangent matrix is symmetric and provides excellent convergence during the iteration process. Numerical examples are presented to demonstrate the accuracy, efficiency, convergence characteristics and the overall superiority of the present

A. Boss and K. S. SURANA

4

formulation. The results are compared with analytical solutions, h-models and those reported in the literature. WEAK FORMULATION

The Fourier heat conduction equation for axisymmetric heat conduction with globally orthotropic material properties can be written in the cylindrical coordinate system as -;

($$+g>

other parameters are known numerically or by explicit expressions. The weak formulation of the Fourier heat conduction equation (1) is constructed by taking the scalar product of the differential equation (1) with a test function (say V) over a typical element P. Through integration by parts, one order of differentiation can be transferred from the dependent variable T to the test function v. After substituting the boundary conditions (2) and rearranging the terms we obtain the following

(1)

+&IQ(T)=0

NO)

(8)

= Qv),

where

with g,n,+g,n,+q+q,(T)+q,(T)=O

on r,

(2)

and onr,,

T=T,

where VT@ + a,)2nr dr

+

(9)

f l-1

g,=rk$+rk,g

and (3)

,=rk,,;+rk,,$ n, and n, are the r and z component

of the outward normal A to the boundary r, . Here we assume that the matrix k[T’j of thermal conductivities given by [k,,,] =

I, 12 1

:_ [

>

(4)

is symmetric, i.e., k,, = k,s. These equations describe steady-state axisymmetric heat conduction in axisymmetric bodies with thermophysical properties dependent on temperature and Q is the rate of internal heat generation. The variables q, qE and q1 represent imposed heat flux, the heat flux due to convection and radiation, respectively. The following relationships are well established in finite element formulation literature

qc=B(T)(T- Tm) 4, =

MT4 - T:,) = a,V)(T - To,),

(9 (6)

where /I is the convective heat transfer coefficient, c is the emissivity, o is the Stefan constant, a, is a parameter related to the effects of radiation by the expression a,(T)=m(T2+

T;,)(T+T,).

(7)

T, and T,, are equilibrium temperatures for which no convection and radiation occurs. In general, we assume that [k], Q, /I and a, are temperature-dependent functions with a, expressed explicitly by (7) and the

veQ2nr dr dz +

I(v) =

vT,fi2nr

dT

I rl +

s

VT,a,Zxr

r1

dT -

vq2nr dr.

(10)

s rI

Equation (5) constitutes the weak formulation of (1) with boundary conditions (2). The quadratic functional associated with this weak formulation can be obtained using the following if [k], Q and fl are not dependent on temperature and if the radiation effects are absent I’(T) = fB(T, T) - I(T).

(11)

ELEMENT GEOMETRY AND MAPPINGS

Consider the shell element to be composed of an arbitrary number of generally orthotropic laminas with arbitrary lay-up pattern bonded together perfectly at the inter-lamina boundaries. Further assume the material directions and each layer thickness may vary from point to point in the longitudinal direction of the element. Figure l(a) shows a map of the three-node axisymmetric shell element in the r-q natural coordinate system under consideration here. The element nodes are located at the middle surface (q = 0) and the bottom and the top surfaces of the element (q = - 1, q = + 1) are defined by these coordinates and the node point vectors. Thus the r and z coordinates of a point P(& r)) can be expressed as

p-Version axisyrnmetric shell element where R,(c); i = 1,2,3 are the parabolic shape functions. Vi(r); i = 1,2,3 are the nodal vectors at the element nodes defining the bottom and top surfaces of the laminate. In order to facilitate the development of the element properties two different natural coordinate systems are utilized. The first coordinate system is c-q in which the entire element is mapped into a two-unit square with the origin of the coordinate system located at the center of the square. Figure l(b) shows the three-node shell element and the kth lamina details. Figure l(c) shows the map of the element and the kth lamina in the c-r) coordinate system. The C coordinate is in the longitudinal direction of the element whereas the tl coordinate is in the transverse or the lay-up direction of the laminas. Since the laminas are considered to be of variable thickness, the map of the kth lamina in the t-rl coordinate space

5

will be of variable thickness also. The second type of natural coordinate system is the lamina natural coordinate system c-kq in which the kth lamina is mapped into a two unit square [Fig. l(d)]. The following simple but useful relations can be easily established from purely geometric considerations. For the kth lamina the lamina thickness k/z=%(r) at a point r within the lamina can be determined using

kh = kh(() =

5 &+i((5)khi,

(13)

i=l

where khi; i = 1,2,3 are the thicknesses of the kth lamina at the locations corresponding to { = - 1, 0 and + 1 of the laminate nodes.

--------------

r

I

kth lamina node

(4

@)

k-

Vl

kv2

kv3

k2

Fig. 1. Three-node axisymmetric laminated shell element. (a) Element map in the T-q space and the laminate nodal vectors. (b) kth Lamina details. (c) kth Lamina map in the l-q coordinate space. (d) kth Lamina map in the Tq coordinate space.

A. Bow

6

and K. S.

The relationship between the q and kq coordinate system can be described as

q-1+$ -‘h(l-‘q)+2i’h.

(

j=l

>

(14)

The inverse of the mapping defined by (14) can be written as t ‘9=-l+F(‘+‘1)-$,-$Ij’~

2 k

(15)

where r= i+r;

m being the total number of laminas.

j-l (16)

Using eqns (13) and (14) inter-lamina values of q coordinate can be easily established for different values of t and q. GENERAL

CONSIDERATIONS

AND PROCEDURE

For each lamina a hierarchical temperature approximation is constructed in the lamina coordinate system (i.e., e-kt) coordinate system for the kth lamina). The temperature approximation for the kth lamina can be of arbitrary polynomial orders pI; and “p, in the longitudinal (5) and transverse (kq) directions of the lamina. It is worth noting that the order of temperature approximation in the l direction is the same for all laminas (necessary for inter-lamina continuity of temperature) but in the transverse direction of the laminas the order of approximation can be different (jp,;j = 1,2, ...,m) for each lamina. The approximation functions and the nodal variables for the kth lamina are derived directly from the Lagrange family of interpolation functions of order pt and "p, . The nodal variables for the entire laminate are determined using the nodal variables of the laminas and the inter-lamina continuity conditions of the temperature. The following is a step-by-step outline of the general procedure for establishing lamina temperature approximation, consideration of degrees of freedom, inter-lamina continuity conditions of temperature and computations of lamina and laminate matrices and vectors. First using the geometric description of the entire laminate given by (12) and eqns (14)-( 15) the lamina nodes are located for kth lamina (k = 1, . . . , m). Figure l(a) shows lamina nodes kl, k2 and k3 for the k th lamina. Next we establish nodal vectors kv, , kV,, kVJ at these nodes of the k th lamina. The k th lamina is mapped into its own natural coordinate spaces { -kq (- 1 < r < 1, - 1 f ‘r~ < 1; Fig. l(d)) and thus an equation similar to (12) can be used to describe the geometry of the kth lamina. Now we construct hierarchical temperature approximation kT(& ‘q) for the kth lamina in the r-kq

SURANA

coordinate system. The heat conduction matrix and the equivalent nodal heat vectors for the kth lamina can now be derived using the hierarchical temperature approximation and the weak form of the Fourier heat conduction equation given by (8). From the individual lamina temperature approximation; we now construct inter-lamina continuity conditions. That is; the temperatures k-‘r(& k- ‘q) and kT(& kq) are equated at (pr + 1) equally spaced points at the interface between the (k - 1)st and the kth laminas. These conditions are arranged in matrix form which provide the essential and needed ingredients for eliminating some degrees of freedom from the composite laminate. We also note that equating the temperatures between (k - 1)st and kth lamina at (pt + 1) equally spaced discrete points for the ptth degree polynomial temperature approximation ensure continuity of T everywhere along the entire interface between (Ic - 1)st and kth laminas. For the first lamina all degrees of freedom are retained in the laminate degrees of freedom. For the second and subsequent lamina (pe + 1) degrees of freedom are eliminated using inter-lamina continuity relations and the remaining degrees of freedom are added to the laminate degrees of freedom. To make this process systematic and recursive, inter-lamina continuity equations are used to construct transformation matrices which transform lamina degrees of freedom to the laminate degrees of freedom. The individual lamina heat conduction matrices and equivalent heat vectors are transformed using these transformation matrices so that each lamina equilibrium equation only reflect the degrees of freedom present for the laminate. The transformed lamina matrices and vectors are then summed to give the matrices and the vectors for the entire laminate. The numerical values of the lamina matrices and equivalent heat vectors are calculated using Gaussian quadrature with pc + 1 and “p,,+ 1 rule for c and kr) directions.

LAMINA

TEMPERATURE

APPROXIMATION

Consider a typical, say k th lamina of the laminate. The geometry of this lamina can be easily described in the lamina coordinate system t-k~

The approximation functions and the nodal variables for the k th lamina are derived using one-dimensional Lagrange interpolation functions, the details of which can be found in [42]. A brief outline is presented in the following. The pc + 1 equally spaced nodes in t direction between - 1 and + 1 with their corresponding Lagrange interpolation functions and

p-Version axisymmetric shell element nodal temperature values can be reduced to an equivalent three-node configuration for which we can write

[jq]=

[

1

r2-1 r3-r

[N21_

c -

!$

[

r4-1

rP( -a

--jj--‘F’-$--‘“.‘------

(Pt)!

1 (24)

where a =

1 r

if i is even, if i is odd.

The “p,,+ 1 equally spaced nodes in the kq direction between - 1 and + 1 with their corresponding Lagrange interpolation functions and nodal temperature values can be reduced to an equivalent one-node configuration for which we can write the following for a typical node “j (9 = 1,2,3) k.-f(b$

=

[Af,,,43, . . . , A&J = ($)‘pv

01 kti h

. . . . (“P,>! _Z

;

i=l,2,3.

(26)

The temperature approximation for the kth lamina given by (20) is hierarchical. We note that the elements of [N,] are functions of r only.

2 tk’f)’‘tj ’ a’tkT) ,_o&)(7s-)k@‘~ (19)

where k$ is the nodal vector norm at node “j and ky is the same coordinate direction as lrq with

ktj

LAMINA PROPERTIES

Substituting the temperature approximation for the kth lamina given by (20) into the weak form given by (8) and choosing v = Ni gives the lamina properties which can be arranged in the following form

kt, 2

-PLY<‘. Equations (18) and (19) provide the required onedimensional approximation functions and the nodal variable operators whose products give the desired approximation functions (details are omitted for the sake of brevity) and the corresponding nodal variables for the kth lamina with three nodes ‘1, k2 and k3. Thus we can write the following for the temperature approximation for the k th lamina

[kmkS) =

(27)

It/-~9

where [kHij] of [kH] is given by kHij =

[kCi]T[k][kC,] det[J]2nr dr dkq

[kNIT[kN](j? + a,)2nr dr

+

kT = kT(& ktl) = fkN(& ‘rl)l{“S}

(28)

I r!

(20)

[kN]T[kN]Q~2nr det[J] dc dkq

where

-

df

PNlTq2nr

s r1 {k6JT=FT,,(T);(?$)z [kN]T(BTz + T,ct,)k

+

dT

(29)

s rt ...,(q);kT,]

(21)

and akT ar

{~a,)‘=[~~],[!],

akT . . ..pw]]

we1 = [AWN ml9

Nil

= [kC](kS).

(30)

ii aZ

(22) Numerical values of (23)

using Gaussian rule.

(“f} are calculated with pe + 1 and “p,,+ 1

[kH] and

quadrature

A. bosx and K. S. thtANA

8

INTER-LAhUNA CONTINUITY CONDITIONS, LAMINA TRANSFORMATION MATRICES AND LAMINATE PROPERTIES

For each lamina the hierarchical temperature approximation is defined by eqn (20). The p-level in the { direction is the same for all laminas whereas the p-level in the ktf direction can be different for each lamina. Consider the inter-lamina continuity of temperature between lamina 1 and 2, i.e.

*T(~i,12tli)=1T(ri,‘2tli);

i=l,2

)...,

p,+l, (31)

where ti and 12rliare the 5 and rl coordinates of the pc + 1 equally spaced points at the interface between lamina 1 and 2. Equation (3 1) can be expanded using (20) and arranged in the following form

The matrix tT] transforms the degrees of freedom for lamina 2 to the degrees of freedom for the laminate consisting of laminas 1 and 2. Now we consider inter-lamina continuity condition of temperature between lamina 2 and 3, i.e. ‘T(&, 23Q)= V(&, %li);

i = 1,2, . . . ,pc + 1, (37)

where & and 23qiare the 5 and q coordinates of the pc + 1 equally spaced points between lamina 2 and 3. Equation (37) can be expanded using (20) and we can write the following [[%I, [3~Jl{;:;;;}

= ti2pel, i2prll{ ;:“6:;) .

(38)

From (37) we can write [t2~el,i2~,ll{ ;;;)

= t[‘p,lv 12pJl{ I:;;)

.

(32) I’S,) = rrfl, rm1 I’mI

(39)

If we define {IX} = {il:i}

and

{‘x} = { iis:),)

where (33) [‘El = [3Pel-‘[2Pel= vl

We can then write the following using (32) [28] = [3pe]-1[2p,] and

[‘k] = -[3pp,]-1[3p,]. w

Substituting for {‘S,) from (34) into (39) we obtain the following =[[‘a

rm1 jg 1

3 1

(34)

where [‘El = [2Prl-‘[‘Pel = VI,

(41)

[‘WI = [2Pel-‘[‘P,

or

and

Mn

{‘h) = N2Rl, I221 = -[2Pel-‘[2Prl.

(35)

Here we note that {Ix} are the total laminate degrees of freedom for lamina 1 and {‘x} are the total degrees of freedom for the laminate consisting of lamina 1 and 2. Thus the inter-lamina continuity relations (32) between lamina 1 and 2 allow us to eliminate {26,} degrees of freedom from lamina 2. Obviously all degrees of freedom of lamina 1 are retained in the laminate degrees of freedom. Using (34) we can derive the following transformation matrix for lamina 2

(36)

(42)

If we define (43) Then using (41) we derive the following transfonnation matrix for lamina 3

The matrix [‘T] transforms degrees of freedom for lamina 3 to the degrees of freedom for the laminate consisting of lamina l-3. Similar transformation matrices can be established for lamina 4, 5, etc. In

p-Version axisymmetric shell element general between lamina k - 1 and k we can write the following

=

PTl{‘xl,

using inter-lamina continuity relations. Inter-lamina continuity equations are conveniently arranged in the form given by (48) which permits transformation of lamina degrees of freedom to the laminate degrees of freedom easily. For each lamina the conductivities can be defined in a material coordinate system which may be different for each lamina and may also vary in the < direction for each lamina. The details of transforming [k] are omitted for the sake of brevity.

(46)

SOLUTION OF DISCRETIZED EQUATIONS OF EQUILIBRIUM

where [k-‘R] = [[k-‘Jil[k-*R], [[“-‘JQ[“-‘J]

+p-‘WJJ]. (47)

Thus we have (G} = [‘T]{‘x}; j = 2, 3, . . . , m.

9

When the thermal conductivities of the material, convection coefficients and the internal heat generation are not dependent on temperature and when the radiation effects are absent the assembled equations of equilibrium for the entire mesh

(48)

[‘H](G) = (‘s);

j = 1,2,3, . . . , m.

(49)

Substituting from (48) and (49) and then pre-multiplying the entire equation by [jT]r we obtain the following [?f]r~~][jT](‘x}

= [jT]r{‘f};

j = 2,. . . , m, (50)

where {“x} represents the total degrees of freedom for the laminate. The matrices [jT]r[jH][jT] and the vector [‘Tlr(‘f} provide the non-zero contribution of each lamina to the laminate conductivity matrix and equivalent heat vector. Thus for the entire laminate we can symbolically write the following (only the non-zero contribution of each lamina shown)

(a} = {“x}.

represent a system of linear simultaneous algebraic equations which can be easily solved using elimination methods to obtain the unknowns (6). In the present work these equations are assembled and solved using wavefront technique. When the thermal properties and the internal heat generation are temperature dependent and when the temperature dependent convective and radiation boundary conditions are present, the coefficients of the assembled matrix [H] and the right-hand side vector {f } in (5) are functions of temperature which is unknown. Thus, eqn (54) then represents a system of simultaneous non-linear algebraic equations in the unknowns (6). Hence their solution requires the use of iterative procedures. In the present work the modified Newton-Raphson procedure is employed. Since [H] and {f} in (54) are functions of T (which in turn is a function of (6)) which is unknown, an initial guess given by (6) will not satisfy (54) and thus we can write

(52)

At this point we make several comments. First we note that [%,I = [%%I= . . . = rhl.

(54)

[HI@) = If1

Also from (27) for each lamina we have

(53)

Secondly if each lamina has constant thickness (but the thickness may vary from lamina to lamina) that is if kh is not a function of < then it can be shown that inverse of [jp,] is not necessary (details omitted for the sake of brevity). This can be confirmed by closely examining the approximation function matrices defined by eqns (24)-(26). We also note that all degrees of freedom for lamina one are retained in the laminate degrees of freedom whereas for lamina 2,3, . . . , m: only {/s,};j = 2,3, . . . , m degrees of freedom contribute to the laminate degrees of freedom and (‘6,); j=2,3,..., m degrees of freedom are eliminated

{H&l)

- U1=

P%

(55)

where {Y } are called residuals (or unbalanced generalized nodal heat vector). The procedure is to find a correction to (6) represented by {Ad} such that {Y({61+ {Aa})) = (0). Expanding {Y({6} + {AS})} in a Taylor series about (6) and retaining only the linear terms gives P’(@1+

{As))} = {Y({6))1

wu)

+

[

m

161{A61+...=0.

(56)

1

Thus (A6) =

(57)

A. ROSEand K. S. SURANA

10

The improved value of (6) is given by ({S} + {Aa)). Similarly In eqn (57) Y((6j) are the residuals calculated at a{f>la@) (~51 using (55). [~{Y}/~{~}] is called the tangent matrix which is to be determined using (55). Noting that [H]=[H(T)], {f} ={f(T)}and T= T((d}), then

using

(58)

+

s r1

(B + ~,Wl%V dr

following

results

for

s

+v

aT

TmWl E.mdr

r1

Recall that

IWIklIN dQ +

the

TaQzdR [N] Ea{s)

+

f#=[Hlf$+}-#.

IHI =

(60)

Tau,aT TJNI aTaIs)dr. s r1

(W

Substituting from (63) and (64) into (58) and noting that aT/a{tS I= [N], the tangent matrix becomes

(59) z

and

= [HI +

[Nlr[N] !$ + 2

1.,

{f>= JWQ s

da+

+

s

(

)

([Nl{6}) dr

~,BTm[NITd~

s

u,T,,[N]~ dJ- -

r1

q dT.

-

(60)

s r1 -

From (59) by defining [a = [R]qk][R],

Tm;T[N]r[N] dT s rl

a[Hl -= am

-a ([Rb dQ s ne w

+

s

-?r,a{s)

-

T$[N]=[N] dT.

If we substitute [N]{6} = T in (65) then

(WlrbWB + a,)) dr.

(61) (66)

But

aim

-a@)

(65)

a[Bj - az-aTa{s}‘aT

arm

where

rN1 t4]=IH]+~~,(~+~)T[Nl~~NldT

and

or

Tm;T[N]'[N] dT ?-

([Nl’[Nl(P

ais1

+ a,)) = [N]‘[Nl

-

[Nl.

To,% [NITIN] dl-.

(62) and Thus

@“A= ,g s

+ Jr, WITINl($ + $Wl

dr.

(63)

(T)dR.

(68)

Note that [fiJ is the symmetric part of the tangent matrix whereas [fi,J is the non-symmetric part of the tangent matrix. In the present work [fi,,J is neglected to avoid the additional complications arising in the

p-Version axisymmetric shell element solution of algebraic equations due to a non-symmetric tangent matrix. Thus the modified tangent matrix can be written as

11

FPS or metric system can be selected. For this reason units are not given for the examples presented in this paper. Computations were performed on an Apollo DN3500 computer using double precision arithmetic. Example 1

and the incremental

correction becomes

{Ad}= -

g [

The main computational tive procedure are: Step

Step Step Step Step Step Step

1

{Y}. nt

(70)

steps involved in the itera-

1. Assume that [k], Q, c(, and p are not functions of temperature. Calculate [H] matrix and {fj vector using (59) and (60) assuming that (6) = (0). Now solve the assembled equations [H](6), = {f}. The solution {a}, constitutes a good starting vector for the Newton-Raphson iterative procedure. 2. Calculate the residuals {(Y({6,})} using (55). 3. Calculate [aJ] = [a {Y}/a {S}],,, using (67) and {S},. 4. Solve for {AS}, using (70). 5. Form improved solution vector using @j>,+, = G%+ IAh),. 6. Calculate residuals {Y((6}, + ,)}. 7. Check for convergence, i.e.

Step 8.

where A1 and A2 are preset tolerances and N* is the total number of degrees of freedom. If not converged then set (6)” = (6)” + 1and go to Step 3. If converged then stop the iteration process.

The numerical examples presented in the subsequent sections for non-linear steady state heat conduction are solved using this procedure. NUMERICAL. EXAMPLES

Here we present three numerical examples to establish the accuracy of the present formulation and to demonstrate its overall superiority for heat conduction in laminated composites. The results of the present formulation are compared with those available in the literature and that of Surana and Orth [42]. During the implementation or during the processing of data for the examples, conversion of units is not done. Hence any set of consistent units either in

In this example, we consider an infinitely long laminated circular cylindrical tube of mean radius 10.0 and thickness 10.0 (Fig. 2a). The cross-section of the tube consists of three different layers of materials with temperature-dependent conductivities. The inside surface of the tube (r = 5) is subjected to a uniform temperature field of 2000 while the outer surface (r = 15) is maintained at zero. The thermal conductivities for the inside, middle and outside laminas are given by the following polynomial expressions ‘k, = ‘k, = 7 - 5 x lo-3T + 3 x 10-62-2 -5 x lo-‘?) ‘k,, = ‘kzz = 2 + 5 x 10_3T - lo-Y2 3kr,=3kZI= 10+4x

10-2.

This problem has been analyzed in detail by Azevedo and Wrobel[13] using the boundary element method. Because of point symmetry along the length of the tube, the tube can be modelled with only one threenode p-version axisymmetric shell element whose nodes are located at the middle surface of the tube; i.e., at r = 10 (Fig. 2b). Since the temperature varies only in the radial (‘r,t) direction, the p-level in the axial (l) direction, i.e., pc is chosen to be the minimum value 2. p-levels in the kn direction of the laminas are progressively increased from two to six. The results were also computed using the formulation of Surana and Orth [42]. Also the tube is modelled using only one three-node p-version shell element. While pc was chosen to be 2, the p-level in the q direction (p,,) for the entire laminate was increased from 1 to 11. It is rather inconvenient to specify the temperatures at the inside and outside surfaces for both formulation in their present forms. So an equivalent set of heat flux and convection parameters (q = 1774.7822, B = - 118.3188 and T, = 5.0) was chosen as boundary conditions. Since the purpose of this example is primarily to demonstrate the numerical accuracy of the present formulation, specifying heat flux and convection coefficients at the boundaries poses no theoretical discrepancy in simulating the original problem studied by Azevedo and Wrobe1 [13]. The results from both formulations are compared with the analytical solution presented by Azevedo and Wrobel [131. Figure 3 shows plots of radial temperature distribution in the tube for both formulations. We note that the present formulation gives results very close

12

A. BOSEand K. S. SURANA

Section A-A r,=lO,

t=lO,

h,=5,

hp3,

h,=2

rm

(b) Fig. 2. An infinitely long composite tube and its finite element model (Example 1). (a) Cylindrical tube details. (b) Finite element model.

to the converged values even for “p,,= 2. For “p, = 3-6 the results are virtually indistinguishable. The formulation shows inter-lamina continuity of temperature but discontinuity in temperature slopes (due to different thermal conductivities of the laminas) which is in agreement with the physics of heat conduction. Since the temperature-dependent conductivities chosen for this example do not have sharp variation at inter-lamina boundaries, the discontinuities in the temperature slopes are relatively mild. And so, the solution obtained from the formulation of Surana and Orth [42] matches closely with the analytical solution. Even then, it is worth noting that the formulation based on single higher order polynomial theory requires much higher p,, levels than the present formulation for optimum convergence. For both formulations, it should be noted that with an increase in p-levels, it takes less iterations for the modified Newton-Raphson procedure to obtain a converged solution. This is due to the fact that at higher p-levels, the higher order derivatives of temperatures are calculated more accurately, thereby reducing the relative magnitudes of the unbalanced generalized nodal heat vectors or residuals. This accelerates the convergence rate at higher p-levels.

Formulation of Surono k Orth I421

Analytical Solution ’ Azevedo k Wrobel I131 1200.0

-

800.0

-

400.0

-

a.01 ’ ’ 5

7

’

’ 8

’

’

’

’

10 12 Radial Distance

’

I-Y-= 13

15

Fig. 3. Radial temperature distribution (Example 1).

p-version axisymmetric shell element

Example 2 For this example, we choose the same infinitely long laminated composite tube geometry as in Example 1 but with different temperature-dependent thermophysical properties. Also, the inside surface of the tube is subjected to a uniform heat flux, q of 200.0 and the outer surface has convection coefficients as @ = 5.0 and T, = 30.0. The thermal conductivities for each lamina are chosen to be such that there is sharp variation of these thermo-physical properties in-between successive laminas across the thickness of the tube. The thermal conductivities for the inside, middle and outside laminas are expressed as fourthorder polynomials in temperature ‘k,, = ‘k, = 24.0 - 2 x 10-3T + 23 x 10-6TZ -2 x lo-“T3 + 25 x 10-14T4 *k,, = *k,, = 1.0 + 0.9 x IO-*T - 0.3 x 10-6T2 ‘k,, = ‘k, = 15.0 + 1 x lo-*T. Also the laminas are isotropic as in the first example. The element models and p-level controls are maintained exactly the same as for Example 1. Figures 4 and 5 show linear and non-linear radial temperature distributions for the present formulation and also for the formulation of Surana and Orth [42]. For both linear and non-linear cases, the present formulation based on piecewise hierarchical ap-

. P-7

360.0

preach converges at relatively low p-levels. For “p, = 3-5, the results are virtually indistinguishable. Moreover, the formulation shows inter-lamina continuity of temperature as well as distinct discontinuity in the temperature slopes as required by the physics of heat conduction. On the contrary, the formulation of Surana and Orth [42] based on single higher order approximation theory completely fails to simulate the linear and non-linear temperature profiles even at such high p,, levels as 12. Also worth noting is the fact that the formulation of Surana and Orth [42] is unable to simulate the ‘lack of C’ continuity’ at the inter-lamina boundaries. Example 3 Here we consider a cylindrical tube of finite length with a three-lamina generally orthotropic composite construction as shown in Fig. 6(a). Each lamina has generally orthotropic temperature-dependent thermal conductivities. The inside surface of the tube is subjected to convection with p = 5.0 and T, = 30; and the outer surface has a heat flux of q = 50.0 acting over a small axial length of the tube as shown in Fig. 6(a). The polynomial expressions for the lamina conductivities are chosen to be Lamina 1: k ,, = 20.0 - 2 x 10-3T + 23 x lO+‘T* (inside) -2 x IO-“‘T’ + 25 x 10-14T4 kzz = 1.0 + 2 x 10-3T + 23 x lO+T* -2 x 10-‘0T3 + 25 x 10-14T4

r

292.0

13

250.0 Present formulation

Present formulation ‘....p,=2 ‘p,=j p,=2 ‘p,=3-:

Formulation of Surano & Orth I421

x---p,= 2 p,= 1 --p,= * - -PI= *---PI=

204.0

l

2 p,= 4 2 p,= 6 2 p,=12

158.0

224.0 e 3 K !

112.0

156.0

68.0

20.0

5

Fig. 4. Radial

20.0

a

10 Radial Distance

13

15

temperature distribution (linear case; Example 2).

k 5

8

10 Radial Distance

13

15

Fig. 5. Radial temperature distribution (non-linear case; Example 2).

A. k3SE and K. S. SIJRANA

14

r,=lO,t=l8, h,=6, h,=6. h,=6

(a> r _““‘_“_“_f’_‘T’_.

=

I

1

I I I

I I

b-t--

2

:3

4 ,5

6 17

I I -1___L__J

---__-_-__-__

(b) Fig. 6. Finite length composite tube and its finite element model (Example 3). (a) Cylindrical tube details. (b) Finite element model.

Lamina 2: (middle)

k, = 0.1 + 0.9 x IO-*T - 0.3 x 10e6T2 klz = 1.0 + 0.9 x lo-*T - 0.3 x 10-6i”2

Lamina 3: (outside)

k,, = 20.0 + 1.0 x lo-*T k,, = 1.0 + 1.0 x lo-*T.

The tube is modelled using a graded mesh of three elements as shown in Figure 6(b). The p-level in the axial direction (pt) of the elements is increased from 2-5. The p-levels in the radial direction (“11)of the three laminas are chosen to be the same and are progressively increased from 2-5. The tube is also modelled using p-version formulation of Surana and Orth [42]. For comparison purposes, the same graded mesh is also used for this formulation. Figures 7 and 8 show plots of linear and non-linear axial temperature distributions at the inside surface of the tube. Similar responses for the outside surface are shown in Figs 9 and 10. It is to be noted here that the linear and non-linear temperature response of the composite tube are quite different in nature and magnitude. In both cases, the present formulation yields to the ‘peak’ values at the center of the heat flux (i.e., at 2 = -12) at relatively lower pq levels. The results for “p, = 3-5 are almost indistinguishable. The axial temperature distribution simulated by the

I

1

46.0

f/-~‘“.,,,~ Formulation

42.0

of Surono

and Orth

1421

40.0 0

-2

-4

-6

Distance

along

-0 the

-10

-12

tube

Fig. 7. Axial temperature distribution (linear case) at the

inside surface (r = 1.0) (Example 3).

p-Version axisymmetric shell element

50.0 ~ -

Finite

Present

46.0

-

44.0

v--c

Element

Solutton

‘;;~I E Finite

,

(r=l)

15

formulation

Element

Solution

(r=19)

115.0

e z e $ ;

,,’ ,’

-

,’ 42.0

I

,’

’

Formulation

d’

of Surano

x---PI=

and Orth

(421

,

2 p,=

A- - pt=5 p,=5 0 - - pt= l---pt=ll

z’

a p.=a p*=ll

i

65.0

J

0

-2

-4

-6

-8

Distance

along

the

-10

40.0

-12

0

-2

tube

-4

-6

Distance

along

-8 the

-10

-12

tube

Fig. 8. Axial temperature distribution (non-linear case) at the inside surface (r = 1.0) (Example 3).

10. Axial temperature distribution (non-linear case) at the outside surface (r = 19.0) (Example 3).

formulation of Surana and Orth [42] completely fails to capture these peak values immediately beneath the heat flux. Figures 11 and 12 show plots of linear and non-linear radial temperature profiles for both for-

mulations. From Fig. 11, it can be seen that the highest temperature (for pc =p,, = 11) calculated by the formulation of Surana and Orth [42] at the outside surface is nearly half the converged values

600.0

,

466.0

_

376.0

-

600.0

Present

formulation

466.0

c

Finite

Present

Element

Solution

(z=-12)

_

formulation

376.0 ? a e $ : 264.0 Formulation

of Surono

and

Orth

[421

x---pt= 2 p,=2 h - - PI= 5 p,= 5 152.0

0

a p.=a

Q - -pr= l---p,=ll

-

-2

p-=11

-4

-6

Distance

along

152.0

-a the

-10

-12

tube

Fig. 9. Axial temperature distribution (linear case) at the outside surface (r = 19.0) (Example 3). CAS 47,1-s

10 Radial

Distance

Fig. 11 Radial temperature distribution (linear case) at the middle surface (z = - 12.0) (Example 3).

A.

16

BOSE and

reached by the present formulation. These profiles and their ‘lack of C’continuity at the inter-lamina boundaries, as seen from these graphs are a testimony to the performance of the present formulation in accurately simulating the physics that governs the linear and non-linear heat conduction in laminated composites. This non-differential nature of the radial temperature distribution cannot be simulated by the formulation of Surana and Orth [42]. Figure 13 shows a graph of the quadratic functional (I) versus degrees of freedom for the linear case, i.e., using the non-temperature dependent coefficients of the polynomials used for lamina conductivities. The quadratic functional converges rapidly for the present formulation as “p,, is increased for the lamina. The quadratic functional values obtained using the formulation of Surana and Orth [42] are far from being converged, even for ps as high as 11. Figure 14 shows a convergence plot of the ‘peak temperature at the center of the heat flux (i.e., at r = 19 and z = - 12) with number of equilibrium iterations using the modified Newton-Raphson procedure. As expected, the presented formulation converges quickly to the ‘peak’ value at a relatively lower p-level in the < and 11directions. On the contrary, the formulation based on single higher order approximation theory converges to a smaller peak value even at high p-levels. SUMMARY

AND CONCLUSIONS

A finite element formulation has been presented for the laminate three-node axisymmetric shell element

K. S.

SURANA -20.0

o-O-0.

-66.7

-s

-

‘a

-*

\

\

q . bo\

s

-113.3

x 3

\

-

\ b

c

2 %

z

0

-206.7

x-Present formulation 0 - -Formulation of Surano

-253.3

and Orth

I421

t t -300.0

1 0

’

’

’

100 Number

’ 200

’

’ 300

of Degrees

’

’ 400

’

of Freedom

’ 500 (N)

’ 600

Fig. 13. Quadratic functional vs degrees of freedom (linear case) (Example 3).

based on piecewise hierarchical p-version approach for steady-state heat conduction. The order of approximation in the longitudinal direction (e) of the element for all lamina is the same (pc) and can be chosen arbitrarily. The order of the temperature 600.0

135.0

Finite

Present

116.0

Element

Solution

(z=-12)

formulation resent

-

Formulation of Surona and Orth x---pr=

2 p.=

- A- - pt=5 .---p,=

l-pt=ll

p.= % p,=

Formulation of Surano and Orth

formulation

[421

x---p,=

2

p,=

.--p,=% .-.-p,=

,,

p,=F p,=

,

[421

, 5 %

p,=ll

160.0

-I

40.0

0

10

5 Radial

15

20

Distance

Fig. 12. Radial temperature distribution (non-linear case) at the middle surface (z = - 12.0) (Example 3).

50.0 1

5

12

8 Iteration

15

Number

Fig. 14. Convergence of temperature (non-linear case) at the center of heat flux (Example 3).

p-Version axisymmetric shell element approximation (“p,,) in the transverse direction (%l) of the laminas can be different for each lamina. The temperature approximation for the element is developed first by establishing a hierarchical temperature approximation for each lamina of the laminate and then by imposing inter-lamina continuity conditions of temperature at the interface between the laminas. The approximation functions and the nodal variables for each lamina are derived directly from the Lagrange family of interpolation functions of order p( and “p,,. This is accomplished by first establishing one-dimensional hierarchical approximation functions and the corresponding nodal variable operators in the [ and krl directions for the three- and one-node equivalent configurations that correspond to pe + 1 and “p,,+ 1 equally spaced nodes in the { and krl directions and then taking their products. The nodal variables for the entire laminate are derived from the nodal variables of the laminas and the inter-lamina continuity conditions of temperature. The element formulation ensures Co continuity or smoothness of temperature across inter-element as well as the interlamina boundaries. Weak formulation of the Fourier heat conduction equation in the cylindrical coordinate system is constructed and the individual lamina matrices and equivalent nodal heat vectors are derived using this weak formulation and the hierarchical temperature approximation for the laminas. Inter-lamina continuity conditions are used to construct the transformation matrices for the laminas. These matrices permit transformation of the lamina degrees of freedom to the laminate degrees of freedom. Using these transformation matrices individual lamina matrices and the equivalent heat vectors are transformed and then summed to obtain the laminate element matrix and equivalent heat vectors. There is no restriction on the number of laminas and their lay-up pattern. Each layer can be generally orthotropic. The material directions and the layer thickness may vary from point to point within each lamina. The geometry of the axisymmetric shell element is defined by the nodes located at the middle surface of the element and the lamina thicknesses. The formulation allows temperature dependent thermal conductivities, internal heat generation, film coefficients and the radiation effects. The numerical examples presented for the laminated composite tubes demonstrate extremely good performance of the present formulation in terms of the accuracy of the calculated temperature. Since the physics of heat conduction in laminated composite requires a piecewise behavior, the present formulation is ideally suited to analyze such problems. The numerical examples clearly demonstrate the ability of the formulation to yield inter-lamina continuity of temperature and discontinuous first derivatives of the temperature in the lamina lay-up direction. As the thermal conductivities of the adjacent laminas differ more and more, this effect becomes more pronounced (Examples 1 and 2) and cannot be simulated accu-

17

rately without the present formulation (or extremely elaborate h-models). The formulation of Surana and Orth [42] where the temperature distribution for the laminate in the direction transverse to the lamina lay-ups is approximated by a single polynomial of higher order totally fails to simulate the physics of heat conduction (which requires only Co continuity in the lay-up direction) in the vicinity and at the interlamina boundaries. Such formulation may produce highly inaccurate temperature values (depending on lamina conductivities) and wrong slopes of the temperature at the inter-lamina boundaries (Examples 2 and 3). Relatively low p-levels in the krl direction produce extremely good results. It is needless to overemphasize the extremely simple models used for all the both examples to obtain high degree of accuracy and faster convergence of results. The polynomial orders in the < direction (PC) and the krl direction (“p,) for each lamina can be individually selected to obtain a faster convergence rate. The formulation is equally effective for extremely thin as well as thick axisymmetric shells and provides a highly accurate and extremely simple (in terms of finite element models) means of numerical simulation of temperature distribution in laminated composite axisymmetric shells. The modified Newton-Raphson procedure yields a symmetric system and provides good convergence during equilibrium iterations. The hierarchical properties of the formulation yield a computationally efficient system at the element level as well as at the assembly and solution level. The computational effort devoted to lower p-levels need not be repeated when computations are performed for higher p-levels. Acknowledgements-Part of this paper (linear heat conduction) was presented at the 1991 ASME National Heat Transfer Conference, Minneapolis, MN, July 28-31, 1991. The computing facilities provided by the Computational Mechanics Laboratory (CML) of the Department of Mechanical Engineering are gratefully acknowledged. This paper in whole or part is accepted for publication at the Energy Technology Conference and Exhibition (ETCE), Composite Mechanics Symposium, Houston, TX, January 26-30, 1992. REFERENCES

1. 0. C. Zienkiewicz and Y. K. Cheung, Finite elements in the solutions of field problems. The Engineer 220, 507-510 (1965).

2. E. L. Wilson and R. E. Nickell, Application of the finite element method to heat conduction analysis. Nuclear Engng Des. 4, 276-286 (1966).

3. 0. C. Zienkiewicz and C. J. Parekh, Transient field problems: two dimensional and three dimensional analysis by isoparametric finite elements. Inr. J. Numer. Meth. Engng 2, 61-71 (1970).

4. A. F. Emery and W. W. Carson, An evaluation of the use of the finite element method in the computation of temperature. Trans. ASME J. Heat Transfer 93, 13&145 (1971).

5. J. Donea, On the accuracy of finite element solutions to the transient heat conduction equation. Int. J. Numer. Meth. Engng 8, 103-I 10 (1974).

A.

18

J~OSE

and K. S.

6. J. C. Bruch and G. Zyvoloski, Transient two dimensional heat conduction problems solved by the finite element method. Int. J. Numer. Meth. Engng 8,481-494 (1974).

1. J. Padovan, Semi-analytical finite element procedures for conduction in anisotropic axisymmetric solids. Int. J. Numer. Meth. Engng 8, 295-310 (1974).

8. J. Padovan, Steady conduction of heat in linear and nonlinear fully anisotropic media by finite elements. Trans. ASME J. Heat Transfer 96, 313-318 (1974). 9. K. A. Rathjen and L. M. JiJi, Heat conduction with melting or freezing in a corner. Trans. ASME J. Heat Transfer 93, 101-109 (1971).

10. G. Comini, S. Delguidice, R. W. Lewis and 0. C. Zienkiewicz, Finite element solution of nonlinear heat conduction problems with special reference to phase change. Int. J. Namer. Meth. Engng 8, 613624 (1974).

Il. L. G. Tham and Y. K. Cheung, Numerical solution of heat conduction problems by parabolic time space element. ht. J. Numer. Meth. Engng 18,461-414 (1982). 12. J. Banasek, A conservative finite element method for heat conduction problems. Int. J. Numer. Meth. Engng 20, 2033-2050 (1984). 13. J. P. S. Azevedo and L. C. Wrobel, Nonlinear heat conduction in composite bodies: a boundary element formulation. Int. J. Numer. Meth. Engng 26, 19-38 (1988).

14. S. J. Fan and C. I. Tsai, Finite element analysis of welding thermal behavior in transient conditions. American Society of Mechanical Papers, 84-HT-80 (1984).

Engineering

Technical

15. J. Goldak, A. Chakravarti and M. Bibby. A new finite element model for welding heat sources. Metallurgical Transactions BlSB, 299-305 (1984). 16. M. Kleiber and A. Sluzalec, Jr, Finite element analysis of heat flow in friction welding. Rozprawy Inzynierski, Polish Academy of Sciences, Institute of Fundamental Technological Research, Warsaw, Poland, 32, 107-l 13 (1984). 17 I. Babuska and M. R. Dorr, Error estimates for the combined h and p-versions of the finite element method. Numersche Mathematik 31, 251-211 (1984). 18. I. Babuska, I. N. Katz, B. A. Szabo and A. P. Greens-

felder, Hierarchic families for the p-version of the finite element method. Proceedings of the Third IMACS International Symposium, pp. 278-286, Lehigh University, Bethlehem, PA (1979). 19. A. G. Peano. Hierarchies of conforming finite elements. Doctoral dissertation, Washington University, St Louis, MO (1979). 20. 1. Babuska, B. A. Szabo and I. N. Katz, The p-version of the finite element method. SIAM Numer. Anal. 18, 515-545 (1981).

21. P. K. Basu and R. M. Lamprecht, Some trends in computerized stress analysis. Proceedings of 7th ASCE Conference Electronic Comput., pp. 312-330, St. Louis (1979). 22. P. K. Basu, B. A. Szabo and M. P. Rossow, Theoretical manual and user’s guide to COMET-X, an advanced computer program for stress analysis. Federal Railroad Administration, FRA/ORD-II/60 (1977). 23. L. J. Brombolich and P. L. Gould, Finite element analysis of shells of revolution by minimization of the potential energy functional. Proceedings of Symposium on Applied Finite Element Method, pp. 219-307, Civil Engineering Department, Vanderbilt University (lY6Y).

SURANA

24. J. Jirousek, Hybrid-Trefftz plate bending elements with p-method capabilities. Int. J. Numer. Meth. Engng 24, 1367-1393 (1987).

25. B. A. Szabo and K. A. Mehta, p-Convergent finite element approximations in fracture mechanics. Int. J. Numer. Meth. Engng 12, 551-560 (1978).

26. B. A. Szabo, P. K. Basu and M. P. Rossow, Adaptive finite element analysis based on p-convergence. NASA Publication No. 2059, pp. 43-50 (1978). 21. D. W. Wang, I. N. Katz and B. A. Sabo. Implementation of a r?’ triangular element based on the pIversion of the finite element method. NASA Publication No. 2245, 153-170 (1982). 28. K. S. Woo and P. K. Basu, Analysis of cylindrical shells by p-version of finite element method. Int. J. Solids Struct. 25, 151-165 (1989).

29. I. Babuska, 0. C. Zienkiewicz, J. Gago and E. R. de A. Oliveira, Accuracy Estimates and Adaptive Refinements in Finite Element Computations. John Wiley, New York (1986). 30. K. S. Surana and P. Kalim, Finite element formulation for axisymmetric shell heat conduction with temperature gradients. Presented at The Winter Annual Meeting of ASME, pp. 17-21, Miami Beach, Florida (1985). 31. P. Kalim, Isoparametric line and transition finite element with temperature and temperature gradients for axisymmetric and planar two dimensional heat conduction. Ph.D. thesis, University of Kansas, Lawrence, KS (1985). 32. G. J. Desalvo and J. A. Swanson, ANSYS (Engineering Analysis System) User’s Manual. Swanson Analysis System, Houston, PA (1985). 33. The MacNeal Schwendler Corporation, NASTRAN User’s Manual (1974).

34. K. S. Surana, J. E. Basas, C. J. Parekh and R. Prabhakaran, NISA (Numerically Integrated Elements for System Analysis) User’s Manual. McAuto in Association with Engineering Mechanics Research Corporation (1976). 35. K. S. Surana and R. K. Phillips, Three dimensional curved shell finite elements for heat conduction. Comput. Struct. 25, 775-785 (1987).

36. R. K. Phillips, Three dimensional curved shell and solid-shell transition finite elements with temperature and temperature gradients as primary variables for steady state heat conduction. Ph.D. thesis, University of Kansas, Lawrence, KS (1986). 31. K. S. Surana and N. J. Orth, Axisymmetric shell elements for heat conduction with p-approximation in the thickness direction. Comput. Struct. 33, 689-705 (1989).

38. K. S. Surana and N. J. Orth, p-Approximation axisymmetric shell elements for heat conduction in laminated composites. Compur. Struct. 33, 1251-1265 (1989). 39. K. S. Surana and N. J. Orth, p-Approximation axisymmetric shell elements for heat conduction in laminated composites. 26th National Heat Transfer Conference, Phi!adelphia, PA (1989). 40. K. S. Surana and G. Abusaleh, Curved shell elements for heat conduction with p-approximation in the shell thickness direction. Comvut. Struct. 34.861-880 (1990). 41. K. S. Surana and N. JI Orth, p-Version hierarchical three-dimensional curved shell element for heat conduction. Comput. Mech. 7, 341-353 (1991). 42. K. S. Surana and N. J. Orth, Completely hierarchical axisymmetric shell element based on p-version for heat conduction in laminated composites. Comput. Struct. 46, 777-789 (1993).