- Email: [email protected]
The impact modeling of structures using solid, shell, and membrane finite elements
Nuclear Engineering and Design 116 (1989) 101-116 North-Holland, Amsterdam
101
T H E I M P A C T M O D E L I N G OF S T R U C T U R E S U S I N G SOLID, SHELL, AND M E M B R A N E FINITE ELEMENTS Samuel W. KEY *, Sharon V. PETNEY, and Robert M. C L A N C Y R E / S P E C lnc., Post Office Box 14984 Albuquerque, New Mexico 87191, USA
Received 21 January 1988
A comparisonof modeling results and efficiencyis made for a 4-node, pressure corrected, plane strain quadrilateral element and the 4-node, constant-stress bending quadrilateral of Belytschko, Lin, and Tsay for modeling the transient dynamic response of a complex buried structure. A detailed derivation of the 4-node, constant-stress bending quadrilateral and associated hourglass control is presented.
1. Introduction For reasons of utility, efficiency, and ease of fabrication, many structural systems are a combination of components which, from a structural engineering perspective, are identified as solids, bars, beams, plates, shells, and membranes. Either by design intent or from likely accident events, these structures can be subject to impact loads. In either case, engineering analyses are required to support the design and establish that the required level of performance will be achieved. Historically, the tools and techniques for large deformation, finite strain, inelastic, transient dynamic analyses have developed in the context of analyzing solids. Recent work by Hughes and Liu [9,10]; Hughes, Liu, and Levit [11]; and Belytschko, Lin, and Tsay [4] in the development of low order shell and membrane elements has led to the possibility of combining these finite elements with those traditionally used in solids modeling. The result is an increase in the quality of impact calculations for structures built from these components. The improvements are in both numerical efficiency and modeling structural features.
2. Structural example Fig. 1 is an example of a buried structure intended to withstand the effects of an applied surface pressure. The pressure pulse is characterized by a high initial amplitude which decays exponentially. The structure is modeled with a vertical symmetry plane so that only the right half is pictured. The structure in question is a length of cylindrical pipe with periodically spaced, external flanges for stiffening, figs. 2 and 3. The ring stiffened cylinder rests on a rectangular frame whose main structural component is a pair of symmetrically placed, vertical plate girders with periodically spaced vertical stiffening flanges. Vertical loads in the cylinder are transmitted to the frame through a vertical plate which is continuous along the length of the cylinder. The assembly is fabricated by welding together steel components bent or cut from plate stock. The steel structure is isolated from the sides and top of the tunnel with additional structural components designed to have limited load-bearing capability. The frame rests firmly on the bottom of the * Present address: The MacNeal-SchwendlerCorporation, 815 Colorado Boulevard, Los Angeles, CA 90041, USA. 0029-5493/89/$03.50 © Elsevier Science Publishers B.V.
S.W. Key et aL / The impact modefing of structures
102
-,~}
-,{~
-300
-300
-¢00
-4100
-500
-50~
-600
-600
-700
-700
-800
-800
-900
. . . . . . :~ : : : T : ~ ? ii j i , i iil i i i
i :
T :
/
-900 I
I
I
0
100
200
I
I
3(]0
¢00
R-AXIS
I 500
I $00
i
0
I 100
I 200
I 300
I t00
1
I
500
600
R-RXI$
Fig. la. Plane strain cross section and mesh of a steel structure in a buried concrete tunnel showing the response at 0 ms to an exponentially decaying pressure pulse applied on the top of the mesh. Only the right half is required due to symmetry. Fig. l b . Plane strain cross section a n d mesh of a steel structure in a b u r i e d concrete tunnel s h o w i n g the response at 50 ms to an e x p o n e n t i a l l y d e c a y i n g pressure pulse applied on the top of the mesh. O n l y the fight half is required due to s y m m e t r y .
tunnel, fig. la. F r o m a modeling standpoint, the interaction between the steel s t r u c t u r e / l o a d limiting c o m p o n e n t s and the tunnel liner is through a sliding interface. W h e n a large deformation response calculation of two or more independently meshed bodies is undertaken, provisions for their physical interaction through impact and the sliding of one b o d y past another is a requisite part of the analysis. In general, the separately meshed bodies can be either in contact at the start of the calculation or come into contact during the course of their individual motions. Here, there is an initial gap between the load limiting c o m p o n e n t s and the tunnel liner. There is a sliding interface between the tunnel liner and the steel structure at the b o t t o m in the event the tunnel floor pulls away during the dynamic response.
3. Loading sequence The cylinder and frame combination is of sufficient length to permit obtaining useful results from a plane strain analysis, fig. 1. The pressure pulse after traveling through the overburden causes the tunnel liner to "cave in". After the gap at the tunnel crown closes, the steel cylinder experiences a nominally constant load transmitted by the load component, fig. 4. This load is applied directly to the cylinder crown, bypassing the external reinforcing flanges. The load is transferred around the cylinder and into the plate girder frame. The resulting velocity history of the b o t t o m horizontal plate tying the left and right plate girders together is delayed by about three milliseconds while the load travels through the cylinder
S.W. Key et aL / The impact modeling of structures
I
I
I
I
103
I
-240
-280
-320
N -'100
-440
-4BO
-520 I
I
I
I
-40
"10
80
120
R RXIS
Fig. 2. Plane strain cross section and mesh of the steel structure. The structure is composed of an externally stiffened cylinder attached to a frame which is composed of two vertical plate girders tied together laterally with continuous plates. Both the continuous components and the periodically spaced flanges are represented by 4-node plane strain quadrilaterals.
Fig. 3. Perspective view and mesh of the steel structure represented by plate elements. The periodically spaced stiffening flanges allow the use of symmetry boundary conditions to isolate a representative length for analysis.
100.0
I ' ' ' ' 1 ' ' ' ' 1 ' ' ' ' 1 ' ' ' '
0.0
0.0 -200.0
*
' - 100.0
-400.0 m
C
i
400.0
|
-200,0
>
3L.
~ -~0.0
"0 u = >~ -12oo.0
I~ -400.0
-1400.0 -500.0 -1600.0 -400.0
0.0
10.0
20.0
30.0
40.0
50,0
r n (,,~:) Fig. 4. History of the vertical stress element, acting on the crown of the cylinder. The m a x i m u m amplitude is at the top; the stress amplitude decays away from the top.
,
0.0
,
,
,
I
10,0
,
,
,
,
I
,
,
,
,
20.0
I
30.0
,
,
,
,
I
40.0
,
,
,
t
50.0
r , , ~ (m,*¢)
Fig. 5. History of the vertical velocity at the base of the steel structure.
104
S.W. Key et al. / The impact modeling of structures
and frame, fig. 5. Another 14 ms elapses before the cylinder and frame reach the nominal free-field velocity of 500 in./s, fig. 5. The deformed mesh at 50 ms shows extensive compaction of the overburden, extensive deformation of the tunnel liner, major displacements in the load limiting component at the top of the cylinder, the tunnel liner collapsed about the cylinder and frame, and a largely undeformed cylinder and frame, fig. lb.
4. Cylinder and frame modeling with solid elements Pressure corrected 4-node quadrilaterals with 2 x 2 Gaussian integration were used throughout the mesh, figs. 1 and 2. While entirely adequate as models of the continuous components (overburden, tunnel liner, load limiting components, cylinder, load bearing strut, and plate girder), special modifications are necessary to model the flange stiffeners of the cylinder and frame. This was achieved by reducing the density, Young's modulus, yield stress, and hardening modulus in those quadrilaterals representing the flanges. The reduction is based on the ratio of their thickness to their spacing. This results in two limitations; one affecting numerical efficiency and the other affectin~ modeling accuracy. Since the calculations were carried out with a conditionally stable, explicit central difference time integrator, the time step is controlled by the element with the smallest elastic wave transit time. In this case, the small square element where the steel load transfer strut of the cylinder meets the top of the plate girder of the frame, fig. 2. For comparison purposes, the ratio of the critical time step in just the overburden and the the tunnel liner alone to the full model is 10 to 1 (3.67 x 10-6/3.53 × 10 7). Second, the potential of "our-of-plane" twisting and, in extreme circumstances, crippling in the flanges is precluded by the modeling with plane strain quadrilaterals. However, in the present analysis these were not significant factors.
5. Four-node bending quadrilateral Since this is a structure built up from components rolled and cut from plate stock, a more nearly equivalent finite element model would use plate bending finite elements to model the separate structural components in their various orientations and intersections. In addition, plate and shell finite elements, with some exceptions, do not consider wave motion in their transverse or thickness direction. The result is a less restrictive bound on the time step below which conditionally stable time intergration schemes must operate. While bending behavior can, in general, be quite complex, the goal is to use an element formulation which is as elementary as possible. One is only interested in the constant membrane and bending response in a simple element. The computational effort in such an element approaches that of solid elements integrated with a "one-point" quadrature and leads to a like scale with the solid element for modeling stress gradients. Such an element is one adopted by Belytschko, Lin, and Tsay [4]. The four-node bending quadrilateral is based on conventional plate bending concepts in which the functional form of the deformation is made explicit in a coordinate normal to the reference surface. Consequently, the deformation is described in terms of displacement variables which only depend on the coordinates of the reference surface. Essentially, the unknown deformation is expanded in a Taylor series in a coordinate normal to the reference surface. Since the plate is thin, only a " f e w " terms are retained in the expansion.
105
S.W.. Key et al. / The impact modeling of structures
Xa
~2 14 la
1
"2 2
Fig. 6. Isoparametric coordinate representation of a flat twodimensional surface arbitrarily oriented in a fixed spatial coordinate system x i. 5.1. Kinematics
Consider a flat two-dimensional surface defined by the coordinates ~a and ~2 fig. 6. Normal to the surface is the coordinate ~3. The surface is a reference surface for the plate: the plate being defined by a thin layer of material of dimension h existing evenly distributed on both sides of the reference surface. The deformation of the plate is characterized by the deformation of the reference surface. The motion is defined in terms of a velocity vi of the reference surface and the region within a thin layer on either side defined by _+h / 2 , v,(x,)
= v , ( ~ ° ) - n k ( x * - x k ~~,, 0 ,m. n mto"~°~,~ ,,.
(5.1)
The greek indices have a range of two, the dimension of the plate reference surface. The spatial gradient of the velocity is given by Vi,j(X k) = Oi,j( ~ a) -- (~.imnnmton( ~a) ) n j -- I'lk( X k -- XkO)( ~.imnnrnto,nj( ~a) ).
(5.2)
In the special case when in isoparametric coordinates ~' coincide with the spatial coordinates x ~, ( x i ) = ( x , y, z } , the velocity of the plate is given by { oi } = ( Ux~ Ztoy, U y - Ztox, Oz },
(5.3)
the stretching, the symmetric part of the velocity gradient, is given by dxx = Vx,x + Ztoy,x,
2 d x y = Vx,y + Vy.x + z ( toy,y - W x,x),
dyy = Oy,y -- Ztox,y ,
2dxz = v,, x + toy,
d= = O,
2dyz = Uz,y - tox,
(5.4)
and the spin, the skew-symmetric part of the velocity gradient, is given by 2toxy = V~,y - Oy,x 'of-z ( toy,y -~ tox,x), 2tox~ = -V~ , x + toy,
(5.5)
2toy~ = - V~,y + % .
A number of observations can be made at this point. The assumptions about the deformation of the plate are those of Mindlin plate theory [12], and are the same assumptions as employed by Belytschko, Lin, and
106
S. W. Key et al. / The impact modeling of structures
Tsay [4]. In particular, transverse shear deformations are modeled. Fibers originally perpendicular to the reference surface and straight remain straight during the motion though not necessarily perpendicular. No extensional strains in the thickness direction are modeled or calculated. Using Mindlin plate theory allows the rotation rates wx and % to be described separately rather than defined by Vz,y and -vz,x, respectively, as in Kirchhoff plate theory. Consequently, the separate finite element assumptions on the velocity and rotation rates are required to be no more than continuous. In the event that the thin plate limit of vanishing transverse shear strains dx= and dy z holds, a Kirchhoff plate theory is recovered. This requires special treatment in the finite elements introduced below to prevent locking. Absent from the deformation is an explicit use of the rotation rate % about the normal vector to the reference surface. It is, however, retained in the development to follow. It is required for two applications: folded plate structures and shells modeled with flat plate finite element facets. When shells are modeled with flat plates, the global curvature properties of the shell are represented by the change in orientation from one quadrilateral element to the next quadrilateral element. In effect, the smooth variation in curvature of the original reference surface is approximated by discrete changes in orientation occurring at element edges; the elements being "chord" approximations to the original shell. This approximation is the same order as the constant membrane and bending stress approximations introduced below for the element integration. The four-node bending quadrilateral element relates the spatial coordinates x ~ to the nodal coordinates x~ through the isoparametric shape functions N 1 as follows: (5.6)
x ' = x~NI(l~'~).
For the four-node bending quadrilateral, the upper case indices have a range of four, corresponding to the four element nodes. The same shape functions are used to define the reference surface displacement field in terms of nodal displacements ui/, (5.7)
u, = u i z U ' ( l~'~).
Since these shape functions apply to spatial coordinates and displacements, their material derivative (represented by a superposed dot) must vanish. Hence, the velocity field and the rotational rate are given by v, = v , , N ' ( ~ ' ) ,
0:, = w , , N ' ( ~ " ) .
(5.8)
The velocity gradient v,,j and the gradient of the rotational rate ~ . s are defined as follows: v,.j = VilNtj,
w,,; = % l N I j .
(5.9)
The shape functions map a unit square in the isoparametric coordinates ~" to a general quadrilateral in the spatial coordinates x ~. The unit square is centered at the origin in the ~"-space so that the shape functions may be conveniently expanded in terms of an orthogonal set of base vectors, given in Table 1, as follows: =
1¢'~A1 + ~'~2F31.
(5.10)
The above vectors represent the displacement modes of a unit square, fig. 7. Though the velocity gradient of a bending quadrilateral is quite complex in description when using a Mindlin plate theory, eq. (5.2), the modes Z t and A~ combine to represent rates of rigid body translation and rotation, and uniform strain rates. The last vector /'] gives rise to modes with linear strain rate variations which are neglected by a mean strain rate quadrature. This vector defines the hourglass patterns for a unit square. The mode /'3/ is referred to as the hourglass base vector.
107
S.W. Key et aL / The impact modeling of structures
Table 1 Element deformation modes Node
~l
1
_! 2
2
1_
3
1
/~2
Zq
A~
_! 2
1
-1
_!
1
1
~
1
1
2~
1
2
2 -12
4
2 2
A½
F~
-1
-1
-1
1 1
-1
1
-1 1
5.2. Mean strain rate quadrature Prior to introducing the assumption of a constant stress quadrilateral in the reference surface it is convenient to deal with the explicit dependence of the velocity on the coordinate ~3 normal to the reference surface. The virtual energy in the Principal of Virtual Work is expanded as follows: ..
/'+1/2 f+l/2
/ ' + 1 / 2 i"
/~ax ,--~1
ft'd,:dv=J_i/2I_1/2I_1/2l'/Oi'jtl;)dCI~ :+1/2
:+1/2
:+1/2
i~:
d~ 2 d~ 3
m nz~ax\
- ,-1
- J - 1 / 2 1 - 1 / 2 J-l~2 t~[ei'~"n to (~ ))njdtl~
-J_1/2J-1/2 J-1/2 t •' + 1 / 2
; + 1 / 2
."+1/2
where h ~ 3 has been introduced for ~xr
ijt
m ('imnn
n "
d~ 2 d~ 3
a "x
~'J(~ ))h~3Jd~1 d~2 d~3,
nk(x k - X*O).The
(5.11)
dependence on ~3 is explicit since J specializes to
(5,12)
Ox-----~Snt"
J = h~rst ~1 ~2
The classical membrane and bending stress resultants, JV"ij and ..giJ, respectively, are introduced j~pij
= f+h/2tij
f+h/2
d ~ 3,
d~ = f+l/2htiJ
•" -- h / 2
" - 1/2
ij'"
(5.13)
-- ~f ) l / 2 h 2 t i j ~ 3
,.it '~j = J_/h/2 t S- d~" -
a/2
d~3"
i
i
I 1
' 'i'
i
J
, -J
]~1
t
i j
,
L_
_
1
A1
tl /' t ,,2't t'),
,,2'
Az1
FaI
Fig. 7. Unit quadrilateral and its displacement modes.
S. IV. Key et aL / The impact modeling of structures
108
A membrane stress r is and a moment stress bds given by
r , j = l s v . , j = f + l/2tu - 1/2
d~3, (5.14)
]LiJ= 1/~¢iJ= f + 1/2hti)~3 d~ 3, - 1/2 provide a reduced virtual internal energy;
+1/2 +1/2 ij ~a x f v tiJdijdO= f - l / 2 f l / 2 ~ vi'/(~ ) d d ~ ] d~2 f+l/2/,+1/2 .. - - L 1 / 2 J 1 / 2 "rtJ(C'imnnWo)n(~al)nJ J d~ 1 dl~ 2 r+l/2 r+l/2
-
...
J-,/2 J-,/= ~"("m°nm~°5(¢°))S d~a d~2"
(5.15)
The integrals in the reduced virtual internal energy, eq. (5.15), are evaluated using a mean stress, thereby considering only a state of constant membrane and bending stress within the element (Flanagan [6]; Flanagan and Belytschko [7[). The preceding expression is approximated as
fvt"d,,dv=V(÷%,+~
,
,,)).
-iJ - (,imnn m--n ~0 )nj-~L--ij ((imnn m~o~
(5.16)
The assumed constant stress fields are denoted by rris and ~'J, which will be referred to as mean stresses. It is assumed that the mean stresses depend only on the mean strains. Mean kinematic quantities are defined by integrating over the element as follows: ~, j
= (1/v) ff,,,
do,
~" = (1/v) f y ° do,
(5.17)
~", = (1/V) f ,4j dr. The gradient operator B/ is defined by
+1/2
B/= f~-l/2
+1/2. 1 - • ~2 L 1 / 2 N , j J CI~ d~ 2,
(5.18)
and an averaging operator A 1 is defined by
AI
f+l/2 f+l/21 o-a/2 N J d ~ a
='-1/2
d~2"
(5.19)
The mean velocity gradient is then given by
~,,j = ( 1 / V ) v , , a / ,
(5.20)
the mean rotational rate gradient by
~7, = ( I l V ) ~ T B / '
(5.21)
S.W.. Key et al. / The impact modeling of structures
109
and the mean rotational rate by (5.22)
~, = (1/V)~TA'. Thus, the virtual internal energy becomes,
fg'Ja, d o =
+
(5.23)
The nodal forces are then given by f i t = rrUB/'
(5.24)
and the nodal torques by
gin= --~.imnnm(~iJnjA' -~ ~_LiJo/).
(5.25)
The development of the gradient operator B/ follows the same steps as detailed in Flanagan [6]; the averaging operator A 1 follows similar arguments. 5.3. Corotational constitutive modeling By its very nature, a shell is a structural component whose thickness is very small compared to the lateral extent or radius of curvature at any point on the shell. As a result, the stress normal to the shell, pressure carrying capacities not withstanding, vanishes compared to the in-plane stresses. Thus, the constitutive models used with shells invariably incorporate a plane stress assumption. That is the case here. As a result, it is far more convenient to implement the constitutive model in a coordinate system local to the shell. The orientation of the normal vector n i is defined by using the cross product of two vectors in the reference surface defining the quadrilateral. If r i and s / are defined as unit vectors between opposite midsides of the quadrilateral, fig. 6, their cross product defines a direction n i, n i = ciJkrjsk,
(5.26)
which is perpendicular to the reference surface along two bisectors. Since, in general, the vectors r i and s i are not orthogonal, a further step is needed to obtain a local orthogonal Cartesian system " a t t a c h e d " to the shell quadrilateral. The system { r i, s i, n ~} is orthogonalized by redefining s ~ to be the cross product of n i and r i. This is a local coordinate system rotating with the shell quadrilateral in which the constitutive model is evaluated [Belytschko and Hsieh [2]; Belytschko and Marchertas [3]. With this procedure, only two-dimensional, plane stress constitutive modeling is required. To accommodate large in-plane shearing which can occur in very energetic applications, a large strain constitutive formulation is retained. However, only the in-plane spin component wrs must be retained, cf. Belytschko, Lin, and Tsay [4], since the corotational coordinate system automatically tracks arbitrary rigid body rotations. The spin components associated with transverse shearing, ~r, and % , , are, however, neglected in the finite strain constitutive modeling. 5. 4. Stress resultant evaluation While the constant stress assumptions result in a mean gradient operator B/ and averaging operator A ~ which act to produce mean strain rates linear through the thickness of the plate, eq. (5.4), material nonlinearities easily result in stress variations through the thickness which are anything but linear. As a result, the integrals which define the membrane and bending stress resultants, eq. (5.13), are implemented
110
S. W. Key et al. / The impact modeling of structures
using a numerical quadrature. Any number of quadrature rules are available, however, the "extreme fiber" stresses occurring at ~"= + h / 2 are generally of the greatest interest and the place where inelastic behavior starts in the presence of bending. Thus, the trapezoidal rule which includes the end points of the interval in the evaluation of the integrand is used (Abramowitz and Stegun [1]). The trapezoidal rule with uniform intervals provides
h/:
= m
LY
'"
which, when applied to the membrane stress resultant, eq. (5.13a), gives jv.iJ = h
t(0) ij + ij t~"n) --m 2 --+tO) "'' +t~m-l~ + 2 j'
(5.28)
and, when applied to the bending stress resultant, eq. (5.13b), gives J t '~j = -m-
~
t~0) -~1)-~1) + ...
l)t~,~ 1) + \
1--i~--m }t~m) .
(5.29)
Thus, the integrals are reduced to evaluating the constitutive equations at a number of stations through the thickness. If the response is elastic, only two stations through the thickness are required without twisting and three stations through the thickness with twisting. If the response is inelastic, no more than five or six stations through the thickness are normally required.
5.5. Bending performance The evaluation of the strain energy in the bending of thick plates reveals an overestimation of the transverse shear contribution. This is corrected by scaling the transverse shear energy back by x~l ) = 4 / 5 (Cowper, [5]). The "correction" is related to the discrepancy between the constant distribution of shear strain predicted by the displacement assumptions and the more nearly parabolic true distribution of shear strain. In the limit of reducing thickness, a plate bending theory with transverse shear becomes arbitrarily stiff in transverse shear response and the transverse shear strains vanish. The result is a Kirchhoff bending theory. Computationally, the result is a 1 / h 2 growth in transverse shear strain energy and an assortment of computational ills generally referred to as shear locking, (Fried, Johnson, and Tessler, [8]). The benefits of C ° interpolation for the velocity vi and the rotational rate ~' disappear. However, if the transverse shear strain energy is viewed as a penalty function approach to approximately satisfying the thin plate limit, the 1 / h 2 growth may be modified and the associated numerical ills avoided. The work of Fried et al. [8] quantifies this approach and is used here. Fried et al., demonstrate conclusively that the minimal beam bending element with quadratic displacement convergence and the strain energy integrated exactly requires a quadratic interpolation of the velocity vi coupled to a linear interpolation of the rotation rate ~i. However, the constant membrane, constant bending stress approximation introduced by Belytschko, Lin, and Tsay [4] results in the same stiffness representation. Thus, the work of Fried et al., will be applied to the present bending element. If A is the area of the element, a transverse shear strain energy scaling of K~2) = 6h2/A provides, in the limiting case of thin plate behavior, quadratic displacement convergence to the Kirchhoff bending result without shear locking in the element. In effect, the transverse shear stiffness is converted to a penalty condition at the element level to enforce, in an approximate fashion, the vanishing of the transverse shear strain dzz ( = d , , ) and d y z ( = ds, ) implying the necessary equivalence between vz,x and a~y, and -Vz,y and o~x.
S. W. Key et al. / The impact modeling of structures
111
A joint shear correction factor r z, which implements switching between thick and thin behavior, is given by x z = m i n ( 0 . 8 / h 2, 6 / A ) h 2.
(5.3)
To implement the transverse shear correction, the transverse shear strain rates are scaled black by before being used in the constitutive evaluation and the resulting transverse shear stresses are scaled by K before the divergence calculation is made.
5.6. Orthogonal hourglass control The mean stress-mean strain rate formulation considers only the linear part of the velocity field. The remaining portion of the nodal velocity is the so called hourglass field. Excitation of these modes may lead to severe, unresisted mesh distortion. A method for isolating the hourglass modes so that they may be treated independently of the rigid body and uniform strain rate modes is required. This is accon~plished by developing an hourglass "gradient operator" connected with hourglass restoring forces. The linear velocity field on which the mean strain rates are based is given by vLIN = Uio -t- ( ~ -- nJ~k )( X k - Xko )Vij, _
1
1
.,,j = -¢ (.,,Bj - ( ,,m.n mco°-,,,A)n, - n~(x~ - x ~ ) ( , , . ° : c o ; ~ ' ) ) , COLIN=
(5.3a)
1 co,o+ 5(,<,x~)<~,,~/. .
The hourglass velocity and rotational rate fields v~ G and coHo may be defined by removing the linear portion of the velocity and rotational rate fields. Thus, v~ G = v , - v~'Y,
CO~G= co,- CO,iN,
(5.32)
or in terms of nodal velocities and nodal rotational rates, H G
.,z
1 j -,,,,-.,o~,--:(,<,-x~,)(,,,:/-(,..°n 1
./
m
n
co~A
J
)",),
(5.33)
,
coi",° = co,z - co,o.~, - --¢ (x, - x~,)co,.,Bl, where
Vio
=
1 I zviz,~,
1 Z ¢0i0= zcoiz2~,
and
x 0i = z1x zi-~z ~ •
(5.34)
Advantage has been taken of nk(x~ - x~,~t) being zero or near zero at the nodes. The hourglass velocity fields, eq. (5.33), are in the improper null space of the gradient operation. (For the rotational rate the gradient operator is B/, but for the velocity the operator is slightly more complex incorporating the average value of the rate of rotation.) Since the rotational rate normal to the quadrilateral co, = coin~ does not participate, the four degrees of freedom it represents may be considered in the improper null space of the mean gradient operator evaluation. The linear velocity and rotational rate fields, eq. (5.31), span 15 degrees of freedom: 3 rates of rigi d body translation, 8 uniform rates of the velocity gradient (the velocity gradient normal to the reference surface v,,, = v~.jnin ~ is not included), and 4 linear rates of the velocity gradient normal to the reference surface, which means that the hourglass subspace is the remaining 5 degrees of freedom - exclusive of the four degrees of freedom represented by co, = co~n(
S. W. Key et al. / The impact modeling of structures
112
Two of the hourglass modes are the familiar inplane " k e y s t o n e " or " b e n d i n g " modes such as the one pictured in fig. 7 for F3t. The third m o d e is an out-of-plane warping mode. The fourth and fifth hourglass modes are related to complex twists involving linear variations of the gradient of the rate of rotation coi. An hourglass gradient operator is constructed from the hourglass basis vector /~1 s as follows:
G I= (V/~)I~S,
(5.35)
where (~ is a generalized element dimension developed in connection with the stability criterion for the explicit central difference time integration (Flanagan, [6]), V is the volume of the element. This scaling provides the hourglass gradient operator with the same dimensional characteristic as the uniform gradient o p e r a t o r B/. The hourglass strain rates ~, and ~/i are developed with G 1 operating on only the hourglass velocity viiH G and rotational rate w,lHG,
q, = (1/V)v,~ oG',
//i = ( l / V ) ~0,,.o,-,lc,.
(5.36)
Alternatively, unrestricted operators m a y be developed by requiring them to satisfy the following condition: Vii71
/ n + ~inb~l U G .o ~l , = v,#
lois71 = ~0HGG/.
(5.37)
Using the hourglass velocity and rotational rate, eq. (5.33), provides
,,,7' +n/.
<0.71=[
(xS- X Zl) ¢ojBJ]G,,
(5.38)
which, when rearranged and using the orthogonality of the mode shapes 521 and Fss, [ I F 3 # = 0, gives 1 j v,,71 + B/.w 7 = v,, G s - -~xJG B)~\) + %,,,,n " ~,t~xjG'A " / 1 j J i ,j), wiD, z = w,l ( G I -
1
j
j
(5.39)
1~
The conditions for unrestricted operators, eq. (5.37), is satisfied if the hourglass operators 7 / and fl[, are defined as
71- (V/,~)(-r'3s- (1/V)xjl"3XB/),
flit = (1/,S)(,,.,.n"xJUA'nj),
(5.40)
from which the hourglass rates are c o m p u t e d as
4, = (1/V)(v, sT'+ ~°Tfl/,),
/1, = ( 1 / V ) o ~ , , 7 ' .
(5.41)
The hourglass strain rate /l~ selects the F / m o d e from ~0, = ~0,n*. T o control the hourglass modes, generalized forces Q~ and torques T s are defined which are conjugate to O~ and /l~, respectively, so that the work rate is given by V tell O.~n 1 __ i l J H G -I- I g n H G - -
VQ'dL +
VT.il".
(5.42)
Utilizing eq. (5.38), the contribution to the nodal forces due to hourglass resistance is given by
f~
=
Q'7',
g .'. o
, i + LT'. = Q/~,,
(5.43)
S. 14:. K e y et aL / The impact modeling o f structures
113
The hourglass restoring forces are calculated from ~i
= ~2/'ttan~
ij
•
•
n
Tn = C2/'ttanSnrn'O '
qj,
(5.44)
where 2~tan is the tangent shear stiffness obtained from the deviatoric constitutive behavior of the mean stress and mean strain state in the element and c is a scaling. The scaling c assures the level of the hourglass restoring forces remains below that of the mean deviatoric stress state. The deviatoric behavior is used since the hourglass modes are constant volume, higher order straining modes of the element. The invariant time derivative of the generalized forces Q+ and torques T, accounts for the finite rotations expected in the appfication of the quadrilateral element in analyzing large deformation, transient dynamic phenomena. The derivative is given by (5.45) where ~0ij is the spin, (vi, / - vj,i)/2. The hourglass restoring forces f ~ and torques g,Hcl are added to those obtained from the divergence of the mean stress state so that the complete result is
:"=(+'+"+e',/), 1~
m
-ij
1
""
i 1
(5.46)
6. Cylinder and frame modeling with plate elements To improve both numerical efficiency and the representation of the structural components, the cylinder and frame combination was remeshed using 4-node plate bending elements, fig. 3. The periodic spacing of the stiffening flanges along the cylinder and the plate girder of the frame permits a " u n i t " length of the cylinder and frame combination to be isolated with symmetry boundary conditions for the purposes of analysis, fig. 3. Since the thickness of the various structural components is carried as a parameter in the plate element, the individual elements appear without any thickness in perspective mesh plots, fig. 3. The spacing for the exterior, tapered flange stiffeners on the plate girder was several times that of the other flange stiffeners and, therefore, was dropped from the plate element model, fig. 3. To produce an isolated analysis of the cylinder and frame modeled with plate elements, the loads observed at the top of the cylinder in the solid element model, fig. 4, were used, along with the velocity history observed at the base of the solid element model, fig. 5. This produces comparable results as long as the structural models closely resemble each other in mass and stiffness distribution.
7. Comparison of modeling results There are four major features of the response on which a comparison can be based: the history of the relative displacement between the top and bottom of the cylinder, the effective stress history in the plate girder of the frame, the distribution of effective stress compared to the yield stress in the stiffeners on the cylinder at the time of maximum ovalization, and the deformed shape at the time of maximum ovalization. The history of the relative displacement between the top and bottom of the cylinder, fig. 8, shows an oscillation period in ovalization for the two structural models which are within a few percent of each other. The long wave length, low frequency response of the two models is essentially the same. The vertical load acting against the lower kinematic boundary in the plate element model shows a static ovaling deflection
114
S. 14( K e y et aL / The impact modeling of structures
....
I
. . . .
I ' ' ' ' 1 ' ' ' ' 1
....
I ....
I ....
I ....
1 ....
1''''
m~
....
0.0
I ....
P....
]
"0.2
l- J---l -l ~~Ki'k~l 1~ "
40000,0
jl ! /
I It
I
I - -
-ff~l~,~mm~
--'.,'l,lll@mlr~
I
I I
-12
~X
I
i
~ i I I
/
tt\]
I I
I X
I I I
/
~
I
t
/ /
I I I
r
I \x/ 0,0
-IJl 0.0
o~
Fig. 8. History of vertical diametral change in the cylinder in the full model based on solid elements and in the isolated model based on plate elements.
,,,I,,,,] Lo
l. 000[-01
=o~
=o.o
3s,o
~0,0
4513
50.0
5A 9 YIO 1. 000([-01 Z. (:X3OE-Ol 3.000£-01 q. O00E-OJ
6. O00E-Ol
E 5. O00E-01 F 6. O00E-OI
Fig. 10. The d i s t r i b u t i o n of effective stress in the m o d e l b a s e d on solid elements at the time of m a x i m u m d i a m e t r a l deflection in the cylinder (22 ms) expressed as a p r o p o r t i o n of the yield stress.
zs.o r,,. (,,~)
I,, ,~o
Fig. 9. History of vertical stress in the vertical plate girder of
A B C O
7,000£-0]
Itt
.... ,o~
2. O00~-OI 3. O00E-O] q. O00E-O]
5. O00~-Ol
\\1~/
the frame in the full m o d e l b a s e d o n solid e l e m e n t s a n d in the i s o l a t e d m o d e l b a s e d o n p l a t e elements.
5IG YLD R B C 0 E F G
if
X\
G 7. O00E-Ol
0 C ~0 C Fig. 11. The distribution of effective stress in the model based on plate elements at the time of m a x i m u m diametral deflection in the cylinder (22 ms) e x p r e s s e d as a p r o p o r t i o n of the yield stress.
S. I,E. Key et a L / The impact modeling of structures
115
r~
Fig. 12. Distorted cross section based on plate elements at the time of m a x i m u m diametral deflection in the cylinder (22 ms). Displacements are magnified by 5.
Llrldeformed
Idogrtificatlon=5
of 0.7 in. while the actual compliant boundaries in the solid element model appear to allow the cylinder and frame to either unload at late times or switch deformation modes. In either case insignificant deflections occur in the flange stiffened cylinder considering its 160 in. vertical inside diameter. The effective stress history in the plate girder of the frame shows a nominally similar behavior, fig. 9. Between 15 and 30 ms, there is an oscillation about an effective stress of 15 000 psi. The effective stress then rises to approximately 46 000 psi at 35 ms., after which it decays to 30 000 psi at 50 ms. The model based on plate elements exhibits an unloading of the frame between 35 and 40 ms at the same time the cylinder exhibits a second peak in its ovalization behavior. While this makes sense from a structural vibration point of view, a similar behavior is absent, or at least greatly diminished, in the model based on solid elements. Perhaps, in the more complete model based on solid elements, the compliance of the tunnel invert masks this behavior. The distribution of effective stress compared to the yield stress in the stiffeners on the cylinder shows the m a x i m u m stresses occur in the flange in compression at the crown at the time of m a x i m u m owli7ation, Figs. 10 and 11. There is an opposite but lower magnitude tensile stress in the flange at the midheight of the cylinder as would be expected in the vertical shortening phase of the ovalization. It would appear that the extreme fiber of the stiffening flange at the top of the cylinder may exhibit plastic yielding, but the constant membrane stress approximation in the plate elements does not allow it to be identified explicitly. The deformed shape at the time of m a x i m u m ovalization is shown in fig. 12. Fig. 12 clearly illustrates the deformation m o d e of the shell model compared to the undeformed shape. These deformed shapes are consistent with the diametral load on the cylinder, except perhaps the rolling of the upper flange on the plate girder. This is a bolted connection which was modeled as perfectly bonded. This result suggests a more detailed model of this connection should be constructed. The outcome could be a recommendation for a more compliant connection or the placement of the inner bolts midway between flanges.
116
S. 14( Key et aL / The impact modeling of structures
8. Evaluation of numerical efficiency Even though the plate element model of the cylinder a n d frame was analyzed as a n isolated structure using loads a n d k i n e m a t i c b o u n d a r y c o n d i t i o n s observed in the full model m a d e up exclusively of p l a n e strain solid elements, it is possible to estimate r e a s o n a b l y closely the savings to be gained from using plate elements in place of solid elements for the steel cylinder a n d frame. The increase in critical time step from 0.35 ms to 1.40 ms w o u l d represent a factor of 4 reduction in effort for the full model.
9. Conclusion Overall, the use of plate elements did not reveal a n y significantly different structural behavior in this analysis. F o r this structure, plate elements do offer the o p p o r t u n i t y to model local detail n o t afforded by the use of p l a n e strain solid elements. F r o m the s t a n d p o i n t of n u m e r i c a l efficiency, there is a factor of 4 r e d u c t i o n in c o m p u t a t i o n effort to be seen by using plate elements.
Acknowledgements S u p p o r t for the initial phase of this work b y the U S A F Ballistic Missile Office, N T S Engineering, a n d A g b a b i a n Associates is gratefully acknowledged. The design a n d effective p e r f o r m a n c e of the steel cylinder a n d frame exhibited in these calculations is the p r o d u c t of the skill a n d experience of H a r o l d F u t t r u p , A g b a b i a n Associates.
References [1] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Applied Mathematics Series, No. 55, Tenth Printing, National Bureau of Standards, United States Department of Commerce (1972). [2] T. Belytschko and B.J. Hsieh, Nonlinear transient finite element analysis with convected coordinates, Internat. J. for Numer. Meths. in Engrg. 7 (1973) 255-271. [3] T. Belytschko and A.H. Marchertas, Nonlinear finite element method for plates and its application to the dynamic response of reactor fuel subassemblies, Trans. ASME, J. of Press. Vessel Teclmol. (1974) 251-257. [4] T. Belytschko, J.I. Lin, and C.S. Tsay, Explicit algorithms for the nonlinear dynamics of shells, Computer Meths. in Appl. Mechs. Engrg. 42 (1984) 225-251. [5] G.R. Cowper, On the accuracy of Timoshenko's beam theory, J. of the Engineering Mechanics Division, EM6, Proceedings of ASCE, pp. 1447-1453 (1968). [6] D.P. Flanagan, Numerical Techniques with One-Point Integrated Hexahedron and Quadrilateral Elements, Ph.D. Dissertation, Northwestern University, Evanston, IL (1981). [7] D.P. Flanagan and T. Belytschko, A uniform strain hexahedron and quadrilateral with orthogonal hourglass control, Internat. J. for Numer. Meths. in Engrg. 17 (1981) 679-706. [8] I. Fried, A. Johnson, and A. Tessler, Minimal-degree thin triangular plate and shell bending finite elements of order two and four, Computer Meths. in Appl. Mechs. Engrg. 56 (1986) 283-307. [9] T.J.R. Hughes and W.K. Liu, Nonlinear finite element analysis of shells: Part I. Three-dimensions shells, Computer Meths. in Appl. Mechs. Engrg. 26 (1981) 331-362. [10] T.J.R. Hughes and W.K. Liu, Nonlinear finite element analysis of shells: Part II. Two-dimensional shells, Computer Meths. in Appl. Mechs. Engrg. 27 (1981) 167-181. [11] T.J.R. Hughes, W.K. Liu, and I. Levit, Nonlinear dynamic finite element analysis of shells, in: Nonlinear Finite Element Analysis in Structural Mechanics, Wunderlich, Stein, and Bathe (Eds.) (Springer-Verlag, 1981) pp. 151-168. [12] R.D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates, J. of Appl. Mechs. 18 (1951) 31-38.
101
T H E I M P A C T M O D E L I N G OF S T R U C T U R E S U S I N G SOLID, SHELL, AND M E M B R A N E FINITE ELEMENTS Samuel W. KEY *, Sharon V. PETNEY, and Robert M. C L A N C Y R E / S P E C lnc., Post Office Box 14984 Albuquerque, New Mexico 87191, USA
Received 21 January 1988
A comparisonof modeling results and efficiencyis made for a 4-node, pressure corrected, plane strain quadrilateral element and the 4-node, constant-stress bending quadrilateral of Belytschko, Lin, and Tsay for modeling the transient dynamic response of a complex buried structure. A detailed derivation of the 4-node, constant-stress bending quadrilateral and associated hourglass control is presented.
1. Introduction For reasons of utility, efficiency, and ease of fabrication, many structural systems are a combination of components which, from a structural engineering perspective, are identified as solids, bars, beams, plates, shells, and membranes. Either by design intent or from likely accident events, these structures can be subject to impact loads. In either case, engineering analyses are required to support the design and establish that the required level of performance will be achieved. Historically, the tools and techniques for large deformation, finite strain, inelastic, transient dynamic analyses have developed in the context of analyzing solids. Recent work by Hughes and Liu [9,10]; Hughes, Liu, and Levit [11]; and Belytschko, Lin, and Tsay [4] in the development of low order shell and membrane elements has led to the possibility of combining these finite elements with those traditionally used in solids modeling. The result is an increase in the quality of impact calculations for structures built from these components. The improvements are in both numerical efficiency and modeling structural features.
2. Structural example Fig. 1 is an example of a buried structure intended to withstand the effects of an applied surface pressure. The pressure pulse is characterized by a high initial amplitude which decays exponentially. The structure is modeled with a vertical symmetry plane so that only the right half is pictured. The structure in question is a length of cylindrical pipe with periodically spaced, external flanges for stiffening, figs. 2 and 3. The ring stiffened cylinder rests on a rectangular frame whose main structural component is a pair of symmetrically placed, vertical plate girders with periodically spaced vertical stiffening flanges. Vertical loads in the cylinder are transmitted to the frame through a vertical plate which is continuous along the length of the cylinder. The assembly is fabricated by welding together steel components bent or cut from plate stock. The steel structure is isolated from the sides and top of the tunnel with additional structural components designed to have limited load-bearing capability. The frame rests firmly on the bottom of the * Present address: The MacNeal-SchwendlerCorporation, 815 Colorado Boulevard, Los Angeles, CA 90041, USA. 0029-5493/89/$03.50 © Elsevier Science Publishers B.V.
S.W. Key et aL / The impact modefing of structures
102
-,~}
-,{~
-300
-300
-¢00
-4100
-500
-50~
-600
-600
-700
-700
-800
-800
-900
. . . . . . :~ : : : T : ~ ? ii j i , i iil i i i
i :
T :
/
-900 I
I
I
0
100
200
I
I
3(]0
¢00
R-AXIS
I 500
I $00
i
0
I 100
I 200
I 300
I t00
1
I
500
600
R-RXI$
Fig. la. Plane strain cross section and mesh of a steel structure in a buried concrete tunnel showing the response at 0 ms to an exponentially decaying pressure pulse applied on the top of the mesh. Only the right half is required due to symmetry. Fig. l b . Plane strain cross section a n d mesh of a steel structure in a b u r i e d concrete tunnel s h o w i n g the response at 50 ms to an e x p o n e n t i a l l y d e c a y i n g pressure pulse applied on the top of the mesh. O n l y the fight half is required due to s y m m e t r y .
tunnel, fig. la. F r o m a modeling standpoint, the interaction between the steel s t r u c t u r e / l o a d limiting c o m p o n e n t s and the tunnel liner is through a sliding interface. W h e n a large deformation response calculation of two or more independently meshed bodies is undertaken, provisions for their physical interaction through impact and the sliding of one b o d y past another is a requisite part of the analysis. In general, the separately meshed bodies can be either in contact at the start of the calculation or come into contact during the course of their individual motions. Here, there is an initial gap between the load limiting c o m p o n e n t s and the tunnel liner. There is a sliding interface between the tunnel liner and the steel structure at the b o t t o m in the event the tunnel floor pulls away during the dynamic response.
3. Loading sequence The cylinder and frame combination is of sufficient length to permit obtaining useful results from a plane strain analysis, fig. 1. The pressure pulse after traveling through the overburden causes the tunnel liner to "cave in". After the gap at the tunnel crown closes, the steel cylinder experiences a nominally constant load transmitted by the load component, fig. 4. This load is applied directly to the cylinder crown, bypassing the external reinforcing flanges. The load is transferred around the cylinder and into the plate girder frame. The resulting velocity history of the b o t t o m horizontal plate tying the left and right plate girders together is delayed by about three milliseconds while the load travels through the cylinder
S.W. Key et aL / The impact modeling of structures
I
I
I
I
103
I
-240
-280
-320
N -'100
-440
-4BO
-520 I
I
I
I
-40
"10
80
120
R RXIS
Fig. 2. Plane strain cross section and mesh of the steel structure. The structure is composed of an externally stiffened cylinder attached to a frame which is composed of two vertical plate girders tied together laterally with continuous plates. Both the continuous components and the periodically spaced flanges are represented by 4-node plane strain quadrilaterals.
Fig. 3. Perspective view and mesh of the steel structure represented by plate elements. The periodically spaced stiffening flanges allow the use of symmetry boundary conditions to isolate a representative length for analysis.
100.0
I ' ' ' ' 1 ' ' ' ' 1 ' ' ' ' 1 ' ' ' '
0.0
0.0 -200.0
*
' - 100.0
-400.0 m
C
i
400.0
|
-200,0
>
3L.
~ -~0.0
"0 u = >~ -12oo.0
I~ -400.0
-1400.0 -500.0 -1600.0 -400.0
0.0
10.0
20.0
30.0
40.0
50,0
r n (,,~:) Fig. 4. History of the vertical stress element, acting on the crown of the cylinder. The m a x i m u m amplitude is at the top; the stress amplitude decays away from the top.
,
0.0
,
,
,
I
10,0
,
,
,
,
I
,
,
,
,
20.0
I
30.0
,
,
,
,
I
40.0
,
,
,
t
50.0
r , , ~ (m,*¢)
Fig. 5. History of the vertical velocity at the base of the steel structure.
104
S.W. Key et al. / The impact modeling of structures
and frame, fig. 5. Another 14 ms elapses before the cylinder and frame reach the nominal free-field velocity of 500 in./s, fig. 5. The deformed mesh at 50 ms shows extensive compaction of the overburden, extensive deformation of the tunnel liner, major displacements in the load limiting component at the top of the cylinder, the tunnel liner collapsed about the cylinder and frame, and a largely undeformed cylinder and frame, fig. lb.
4. Cylinder and frame modeling with solid elements Pressure corrected 4-node quadrilaterals with 2 x 2 Gaussian integration were used throughout the mesh, figs. 1 and 2. While entirely adequate as models of the continuous components (overburden, tunnel liner, load limiting components, cylinder, load bearing strut, and plate girder), special modifications are necessary to model the flange stiffeners of the cylinder and frame. This was achieved by reducing the density, Young's modulus, yield stress, and hardening modulus in those quadrilaterals representing the flanges. The reduction is based on the ratio of their thickness to their spacing. This results in two limitations; one affecting numerical efficiency and the other affectin~ modeling accuracy. Since the calculations were carried out with a conditionally stable, explicit central difference time integrator, the time step is controlled by the element with the smallest elastic wave transit time. In this case, the small square element where the steel load transfer strut of the cylinder meets the top of the plate girder of the frame, fig. 2. For comparison purposes, the ratio of the critical time step in just the overburden and the the tunnel liner alone to the full model is 10 to 1 (3.67 x 10-6/3.53 × 10 7). Second, the potential of "our-of-plane" twisting and, in extreme circumstances, crippling in the flanges is precluded by the modeling with plane strain quadrilaterals. However, in the present analysis these were not significant factors.
5. Four-node bending quadrilateral Since this is a structure built up from components rolled and cut from plate stock, a more nearly equivalent finite element model would use plate bending finite elements to model the separate structural components in their various orientations and intersections. In addition, plate and shell finite elements, with some exceptions, do not consider wave motion in their transverse or thickness direction. The result is a less restrictive bound on the time step below which conditionally stable time intergration schemes must operate. While bending behavior can, in general, be quite complex, the goal is to use an element formulation which is as elementary as possible. One is only interested in the constant membrane and bending response in a simple element. The computational effort in such an element approaches that of solid elements integrated with a "one-point" quadrature and leads to a like scale with the solid element for modeling stress gradients. Such an element is one adopted by Belytschko, Lin, and Tsay [4]. The four-node bending quadrilateral is based on conventional plate bending concepts in which the functional form of the deformation is made explicit in a coordinate normal to the reference surface. Consequently, the deformation is described in terms of displacement variables which only depend on the coordinates of the reference surface. Essentially, the unknown deformation is expanded in a Taylor series in a coordinate normal to the reference surface. Since the plate is thin, only a " f e w " terms are retained in the expansion.
105
S.W.. Key et al. / The impact modeling of structures
Xa
~2 14 la
1
"2 2
Fig. 6. Isoparametric coordinate representation of a flat twodimensional surface arbitrarily oriented in a fixed spatial coordinate system x i. 5.1. Kinematics
Consider a flat two-dimensional surface defined by the coordinates ~a and ~2 fig. 6. Normal to the surface is the coordinate ~3. The surface is a reference surface for the plate: the plate being defined by a thin layer of material of dimension h existing evenly distributed on both sides of the reference surface. The deformation of the plate is characterized by the deformation of the reference surface. The motion is defined in terms of a velocity vi of the reference surface and the region within a thin layer on either side defined by _+h / 2 , v,(x,)
= v , ( ~ ° ) - n k ( x * - x k ~~,, 0 ,m. n mto"~°~,~ ,,.
(5.1)
The greek indices have a range of two, the dimension of the plate reference surface. The spatial gradient of the velocity is given by Vi,j(X k) = Oi,j( ~ a) -- (~.imnnmton( ~a) ) n j -- I'lk( X k -- XkO)( ~.imnnrnto,nj( ~a) ).
(5.2)
In the special case when in isoparametric coordinates ~' coincide with the spatial coordinates x ~, ( x i ) = ( x , y, z } , the velocity of the plate is given by { oi } = ( Ux~ Ztoy, U y - Ztox, Oz },
(5.3)
the stretching, the symmetric part of the velocity gradient, is given by dxx = Vx,x + Ztoy,x,
2 d x y = Vx,y + Vy.x + z ( toy,y - W x,x),
dyy = Oy,y -- Ztox,y ,
2dxz = v,, x + toy,
d= = O,
2dyz = Uz,y - tox,
(5.4)
and the spin, the skew-symmetric part of the velocity gradient, is given by 2toxy = V~,y - Oy,x 'of-z ( toy,y -~ tox,x), 2tox~ = -V~ , x + toy,
(5.5)
2toy~ = - V~,y + % .
A number of observations can be made at this point. The assumptions about the deformation of the plate are those of Mindlin plate theory [12], and are the same assumptions as employed by Belytschko, Lin, and
106
S. W. Key et al. / The impact modeling of structures
Tsay [4]. In particular, transverse shear deformations are modeled. Fibers originally perpendicular to the reference surface and straight remain straight during the motion though not necessarily perpendicular. No extensional strains in the thickness direction are modeled or calculated. Using Mindlin plate theory allows the rotation rates wx and % to be described separately rather than defined by Vz,y and -vz,x, respectively, as in Kirchhoff plate theory. Consequently, the separate finite element assumptions on the velocity and rotation rates are required to be no more than continuous. In the event that the thin plate limit of vanishing transverse shear strains dx= and dy z holds, a Kirchhoff plate theory is recovered. This requires special treatment in the finite elements introduced below to prevent locking. Absent from the deformation is an explicit use of the rotation rate % about the normal vector to the reference surface. It is, however, retained in the development to follow. It is required for two applications: folded plate structures and shells modeled with flat plate finite element facets. When shells are modeled with flat plates, the global curvature properties of the shell are represented by the change in orientation from one quadrilateral element to the next quadrilateral element. In effect, the smooth variation in curvature of the original reference surface is approximated by discrete changes in orientation occurring at element edges; the elements being "chord" approximations to the original shell. This approximation is the same order as the constant membrane and bending stress approximations introduced below for the element integration. The four-node bending quadrilateral element relates the spatial coordinates x ~ to the nodal coordinates x~ through the isoparametric shape functions N 1 as follows: (5.6)
x ' = x~NI(l~'~).
For the four-node bending quadrilateral, the upper case indices have a range of four, corresponding to the four element nodes. The same shape functions are used to define the reference surface displacement field in terms of nodal displacements ui/, (5.7)
u, = u i z U ' ( l~'~).
Since these shape functions apply to spatial coordinates and displacements, their material derivative (represented by a superposed dot) must vanish. Hence, the velocity field and the rotational rate are given by v, = v , , N ' ( ~ ' ) ,
0:, = w , , N ' ( ~ " ) .
(5.8)
The velocity gradient v,,j and the gradient of the rotational rate ~ . s are defined as follows: v,.j = VilNtj,
w,,; = % l N I j .
(5.9)
The shape functions map a unit square in the isoparametric coordinates ~" to a general quadrilateral in the spatial coordinates x ~. The unit square is centered at the origin in the ~"-space so that the shape functions may be conveniently expanded in terms of an orthogonal set of base vectors, given in Table 1, as follows: =
1¢'~A1 + ~'~2F31.
(5.10)
The above vectors represent the displacement modes of a unit square, fig. 7. Though the velocity gradient of a bending quadrilateral is quite complex in description when using a Mindlin plate theory, eq. (5.2), the modes Z t and A~ combine to represent rates of rigid body translation and rotation, and uniform strain rates. The last vector /'] gives rise to modes with linear strain rate variations which are neglected by a mean strain rate quadrature. This vector defines the hourglass patterns for a unit square. The mode /'3/ is referred to as the hourglass base vector.
107
S.W. Key et aL / The impact modeling of structures
Table 1 Element deformation modes Node
~l
1
_! 2
2
1_
3
1
/~2
Zq
A~
_! 2
1
-1
_!
1
1
~
1
1
2~
1
2
2 -12
4
2 2
A½
F~
-1
-1
-1
1 1
-1
1
-1 1
5.2. Mean strain rate quadrature Prior to introducing the assumption of a constant stress quadrilateral in the reference surface it is convenient to deal with the explicit dependence of the velocity on the coordinate ~3 normal to the reference surface. The virtual energy in the Principal of Virtual Work is expanded as follows: ..
/'+1/2 f+l/2
/ ' + 1 / 2 i"
/~ax ,--~1
ft'd,:dv=J_i/2I_1/2I_1/2l'/Oi'jtl;)dCI~ :+1/2
:+1/2
:+1/2
i~:
d~ 2 d~ 3
m nz~ax\
- ,-1
- J - 1 / 2 1 - 1 / 2 J-l~2 t~[ei'~"n to (~ ))njdtl~
-J_1/2J-1/2 J-1/2 t •' + 1 / 2
; + 1 / 2
."+1/2
where h ~ 3 has been introduced for ~xr
ijt
m ('imnn
n "
d~ 2 d~ 3
a "x
~'J(~ ))h~3Jd~1 d~2 d~3,
nk(x k - X*O).The
(5.11)
dependence on ~3 is explicit since J specializes to
(5,12)
Ox-----~Snt"
J = h~rst ~1 ~2
The classical membrane and bending stress resultants, JV"ij and ..giJ, respectively, are introduced j~pij
= f+h/2tij
f+h/2
d ~ 3,
d~ = f+l/2htiJ
•" -- h / 2
" - 1/2
ij'"
(5.13)
-- ~f ) l / 2 h 2 t i j ~ 3
,.it '~j = J_/h/2 t S- d~" -
a/2
d~3"
i
i
I 1
' 'i'
i
J
, -J
]~1
t
i j
,
L_
_
1
A1
tl /' t ,,2't t'),
,,2'
Az1
FaI
Fig. 7. Unit quadrilateral and its displacement modes.
S. IV. Key et aL / The impact modeling of structures
108
A membrane stress r is and a moment stress bds given by
r , j = l s v . , j = f + l/2tu - 1/2
d~3, (5.14)
]LiJ= 1/~¢iJ= f + 1/2hti)~3 d~ 3, - 1/2 provide a reduced virtual internal energy;
+1/2 +1/2 ij ~a x f v tiJdijdO= f - l / 2 f l / 2 ~ vi'/(~ ) d d ~ ] d~2 f+l/2/,+1/2 .. - - L 1 / 2 J 1 / 2 "rtJ(C'imnnWo)n(~al)nJ J d~ 1 dl~ 2 r+l/2 r+l/2
-
...
J-,/2 J-,/= ~"("m°nm~°5(¢°))S d~a d~2"
(5.15)
The integrals in the reduced virtual internal energy, eq. (5.15), are evaluated using a mean stress, thereby considering only a state of constant membrane and bending stress within the element (Flanagan [6]; Flanagan and Belytschko [7[). The preceding expression is approximated as
fvt"d,,dv=V(÷%,+~
,
,,)).
-iJ - (,imnn m--n ~0 )nj-~L--ij ((imnn m~o~
(5.16)
The assumed constant stress fields are denoted by rris and ~'J, which will be referred to as mean stresses. It is assumed that the mean stresses depend only on the mean strains. Mean kinematic quantities are defined by integrating over the element as follows: ~, j
= (1/v) ff,,,
do,
~" = (1/v) f y ° do,
(5.17)
~", = (1/V) f ,4j dr. The gradient operator B/ is defined by
+1/2
B/= f~-l/2
+1/2. 1 - • ~2 L 1 / 2 N , j J CI~ d~ 2,
(5.18)
and an averaging operator A 1 is defined by
AI
f+l/2 f+l/21 o-a/2 N J d ~ a
='-1/2
d~2"
(5.19)
The mean velocity gradient is then given by
~,,j = ( 1 / V ) v , , a / ,
(5.20)
the mean rotational rate gradient by
~7, = ( I l V ) ~ T B / '
(5.21)
S.W.. Key et al. / The impact modeling of structures
109
and the mean rotational rate by (5.22)
~, = (1/V)~TA'. Thus, the virtual internal energy becomes,
fg'Ja, d o =
+
(5.23)
The nodal forces are then given by f i t = rrUB/'
(5.24)
and the nodal torques by
gin= --~.imnnm(~iJnjA' -~ ~_LiJo/).
(5.25)
The development of the gradient operator B/ follows the same steps as detailed in Flanagan [6]; the averaging operator A 1 follows similar arguments. 5.3. Corotational constitutive modeling By its very nature, a shell is a structural component whose thickness is very small compared to the lateral extent or radius of curvature at any point on the shell. As a result, the stress normal to the shell, pressure carrying capacities not withstanding, vanishes compared to the in-plane stresses. Thus, the constitutive models used with shells invariably incorporate a plane stress assumption. That is the case here. As a result, it is far more convenient to implement the constitutive model in a coordinate system local to the shell. The orientation of the normal vector n i is defined by using the cross product of two vectors in the reference surface defining the quadrilateral. If r i and s / are defined as unit vectors between opposite midsides of the quadrilateral, fig. 6, their cross product defines a direction n i, n i = ciJkrjsk,
(5.26)
which is perpendicular to the reference surface along two bisectors. Since, in general, the vectors r i and s i are not orthogonal, a further step is needed to obtain a local orthogonal Cartesian system " a t t a c h e d " to the shell quadrilateral. The system { r i, s i, n ~} is orthogonalized by redefining s ~ to be the cross product of n i and r i. This is a local coordinate system rotating with the shell quadrilateral in which the constitutive model is evaluated [Belytschko and Hsieh [2]; Belytschko and Marchertas [3]. With this procedure, only two-dimensional, plane stress constitutive modeling is required. To accommodate large in-plane shearing which can occur in very energetic applications, a large strain constitutive formulation is retained. However, only the in-plane spin component wrs must be retained, cf. Belytschko, Lin, and Tsay [4], since the corotational coordinate system automatically tracks arbitrary rigid body rotations. The spin components associated with transverse shearing, ~r, and % , , are, however, neglected in the finite strain constitutive modeling. 5. 4. Stress resultant evaluation While the constant stress assumptions result in a mean gradient operator B/ and averaging operator A ~ which act to produce mean strain rates linear through the thickness of the plate, eq. (5.4), material nonlinearities easily result in stress variations through the thickness which are anything but linear. As a result, the integrals which define the membrane and bending stress resultants, eq. (5.13), are implemented
110
S. W. Key et al. / The impact modeling of structures
using a numerical quadrature. Any number of quadrature rules are available, however, the "extreme fiber" stresses occurring at ~"= + h / 2 are generally of the greatest interest and the place where inelastic behavior starts in the presence of bending. Thus, the trapezoidal rule which includes the end points of the interval in the evaluation of the integrand is used (Abramowitz and Stegun [1]). The trapezoidal rule with uniform intervals provides
h/:
= m
LY
'"
which, when applied to the membrane stress resultant, eq. (5.13a), gives jv.iJ = h
t(0) ij + ij t~"n) --m 2 --+tO) "'' +t~m-l~ + 2 j'
(5.28)
and, when applied to the bending stress resultant, eq. (5.13b), gives J t '~j = -m-
~
t~0) -~1)-~1) + ...
l)t~,~ 1) + \
1--i~--m }t~m) .
(5.29)
Thus, the integrals are reduced to evaluating the constitutive equations at a number of stations through the thickness. If the response is elastic, only two stations through the thickness are required without twisting and three stations through the thickness with twisting. If the response is inelastic, no more than five or six stations through the thickness are normally required.
5.5. Bending performance The evaluation of the strain energy in the bending of thick plates reveals an overestimation of the transverse shear contribution. This is corrected by scaling the transverse shear energy back by x~l ) = 4 / 5 (Cowper, [5]). The "correction" is related to the discrepancy between the constant distribution of shear strain predicted by the displacement assumptions and the more nearly parabolic true distribution of shear strain. In the limit of reducing thickness, a plate bending theory with transverse shear becomes arbitrarily stiff in transverse shear response and the transverse shear strains vanish. The result is a Kirchhoff bending theory. Computationally, the result is a 1 / h 2 growth in transverse shear strain energy and an assortment of computational ills generally referred to as shear locking, (Fried, Johnson, and Tessler, [8]). The benefits of C ° interpolation for the velocity vi and the rotational rate ~' disappear. However, if the transverse shear strain energy is viewed as a penalty function approach to approximately satisfying the thin plate limit, the 1 / h 2 growth may be modified and the associated numerical ills avoided. The work of Fried et al. [8] quantifies this approach and is used here. Fried et al., demonstrate conclusively that the minimal beam bending element with quadratic displacement convergence and the strain energy integrated exactly requires a quadratic interpolation of the velocity vi coupled to a linear interpolation of the rotation rate ~i. However, the constant membrane, constant bending stress approximation introduced by Belytschko, Lin, and Tsay [4] results in the same stiffness representation. Thus, the work of Fried et al., will be applied to the present bending element. If A is the area of the element, a transverse shear strain energy scaling of K~2) = 6h2/A provides, in the limiting case of thin plate behavior, quadratic displacement convergence to the Kirchhoff bending result without shear locking in the element. In effect, the transverse shear stiffness is converted to a penalty condition at the element level to enforce, in an approximate fashion, the vanishing of the transverse shear strain dzz ( = d , , ) and d y z ( = ds, ) implying the necessary equivalence between vz,x and a~y, and -Vz,y and o~x.
S. W. Key et al. / The impact modeling of structures
111
A joint shear correction factor r z, which implements switching between thick and thin behavior, is given by x z = m i n ( 0 . 8 / h 2, 6 / A ) h 2.
(5.3)
To implement the transverse shear correction, the transverse shear strain rates are scaled black by before being used in the constitutive evaluation and the resulting transverse shear stresses are scaled by K before the divergence calculation is made.
5.6. Orthogonal hourglass control The mean stress-mean strain rate formulation considers only the linear part of the velocity field. The remaining portion of the nodal velocity is the so called hourglass field. Excitation of these modes may lead to severe, unresisted mesh distortion. A method for isolating the hourglass modes so that they may be treated independently of the rigid body and uniform strain rate modes is required. This is accon~plished by developing an hourglass "gradient operator" connected with hourglass restoring forces. The linear velocity field on which the mean strain rates are based is given by vLIN = Uio -t- ( ~ -- nJ~k )( X k - Xko )Vij, _
1
1
.,,j = -¢ (.,,Bj - ( ,,m.n mco°-,,,A)n, - n~(x~ - x ~ ) ( , , . ° : c o ; ~ ' ) ) , COLIN=
(5.3a)
1 co,o+ 5(,<,x~)<~,,~/. .
The hourglass velocity and rotational rate fields v~ G and coHo may be defined by removing the linear portion of the velocity and rotational rate fields. Thus, v~ G = v , - v~'Y,
CO~G= co,- CO,iN,
(5.32)
or in terms of nodal velocities and nodal rotational rates, H G
.,z
1 j -,,,,-.,o~,--:(,<,-x~,)(,,,:/-(,..°n 1
./
m
n
co~A
J
)",),
(5.33)
,
coi",° = co,z - co,o.~, - --¢ (x, - x~,)co,.,Bl, where
Vio
=
1 I zviz,~,
1 Z ¢0i0= zcoiz2~,
and
x 0i = z1x zi-~z ~ •
(5.34)
Advantage has been taken of nk(x~ - x~,~t) being zero or near zero at the nodes. The hourglass velocity fields, eq. (5.33), are in the improper null space of the gradient operation. (For the rotational rate the gradient operator is B/, but for the velocity the operator is slightly more complex incorporating the average value of the rate of rotation.) Since the rotational rate normal to the quadrilateral co, = coin~ does not participate, the four degrees of freedom it represents may be considered in the improper null space of the mean gradient operator evaluation. The linear velocity and rotational rate fields, eq. (5.31), span 15 degrees of freedom: 3 rates of rigi d body translation, 8 uniform rates of the velocity gradient (the velocity gradient normal to the reference surface v,,, = v~.jnin ~ is not included), and 4 linear rates of the velocity gradient normal to the reference surface, which means that the hourglass subspace is the remaining 5 degrees of freedom - exclusive of the four degrees of freedom represented by co, = co~n(
S. W. Key et al. / The impact modeling of structures
112
Two of the hourglass modes are the familiar inplane " k e y s t o n e " or " b e n d i n g " modes such as the one pictured in fig. 7 for F3t. The third m o d e is an out-of-plane warping mode. The fourth and fifth hourglass modes are related to complex twists involving linear variations of the gradient of the rate of rotation coi. An hourglass gradient operator is constructed from the hourglass basis vector /~1 s as follows:
G I= (V/~)I~S,
(5.35)
where (~ is a generalized element dimension developed in connection with the stability criterion for the explicit central difference time integration (Flanagan, [6]), V is the volume of the element. This scaling provides the hourglass gradient operator with the same dimensional characteristic as the uniform gradient o p e r a t o r B/. The hourglass strain rates ~, and ~/i are developed with G 1 operating on only the hourglass velocity viiH G and rotational rate w,lHG,
q, = (1/V)v,~ oG',
//i = ( l / V ) ~0,,.o,-,lc,.
(5.36)
Alternatively, unrestricted operators m a y be developed by requiring them to satisfy the following condition: Vii71
/ n + ~inb~l U G .o ~l , = v,#
lois71 = ~0HGG/.
(5.37)
Using the hourglass velocity and rotational rate, eq. (5.33), provides
,,,7' +n/.
<0.71=[
(xS- X Zl) ¢ojBJ]G,,
(5.38)
which, when rearranged and using the orthogonality of the mode shapes 521 and Fss, [ I F 3 # = 0, gives 1 j v,,71 + B/.w 7 = v,, G s - -~xJG B)~\) + %,,,,n " ~,t~xjG'A " / 1 j J i ,j), wiD, z = w,l ( G I -
1
j
j
(5.39)
1~
The conditions for unrestricted operators, eq. (5.37), is satisfied if the hourglass operators 7 / and fl[, are defined as
71- (V/,~)(-r'3s- (1/V)xjl"3XB/),
flit = (1/,S)(,,.,.n"xJUA'nj),
(5.40)
from which the hourglass rates are c o m p u t e d as
4, = (1/V)(v, sT'+ ~°Tfl/,),
/1, = ( 1 / V ) o ~ , , 7 ' .
(5.41)
The hourglass strain rate /l~ selects the F / m o d e from ~0, = ~0,n*. T o control the hourglass modes, generalized forces Q~ and torques T s are defined which are conjugate to O~ and /l~, respectively, so that the work rate is given by V tell O.~n 1 __ i l J H G -I- I g n H G - -
VQ'dL +
VT.il".
(5.42)
Utilizing eq. (5.38), the contribution to the nodal forces due to hourglass resistance is given by
f~
=
Q'7',
g .'. o
, i + LT'. = Q/~,,
(5.43)
S. 14:. K e y et aL / The impact modeling o f structures
113
The hourglass restoring forces are calculated from ~i
= ~2/'ttan~
ij
•
•
n
Tn = C2/'ttanSnrn'O '
qj,
(5.44)
where 2~tan is the tangent shear stiffness obtained from the deviatoric constitutive behavior of the mean stress and mean strain state in the element and c is a scaling. The scaling c assures the level of the hourglass restoring forces remains below that of the mean deviatoric stress state. The deviatoric behavior is used since the hourglass modes are constant volume, higher order straining modes of the element. The invariant time derivative of the generalized forces Q+ and torques T, accounts for the finite rotations expected in the appfication of the quadrilateral element in analyzing large deformation, transient dynamic phenomena. The derivative is given by (5.45) where ~0ij is the spin, (vi, / - vj,i)/2. The hourglass restoring forces f ~ and torques g,Hcl are added to those obtained from the divergence of the mean stress state so that the complete result is
:"=(+'+"+e',/), 1~
m
-ij
1
""
i 1
(5.46)
6. Cylinder and frame modeling with plate elements To improve both numerical efficiency and the representation of the structural components, the cylinder and frame combination was remeshed using 4-node plate bending elements, fig. 3. The periodic spacing of the stiffening flanges along the cylinder and the plate girder of the frame permits a " u n i t " length of the cylinder and frame combination to be isolated with symmetry boundary conditions for the purposes of analysis, fig. 3. Since the thickness of the various structural components is carried as a parameter in the plate element, the individual elements appear without any thickness in perspective mesh plots, fig. 3. The spacing for the exterior, tapered flange stiffeners on the plate girder was several times that of the other flange stiffeners and, therefore, was dropped from the plate element model, fig. 3. To produce an isolated analysis of the cylinder and frame modeled with plate elements, the loads observed at the top of the cylinder in the solid element model, fig. 4, were used, along with the velocity history observed at the base of the solid element model, fig. 5. This produces comparable results as long as the structural models closely resemble each other in mass and stiffness distribution.
7. Comparison of modeling results There are four major features of the response on which a comparison can be based: the history of the relative displacement between the top and bottom of the cylinder, the effective stress history in the plate girder of the frame, the distribution of effective stress compared to the yield stress in the stiffeners on the cylinder at the time of maximum ovalization, and the deformed shape at the time of maximum ovalization. The history of the relative displacement between the top and bottom of the cylinder, fig. 8, shows an oscillation period in ovalization for the two structural models which are within a few percent of each other. The long wave length, low frequency response of the two models is essentially the same. The vertical load acting against the lower kinematic boundary in the plate element model shows a static ovaling deflection
114
S. 14( K e y et aL / The impact modeling of structures
....
I
. . . .
I ' ' ' ' 1 ' ' ' ' 1
....
I ....
I ....
I ....
1 ....
1''''
m~
....
0.0
I ....
P....
]
"0.2
l- J---l -l ~~Ki'k~l 1~ "
40000,0
jl ! /
I It
I
I - -
-ff~l~,~mm~
--'.,'l,lll@mlr~
I
I I
-12
~X
I
i
~ i I I
/
tt\]
I I
I X
I I I
/
~
I
t
/ /
I I I
r
I \x/ 0,0
-IJl 0.0
o~
Fig. 8. History of vertical diametral change in the cylinder in the full model based on solid elements and in the isolated model based on plate elements.
,,,I,,,,] Lo
l. 000[-01
=o~
=o.o
3s,o
~0,0
4513
50.0
5A 9 YIO 1. 000([-01 Z. (:X3OE-Ol 3.000£-01 q. O00E-OJ
6. O00E-Ol
E 5. O00E-01 F 6. O00E-OI
Fig. 10. The d i s t r i b u t i o n of effective stress in the m o d e l b a s e d on solid elements at the time of m a x i m u m d i a m e t r a l deflection in the cylinder (22 ms) expressed as a p r o p o r t i o n of the yield stress.
zs.o r,,. (,,~)
I,, ,~o
Fig. 9. History of vertical stress in the vertical plate girder of
A B C O
7,000£-0]
Itt
.... ,o~
2. O00~-OI 3. O00E-O] q. O00E-O]
5. O00~-Ol
\\1~/
the frame in the full m o d e l b a s e d o n solid e l e m e n t s a n d in the i s o l a t e d m o d e l b a s e d o n p l a t e elements.
5IG YLD R B C 0 E F G
if
X\
G 7. O00E-Ol
0 C ~0 C Fig. 11. The distribution of effective stress in the model based on plate elements at the time of m a x i m u m diametral deflection in the cylinder (22 ms) e x p r e s s e d as a p r o p o r t i o n of the yield stress.
S. I,E. Key et a L / The impact modeling of structures
115
r~
Fig. 12. Distorted cross section based on plate elements at the time of m a x i m u m diametral deflection in the cylinder (22 ms). Displacements are magnified by 5.
Llrldeformed
Idogrtificatlon=5
of 0.7 in. while the actual compliant boundaries in the solid element model appear to allow the cylinder and frame to either unload at late times or switch deformation modes. In either case insignificant deflections occur in the flange stiffened cylinder considering its 160 in. vertical inside diameter. The effective stress history in the plate girder of the frame shows a nominally similar behavior, fig. 9. Between 15 and 30 ms, there is an oscillation about an effective stress of 15 000 psi. The effective stress then rises to approximately 46 000 psi at 35 ms., after which it decays to 30 000 psi at 50 ms. The model based on plate elements exhibits an unloading of the frame between 35 and 40 ms at the same time the cylinder exhibits a second peak in its ovalization behavior. While this makes sense from a structural vibration point of view, a similar behavior is absent, or at least greatly diminished, in the model based on solid elements. Perhaps, in the more complete model based on solid elements, the compliance of the tunnel invert masks this behavior. The distribution of effective stress compared to the yield stress in the stiffeners on the cylinder shows the m a x i m u m stresses occur in the flange in compression at the crown at the time of m a x i m u m owli7ation, Figs. 10 and 11. There is an opposite but lower magnitude tensile stress in the flange at the midheight of the cylinder as would be expected in the vertical shortening phase of the ovalization. It would appear that the extreme fiber of the stiffening flange at the top of the cylinder may exhibit plastic yielding, but the constant membrane stress approximation in the plate elements does not allow it to be identified explicitly. The deformed shape at the time of m a x i m u m ovalization is shown in fig. 12. Fig. 12 clearly illustrates the deformation m o d e of the shell model compared to the undeformed shape. These deformed shapes are consistent with the diametral load on the cylinder, except perhaps the rolling of the upper flange on the plate girder. This is a bolted connection which was modeled as perfectly bonded. This result suggests a more detailed model of this connection should be constructed. The outcome could be a recommendation for a more compliant connection or the placement of the inner bolts midway between flanges.
116
S. 14( Key et aL / The impact modeling of structures
8. Evaluation of numerical efficiency Even though the plate element model of the cylinder a n d frame was analyzed as a n isolated structure using loads a n d k i n e m a t i c b o u n d a r y c o n d i t i o n s observed in the full model m a d e up exclusively of p l a n e strain solid elements, it is possible to estimate r e a s o n a b l y closely the savings to be gained from using plate elements in place of solid elements for the steel cylinder a n d frame. The increase in critical time step from 0.35 ms to 1.40 ms w o u l d represent a factor of 4 reduction in effort for the full model.
9. Conclusion Overall, the use of plate elements did not reveal a n y significantly different structural behavior in this analysis. F o r this structure, plate elements do offer the o p p o r t u n i t y to model local detail n o t afforded by the use of p l a n e strain solid elements. F r o m the s t a n d p o i n t of n u m e r i c a l efficiency, there is a factor of 4 r e d u c t i o n in c o m p u t a t i o n effort to be seen by using plate elements.
Acknowledgements S u p p o r t for the initial phase of this work b y the U S A F Ballistic Missile Office, N T S Engineering, a n d A g b a b i a n Associates is gratefully acknowledged. The design a n d effective p e r f o r m a n c e of the steel cylinder a n d frame exhibited in these calculations is the p r o d u c t of the skill a n d experience of H a r o l d F u t t r u p , A g b a b i a n Associates.
References [1] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Applied Mathematics Series, No. 55, Tenth Printing, National Bureau of Standards, United States Department of Commerce (1972). [2] T. Belytschko and B.J. Hsieh, Nonlinear transient finite element analysis with convected coordinates, Internat. J. for Numer. Meths. in Engrg. 7 (1973) 255-271. [3] T. Belytschko and A.H. Marchertas, Nonlinear finite element method for plates and its application to the dynamic response of reactor fuel subassemblies, Trans. ASME, J. of Press. Vessel Teclmol. (1974) 251-257. [4] T. Belytschko, J.I. Lin, and C.S. Tsay, Explicit algorithms for the nonlinear dynamics of shells, Computer Meths. in Appl. Mechs. Engrg. 42 (1984) 225-251. [5] G.R. Cowper, On the accuracy of Timoshenko's beam theory, J. of the Engineering Mechanics Division, EM6, Proceedings of ASCE, pp. 1447-1453 (1968). [6] D.P. Flanagan, Numerical Techniques with One-Point Integrated Hexahedron and Quadrilateral Elements, Ph.D. Dissertation, Northwestern University, Evanston, IL (1981). [7] D.P. Flanagan and T. Belytschko, A uniform strain hexahedron and quadrilateral with orthogonal hourglass control, Internat. J. for Numer. Meths. in Engrg. 17 (1981) 679-706. [8] I. Fried, A. Johnson, and A. Tessler, Minimal-degree thin triangular plate and shell bending finite elements of order two and four, Computer Meths. in Appl. Mechs. Engrg. 56 (1986) 283-307. [9] T.J.R. Hughes and W.K. Liu, Nonlinear finite element analysis of shells: Part I. Three-dimensions shells, Computer Meths. in Appl. Mechs. Engrg. 26 (1981) 331-362. [10] T.J.R. Hughes and W.K. Liu, Nonlinear finite element analysis of shells: Part II. Two-dimensional shells, Computer Meths. in Appl. Mechs. Engrg. 27 (1981) 167-181. [11] T.J.R. Hughes, W.K. Liu, and I. Levit, Nonlinear dynamic finite element analysis of shells, in: Nonlinear Finite Element Analysis in Structural Mechanics, Wunderlich, Stein, and Bathe (Eds.) (Springer-Verlag, 1981) pp. 151-168. [12] R.D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates, J. of Appl. Mechs. 18 (1951) 31-38.