Curved shell elements for heat conduction with p-approximation in the shell thickness direction

CURVED SHELL ELEMENTS FOR HEAT CONDUCTION WITH p-APPROXIMATION IN THE SHELL THICKNESS DIRECTION K. S. SURANAand G. ABUSALEH The Unive~ity of Kansas, Department of Mechanical Engineering, 3013 Learned Hall, Lawrence, KS 66045, U.S.A. (Received 7 March 1989) Ah&act-A finite element formulation is presented for the curved shell elements for heat conduction where the element temperature approximation in the shell thickness direction can be of an arbitrary polynomial order p. This is accomplished by introducing additional nodal variables in the element approximation corresponding to the complete Lagrange interpolating polynomials in the shell thickness direction. This family of elements has the important hierarchical property, i.e. the element properties corresponding to an approximation order p are a subset of the element properties corresponding to an approximation order p + 1. The formulation also enforces continuity or smoothness of temperature across the inter-element boundaries, i.e. C? continuity is guaranteed. The curved shell geometry is constructed using the co-ordinates of the nodes lying on the middle surface of the shell and the nodal point normals to the middle surface. The element temperature field is defined in terms of hiera~hical element appro~mation munitions, nodal tem~ratu~s and the derivatives of the nodal tem~ratures in the element thickness direction ~or~s~n~ng to the complete Lagrange interpolating polynomials. The weak formulation [or the quadratic functional) of the three-dimensional Fourier heat conduction equation is constructed in the Cartesian co-ordinate space. The element properties of the curved shell elements are then derived using the weak formulation (or the quadratic functional) and the hierarchical element approximation. The element matrices and the equivalent heat vectors (resulting from distributed heat flux, convective boundaries and internal heat generation) are all of hierarchical nature. The element formulation permits any desired order of temperature distribution through the shell thickness. A number of numerical examples are presented to demonstrate the superiority, efficiency and accuracy of the present formulation and the results are also compared with the analytical solutions. For the first three examples, the h-approximation results are also presented for comparison purposes.

NOTATION

bilinear functional of a test function e and temperature T represents a typical element vector of equivalent nodal quantities and the secondary variables element matrix quadratic functional for an element e quadratic functional for the entire model quadratic functional per unit volume for the entire model Jacobian matrix matrix of material thermal conductivities thermal conductivities in the global x, y, z co-ordinate space linear functional of test function 0 total number of nodes for the element total number of degrees of freedoms for the entire model, i.e. total number of nodal variables for the entire model direction cosines of an outward normal to the element surface or boundary two-dimensional approximation function for a node i 861

Q 4

[RI

one-dimensional (in c) approximation for a node i approximation function matrix for a node i internal heat generation per unit mass distributed heat flux per unit area matrix relating generalized nodal variables to the derivatives of temperature (T) with respect to x, y and z temperature temperature at a node i hierarchical nodal variables at a node i specified temperature ambient temperature shell thickness at a node i test function normal to the middle surface of the shell at a node i such that //Fri 11= r, with respect to Cartesian co-ordinate directions global co-ordinates of a node i film coefficient or heat transfer coefficient domain of an element (volume in this case) element boundaries on which distributed heat flux and convection act boundaries on which temperature T is specified

K. S. SURANA and G. ABUSALEH

862

vectors of generalized nodal variables a vector of generalized variables at a

node i natural co-ordinates co-ordinate in the thickness direction incremental surface area on an element boundary angular co-ordinate (of ROZ cylindrical system)

INTRODUCTION

The applications of the finite element method to boundary value and initial value problems have been increasing steadily in recent years since the method has been given firm mathematical foundation. The origin of the modern finite element method and its formalization in terms of variational methods can be found in [l-8]. The finite element method for heat conduction studies was first proposed by Zienkiewicz and Cheung [9]. They proposed a finite element solution for steady state heat conduction problems governed by the quasi-harmonic equation and they demonstrated the power of their method by analyzing temperature distribution within complex geometries. A lot of research was done in the following years. Wilson and Nickel1 [lo] extended the stationary formulation to transient heat conduction problems. They formulated time dependent finite elements for two-dimensional axisymmetric geometries. Zienkiewicz and Parekh [ 1l] proposed a finite element solution for two-dimensional and three-dimensional transient heat conduction problems. They used isoparametric elements, where approximation functions for the primary variables and shape functions for geometric description are the same, to approximate the solution in space and then developed a recurressive process to approximate the solution in time. Zienkiewicz and Parekh were the first to use isoparametric finite elements to solve the heat conduction problem. Bruch and Zyvoloski [ 121 proposed a finite element solution for the twodimensional transient heat conduction problem. Subsequent work by many different researchers came about during the same period of time. Padovan [13] reported work on heat conduction for linear and non-linear anisotropic materials, where material properties are functions of spatial variables and temperature. Comini and Delguidice [14] also reported a finite element solution for non-linear heat conduction problems. The conservative finite element method in heat conduction was first reported by Banasek [15], who used the concept of energy balance to formulate the element properties. His proposed solution guarantees local and global energy conservation but may be more difficult to program. Other researchers applied the finite element method to solve field problems such as welding. Work in this area was reported by Goldak et al. [ 161, Kleiber and Shtzalec [ 171 and by Fan and Tsai [18].

The p-approximation technique in the finite element method is a relatively new approach that came about in the 1970s. “p-Approximation” refers to polynomial approximation. In this approach the number of elements and their sizes remain fixed and the accuracy of the solution is improved by increasing the order of approximation for the elements (thereby increasing the total number of degrees of freedom for the entire model). In the h-approximation approach, the solution accuracy is improved by reducing the element sizes, thereby increasing the number of degrees of freedom for the entire model (mesh refinement). Thus, one obvious advantage of the p-approximation approach is that results with progressively increased accuracy can be obtained from only one finite element mesh. The element approximation functions and the element properties resulting from the p-approximation approach are hierarchical in nature, i.e. the properties for p th order approximation are a subset of (p + 1)th order approximation. This important property results in computational efficiency not only at the element level but also at the assembly and solution of the resulting equations for the entire model. In the following, a brief survey of the p-approximation research literature is presented. In 1978 Szabo and Mehta[19] proposed a papproximation finite element procedure to estimate the strain energy release rate in fracture mechanics applications. They increased the order of approximation in every element while keeping the finite element mesh fixed. Their conclusions indicated that the errors in strain energy and strain energy release rates are inversely proportional to the number of degrees of freedom when the order of approximation is uniformly increased for the entire mesh, while they are inversely proportional to the square root of the number of degrees of freedom if the h-approximation concept is employed. Katz [20] implemented hierarchical finite elements using the p-approximation approach and he demonstrated that the convergence rate for singularity problems using the p-approximation is twice that of the h-approximation. Shortly after that, Katz and Wang[21] reported similar results on problems with singular solutions caused by corners. Anderson and Falk [22] performed linear elastic and small strain analysis on threedimensional domains where the basis functions are complete polynomials of the order 12 or less. Their results obtained by adaptive schemes with p-approximation show a high rate of convergence in the energy norm compared to other schemes with the h-approximation. Bennighof and Meirovitch[23] reported that the convergence rates of the p-version are superior to those of the h-version. They also demonstrated the superiority of the p-approximation in obtaining lower eigenvalues when used in structural dynamics applications. A complete list of references on p-

863

Curved shell elements for heat conduction approximation technique can be found in a recent publication edited by Babuska et al. [24]. Though the p-approximation finite element concept has been applied to the two-dimensional field problems, there is no reported work for threedimensional curved shell heat conduction using this approach. In a recent paper Surana and Kalim [25] (also Kalim [26]) formulated axisymmetric shell elements for heat conduction where temperature as well as temperature gradients were retained as nodal variables at the element nodes. Many commercial finite element systems [27-291 report two-node axisymmetric shell and curved shell elements with constant temperature through the shell thickness. Three-dimensional curved shell elements constitute an extremely important class of finite elements for analyzing pressure vessels, piping system and many other sheet metal fabrications. In another recent paper Surana and Phillips [30] (also Phillips [31]) presented a formulation for curved shell elements for three-dimensional heat conduction where temperature and temperature gradient in the shell thickness direction were considered as nodal variables at the element nodes, thereby permitting linear temperature distribution through the shell thickness. In many thick shell applications where shell surfaces are subjected to extreme conditions, the linear temperature distribution through the shell thickness is hardly adequate. For efficient modeling and accurate heat conduction studies in shell structures of this type, curved shell elements with higher order temperature distribution through the shell thickness are essential. In this paper, three-dimensional curved shell elements are presented where the element temperature approximation in the shell thickness direction (0 can be of an arbitrary order p. This is accomplished by introducing additional nodal variables at the middle surface of the element corresponding to the complete Lagrange interpolating polynomials in the direction of the shell thickness. The element approximation functions and the nodal variables resulting from this formulation have a hierarchical property, i.e. the approximation functions and the nodal variables corresponding to an approximation order p are a subset of those corresponding to order p + 1. The weak formulation (or the quadratic functional) of the Fourier heat conduction equation is constructed in the Cartesian co-ordinate space. The element properties are then derived using the quadratic functional and the hierarchical element approximation, The element matrices and the equivalent heat vectors resulting from such process are also hierarchical. The order of approximation in the plane of the element may be selected according to the need (linear, parabolic, cubic, etc.) and is kept fixed for a given application. Numerical examples are presented to demonstrate the accuracy, effectiveness and superiority of these elements over existing shell formulations and even

three-dimensional solid elements. Comparisons are also made with the analytical results. The results obtained from the p-approximation models using these elements are also compared with those obtained using h-approximation and the convergence graph showing the rates at which the errors are decreasing is presented for both p and h approximations. The superiority of the present formulation over existing curved shell elements and h-approximation models using three-dimensional solid elements is clearly illustrated. WEAK FORMULATION

The Fourier heat conduction equation in the Cartesian co-ordinate system x, y, z can be written as

3c+%+~+Qp=o

(1)

ay

with gl%+W,+g,~,+q

+BV-

T,)=O

on rl

(2)

and T=T,

on TZ,

(3)

where

aT

&=K&K,--+Kzz;.

The above equations compact form

(4)

ay

can be written

in a more

with it.{[Kl{VT}}+q+/3(T-T,)=O T=T,

onr,

onr,,

(6) (7)

where {V}T= and

[

&g

Kx Kxy Kxz 1 1

WI = K,x

K,

K,rz .

[ K,,

K,

K,

(8)

(9)

K.

844

S. SURANA

and G.

Here we assume that the matrix [K] of thermal conductivities is symmetric, i.e. Ky.x= Kxy , Kz.r= K,, and KY = KYr. The weak formulation of the Fourier heat conduction eqn (1) [or eqn (5)] is constructed by taking the scalar product of the differential eqn (1) with a test function (say v) over a typical element domain R’. dfZ=Q,

w

ABUSALEH

and

It is rather obvious that B(v, T) is a bilinear and symmetric functional and Z(v) is a linear functional associated with the weak formulation given by eqn (15). The quadratic functional associated with this weak fo~ulation can be easily written as I’(T) =

where df2 = dx dy dz.

(11)

Through integration by parts, one order of differentiation can be transferred from the functions g, , g, and g, to the test function v.

$3(T, T) - Z(T)

(18)

or I’(T)=;

{VT.Jr[K]{VT}dR+; s*

The weak formulation given by eqn (15) or the quadratic functional given by eqn (19) forms the basis for deriving the properties of the curved shell elements presented here. Here we clearly see that Tis the primary variable and (g,nX + g2nY+ g,n,) is the secondary variable. We also note that eqn (2) represents natural boundary conditions whereas essential boundary conditions are given by eqn (3). Substituting from (2) into (12), utilizing (3) and then rearranging the terms, we can write the following: VT@df

The first term of eqn (13) can also be written in a slightly different form, giving

wwuV~~

s ne

d&-B f ~~,~T~df=~~v~Q~ +

I Fl

UT,/3 dl- -

uq dl-

(14)

s ri

or B(v, T) = l(v),

(15)

where E(Y, T) =

(Vvj~KffVT) s i%

dR +

VT/~dT

s i-1

(16)

ELEMENT

APPROXIMATION

Geometry

The geometry of the curved shell elements is derived from the three-dimensional solid elements in which two higher order surfaces each containing an equal number of nodes are connected by linear edges (Fig. la and b). In Fig. la, a l&node threedimensional solid element is shown where two eightnode parabolic faces are connected by linear edges. Likewise, Fig. lb shows a 24-node cubic-linear element. The external faces of these elements are curved, while the sections across the thickness are generated by straight lines. For such solid elements each pair of nodes i, and ib with given Cartesian co-ordinates completely describes the element geometry. In describing the configuration of the curved shell elements this geometry (of Fig. la and b) is replaced by a middle surface containing the nodes and the normals to these nodes at the middIe surface such that if pXi is such a normal at node i then /Ip3i /I = r, represents the nodal thickness of the shell element at node i. By using such procedures the equations describing the geometry of the curved shell elements shown in Fig. lc and d can be derived from those of Fig. la and b. For a curved shell element with n nodes the coordinates x, y, z of a point P within the element can be determined using

Curved shell elements for heat conduction

3

(b)

(a)

X

3 NormalT (cf

at node i

(d)

Fig. 1. Thick shell and curved shell finite elements. (a) Parabolic thick shell element. (b) Cubic thick shell element. (c) Parabolic curved shell element. (d) Cubic curved shell element.

where 5 and y are the natural co-ordinates in the plane of the element and [ is the third natural co-ordinate in the thickness direction. Temperature field Figure 2a-c show a special class of threedimensional solid elements in which the temperature field is parabolic (serendipity family) in the plane of the element but linear, parabolic and cubic in the thickness direction. The temperature approximation for such elements can be easily derived by using nodal temperatures and the nodal approximation functions. Furthermore the nodal approximation functions can be simply obtained by taking the tensor product of the inpiane approximation functions N,({, q), i=I,...,n, with the one-dimensional approximation functions iir,([), k = 1, . , . , nyr where E<,,are the number of nodes in the plane of the element and ni are the number of nodes in the element thickness direction. By using such procedure higher order temperature approximation in the element thickness direction can be easily established. The number of nodes in the plane of the element can be selected according to the requirement (four nodes for linear temperature, 12 nodes for serendipity cubic element

etc.). As is obvious from Fig. 2a-c, as the order of the temperature approximation in thickness direction increases, the number of nodes required also increases dramatic~ly, which makes these elements completely impractical for applications. Within the assumptions stated here the geometry of the elements shown in Fig. 2a-c can be replaced by the element shown in Fig. 2d where the nodes are located at the middle surface of the element and the magnitudes of the nodal normals define the nodal thicknesses, thereby establishing the threedimensional geometry of the element. If the element of Fig. 2d is to have linear, parabolic, cubic, etc., temperature approximation in the thickness direction (corresponding to Fig. 2a-c) then the nodes of the element in Fig. 2d must have nodal temperatures and additional nodal variables as degrees of freedom corresponding to the complete Lagrange interpolating polynomials in the thickness direction. To illustrate the process of deriving hierarchical approximation functions and the corresponding nodal variables, consider the elements shown in Fig. 3a and b. For the element of Fig. 3a, we can write

866

K. S. SURANA and G. ABUSALEH

e:z$G??, Direction

(4

(b)

e-%g$@ Direction

Cd)

Cc)

Fig. 2. Solids and shells with higher order approximation in the thickness direction. (a) Parabolic thick shell-linear temperature approximation in thickness direction. (b) Parabolic thick shell-parabolic temperature approximation in thickness direction. (c) Parabolic thick shell-cubic temperature approximation in the thickness direction. (d) curved shell element with nodes on middle surface.

Substituting

where

from eqn (25) into (24) we obtain ni4

T(t,~,i)=

C Ni(5,~) i= I

TI+i(T,-T,,)

(22) +‘(Ti,-2Tn+Tt3)

and nc,, is the number of nodes in a & plane. Substituting from eqn (22) into (21) we obtain:

2

>

. (26)

Similarly for the element of Fig. 3c we can write

(( >

XFor the element

T, + T/Z 2

+;(T,-

Ti,) . (23)

of Fig. 3b we can write +27T,-TM)+$(9Ti,

*tn

- 27T,2 - 27T,, + 9T&

+ fikZ)r,

+ %(Wi&

(24) +;(-9Ti,

where

1

+ 27T, - 27T,3 + 9Ti4) .

(27)

fi,(r;)=ir(r fi#)= &(i)

- 1)

l-i2 = trci + 1).

(25)

For the elements of Fig. 3a-c at a location i, the linear, parabolic and cubic Lagrange polynomial can be established in the y coordinate system (- t,/2 < y < ti/2) using the temperatures (c, , T,&,

Curved shell elements for heat conduction

867

iocation

i,

Node i

W Fig. 3. Thick shells with nodal variables in y space. (a) Linear temperature approximation in thickness direction. (b) Parabolic temperature approximation in the thickness direction. (b) Parabolic t~rn~mture approximation in the thickness direction. (c) Cubic temperature approximation in the thickness direction. (d) Curved shell with nodes on middle surface.

These are given by

and

Using eqns (2&)-(30), the temperature YT~and its derivatives with respect to y can be determined at y = 0 (corresponding to the midplane element of Fig. 3d). These are given in Table I. Substituting from Table 1 for 1;:, iWJt?y, . . . , etc. into eqns (23), (26) and (27) we obtain

(32)

Equations (31x33) represent the temperature approximations (linear, parabolic and cubic in the

K.S. SIJRANA and G. ABUSALEH

868

Table 1. Derivatives of temperature w.r.t. 7 at a shell node f for linear, parabolic and cubic p-levels. Configuration

For location aty=O

Equation

i

Fig. 3a

Fig. 36

(32)

Fig. 3c

(33)

$&-277;,+27q3-Tr,) , $=$“Ti,-27T,z-27T,3+9T,) I

aq_3

,,,-$-9T,r I

thickness direction) for the curved shell element of Fig, 3d. Thus if p is the polynomial or&r of temperature approximation in the thickness direction then eqns (31)-(33) correspond to p = 1, 2 and 3, respectively. Using the process of induction, for an element with p order of approximation in the thickness direction we can write (dropping the [q subscript for n)

+27T,z-27&+9Ts,)

Equation (36) describes the hierarchical approximation functions for a node i of the shell element. The corresponding hierarchical nodal variables for the node i are given by eqn (37).

ELEMENTPROPERTIES The derivatives of the temperature with respect to x, y and z (required in the weak formulation (14) or the quadratic functional (19)) can be obtained by using the Jacobian of transfo~ation (between 5, q, 6 and x, y, z co-ordinate spaces) and we can write or

(38)

where a typical [nil and {ai}, for a node i (with p order of approximation in thickness direction) can be written as The partial derivatives of the temperature T(& ‘I, [) with respect to 5, q and i can be easily obtained using the element temperature approximation (34). Upon

Curved shell elements for heat conduction substituting

869

from these into eqn (38) we obtain

= YIP

i=l,2,...,n

ary T) -=O;

ter

i=t2,...,n.

(39)

(43)

These equations can be arranged in the matrix form to give where [RI] for a typical node i is given by

or

A,, 4

A,

[RI = Bi

E,,

&

Bn

G

Gl

G

cl3

4

Equation

A,

1

. . . Bip . C,

WW%

= UL

(44)

(41)

(39) implies tbat

(42)

The derivatives of the temperature with respect to x, y and z can now be substituted from eqn (42) into the quadratic functional P(T) given by eqn (19). Minimization of I’(T) requires that aZ’(T) -=O; al;: azy T) -=oo; a!$ c >

i=l,2,...,n

i=1,2,...,n

The integrals over Zl refer to the surface integrals. These integrals are to be evaluated only for those surfaces that have convection and distributed heat flux q acting on them and thus appropriate elements of [q must be selected. Also, dZ must be transformed to e, Q and i space. Since both the approximation function matrix [@ and the [R] matrix are hierarchical, the hierarchical nature of [EZq and (f Je is rather obvious. It is also clear that the element matrices are symmetric (provided that Kyx= K,,, K, = Kxx and Kr,,= K,,,). The numerical values of the element matrices and the equivalent heat vectors are computed using Gaussian quadrature.

K. S.

870

SURANA and

NUMERICAL EXAMPLES

Numerical examples are presented here to establish the accuracy of the formulation and to demonstrate its applications. The examples are selected in such a way that their solutions can also be obtained using two-dimensional formulation as well as analytical methods. For each example the results obtained from the present formulation are also compared with those obtained from the comparable finite element models using mesh refinement (h-approximation). The convergence graphs showing the rate at which the errors are decreasing are also presented and results are compared with the h-approximation finite element models. The element formulation presented here has been implemented in a separate subsystem of FINESSE (Finite Element System [32]). During the implementation and processing of data for the examples,

x-

portion

G.

ABUSALEH

no conversion of units is done. Hence any set of consistent units either in FPS or the metric system can be selected. For this reason units are not given for the examples presented here. Example 1

Consider an infinitely long hollow circular tube of mean radius r,,, = 10 as shown in Fig. 4a. The inside surface of the tube has a distributed heat flux of q = 200 acting on it. The outer surface of the tube has convection with /l = 5 and T, = 30. Since the temperature does not vary axially and circumferentially, it suffices to model a radial sector of the tube (3” in this case) with an arbitrary axial dimensional (0.5 in this case). A finite element model of the tube using a single 6 node curved shell is shown in Fig. 4b. A comparable finite element model of the tube sector using three-dimensional solid elements is shown

modeled

4-+

---._--_--M __-------

Casea Case b

(4

t=16 t=O.l

rm = IO K,,=

KY,= K,,=

3.0

0.0

I

(b)

1, Id-node thick shell elements with equal thickness ( 96 nodes )

Fig. 4. An infinitely long cylindrical tube and its finite element models (example 1). (a) Cylindrical tube details. (b) Finite element model of the tube using curved shell element (one element, eight nodes). (c) Finite element model of the tube using solid elements (11 16-node solid elements, 96 nodes; eight 20-node solid elements, 104 nodes).

Curved shell elements for heat conduction

871

radial heat conduction

equation

(T),~~=~[Ln~)+~]+

0

2

4

6

10 12 14 16 18 20

8

RADIAL DISTANCE

Fig. 5. ~em~rature distribution along the radius of the tube for different p-levels (example 1).

in Fig. 4c. Here we consider two different tube thicknesses: t = 18 in. (r,,,/f = 0.55556) and d = 0.1 in, (r,/t = 100). Case a: (t = 18). Figure 5 shows the tem~rature distribution across the thickness of the tube for different p-levels. Results for p-levels of 5 through 11 are almost the same as the theoretical solution obtained from the solution of the one-dimensional 240

T,.

(47)

The temperatures at the inside, middle and outer surface of the tube as a function of p-level are shown in Fig. 6. It is interesting to note the non-monotonic convergence of the inside and the middle surface tem~ratures as the p-level is increased. The graphs of the quadratic functional per unit volume (1) versus the total number of degrees of freedom for p-approximation curved shells, 16-node thick shells (h-approximation) and 20-node threedimensional solids (h-approximation) are shown in Fig. 7. Figure 8 shows plots of the log of percentage error in 1 [theoretical value of f is computed using eqns (19) and (47)] versus log of the total number of degrees of freedom N for the curved shells, thick shells and the three-dimensional parabolic solid elements. The extremely high accuracy of i and the faster rate of convergence for the ~-approximation curved shells is rather obvious from Figs 7 and 8. It is also worth noting that the p-approximation results are obtained from a single finite element model (with element nodes at the middle surface of the tube) whereas each data point for the other two elements requires a new finite element model. Cafe 6: (t = 0.1). Figure 9 shows temperature distribution as a function of radial distance for different p-levels. The graphs of temperature versus p-level for the inside, middle and outer surface of the tube are shown in Fig. 10. Here, since the tube is extremely thin, as expected, a11p-levels give accurate results. It is worth noting that even with such small

220

-335I

200

i

180

:

I 160 i

ii

THICKNESS - 1s AMLYMAL SOLUTION FINITE ELEMENT SOLVIION ‘_____ INSIDE SURFACE (1.._... __.__ MIDDLE SURFACE .--___ OUTSIDE SURFACE

2 140 d li;i 5 120 t

<-

-341

% 3

-346

9 t zj

-354

v

iz

100 a0

$

-360

y

-366

2 g

-373

60 5: i9

-379

-385 0

12

3

4

5

6

7

8

9

1011

p-LEVEL

Fig. 6. Temperature versus p-level for inside, middle and outside surface of the tube (example 1).

1 1 I 0 N"Mf&

1 I 0F4iE&:

1 I

1 I

1 I

OF it&ED;i”(

1 120 N )

Fig. 7. Quadratic functional per unit volume (1) versus number of degrees of freedom (N) (example I).

K. S.

872

SURANA and

ABUSALEH

G.

t OY’JNDRICU TUBE I 0 0

S-noda

cumd

ahoIls

lbnode thkktill8 A ZO-nodeAid.

10' lo5

10'

10'

Lo9 ( N )

Fig. 8. Log of percentage error in i versus log of N (example I).

Example 2

segment using three-dimensional solid elements is shown in Fig. llc. The temperature distribution across the thickness of the sphere for different p-levels and the theoretical solution are shown in Fig. 12. The theoretical solution is given by

Here we consider a hollow sphere of mean radius 10 with a thickness of 18 as shown in Fig. lla. The inside surface of the sphere is subjected to a distributed heat flux of q = 1000 whereas the outer surface has convection with B = 5 and T, = 30. Due to point symmetry in meridinal and circumferential directions it suffices to model the sphere with one curved shell element as shown in Fig. 11b. A comparable finite element model of the sphere

Figure 13 shows plots of temperature versus p-level for inside, middle and outside surfaces of the sphere. The graphs of 1 versus N (for curved shell and ldnode thick shell models) are shown in Fig. 14. Plots of percentage error in i versus N (on a log-log

shell thickness there are no numerical difficulties encountered when higher p-levels of 8 through 11 are used and there is no deterioration of the computed results.

80

, [ CYUNDRICAL TUBE I MEAN RADIUS - 10 THICKNESS 0.1 D- - - - F.E. SOLUTION ( Au p-LEVELS BAN&YW.Xl SOLUTION

78 f )

76

75 g 2 d 5 P

74 -

;

74

73 _ 72 -

2 B k!

73

P

71 -

72

70 69 68 67 -

69

66 r 9.94

Fig.

9.

68

9.96

Temperature

10.00 10.02 RADIAL DISTANCE

9.98

versus radial distance (example 1).

10.04

10.06

I

I

I

I

I

I

I

I

I

’

01234567691011 p-LEVEL

for f = 0.1

Fig. 10. Temperature versus p-level for inside, middle and outside surface of the tube (example 1 - t = 0.1).

Curved

a13

shell elements for heat conduction

IIll = 10 KXX= Kyy= K,i K xF

K,,=

3.0

KY,= 0.0

(a)

lb)

11, 16-node thick shell elements with equal thickness ( 96 nodes )

Fig. 11. Hollow sphere and its finite element models (example 2). (a) Sphere details and material properties. (b) Finite element model of the sphere using curved shell element (I element, eight nodes). (c) Finite element model of the sphere using solid elements (11 16-node elements, 96 nodes).

400

I

360 320 260 2 240 z d 200 k! P

320

C SPHERE I MEAN RADIUS - 10 - 18 THICKNESS OANALYllCAL SOLUTlON FlNlTE ELEMENT SOLUTION a__.______,p-1 ~--_--_p-2 e_____ p-3 Op - 6-11

?

260 240

i

FINITE ELEMENT SOLUllON INSIDE SURFACE MIDDLE SURFACE OUTSIDE SURFACE

160 120 80

0

2

4

6

6 10 12 14 RADIAL DISTANCE

16

16

20

Fig. 12. Temperature distribution along the radius of the sphere for different p-levels (example 2). C.A.S. )4X-F

0

12

3

4

5

6

7

6

91011

p-LEVEL

Fig. 13. Temperature versus p-level for inside, middle and outside surface of the sphere (example 2).

K. S. SURANA and G. ABUSALEH

814 -370

<-

,

using three-dimensional solid elements is shown in Fig. 16~. Figures 17-19 show the plots of temperature versus angle 0 for inside, middle and outside surfaces of the torus tube obtained using different p-levels. Results for p = 6 through 11 are almost identical. The graphs of 1 versus N for p-approximation curved shell and ldnode thick she11 (h-approximation) element models are shown in Fig. 20. Figure 21 shows plots of log of percentage error in i versus the log of N for both models (1 for p = 11 is treated as the theoretical value of i). Again the faster convergence rate of i for the curved she11 element is rather obvious from Figs 20 and 21. For the curved shell model 50 elements are used along the circumference of the torus tube to reduce the modeling error in the geometry. Cubic or quartic elements (p = 3 and 4 in [ and 9 directions) can be easily used to completely eliminate the geometry error. Models using these elements will also result in a lower number of elements as well as a lower total number of nodes for the entire model.

I 0 0

-380

[ SPHERE 1 S-node cm/ad shsllo l&node thick ahella

Y 2 > t f

-390

-400

i J

-410

P P 4

-420

z

0

NUMB::

OF&RE;:

100 125 150 OF FREEDOM ( N )

Fig. 14. Functional f versus N ( example 2).

Example 4

scale) are presented in Fig. 15. Faster convergence of i for the curved shell model is quite clear from Figs 14 and 15. Example 3

Consider a hollow torus of circular cross-section and of mean radius (It,) 20 as shown in Fig. 16a. The circular cross-section of the torus has inside and outside radii of 1 and 19 (i.e. mean radius r, of 10 and a thickness of 18). The inside of the torus tube is subjected to a distributed heat flux q of 200 whereas the outer surface has convection with /I = 5 and T, = 30. Due to symmetry we need to consider half of the tube cross-section over a typical circumferential sector 4 (4 = 10” used in this case). Figure 16b shows the finite element mode1 of the torus sector using eight-node p-approximation curved shell elements. A typical finite element model of the torus

Consider a cylindrical tube with step change in the thickness, the details of which are shown in Fig. 22a. The inside surface of the tube is subjected to a distributed heat flux of 200 whereas the outside surface has convection with B = 5 and T, = 30. Since the temperature distribution in the tube is invariant of circumferential direction (e), it is necessary to mode1 only a sector of the tube (say 3”). Figure 22b shows a finite element mode1 of the tube using eight-node p-approximation curved shell elements. Figures 23-25 show graphs of temperature versus axial distance along the length of the tube at r = 8, 10 and 12 respectively as obtained using different p-levels. A plot of 1 versus N corresponding to different p-levels is shown in Fig. 26. Figure 27 shows a graph of the log of percentage error in 1 (obtained by treating j corresponding to p = 11 as an exact

[ SPHERE I

0 B-node q l&nods

curved shells thick shells

Fig. 15. Log of percentage error in 1 versus log of N (example 2).

875

Curved shell elements for heat conduction

A--’

Section

A-A

K,,= KY,=K,, 3.0 KXp K,,= KY,=0.0 (a) 50 equally spaced 16-node thick shell

(

3036 nodes)

50 equally soaced shell

Y

A

e

(b)

I I equally spaced t6-node

thick

S'ISllS CC)

Fig. 16. Hollow torus of circular cross-section and its finite element models (example 3). (a) Torus details and properties. (b) Finite element model of the torus using curved shell element (50 equally spaced

elements, 253 nodes). (c) Finite element model of the torus using solid elements (550 ldnode elements, 3036 nodes). value of i) versus log of N. Rapid convergence of the

results (temperature as well as i) is rather obvious from these graphs as the p-level is increased. This example clearly demonstrates that even when the region modeled has a sharp change in the thickness, the formulation is quite accurate and effective in predicting the temperature distribution. Furthermore, the temperature distribution across the stepped face is continuous. From the graph presented in Fig. 27 we note that even though the convergence is monotonic the rate of convergence is not monotonically increasing. This is due to the fact that the section of the tube with a step change in the thickness represents the singular region in the vicinity of which a combination of mesh refinement and increase in p-level is necessary to achieve the fastest rate of convergence.

SUMMARY

A finite element formulation for curved shell elements for heat conduction has been presented

where in the element temperature approximation in the shell thickness direction is hierarchical in nature and can be of an arbitrary polynomial order p. The element matrices and vectors are also hierarchical in nature, i.e. the element properties corresponding to an approximation order p are a subset of the element properties corresponding to an approximation order p + 1. The element formulation ensures continuity of temperature across the inter-element boundaries, i.e. C” continuity is guaranteed. All four numerical examples presented here demonstrate the superiority of the p-approximation curved shell element formulation in terms of the accuracy of the temperature results and the rate of convergence when compared with corresponding h-approximation models utilizing thick shells or even 20-node parabolic threedimensional solid elements. Retaining true temperature derivatives at the nodes, i.e. temperature derivatives in y space (i.e. element thickness space and not the Cspace), enables modeling of structures which exhibit sharp changes in thickness (example 4). The element gives continuous temperature distribution along a stepped face of example 4 (2 = -6).

43 -.._a-

. ..f. _.._*_

____ -_____--___ -

44

a. 1 p-3 p-4 p-5

p

p-2

p -

43

SF*“_.

.f.-..“..t-

42

8-11

--_c *_____”

41

k-2 P’f P’4

p = B-11

40 g

39

$38 g P

3J 3s 35 34 33 32 31 30

Fig. 1% ‘&?mperatureversusan& @fan the outside surface of the torus at differmt p&v& fexarnple 3).

The simplicity and the economy of the curved shell models when compared with three-dimensional solid models is quite obvious even for the relatively simple geometries used here. The polynomial approximation in the shell thickness direction can be

-

p-elf

increased to any desired level without changing the mesh (thus only a single finite element model is needed when curved shell elements are used) to obtain the required accuracy of temperature or the quadratic functional as seen in the numerical examples and the comparisons with the analytical solutions. As expected ~e~-a~~oxima~o~ shell elements yield the faster rate of convergence in aff cases (i.e. more accurate solution for a @en rmmber of degrees of freedom when compared with ~-approximation

<3Y -350 B t 4

-360

% s! g

-3711

f p s

-380

Curved shell elements for heat conduction

c4

<-

10’

.c E

loa

9)

2

10”

-I8 10-l

Fig. 21. Log of percentage error in P versus log of N (example 3).

Lick t2

+-

24

-

y=-

:_

L

_--_--_--_-__ _-

rmp lo ; t,=4

K xX= KY,= K,=3 K ry” K,,=

__._

_--

;

t,=e

K,=o

(a)

Fig. 22. Stepped cylindrical tube and its finite element model (example 4). (a) Stepped cyIindrica1 tube details. (b) Finite element model using curved shell elements.

K. S.

818 260

[ STEPPED

TUBE 1 F.E. SOLUTION p-1 p-2 p-3

. . . .._. ._ ______ -_-

r

-12

-11

-10

SURANA and

ABUSALEH

G.

( at r -

-7 -6 -6 DISTANCE ALONG

-9

S ) ________

-3

-5

-4

THE

Z-AXIS

-2

-1

0

Fig. 23. Temperature versus axial distance along the tube (at r = 8) for different p-levels (example 4).

STEPPED TUBE 1 F.E. .SOUJllOt’d( at r - 10 ) p-1 p-2 P’S

I

em.. . ..____ . . . .._.

190

L+

_I...

. . . . . _. . . . . __---_-

*...

p-5-11

160

z it P

170

160

I

150 -12

-10

-11

I

I

I -9

I

I

I

I

I

-5

-4

-3

-2

-1

I

-6 -8 -7 DISTANCE ALONG

0

THE Z-AXIS

Fig. 24. Temperature versus axial distance along the tube (at r = 10) for different p-levels (example 4).

130

:

1.20

C----..--...-

[ STEPPED

_......

110

_____..... _

.-....._______.____~~~ __

;

. ...

TUBEI 1 F.E.sownON( at r - 12 ) p-1 p-2 p-3

__-----

%* ._

p-5-11

70 y 60 = 50-12

I -11

I

I

-10

-9

I

I

I

-7 -6 -6 DISTANCE ALONG

I

I

I

I

I

-5

-4

-3

-2

-1

0

THE Z-AXIS

Fig. 25. Temperature versus axial distance along the tube (at r = 12) for different p-levels (example 4).

879

Curved shell eiements for heat conduction -1110

The hierarchical element properties yield computationally efficient systems at the element level as well as at the assembly and solution level. The computational effort devoted for a level p need not be repeated when computations are performed for level p + 1.

c STEPPED 0

B-node

TUBE I curved shells

rs

*c-1120 Y 3 9 +-I130

Acknowledgements--The computing facilities provided by the Comuutational Mechanics Laboratorv KML) of the Mechanidal Engineering Department of tie ‘Univ&sity of Kansas are gratefully acknowledged.

3 B is -1140 0 5

RETRENCH 1.R. Courant, Variational methods for the solution of problems of equilibrium and vibrations. BUN, Am. Math. Sot. 49, 1-23 (1943). 2. R. W. Clough, The finite element method in plane stress analysis. Proc. Am. Sot. Civil Engrs 87, 345-378 (1960). 3. J. F. Besseling, The complete analogy between the matrix equations and continuous field equations of structural analysis. Colloque International des Techniques de Calcul Analagique et Numerique de L’Aeronautique, Liege, pp, 223-242 (1963). 4. R. J. Meiosh, Basis for derivation of matrices for the direct stiffness method. AZAA Jnl 1, 1631-1637 (1963). 5. R. H. Gallagher, A Correlation Study of being of Matrix Structural Analysis. Pergamon, Oxford (1964). 6. L. R. Hemnann, A bending analysis for plates. Proc. 1st Conf. on Matrix Meth. in Struct. Mech. AFFDL-TR66-80, Wright-Patterson APB, pp. 577-604 (1965). 7. J. T. Oden, A general theory of finite elements. I, topological considerations, II, applications. Int. J. Numer. Meth. Engng 1, 205-221, 247-259 (1969). 8. B. A. Szabo and G. C. Lee, Derivation of stiffness matrices for problems in elasticity by Galerkin’s method. ht. J.-Numer. Meth. Engng 1, 301-310 (1969). 9. 0. C. Zienkiewicz and Y. K. Cheune. Finite eiements in the solution of iield problems. se Engineer 220, 507-S 10 (1965). 10, E. L. Wilson and R. E. Nickell, Appli~tion of the finite element method to heat conduction analysis. Ntu$. Engng Des. 4, 276-286 (1966). 11 0. C. Zienkiewicz and C. J. Parekh, Transient field problems: two dimensional and three dimensional analysis by isoparametric finite elements. Znr. .I. Numpr. Meth. Engng 2, 61-71 (1970).

z-1150 P d F2 3-l

160

-1170 loo

700 800 300 400 500 600 200 NUMBER OF DEGREES OF FREEDOM ( N )

Fig. 26. Functional f versus the total number of degrees of freedoms N (example 4).

models). As demonstrated in the numerical examples the formulation is equally effective for thin as well as thick shells. For extremely thick shell, or even three-dimensional solid, applications the p-level can be conveniently increased without changing mesh to obtain the desired accuracy of results. It is worth pointing out that even though in describing the element geometry node point normals have been used (for ~nvenien~), the element temperature approximation is not restricted to these normal definitions. Thus, the element formulation presented here permits non-differentiable geometric connections between the elements (i.e. sharp corners).

t STEPPED TUBE 1 D 8-node curved shells

IO’

10* W

N)

Fig. 27. Log of percentage error in i versus log of N {example 4).

K. S. SURANA and G AEWSALEH

880

12. J. C. Bruch and G. Zyvoloski, Transient two dimensional heat conduction problems solved by the finite element method. Int. J. Numer. Meth. Ennnn 8.481494 II

I

(1974).

13. J. Padovan, Semi-analytical finite element procedures for conduction in anisotropic axisymmetric solids. In?. J. Numer. Meth. Engng 8, 295-310 (1974). 14. 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. Inr. J. Numer. Meth. Engng 8, 613-624 (1974). 15. J. Banasek, A conservative finite element method for heat conduction problems. Inr. J. Numer. Meth. Engng 20, 2033-2050 (1984). 16. J. Goldak, A. Chakravarti and M. Biddy, A new finite element mode1 for welding heat source. Metall. Trans. E lSB, 299-305 (1984). 17. M. Kleiber and A. Sluzalec Jr, Finite element analysis of heat flow in friction welding. Rozprawy Inzynierskie, Polish Academy of Sciences, Inst. of Fundamental Technological Research, Warsaw, Poland, 32, 107-l 13 (1984). 18. S. J. Fan and C. I. Tsai, Finite element analysis of welding thermal behavior in transient conditions. American Society of Mechanical Engineers Technical Papers, 84-HT-80 (1984). 19. B. A. Szabo and A. K. Mehta, P-convergent finite element approximations in fracture mechanics. Inr. J. Numer. Meth. Engng 12, 551-560 (1978). 20. I. Katz, Development and application of the p-version

of the finite element method. 21 Nov. 1985, Annual Report 30 Sept., 1984-29 Sept. (1985). 21. I. N. Katz and D. W. Wang, p-version of the finite element method. SIAM J. Numer. Anal. 22, 1082-l 106 (1985). 22. B. Anderson

and U. Falk, Finite element analysis of three-dimensional structures using adaptive p-extensions. Presented at NATO Advanced Study Inst. on Computer Aided Optimal Design, Structural and

Mechanical Systems, Toria, Portugal, 29 June-1 1 July, (1986). 23. J. K. Bennighof and L. Meirovitch, Eigenvalue convergence in finite element method. Int. J. Numer. Meth. kngng 23, 2151-2165 (1986). 24. I. Babuska. 0. C. Zienkiewicz. J. Gaeo and E. R. de A. Oliveira, Accuracy Estimates and AdGptive Refinements in Finite Element Computations. John Wiley, New York

(1986). 25. K. S. Surana and P. Kalim, Finite element formulation for axisymmetric she11heat conduction with temperature gradients. Presented at the Winter Annual Meeting of ASME, Miami Beach, Florida, 17-21 November (1985). 26. P. Kalim, Isoparametric line and transition finite elements with temperature and temperature gradients for axisymmetric and planer two dimensional heat conduction. Ph.D. thesis, University of Kansas, KS (1985). 27. G. J. Desalvo and J. A. Swanson, ANSYS (Engineering Analysis System) User’s Manual, Swanson Analysis Systems, Houston, PA (1983). Corporation, NASTRAN 28. The Macneal-Schwendler User’s Manual (1974). 29. 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, Oak Park, Michigan (1976). 30. K. S. Surana and R. K. Phillips, Three dimensional curved she11finite elements for heat conduction. Comput. Struct. 25, 775-785 (1987). 31. R. K. Phillips, Three-dimensional

curved she11 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, KS (1986). 32. K. S. Surana, FINESSE (Finite Element System) User’s Manual. University of Kansas, Lawrence, KS (1987).