A preliminary evaluation of the B2000 nonlinear shell element Q8N.SM

(1)

Series 05

Aerospace Structures and

Computional Mechanics 03

A Preliminary Evaluation of the

82000 Nonlinear Shell Element

Q8N.SM

(2)

(3)

A Preliminary Evaluation of the

('

82000

Nonlinear Shell Element Q8N.SM

(4)

Series 05: Aerospace Structures

(5)

A Preliminary Evaluation of the

82000 Nonlinear Shell Element

Q8N.SM

c.

Wohlever

Delft University Press / 1 998

2392

317

(6)

Published and distributed by: Delft University Press Mekelweg 4

2628 CD Delft

The Netherlands Telephone

+31 (0)152783254

Fax

+31 (0)15278 1661

e-mail: DUP@DUP.TUDelft.NL by order of:

Faculty of Aerospace Engineering Delft University of Technology Kluyverweg 1

P.O.

Box

5058

2600 GB

Delft The Netherlands Telephone

+31 (0)152781455

Fax

+31 (0)15278 1822

e-mail: Secretariaat@LR.TUDelft.NL website: http://www.lr.tudelft.nl

Cover: Aerospace Design Studio,

66.5 x

45.5 cm, by:

Fer Hakkaart, Dullenbakkersteeg 3, 2312 HP Leiden, The Netherlands Tel.

+ 31 (0)71 512 67 25

90-407-1585-8

1 998 by Faculty of Aerospace Engineering

No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or

mechanical, including photocopying, recording or by any information storage and retrieval system, without written permission from the publisher: Delft University Press.

(7)

1

2 3 4

5

Introduction

Kinematic description of the shell model "Standard" shell model

Numerical implementation of Q8N.SM into B2000 Numerical examples

3

6 7 8

5.1 Snap-through of a hinged cylindrical panel 8

5.2 Pinched hemisphere 9

5.3 In-plane bending of a beam --- mesh distortion 15 5.4 Buckling and post-buckling behavior of a flat plate 15 5.5 Buckling and post-buckling of an axially-compressed cylindrical shell 20

5.6 Linear computations 31

5.6.1 Beam bending; linear analysis 31

5.6.2 Pinched hemisphere; linear analysis 6 ConcIusions

32

34 35

36

References Acknowledgements

(8)

(9)

1 Introduction

This work is a preliminary evaluation of a new B2000 [1] nonlinear shell element designated

"Q8N .SM". The underlying theory used in the development of Q8N .SM is based upon a geomet-rically exact stress resultant shell model and the reader is referred to [2], [3], [4], [5] for an in-depth

discussion. The element Q8N.SM is an implementation of the "standard" nonlinear 4-node shell element described in [5]. The theoretical background for Q8N.SM presented in this report is lim-ited to a brief description of the kinematic assumptions used in developing the shell model. An unusual feature of Q8N.SM is that it was programmed into B2000 as an 8-node element having three dis placement degrees of freedom per node. The reason for the 8-node implementation is

discussed in section 4.

2 Kinematic description of the shell model

Physically, a shell is a three-dimensional deformable body whose motion is described by the balance laws (equilibrium equations) of continuum mechanics. Shell theory is an at tempt to describe the equilibrium state of a "thin" three-dimensional body by a set of two-dimensional balance laws defined on a "mid-surface" of the body. Thus, shell theory is approximate and can only be expected to give partial information about the stress and strain state of a three-dim en sion al body. However, for many engineering applications the partial information that shell

theory provides is sufficient and furthermore, the two-dimensional shell equations are of ten more tractable than the full three-dimensional equilibrium equations.

The three-dimensional body, which is referred to as shell, (Fig. 1), is parameterized by a general convected curvilinear coordinate system (~o,Ç). Throughout this work Greek indices run from 1 to 2 while Roman indices run from 1 to 3. (~,,) are the surface coordinates which parameterize the mid-surface S of the shell. S is assumed to be a smooth orientable surface with a piecewise smooth boundary {) S. The truckness coordinate

~

E [-

~, ~

1

is a measure of the distance between material points in the shell and the mid-surface. h is the constant thickness of the shell in the reference configuration. A standard assumption in shell theory is that thickness h is small compared to the ot her two dimensions of the shell and the mean curvature of S. Hence the convention of referring to a shell as a "thin" three-dimensional body.

A position vector which locates a material point in the undeformed reference configuration of the shell has the form

(1) X(~o) :

JR2

~

JE"

locates the position of the mid-surface Sof the shell and D(~,,) :

JR2

~

JE"

is a unit vector which defines the direction of a materialline (fiber) in the shell. Associated with R is a natural set of covariant basis 'lectors (Go, G3 ) defined as

G" =

~~

= X,o

+~D

.

"

G3 =

~~

=

D (2)

where 0,,, ==

at·

Another set of basis vectors are the covariant vectors (A", A3) associated with the mid-surface S. (Ao , A3 ) are defined by

A_" _aR _- _{a~Q ~=o}

I

=X,o

(10)

Coordinate Lines---,..:;c

x,

Figure 1: Kinematic shell model.

Reference Configuration

u

Defonned Configuration

X,

(11)

When the shell deforms, (Fig. 2), a material point originally located by a position vector R in the reference state moves to a new point located by the position vector r defined as

(4)

The kinematic assumption built into eq.

(4)

is that an initially straight material line along the

director D in the undeformed configuration is deformed into another straight materialline along

the direct or d in the deformed configuration. Equation (4) describes the kinematics of what is known as an extensible one-director shell model. The displacement U(ço,Ç) of a material point

is given by

U(ço,ç) = r(ç""Ç)-R(ço,Ç)

[x(çol - X(ço)] +([d(çol - D(çol]

' - . . . - ' ' - . . . - '

u v

(5)

where u(ça) is the dis placement of the mid-surface, v(ço) is the director displacement.1

From eq. (4

1

two sets of natural covariant convected basis vectors (go, g3) and (a"" a3) defined with respect to the deformed shell are

go ₌

:~:

= x,o

+

çd,o

g3 ₌

~~

=d (6)

and

a", _-

_

~

a

r

I

~=o = x,o

a3

_

_-

ar

a~ ~=o

I

=d (ï)

The sets (Go, _{G 3) anc}(go, g3) are referred to as the body convected bases while the sets (A"" A 3) and _{(ao ,a3)}are referred to as the surface convected bases. Figure 3 illustrates the surface

con-vected bases in the deformed and undeformed configurations. One may also write the deformed surface convected basis as a function of the displacement fields u and v as

a o

=

X,o

+

u,o

d

=

D+v (8)

It is assumed that for physically meaningful deformations, the following condition must always hold,

(9)

Condition (9) is analogous to the requirement in three-dimensional elasticity [8] that the deter-minant of the deformation gradient is always greater than zero.

1 Note that unlike classica! linear shell theory [6] and some nonlinear shell theories [7], the deformed director d

is not constrained to be of unit length, i.e., the shell is able to change its thickness during the deformation. The director disalso free to rotate with respect to the shell mid-surface 5, thus allowing for shear.

(12)

Mid-Surface in Reference Configuration D

)"'-_ _ _ x.

X, x Mid-Surface in Deformed Configuration

Figure 3: Surface convected bases in the deformed and undeformed configuration.

3 "Standard"

shell model

In most standard FE shell models there are five "true" degrees of freedom per node - three

dis placement variables for the mid-surface and two rotations for the surface normal. The kinematic description of the shell discussed in section 2 affords a rather unique implementation for a nonlinear shell element. By allowing the direct or field d to be extensible, one can assign six "true" degrees of freedom per node - three displacement variables for the mid-surface and three displacement variables for the director d. Referring to eq. (5), the displacement variables u and y can be written as

(10) The set (El, E2 , E3 ) is the standard orthonormal basis in Euclidean three space. The advantage

of such an approach is that the configuration updates, which can occur in a nonlinear analysis, can take place in a linear space. In other words, iterative updates of a displacement field which occur in a Newton-Rhapson process are simple additions and would be of the form,

uk+! uk

+

llu k+!

yk+l

=

yk

+

llyk+! (11)

where uk and yk are the displacement fields from step k, and lluk+! and llyk+l are the current

corrective updates. Having the configuration updates take place in a linear space avoids having to keep track of things such as updated reference frames or having to maintain inextensiblilty-constraints usually found with ot her nonlinear shell theories. In fact, the implementation of this type of element is very similar to that of a full 3-D continuum element. The disadvantage, as discussed in [5), is that this element tends to become ill-conditioned in the thin sheillimit. So at this point in time, the shell element Q8N .SM is restricted to the rat her vague class of

(13)

should be aware of this short-coming and use the proper caution - of course, caution should

always be used in nonlinear FEM analysis.

To avoid common problems encountered with shell elements such as "membrane-locking" and "shear-locking" the authors in [5] used ~ mixed-variational formulation. The mixed-variational approach avoids the need to invoke procedures such as uniform or reduced integration. Thus, full 2 x 2 Gauss integration is used and the element does not suffer from the problem of spurious modes. It should be noted however that the mixed-variational formulation in [5] assumes the existence of an internal potential energy function in which the membrane, bending and shear fields are all decoupled from one another. Thus, the extension of this element to include plastic

analysis is questionable.

4 Numerical implementation of Q8N.SM into B2000

Due to the simp Ie nature of the nonlinear updates, the implementation of the element Q8N.SM was relatively straight forward. The no dal pattern for the element was chosen to match that of the B2000 8-node continuum element HE8.BA. As depicted in Fig. 4, nodes 1-4 are used to locate the mid-surface of the shell and no des 5-8 are used to locate the outer-surface of the shell. Although this way of defining a shell element is unusual it is a convenient way to define both the thickness of the shell and the initial orientations of the directors at the nodes. It is important to keep in mind however, that although the implementation of this element in B2000 looks three-dimensional the

underlying field equations and constitutive laws are strictly two-dimensional.

X₃

i

; i I I

/

XI 8 Top Surface 7 3 ~ •... --~~---~~·6 5 Mid·Surface / " ' i r -h 2 I h 2 _{Bottom Surface} • - Element Node X₂

(14)

Some general comments about the shell element as implemented in B2000 are included here: • At this particular time, the implementation of Q8N.SM is limited to an isotropie linear elastic constitutive law. Thus, the only material parameters defined for the shell are Young's modulus (E) and Poisson's ratio (/I). The capability for a non-isotropic material or a layered material could be included.

• For proper behavior, the shell must use the 2 X 2 Gauss integration discussed above, I.e.,

nint

=

2.

• This shell element allows for the definition of a surface load, as opposed to a more standard load applied to the mid-surface. However, this ability has not been explored and the reader is referred to (5) for more detail on this subject.

• Only a nonlinear version of this particular shell element has been programmed into B2000. • A gradient package to compute internal stresses is not available in B2000 at this time.

• The element Q8N.SM has a prevariational package oflength 65. This is substantially smaller than the prevariational packages for any of the other nonlinear shell elements in B2000. Given the simple update structure of Q8N .SM, it is possible to elintinate the prevariational package altogether - although this would result in a slight increase in computational over-head as some initial computations for quantities such as the element area would have to be repeated at every load step.

• The required element attributes for Q8N.SM are as follows:

mid m specifies element material number, which must identify both Young's modulus (E) and Poisson's ratio (/I);

type Q8N .SM specifies current element type.

5 Numerical examples

In this section, several numerical simulations are presented which demonstrate the abilities and shortcomings ofthe element Q8N.SM. The numeri cal .oxamples include some "standard" nonlinear benchmark problems which test the element's ability to capture large deformation response of thin shell structures with relatively coarse meshes. There are also two buckling and post-buckling studies to demonstrate the element 's performance in more challenging large scale analyses. Finally, two linear benchmark problems are used to discuss the element's sensitivity to mesh distortion and convergence properties.

5.1 Snap-through of a hinged cylindrical panel

The snap-through of a shallow hinged cylindrical panel loaded by a single point load, as shown in Fig. 5, is a popular nonlinear benchmark. As the concentrated load is slowly increased the initia! shape of the panel, characterized by its negative curvature, loses stability and the panel "snaps through" to a new equilibrium shape with positive curvature. Referring to Fig.5, the geometrie and material parameters for this problem are: R=2540 mm, L=254 mm, Young's modulus E= 3120. i5 N /mm2, Poisson's ratio /I = 0.3, 0 = 0.1 radians. Simulations were run for two values of shell thickness(t) giving radius-to-thickness ratios of R/t=200 and R/t=400. In

(15)

F

Figure 5: Hinged cylindrical panel.

addition, each simulation was performed using two different meshes, a 4 x 4 and an 8 x 8 mesh both modeling one-quarter of the panel.

The computations for the case Rjt=200 were compared to results reported by Simo et al. in [4] and are displayed in Fig. 6. The results in [4] were computed using a 4 x 4 mesh of a four-node element very similar to Q8N.sM. The main difference bet ween the two elements is that the element used in [4] is based on an inextensible director, and thus has five DOF per node. The computed results from Q8N.SM compare weil with Simo's results - the only noticeable point of discrepancy is in the vicinity of the second limit point.

The computations for the case Rjt=400, were compared to results reported by Sansour and Bednarczyk [9]. The computations in [9] were performed using a 4 x 4 mesh of rune-node elements. The higher radius-to-thickness ratio is a more challenging problem but the 8 x 8 Q8N.SM mesh compares extremely well to results reported in [9] as is illustrated in Fig. ï.

5.2 Pinched hemisphere

The pinched hemisphere problem was originally developed as a benchmark for linear sheil elements [10]. The finite deformation version of the pinched hemisphere has become a popular benchmark for nonlinear shell elements. The problem consists of a hemisphere, with an 18°-hole at the

(16)

-_. 5 NodesiSide (16 Q8N.SM Elements) - 9 NodesiSide (64 Q8N.SM Elements)

D 5 NodesiSide (16 Four·Node Elements), Simo et al.

2.5 2 Load (KN) 1.5

Figure 6

10 15 Deflection (mm) 25

Snap-through of a hinged cylindrical panel: Rlt=200

.

Plot shows the

mag-nitude of the point load versus the shell deflection at the point of

applicati-on.

(17)

0.8 -_. 5 Nodes/Side (16 Q8N.SM Elements) - 9 Nodes/Side (64 Q8N.SM Elements)

a 9 Nodes/Side (16 Nine·Node Elements), Sansour & Bednarczyk

0.6 0.4 Load (KN) 0.2 ·0.2 -0.4

Figure 7

o

1'0

, , ,

·

I I I I I I , , , ,

·

/0

.

, , , 15 Deflection (mm)

20

25 30

Snap-through of a hinged cylindrical panel: RJt=400. Plot shows the

mag-nitude of the point load versus the shell deflection at the point of

applicati-on.

(18)

Figure 8: Finite deformation version of a pinched hemisphere: Problem configuration. top, which is deformed by two inward and two outward point loads as shown in Fig.8. The

geometrie and material parameters for this problem are: sphere radius R=lC ;urn, shell thickness

t=0.04 mm, E=6.825 x 10ï_N_/mm2 _and_{v = 0.}_{3. Symmetry boundary conditions were used and}

only the first octant of the sp here was modeled. The computed results for a 16 x 16 and a 32 x 32

mesh are shown in Fig. 9 and are compared to results reported in [4] which were computed with

a 16 x 16 mesh. Figure 10 displays the final deformed configuration of the hemisphere without magnification.

The 32 x 32 mesh shows excellent agreement with the results of Simo et al. [4], while the

16 x 16 mesh gives slightly stiffer results. The stiffer response of the 16 x 16 mesh is consistent with observations discussed in [5] and could indicate an ill-conditioning in the thin shell limit. However, a more refined mesh in this particular case seemed to mitigate the stiff response.

(19)

120r---~

100

80

_. 17 Nodes/Side (256 Q8N.SM E1ements) - 33 Nodes/Side (1024 Q8N.SM Elements)

a 17 Nodes/Side (256 FouT-Node Elements), Simo et al.

, , , , , , , , , , , , Load (N) 60 Figure 9 40 20 3 6 Deflection (mm)

Finite defonnation of a pinched hemisphere Rlt=250. Deflection-Ioad curves for different mesh configurations. The curves on the right plot the magnitude of the inward point loads versus the inward deflection of the shell at the point of application (y-direction). The curves on the left show the magnitude of the outward point loads versus the outward deflection of the shell at the point of application (x-direction).

(20)

(21)

j

I~

F regular elements

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

j

\S

/

ZZ

I

\

i

~

F skewed elements

Figure 11: Regular and skewed elements for nonlinear beam bending.

5.3 In-plane bending of a beam - mesh di stort ion

The nonlinear in-plane bending of a beam can be used as a measure of an element's sensitivity to mesh distortion. The goal is to compute the large deformation response of a cantilevered beam loaded by an in-plane shear force. To compare the relative effect of mesh distortion, the computa-tions were performed with two meshes, (Fig. 11): alO-element mesh of regular shaped elements, and alO-element mesh of irregular shaped elements. The material and geometric properties used for this problem are as follows: E= 10 x 106_,-_V _{= 0_3, beam length L= 1.0, beam thickness and}

width h=w=O_l. The two meShes-áii1:r-flfematerial and geometrie proper ties described above we re the same ~sed in [4]. In Fig. 12 the tip displacement versus load for both meshes is shown. The results In Fig. 12 clearly show that for the material and geometric properties chosen for this analysis, the element shows almost no sensitivity to mesh distortion and compares very well with the results reported in [4]. However, as noted in [9], the apparent sensitivity of a shell element to mesh distortion can be dramatically affected by the relative magnitudes of the element's width, thickness and length. The question of mesh sensitivity of Q8N.SM was further investigated in a linear analysis, the results of which are discussed in section 5.6.1.

5.4 Buckling and post-buckling behavior of a flat plate

In this section, the results of a buckling and post-buckling analysis of a flat plate will be presented. These results will demonstrate the performance of the shell element Q8N .SM in a large buckling analysis.

The static nonlinear elastic deformation of a rectangular flat plate with classical boundary conditions and subjected to in-plane compressive forces is a much studied and relatively weil understood problem - at least for "moderate" deformations. Consider the situation in Fig. 13 which illustrates a flat plate made of a homogeneous elastic material of uniform thickness and deformed by a uniform end-displacement. The plate has simple-support boundary conditions along the unloaded edges and is clamped along the loaded edges. As is weil known from classical plate

(22)

"-1J'~ .. t • ";.- ~.::. . .

I~.---~

Regular Mesh (10 Q8N.SM Elements) Skewed Mesh (10 Q8N.SM Elements)

800 D Regular Mesh (10 Four-Node Elements), Simo el al.

600 Load 400 200 0.2 0.4 Tip Deflection

(23)

(y,v) (z,w) cJamped "---- simply supponed boundary conditions boundary conditions

Figure 13: Axially-Ioaded flat plate.

(X,u)

cJamped

boundary conditions

theory [11], the initial equilibrium shape of the plate, starting from zero load, maintains its original

flat shape, i.e., the symmetry of the initial deformed shape is identical to the symmetry of the undeformed plate. As the load is further increased, a stabie symmetrie pitchfork bifurcation point is reached and the fiat shape loses stability and the plate buckles into a new stabie configuration which is characterized by one wave in the transverse direction and one or more waves in the longitudinal direction. The new stabie equilibrium shape initially supports an increasing load, i.e., there is no immediate catastrophic failure. Referring to Fig. 13 and with respect to the

standard XYZ-coordinate system, the variables (u,v,w) represent the displacement of the mid-surface of the plate and the variables (0:,,:3,7) represent the dis placement of the directors. The boundary conditions along the various edges can now be expressed as:

l. Loaded edges w = 0: =

!3

= 0;

2. Simply supported edge w = v = 0: = 0;

3. Symmetry boundary conditions v = ,:3 = O.

The initial analysis performed on the plate consisted of computing the first 4 buckling modes oft" the fiat plate state. The buckling analysis in this study was carried out by computing the

nonlinear primary solution branch with a path-following routine and solving the linear eigenvalue problem in the vicinity of the bifurcation points - for the element Q8N .SM this was carried out using the B2000 group-theoretic nonlinear continuation processor "b2contsym". These results are compared to the results from a similar analysis done with STAGS2 _{using the 4-node nonlinear}

shell element E410 [12].

Referring to Fig. 13, the geometrie and material parameters for the plate are: L=1076.86 mmo

t=2 mm, h=200 mm, E=70000 N/mm2 and v = 0.3. The uniform finite element mesh used for

(24)

the buckling analysis consisted of 70 elements along the length of the plate and 8 elements along the width. The shapes of the first four buckling modes are shown in Fig. 14 and their respective normalized buckling loads (>.) are displayed in Table 1.3 _{The results}_{of the separate buckling}

analysis performed using the elements Q8N.SM and E410 are in very good agreement.

Table 1: Loads for the first .:1 buckling modes of an axially-compressed fiat plate (70 x 8 mesh).

Normalized Buckling Normalized Buckling

Mode Number Load

p.)

Q8N.SM Load

p.

)

E410 % difference

1 0.699 0.696 0.43

2 0.706 0.702 0.57

3 0.ï42 0.739 0.41

4 0.768 0.764 0.52

In the second analysis, the first two primary bifurcation branches for an imperfection-free plate

were computed. The eigenrnodes from the initial buckling analysis were used to define initial step

directions along the bifurcation branches. To investigate the con vergen ce properties of Q8N .SM

the first primary bifurcation branch, (Fig. 15), was computed with two different meshes: the 70 x 8 mesh des cri bed above and a uniform 112 x 10 mesh. For further comparison, the location

of the first six limit points and secondary bifurcation points as computed by the two mes hes are

shown in Table 2. Note that in Fig. 15 and Table 2 the independent parameter is the pres cri bed displacement factor PA instead of the normalized axial-load À. The tot al end-shortening of the

plate can be computed by multiplying PA by 0.8033 mmo As is evident in Fig. 15 and Table 2, the results from the two meshes agree very weil with each other suggesting that the coarser 70 x 8

mesh is adequate for this particular analysis.

Table 2: Displacement factor PA for the first six limit and bifurcation points along the first

primary bifurcation branch as computed by two different meshes.

displacement factor PA displacement factor PA

70 x 8 mesh 112 x 10 mesh % difference

Ist _{limit point} _0.442 0.440 0.45 2nd _{limit point} _0.433 0.431 0.46 3rd _{limit point} _1.20 1.18 1.7 4th _{limit point} _0.925 0.909 1.8 .5th _{limit point} _4.30 ~ .12 4.4 6th _{limit point} 2.15 _.11 1.9 p t bif. point 0.411 0.410 0.24 2nd _bif._point _0.455 0.453 0.44 3rd _{bif. point} _1.15 1.14 0.88 4tn _bif_._point ₀_.965 0.953 1.3 .5tn _{bif. point} _4.15 4.02 3.2 6th _bif_._point _2.22 2.18 1.8

(25)

(26)

In Fig. 16, the stabIe and unstable parts of the first primary bifurcation branch are plotted. The computed results show that the primary bifurcation branch is stabIe up to a load factor of PA ~ 0041. At that point a secondary bifurcation point and a limit point occur almost on top of each other and the branch becomes unstable. However, further along the path another secondary bifurcation point and limit point occur, again in very close proximity to one another, and the primary bifurcation branch regains its stability. This pattern of change-in-stability along the primary branch is repeated a number of times and is clearly illustrated in Fig. 16.

The second primary bifurcation branch was computed with the ,0 x 8 mesh and is shown in Fig. 1 ï. Coming off the primary solution branch the second bifurcation branch is initially unstable but at a normalized load of PA ~ 0.35, a bifurcation point is detected and the branch becomes stabIe. It remains stabIe until PA ~ 0.64 at which point a limit and bifurcation point occur very close to one another and the branch loses stability. However, as was found with the first primary bifurcation branch, another limit and bifurcation point occur at PA ~ 0.612 and PA ~ 0.642 respectively and the branch regains stability.

It is important to keep in mind that all the analyses discussed in this section were done using symmetry boundary conditions along the center of the plate as illustrated in Fig. 13. Therefore, statements about the absolute stability of the above solutions can only refer to stability with respect to small perturbations which are symmetrie about the center of the plate. To determine the absolute stability of these solutions with respect to arbitrary perturbations, the full plate must be modeled. The full plate was modeled in a follow-up analysis with a 70 x 17 mesh and it was found that that the stabie branches shown in Figs. 16 and 17 are indeed stabIe with respect to an)' small perturbation.

5.5 Buckling and post-buckling of an axially-compressed cylindrical shell

In this section, the results of a buckling and post-buckling analysis of an axially-compressed cylindrical shell will be presented. These results will further demonstrate the performance of Q8l\ .SM in a rather large and sophisticated buckling analysis.

Consider the situation in Fig. 18 which illustrates a cylindrical shell made of a homogeneous elastic material of constant thickness and deformed by a uniformly applied compressive load. A buckling and post-buckling analysis of an axially-compressed cylindrical shell is one of the most challenging problems in nonlinear bifurcation analysis. It is weil known that the primary axisym-metric solution branch can be riddled with cJosely spaeed symmetry-breaking bifurcation points.

In a numerical arc-length continuation scheme the close proxirnity of the bifurcation points on the primary path manifests itself in severe ill-conditioning of the tangent stiffness matrix. To cope with the numeri cal ill-conditioning which this problem poses, the B2000 group-theoretic nonlinear continuation processor "b2contsym" was used to carry out the analysis with the Q8N.s:M shell element. The cylinder considered in this analysis was modeled with simple-support (SS-3) bound-ary conditions on both ends and the following geometrie and material parameters: L=200 in .. R=100 in .. t=l in., E=106 psi and v = 0.3. The full cylinder was modeled with a uniform FE mesh which consisted of 50 Q8N .SM elements along the length and 144 elements around the circumference.

As was done in section 5.4, two analyses were performed on the cylindrical shell. The ini-tial analysis consisted of detecting the first six buckling modes along the primary axisymmetric solution branch. The buckling analysis using Q8N.SM was carried out by solving the linear eigen-value problem in the vicinity of the bifurcation points along the primary solution branch. The buckling loads and mode shapes from the Q8N.SM analysis were compared to computed results

(27)

PA Figure 15 8 70x 8 mesh ll3 x 10 mesh 4 2

o

_ó

₃ ₆

Nonnal Displacement of Plate Center

NonnaJ displacement of center of the plate versus the displacement load factor (PA). First primary bifurcation branch for the tlat-plate computed with 2 different meshes.

(28)

PA Figure 16 4 - - - slab Ie branch --- unslable branch 3 2

o

Ó

10

20 30 50

60

70 IDisplacemenl1

EucJidean norm of all plate displacements versus the displacement load factor (PA)_ First primary bifurcation branch for the flat-plate: 70

* 8

mesh_

(29)

PA Figure 17 4 - - - stabJe branch ---. unstable branch 2 ---

~

o

Ó

10

20

30

40

50

70

IDispJacementl

EucIidean nonn of all plate displacements versus the displacement load factor (PA). Second primary bifurcation branch for the fIat-plate: 70

*

8 mesh.

(30)

/axialioad

~..--,...--.---,..-R

L

,

, 1-:r-j---'--t---l-1- '

I ... - - ... I (' .!..

-

_-

'

_-

-Figure 18: Axially compressed cylindrical shell.

based upon a semi-analytic discretization, i.e., Fourier decomposition, of Donnell-type shell equa-tions, referred to as a "Level 2" analysis.4 _{The Fourier decomposition leads to a set of ordinary} differential equations which can be solved by the "Shooting Method" [13].

A classical analysis for the cylindrical shell described above would predict the first buckling mode should occur at a normalized axial-Ioad of approximately 0.85 and have 8 or 9 circumferential waves [14). The Q8N.SM analysis found the first buckling mode to have 8 circumferential waves at a normalized load of 0.843. These results were supported by the Level 2 analysis which predicted the first buckling mode to have 8 circumferential waves and to occur at a normalized load of 0.842. The radial displacement and cylinder cross sections for the first six buckling modes, as computed with Q8N.SM, are shown in Figs.19 and 20. The buckling loads and wave numbers for the first six buckling modes, as computed with both Q8N.SM and the Level 2 analysis, are displayed in Table 3. In Table 3, the column "Type" identifies the symmetry of the buckling mode as being either symmetric (Sym) or anti-symmetric (Asym) with respect to the mid-plane of the cylinder. The results from the buckling analysis done with the element Q8N.SM compare extremely weil with the results from the Level 2 analysis. The only notabie discrepancy occurs for the anti-symmetric buckling mode with ï circumferential waves. The reason for this discrepancy is unclear at this point. Perhaps a more refined mesh in the axial direction with the Q8N.SM element is necessary for this particular buckling mode.

In the second analysis, the primary bifurcation branch associated with the first buckling mode and one secondary bifurcation branch were computed with the Q8N.SM element. As shown in Fig. 19 the eigenmode associated with the first bifurcation point is symmetrie ab out the mid-plane of the cylinder and has 8 circumferential waves. It can be shown using results from group theory [15) that the shape of the cylinder along the first primary bifurcation branch will always maintain an 8-fold symmetry about the circumference and a reflection symmetry about the mid-plane of

(31)

Figure 19

Cross sections and radial displacements of the 1 st three buckling modes for

the cylindrical shell.

(32)

Figure 20 Cross sections and radial displacements of the 2nd three buckling modes for the cylindrical shell.

(33)

Table 3: Loads for the first 6 buckling modes of an axially-compressed cylindrical shell.

Normalized Normalized

Mode Wave Buckling Load Buckling Load

Number Number Type Q8N.SM Level 2

%

difference

1 8 Sym 0.843 0.842 0.12 2 8 Asym 0.849 0.846 0.35 3 9 Sym 0.856 0.851 0.59 4 9 Asym 0.857 0.852 0.59 5 7 Sym 0.861 0.859 0.23 6 7 Asym 0.862 0.84.5 2.01

the cylinder. Therefore, the first primary bifurcation branch was computed by modeling half the cylinder and using symmetry boundary conditions. The fini te element mesh consisted of 20 elements along the length of the cylinder and 144 elements around the circumference. While the bifurcation branch was being computed, a limited secondary bifurcation analysis was performed. The secondary bifurcation analysis was limited in the sense that only bifurcation points associated with eigenmodes which had ol-fold symmetry and were symmetrie about the mid-plane of the

cylinder were diagnosed. Such a secondary bifurcation point was found at a normalized load of

0.443. The secondary bifurcation branch passing through this point was then computed. The

results from post-buckling analysis are shown in Fig. 21. The load and displacement in Fig. 21 are respectively the normalized axial-Ioad and normalized end-displacement. The primary solution

branch in Fig.21 is stabIe until the first bifurcation point is reach at a normalized load of 0.84.

Both the primary bifurcation branch and the secondary bifurcation branch are unstable. It is interesting to note that the secondary bifurcation branch has a limit point which occurs at a normalized axial-load of approximately 0.11. The cross section of the mid-plane of the cylinder and the shape of the cylinder ne ar this limit point are shown in Figs.22 and 23 respectively.

(34)

Primary Solution Branch Primary Bifurcation Branch 0.8 Secondary Bifurcation Branch

0.6 Load Figure 21 0.4 0.2 0:2 0.4 , , ... Displacement ~~~ ~~~ ~~~

----

---0.8

---Post-buckling response of the cylindrical shell. Normalized axial-load versus norrnalized end-displacement for the axially compressed cylindrical shell.

(35)

Figure 22 lOO 50 -50 -lOO -foo -50 50 lOO

Cross section of the axially compressed cylinder along the secondary bifur-cation branch, as displayed in Fig.21, near the limit point which occurs at a normallzed axial-load of 0.11.

(36)

Figure 23 Shape of the axially compressed cylinder along the secondary bifurcation branch, as displayed in Fig.21, near the limit point which occurs at a nor-malized axial-Ioad of 0.11.

(37)

Figure 24

I

IL.-/ _ _ / _ _ _ _ _

~

- extenslOn

V'~'-"""-L

/

.

/ L - - - ) . " in-plane shear

~

Z

7

Z

7 ~

I

_I

I

regular elemeDls

~

Û

rK\\\\

paraJlelogrnm elcmenlS

Cantilevered beam loaded on the free end_ Rectangular and parallelograrn shaped elements_

5.6 Linear computations

Although the element Q8N_SM is a nonlinear element, it is instructive to consider the element's

performance in two standard linear benchmark problems taken from [10]- The linear beam bend-ing problem discussed in section 5_6_1 further investigates the sensitivity of Q8N_SM to mesh distortion, while the linear pinched hemisphere problem presented in section 5_6_2 illustrates some convergence characteristics of the element- The "linear" results presented in this section are

the initia! Euler step predictions computed by the B2000 nonlinear continuation module "b2cont"

[lJ-5.6.1 Beam bending; linear analysis

The linear deformation of straight cantilevered beam loaded at the free end, (Fig_ 24), is a common

linear benchmark for plate and shell elements [10]- A5 noted in section 5_3, by performing the

computations with irregular shaped elements thls test can be used as one measure of an element '5

sensitivity to mesh distortion_ The geometric and materia! parameters for beam ana!yzed here are: beam length L=6~· beam thickness h=O_l; E=107 and v = 0_3_ The dellectlvns of the loaded beam were computed for four different values of width: w=l, 0_5, 0_2, 0_1- The different values of width were chosen to demonstrate how the aspect ratio

ç

of the element length-to-width can

affect its performance. A5 recommended in [10], all the computations were done with a 1 x 6 mesh and the mesh was chosen 50 that all the elements had equa! volumes_ The numerica! results in

Table 4 have been normalized by theoretica! values from classica! linear beam theory_

The normalized values in Table 4 demonstrate the excellent performance of the element

Q8N_SM when used in a rectangular mesh_ Indeed, all the numerically computed results are within three-percent of theoretica! values and are relatively unaffected by the aspect ratio

ç

(38)

Table 4: Normalized defiections for the straight cantilevered beam. rectangular elements parallelogram elements

: out-of-plane in-plane axial out-of-plane in-plane axial

( shear shear extension shear shear extension

1 .981 1.01 .996 1.05 1.02 1.07

2 .978 .996 .996 1.05 .892 1.12

5 .981 .991 .996 1.06 .676 1.32

10 .980 .973 .996 1.06 .590 7.78

However, the same does not hold true for the parallelogram-shaped elements. For an element as-pect ratio of (

=

1, the parallelogram shaped elements still perform weil - but this performance degrades as the aspect ratio increases. For example, the computed values for the in-plane shearing and axial extension show more than ten-percent error for an aspect ratio of 2. This sensitivity to mesh distortion is disquieting and should be taken into account while using this element. It is interesting to note that while discussing their sheil element's sensitivity to mesh distortion in the cantilevered beam problem, the authors in both [4] and [9] used an aspect ratio of ( = 1. 5.6.2 Pinched hemisphere; linear analysis

The pinched hemisphere problem was originally developed as a benchmark for linear sheil elements [10]. As discussed in [4], this test is a good measure of an element's ability to capture inextensible bending behavior and to model rigid body modes. The problem parameters are identical to those used in Section 5.2. The results for the displacements of the sheil at the load points, normalized by 0.094 as recommended in [10], for five mesh configurations are listed in Table 5. The results in Table 5 compare weil with those discussed in [5].

Table 5: Results for the pinched hemisphere: Linear analysis. Mesh Configuration Normalized Displacements

4x4 0.098 8x8 0.862 16 x 16 0.988 32 x 32 0.994 .54 x 54 0.995

6 Conclusions

In this work, a first look at the implementation into B2000 of the finite-deformation small-strain 8-node nonlinear sheil element Q8N.SM has been given. The underlying theory and numerical implementation of Q8N .SM is based on a standard extensible one-director sheil model as discussed in [2]-[5]. The performance of Q8N.SM in the B2000 environment has been demonstrated in a number of challenging problems and the basic conclusions and recommendations for this element are summarized here:

(39)

• Overall the performance of Q8N .SM compares weil with ot her state-of-the-art nonlinear sheil elements and can be recommended for use in a wide range of problems, including buckling and post-buckling analyses.

• With regard to the previous statement, the user should keep in mind that the element's ability to produce accurate answers with coarse mes hes is considered to beonly average. • The linear nature of the update space for Q8N.SM leads to a computationaily attractive

framework by circumventing the need for nonlinear rotational updates.

• Although the current implementation uses a smail prevariational package, the simple struc -ture of the update space obviates the need for one, and it could be eliminated.

• As reported in [5], the element exhibits some ill-conditioning in the thin-sheil limit as com -pared to the inextensible director formulation discussed in [4] and the multiplicative director formulation discussed in [5]. However, it appears that at least in some situations, the iil-conditioning can be overcome by refining the mesh - which of course leads to increased computationaloverhead.

• As discussed in section .5.6.1, the sensitivity to mesh distortion for various aspect ratios ( is disturbing and merits further study.

(40)

References

[1] S. M. Merazzi and A. de Boer. B2000 processors reference manual, version 1.8. SMR Cor-poration, 1995.

[2] J. C. Simo and D. D. Fox. On a stress resultant geometrically exact shell model. part

i: Formulation and optimal parametrization. Comp. Meths. App. Mech. Eng., 72:267-304,

1989.

[3] J. C. Simo, D. D. Fox, and M. S. Rifai. On a stress resultant geometricallY exact shell model. part ii: The linear theorYi computational aspects. Comp. Meths. App. Mech. Eng., 73:53-92,

1989.

[4] J. C. Simo, D. D. Fox, and M. S. Rifai. On a stress resultant geometrically exact shell model. part iii: Computational aspects of the nonlinear theory. Comp. Meths. App. Mech. Eng.,

79:21-70,1990.

[5] J. C. Simo, M. S. Rifai, and D. D. Fox. On a stress resultant geometrieal!y exact shell model. part iv: VariabIe thiekness shells with through-the-thiekness stretching. Comp. Meths. App. Mech. Eng., 81:91-126,1990.

[6] A. E. H. Love. On the smal! free vibrations and deformations of thin elastic shells. Phil.

Trans. Royal Soc. (London), l79A, 1888.

[7] J. L. Sanders Jr. Nonlinear theories for thin shells. Quart. App. Math., 21:21-36, 1963. [8] M. E. Gurtin. An Introduction to Continuum Mechanics. Academie Press, 1981.

[9] C. Sansour and H. Bednarczyk. The cosserat surface as a shell model, theory and fini te-element formulation. Comp. Mahs. App. Mech. Eng., 120:1-32, 1995.

[10] R. H. MacNeal and R.L. Harder. A proposed set of problems to test finite element accuracy. In

Proceedings of the AIAA Conference on Structures and Structural Dynamics, Palm Springs, CA, USA, 1984.

[11] D. O. Brush and B. O. Almroth. Buckling of Bars, Plates and Shells. McGraw-Hill Inc., 1975.

[12] F. A. Brogan, C. C. Rankin, and H. D. Cabiness. Stags user manual, version 2.0. Lockeed Palo Alto Research Laboratory, 1994.

[13] J. Arbocz and J. M. A. M. Hol. Koiter's stability theory in a computer-aided engineering

(CAE) environment. Int. J. Sol. Struct., 26(9/10):945-973, 1990.

[14] N. Yamaki. Elastic Stability of Circular Cylindrical Shells. North-Holland, 1984.

[15] T. J. Healey. A group theoretic approach to computational bifurcation problems with sym-metry. Comp. Meths. App. Mech. Eng., 67:257-295, 1988.

(41)

Acknowledgements

The work described in this report has been in part supported by the "Dwarsverband - Commissie Niet-lineaire Systemen" from NWO ("Nederlandse Organisatie voor Wetenschappelijk Onder-zoek") via a 6 month fellowship for Dr J.C. Wohlever's stay in the Netherlands. This aid is gratefully acknowledged.

(42)

A

Sample input file for Q8N .SM

This appendix contains an input file used for the rectangular mesh beam ben ding problem dis-cussed in section 5.6.1. Figure 25 iL .~trates the FEM beam model described in the input file

"beam.rect.ip." Note that odd-numbered nodes locate points on the mid-surface ofthe beam and even-numbered no des locate points on the outer-surface of the beam.

z 13 17 ₂₁ 2 :: 6 : 10 : 14 : 18: 22 : 26 I' I , , I I I I , ,

/

l\::::-l~::.·:::-i!

:

:::::.~~:::

:

::j~::::::

l?::::::.

/ 4' ( . ' \ 8 (;:;\12 (.;"\16 ( ; \20 (;;\24 r?\28 X

000~00/

i. _ node number i

(0 .

element number i y 27

Figure 25: Node and element definition for cantilevered beam: Rectangular elements.

*****************************************************************

title 'beam.reet.ip: (b=O.OS)

(w=0.50)

adir

analysis collapse

Input file for reetangular mesh'

<=

Define b to be half the beam thickness

<=

Define II to be half the beam width

pas 0.0 dpas 1.0 pamax 1

leb 0 neut 1 nfaet 2 nstrat 0 maxit 15 maxstp 1 epsdis 1.0 epsr 0.000001 end branch 1 bdir deform nonlinear end nodes 1 O. 0 (w) 2 (h) 0 (w) 3 O. 0 ( -w) 4 (h) 0 (-w) 5 O. 1 (w) 6 (h) 1 (w) 7

o

.

1 (-w) 8 (h) 1 (-w) 9 O. 2 (w)

(43)

10 (h) 2 (w) 11 O. 2 (-w) 12 (h) 2 (-w) 13 O. 3 (w) 14 (h) 3 (w) 15

o

.

3 (-w) 16 (h) 3 (-w) 17

o

.

4 (w) 18 (h) 4 (w) 19 O. 20 (h) 21

o.

22 (h) 23 O. 24 (h) 25 O. 26 (h) 27

o.

28 (h) end elem 4 (-w) 4 ( -w) 5 (w) 5 (w) 5 (-ll) 5 (-ll) 6 (w) 6 (w) 6 (-ll) 6 (-ll)

type Q8N.SM mid 1 nint 2

1 1 3 7 5 2 4 8 6 2 5 7 11 9 6 8 12 10 3 9 11 15 13 10 12 16 14 4 13 15 19 17 14 16 20 18 5 17 19 23 21 18 20 24 22 6 21 23 27 25 22 24 28 26 end bound loek LLL 1 3 loek FLL 2 4 end force case 1 dof 1/3/2

P

.5 25 27 end endbranch emat mid 1

type ISOTROPIC e 1.e7 p 0.3

endmid end run

(44)

(45)

Series 01: Aerodynamics

01. F. Motallebi, 'Prediction of Mean Flow Data for Adiabatic 2-D Compressible Turbulent Boundary Layers'

1997 / VI

+

90 pages / ISBN 90-407-1564-5

02. P.E. Skäre, 'Flow Measurements for an Afterbody in a Vertical Wind Tunnel'

1997 / XIV

+

98 pages / ISBN 90-407-1565-3

03. B.W. van Oudheusden, 'Investigation of Large-Amplitude 1-DOF Rotational Galloping'

1998 / IV

+

100 pages / ISBN 90-407-1566-1

04. E.M. Houtman / W.J. Bannink / B.H. Timmerman, 'Experimental and Computational Study of a Blunt Cylinder-Flare Model in High Supersonic Flow'

1998 / VIII

+

40 pages / ISBN 90-407-1567-X

05. G.J.D. Zondervan, 'A Review of Propeller Modelling Techniques Based on Euler Methods'

1998/ IV

+

84 pages / ISBN 90-407-1568-8

06. M.J. Tummers I D.M. Passchier, 'Spectral Analysis of Individual Realization

LDA Data'

1998/ VIII

+

36 pages / ISBN 90-407-1569-6 07. P.J.J. Moeleker, 'Linear Temporal Stability Analysis'

1998/ VI

+

74 pages / ISBN 90-407-1570-X

08. B.W. van Oudheusden, 'Galloping Behaviour of an Aeroelastic Oscillator with Two Degrees of Freedom'

1998/ IV

+

128 pages / ISBN 90-407-1571-8

09. R. Mayer, 'Orientation on Quantitative IR-thermografy in Wall-shear Stress Measurements'

1998 / XII

+

108 pages / ISBN 90-407-1572-6

10. K.J.A. Westin / R.A.W.M. Henkes, 'Prediction of Bypass Transition with Differential Reynolds Stress Modeis'

1998 / VI

+

78 pages / ISBN 90-407-1573-4

11. J.L.M. Nijholt, 'Design of a Michelson Interferometer for Quantitative Refraction Index Profile Measurements'

1998/ 60 pages / ISBN 90-407-1574-2

12. R.A.W.M. Henkes / J.L. van Ingen, 'Overview of Stability and Transition in External Aerodynamics'

1998 / IV

+

48 pages / ISBN 90-407-1575-0

13. R.A.W.M. Henkes, 'Overview of Turbulence Models for External Aerodyna-mics'

(46)

Series 02: Flight Mechanics

01. E. Obert, 'A Method for the Determination of the Effect of Propeller Slip-stream on a Static Longitudinal Stability and Control of Multi-engined Aircraft'

1997/ IV

+

276 pages / ISBN 90-407-1577-7

02. C. Bill / F. van Dalen / A. Rothwell, 'Aircraft Design and Analysis System (ADAS)'

1997 1 X

+

222 pages 1 ISBN 90-407-1578-5

03. E. Torenbeek, 'Optimum Cruise Performance of Subsonic Transport

Air-craft'

1998 1 X

+

66 pages 1 ISBN 90-407-1579-3

Series 03: Control and Simulation

01. J.C. Gibson, 'The Definition, Understanding and Design of Aircraft Handling Oualities'

1997 1 X

+

162 pages 1 ISBN 90-407-1580-7

02. E.A. Lomonova, 'A System Look at Electromechanical Actuation for Primary

Flight Control'

1997 1 XIV

+

110 pages 1 ISBN 90-407-1581-5

03. C.A.A.M. van der Linden, 'DASMAT-Delft University Aircraft Simulation

Model and Analysis TooI. A Matlab/Simulink Environment for Flight Dyna-mics and Control Analysis'

19981 XII

+

220 pages 1 ISBN 90-407-1582-3

Series 05: Aerospace Structures and

Computional Mechanics

01. A.J. van Eekelen, 'Review and Selection of Methods for Structural Reliabili-ty Analysis'

1997 1 XIV

+

50 pages 1 ISBN 90-407-1583-1

02. M.E. Heerschap, 'User's Manual for the Computer Program Cufus. Ouick

Design Procedure for a CUt-out in a FUSelage version 1 .0' 19971 VIII

+

144pages/ISBN90-407-1584-X

03. C. Wohlever, 'A Preliminary Evaluation of the B2000 Nonlinear Shell

Element 08N.SM'

1998 1 IV

+

44 pages 1 ISBN, 90-407 -1585-8

04. L. Gunawan, 'Imperfections Measurements of a Perfect Shell with Specially

Designed Equipment (UNIVIMP)

(47)

Series 07: Aerospace Materials

01. A. Vasek / J. Schijve, 'Residual Strenght of Cracked 7075 T6 AI-alloy Sheets under High Loading Rates'

1997 / VI

+

70 pages / ISBN 90-407-1587-4

02. I. Kunes, 'FEM Modelling of Elastoplastic Stress and Strain Field in Centre-cracked Plate'

1997 / IV

+

32 pages / ISBN 90-407-1588-2

03. K. Verolme, 'The Initial Buckling Behavior of Flat and Curved Fiber Metal Laminate Panels'

1998 / VIII

+

60 pages / ISBN 90-407-1589-0

04. P.W.C. Prov6 Kluit, 'A New Method of Impregnating PEl Sheets for the

/n-Situ Foaming of Sandwiches'

1998 / IV

+

28 pages / ISBN 90-407-1590-4

05. A. Vlot / T. Soerjanto / I. Yeri / J.A. Schelling, 'Residual Thermal Stress es around Bonded Fibre Metal Laminate Repair Patches on an Aircraft Fusela-ge'

1998/ IV

+

24 pages / ISBN 90-407-1591-2

06. A. Vlot, 'High Strain Rate Tests on Fibre Metal Laminates' 1998 / IV

+

44 pages / ISBN 90-407-1592-0

07. S. Fawaz, 'Application of the Virtual Crack Closure Technique to Calculate Stress Intensity Factors for Through Cracks with an Oblique Elliptical Crack Front'

1998 / VIII

+

56 pages / ISBN 90-407-1593-9

08. J. Schijve, 'Fatigue Specimens for Sheet and Plate Material' 1998/ VI

+

18 pages / ISBN 90-407-1594-7

Series 08: Astrodynamics and Satellite Systems

01. E. Mooij, 'The Motion of a Vehicle in a Planetary Atmosphere' 1997 / XVI

+

156 pages / ISBN 90-407-1595-5

02. G.A. Bartels, 'GPS-Antenna Phase Center Measurements Performed in an Anechoic Chamber'

1997/ X

+

70 pages / ISBN 90-407-1596-3

03. E. Mooij, 'Linear Quadratic Regulator Design for an Unpowered, Winged Re-entry Vehicle'

A preliminary evaluation of the B2000 nonlinear shell element Q8N.SM

Series 05

Aerospace Structures and

Computional Mechanics 03