- Email: [email protected]
Degenerated plate-shell elements with refined transverse shear strains
Computerr & Srruciures Vol 60, No 6. pp 1079-1091. 1996 Copyrlghl B 1996 Pubhshed by Elsevw Scmce Ltd Pnnted I” Great Bntam. All nghts reserved W45-7949/96 $I 5 00 + 0 00
Pergamon 0045-7949(95)00449-l
DEGENERATED PLATE-SHELL ELEMENTS WITH REFINED TRANSVERSE SHEAR STRAINS M. Ziyaeifar and A. E. Elwi University
of Alberta,
Edmonton,
Alberta,
Canada
(Receined 18 November 1994) Abstract-Although traditional three-dimensIona plate-shell elements relax the constraint so that normal cross-sectlons remam normal to the neutral plane during transverse shear deformation, the section is still constrained to remain plane. The work reported here relaxes this constraint by introducing shape functions
across the thickness to approximate the transverse shear stram field in the thickness direction. These shape functions are treated in the manner of generalized angles undergoing small deformations. They are added as new degrees of freedom to the ordinary displacement field of the degenerated shell elements. The field is still able to simulate large deformation behavior of the element. Each shape function yields two independent variables, one in each direction. In this work, two types of shape functions are proposed allowing a parabolic transverse shear strain, as well as an unsymmetric transverse shear strain distribution in the tluckness directlon. These two modes of deformation are particularly important in the case of diffused material failure in shell structures. Displacement field representation and finite element formulation based on a total Lagrangian approach are given. Examples are presented demonstrating the applicability of this element in a variety of problems. Copyright 0 1996 Elsevier Science Ltd
1. INTRODUCTION
Early plate-shell elements were usually based on some basic assumptions, such as the Love-Kirchhoff hypothesis in which the cross-section remains plane and perpendicular to the neutral axis after deformation. Such assumptions reduce a three-dimensional continuum medium to a two-dimensional surface that could be a plate or shell. Although this approach reduces the required computational effort to analyze the structure by orders of magnitude, it costs too much accuracy when the hypothesis is no longer applicable. In conventional shell elements the stress picture is close to that of a plane stress field. Whereas, normal stresses in the direction of the thickness are excluded from the formulation, transverse shear forces are included. At best, however, the shear strain distribution is assumed uniform across the thickness. Strain fields based on prehminary assumptions could be very suitable for a large variety of problems dealing with general behavior phenomena (such as deflections, vibration frequencies, etc.). Also they are adequate if the structure is thin and material behavior is elastic. However, where the thickness of the shell increases or the material response depends on the strain level in the element, a more accurate representation for the strain field is required. The first stage of refinement that is currently in common practice, is referred to as first-order shear deformation shell elements. These elements that are formulated using independent interpolation functions for displacements and rotations, are commonly
addressed as Mindlin-Reissner plate elements. This type of formulation relaxes one of the main assumptions of the Love-Kirchhoff hypothesis in bending; the cross-section remains plane, but is not restricted to be perpendicular to the neutral axis. Therefore, a constant shear strain distribution across the thickness will result if a transverse shear force exists in the plate [ 11. In the next phase of refinement, other assumptions in classic bending theory can be relaxed. The crosssection need not remain plane after deformation. This means that, the cross-section can warp to comply with the shear strain distribution across the thickness. To have a parabolic shear strain distribution, a third-order variation of the in-plane displacement is required, which is named as third-order shear deformation theory [2]. Moreover, nonlinear material responses remove symmetry of the shear strain distribution across the thickness. This 1s pronounced if the differences in the material behavior in tension and compression are large. Analysis of diffused failure in shells and plates, particularly in concrete panels, is a case in point. To deal with this problem, the transverse shear strain treatment should allow unsymmetric distributions. This could be done for example, by adding an anti-symmetric function to the parabolic shear strain distribution. Traditionally there are two different approaches in dealing with plate-shell finite elements, shell theory based elements and continuum based elements (the so-called degenerated elements). The relationship between the two approaches was recently discussed by
1079
1080
M. Ziyaeifar and A. E. Elwi
Buechter and Ramm [3]. Several higher order (refined) shell or plate elements were formulated based on shell theory. Among those, the work by Reddy [2,4,5,6] is well known. Reddy’s work, which incorporates the third-order shear deformation theory, is particularly directed to laminated plates. In these problems in order to have a continuous shear stress distribution across the thickness, considering different material properties in each layer, a discontinuous shear strain distribution is required. Along the same lines higher order shell theories were expanded recently to large displacement applications by incorporation of finite-rotation shell theory instead of the usual moderate-rotation model [7]. A continuum based layered refined shell element for laminated shell elements in which a static condensation technique was used to reduce the computational efforts was introduced by Owen and Li [8]. Other than this particular effort, not much work has been done to refine continuum based shell elements to include out of plane shear strain distribution. This paper introduces a continuum based degenerated plate-shell element that allows parabolic and unsymmetric transverse shear strain distributions across the thickness. The displacement field and the finite element formulation are presented. Examples of the application of the element are given and compared with the performance of traditional degenerated plateshell elements. The accuracy of the resulted strain field using refined elements is demonstrated in comparison with a plane stress analysis. 2. FORMULATION
OF THE DISPLACEMENT
FIELD
The main concept of degenerated shell [9-l l] based on a three-dimensional continuum representation of the strain field, in which an isoparametric interpolation imposes kinematic constraints similar to those in the first-order shear deformation theories (i.e. Mindlin-Reissner theory). In other words, independent isoparametric interpolation is carried out for translational and rotational degrees of freedom of each node. Figure 1 shows the deformed elevation of a degenerated plate-shell element. The transverse shear strain at the middle surface, /lo, is defined as
Tii Bo
a
Fig. 1. Deformation of the cross-section considering constant shear strain distribution.
thickness, which is the ideal distribution of the linear elastic shear strain, satisfying the equilibrium conditions of shear stresses at the top and bottom fibers. The counterpart normal displacement field, ug, of such a strain distribution is a third-order function described as
in which fi,, is the shear strain at the middle surface, z is the coordinate along the thickness and h is the thickness of the element. The resulting strain field in the thickness direction is described as
Lx(Z)= eg_!=~o1_J!z2~
( >
(3)
Figure 2b illustrates eqns (2) and (3). The total longitudinal displacement is written as
u(z)=uo+$z -ufl(z).
(4)
Here, u, is the displacement of the middle surface. This is illustrated in Fig. 2a. Substituting for the derivative from eqn (1) eqn (4) takes the form u(z) =
Here, CIis the degree of freedom of rotation of the cross-section, while the first term on the right hand side is the rotation of the normal to the middle surface. Both terms are not functions of z, the coordinate along the normal to the middle surface. The direct consequence is the contribution of a constant shear strain distribution across the thickness. To have a more accurate strain field in the element, the assumption that the cross-section remains plane in bending theory must be relaxed. This should result in a parabolic shear strain distribution across the
tz
dwldx
110 +
ctz +
poz-
up(z).
Substituting for the parabolic shear adjustment term, uB, from eqn (2) yields u(z)=u,+az
(5) strain
+po $z3
( >
U(Z) = u,, + u,(z) + u, _P(z)
(symbolically).
(6b)
The first two components are a translational and a rotational degree of freedom that are already implemented into ordinary degenerated shell elements by
Degenerated plate-shell elements
1081
(4
@I
Fig. 2. Cross-section deformation due to parabolic shear strain distribution across the thickness.
imposing a constraint in which a cross-section normal to the middle surface before deformation remains plane, but not normal after deformation. The third component relaxes the requirement that the crosssection remains plane. Therefore, two additonal degrees of freedom are added, one in each direction. Since the degrees of freedom that generate the u, and u, contributions already exist, the additional degrees
yo=o ya=$+=130
Y,+
of freedom, associated with each out of plane shear strain components can be denoted &, and &,, such that BOX * YXZ and
BOY * r,, .
Figure 3 shows the different components displacement field in one dimension.
%z
Fig. 3. Normal displacement and shear strain distribution across the thickness.
of the
1082
M. Ziyaelfar
and A. E. Elwi
Conversely, the third-order shear deformation theory proposed by Reddy [2], possesses a displacement field of the form
h
in which the additional degree of freedom in each direction with respect to the classic shell elements is only a. In the above representation, the third term imposes a strict constraint on the parabolic shear strain distribution. Therefore, it is always consistent with the parabolic shear stram distribution and satisfies zero shear strain requirements on the top and bottom of the element. Involving the derivative of the normal displacement degree of freedom, W, is numerically inconvenient in finite element formulation. In the degeneration approach, the derivative is eliminated at the expense of introducing additional degrees of freedom. This may be viewed as a potential source of problems as will be discussed below. On the basis of the field described in eqn (6), the structure has two different channels to move into, in addition to axial deformation. Its share from each of those, is inherited from the basic concept of finite element. In other words, the structure would choose a combination from both of those channels of deformations in which the resulting strain distribution satisfies equilibrium and simulates the actual strain field. In the presence of transverse shear forces, the third displacement component in eqn (6) is directly proportional to the shear strain at the neutral axis, pi,. This component and its associated shear strain distribution are shown in Fig. 3. The second component of displacement in eqn (6) is representative of the rotational degree of freedom. Its derivative in the longitudinal direction provides the normal flexural stresses associated with bending moments. It also provides a constant shear strain contribution. In the following, the second and third terms of eqn (6) will be referred to as shape functions, albeit, across the thickness. The shape of the full shear strain distribution is a combination of the two contributions, and is depicted in Fig. 3. It is noted here that there is no strict constraint that imposes a perfect parabola shape on the shear strain distribution across the thickness, similar to that in the third-order shear deformation theory. Since each component of longitudinal displacement across the thickness having a contribution to the out of plane shear strain distribution is a totally independent degree of freedom, there is no assurance of a zero shear strain at the top and bottom fibers. Everything is relegated to the self adjustment capability of the finite element these new components formulation. Moreover, contribute to the longitudinal normal strains as well as the transverse shear strains. This may cause a discrepancy in the results in the case of discontinuity in moment and/or shear force in the structure. This issue will be discussed later.
,
uy (z> Fig. 4. Unsymmetric
r, (4
shear strain and its normal displacement distribution.
In the case of anisotropy, elasto-plastic behavior or material models where the behavior in tension is different from compression, unsymmetric shear strain distribution across the thickness may follow. These strain fields can be better represented by introducing more sophisticated shape functions in the same way as above. In this work a new shape function is added to the longitudinal displacement distribution across the thickness, eqn (6). This function must comply with zero transverse shear strain conditions on the top and bottom surface of the shell. Figure 4 is an example of the desired function and its contribution to the transverse shear strain distribution across the thickness that is anti-symmetric. This new shape function is constructed from spline functions to comply with the requirement of zero slope on top, bottom and middle surface:
iz2-$z’sign(z)
1 @aI (gb)
in which y0 represents the maximum transverse shear strain contribution at one fourth of the height of the cross-section. This shape function has no contribution to the flexural deformation of the cross-section. On the other hand, close attention must be paid to its contribution to the axial deformation of the element, which is normally represented only by the translational degree of freedom uO. This contribution is particularly pronounced at the discontinuity locations of axial forces such as point loads and support reactions. This problem will be addressed later. It should be noted that for each across-thethickness shape function, one degree of freedom is added in each direction at each node. Thus, for the parabolic and antisymmetric shear distributions, two additional degrees of freedom in each dtrection must be added to the ordinary five degrees of freedom at each node. The increase in computational effort is thus substantial. The use of the additional degrees of freedom should be based on the specific characteristics
1083
Degenerated plate-shell elements of each problem. If there is no possibility of having an unsymmetric behavior in the shell thickness, it is better to suppress or remove these degrees of freedom from the formulation. It is convenient to give the new displacement shape functions, the third component of eqn (6) and that of eqn @a), the same extreme value as that of the rotational degree of freedom contribution, the second term of eqn (6), at the top and bottom of the crosssection. This has no effect on our formulation because the additional components will be treated as shape functions. The intention is to provide the ability to compare the relative contributions of each component. Introducing t as the local coordinate in the thickness direction, the final equations for shape functions across thickness are as follows: u(t) = ug + U,(f) + u,_,(t)
+ u;.(t)
Fig. 5. Definition of rotational degrees of freedom. ence. The rotational degrees of freedom were defined as changes in the angles, $, the normal to the middle surface described with respect to the global frame of reference. Thus, “‘tl, and “‘~1~ are defined as
t = 2z/h u,(t) = ahtj2
“‘al = (“‘tj, - ‘$])
and
ma2= (“tiz - 0tj2), (10)
u1 -p(t) =&f(t) u,(t) = ?“g(t) f’(t) = ht’/2 g(t) = 3ht2/2 - ht3 sin(t).
(9)
Here g(t) and f(t) are shape functions. The concept of shape function that was used in the above representation can be easily expanded to the other type of required refinements in the shells, particularly in multi-layered shell elements. This can be achieved by introducing a piecewise smooth shape function across the thickness as a broken line. This function can accommodate the discontinuity of transverse shear strain across the thickness to provide the required continuity of shear force between layers. The substitution of such a function instead off(t) or g(t) does not need any special treatment in the formulation. 3. GEOMETRIC
TRANSFORMATIONS
The model developed above was implemented into Program NISA-80 that was developed at the University of Stuttgart by Stegmuller et al. [12]. This program was written in FORTRAN77 and contains the three-dimensional degenerated plate-shell elements developed by Ramm [ 10, 111. In expanding the formulation of the degenerated shell element, the notation used by Ramm [lo, 111,is used with minor changes. In the original element the degrees of freedom at each node were three translations and two rotations: %, , m~2,“‘u3, “‘a, and “‘u2. The translations take place at the middle surface and are defined with respect to the global axes of refer-
in which m denotes the current configuration and o denotes the original undeformed configuration. Figure 5 shows the definition of these angles. The four new degrees of freedom, “/3,, “f12, “y, and “‘yzintroduced in this work are naturally similar to, and are modeled after, the rotational degrees of freedom m~l and m~2. There are, however, two main differences. The new degrees of freedom have zero values in the original configuration. More importantly, since they have the magnitude of shear strain values they are not likely to have extremely large values. Therefore, they are treated as small displacement changes. This is in contrast to the contribution of formal changes in the angles of normal “$, and m$z (or m~l and Yx2), which existed in the original formulation, and can be as large as necessary. 4. FINITE ELEMENT DISPLACEMENT
FIELD
Ramm [I l] describes the mapping of a vector of length L along the normal to the middle surface in the current configuration that undergoes small changes in the angles “‘Ic/,and “‘tj2. The magnitudes of these small changes are CI,and a,. The quantities c(, and a2 are thus the incremental changes in angles “$I and “‘$*. The translations of the tip of the vector relative to the global frame of reference are thus A”x, A”x, 1 A”x, i 0 =
-sin(“$,) sin(“&) sin(“$,) ~os(~+~)
-sin(“ll/,) cos(?j,) cos(m+2) z: L cos(m$,) sin(“$,) 111 (1 la)
M. Ziyaeifar and A. E. Elwi
1084
{Amx} =[F(“$,,“‘$,]{cc}
(symbolically).
(llb)
If the changes in the original degrees of freedom of rotation are small from one step to the next, the form of eqn (11 b) can be used to obtain the corresponding contribution of the new degrees of freedom to the displacement field. Considering the above, the full displacement field of a point on the cross-section of an element at position t along the normal to the middle surface is described as
followed and material nonlinearities are allowed. Assuming that the response is linear in the increment between configurations m and m + 1, the principle is written as:
s
e,,Ct,kdekl dV +
‘“St,h, d V
Y
=
s s
m+ ‘f&,
m+ ‘fkduk dA +
A
-
“‘S,,e,, d V.
”
+ f
cwv, I- cos(“kh)
1
sin(?j,) cos(m+2) - sin(O$,) cos(Oi+$) sin(“tj,) sin(“&) - sin(“+,) sin(“&) 1
+f(t) [wk,
)
V2)1
1:;II 2
+ g(t)[wv,
7 “ti2)l
;I {
.
dV
sY
(14)
Here e is the linear part of the increment of Green-Lagrange strain tensor, rl is the nonlinear part, S is the second Piola-Kirchhoff stress tensor, C is the constitutive tensor, t is the vector of surface tractions and f is the vector of body forces per unit volume. The tensors e and 1 are defined in terms of the displacement derivatives and its incremental values as
(12)
I
To obtain the increment of displacement field, one notes that the corresponding increments in the rotational degrees of freedom are also small enough to be represented by a transformation based on eqn (11 b). If the same interpolation shape functions are used over the middle surface of the element, the increment of displacement field is written in terms of the nodal degrees of freedom as
In eqn (15) at a point defined by I, q and t in the element the strain tensor, e and increment of displacement gradient, u,,,, are written in matrix form in terms of the nodal displacement increments as
{e> =n&l[m~Lln{~n~
UW
{Vu1= ,;, [m4dn{~n~~
(16t.4
Substituting in eqn (14) and recognizing the arbitrary variation of the displacement increment, the equation of equilibrium is recovered as [m~L+mK~J(~} ={“‘+‘I? -mQ}
(13) Here, /IIt /12,y, and y2 are increments in the degrees of freedom “‘/l,, mjJ2, my, and “‘y2, respectively. The node number n varies from 1 to the total number of nodes in an element, I and [N], is the interpolation submatrix at node n. 5. FINITE
ELEMENT
(17)
in which the first bracket on the left hand side represents the linear (tangent) and geometric stiffness matrices, while the right hand side represents the current level of loads and the internal equilibrating forces. The stiffness matrices are recovered as
V’KJ =
[“&1T[C1[m4_1 dK sY
FORMULATION
In the following the element formulation is presented in the context of the incremental principle of virtual work. A total Lagrangian approach is
In the following, the forms of the matrices on the right hand side of eqns (16) and (18) are defined.
1085
Degenerated plate-shell elements 6. DERIVATION
OF THE STRAIN MATRIX
Next the transformation reference takes the form
DISPLACEMENT
The third and fourth terms of eqn (15) involve products of the increment of displacement gradient and the displacement gradient tensor in configuration”. Thus, the final representation of strain-displacement derivatives can be written as @I = rmfflPI
(19)
in which [“‘HI contains a template of the current displacement gradient, while {Vu} is the gradient of the increment of displacement. Equation (19) is shown in full in Fig. 6. Components of the gradient of the increment of displacement vector are obtained using the form of eqn (16b). Thus, the first step is the calculation of displacement derivatives matrix [I&]. Displacement derivatives must be first interpolated from nodal degrees of freedom in the natural local coordinate system (5, q, t), and then transformed to the global coordinates. For example, in order to evaluate the gradient of the increment of displacement vector, one must first evaluate the gradient with respect to the local coordinates:
i=l,2,3
and
{=l,q,t
to the global frame of
in which [J]-’ is the inverse of the Jacobean of transformation matrix. Combining eqns (20) and (21) and comparing to eqn (16b), the displacement gradient is obtained in the form
Pi
= 1 [Jl-'["N'(I3
9,
t>ln{U>,
“=I
=,i,
[m4-JL1n~~~“~
The vector on the right hand side represents the displacement degrees of freedom of node n, while [“‘&In stands for node n submatrix, and is shown in Fig. 7. In this figure R’; are the components of the transformation matrix defined in eqn (1 l), while the shape function derivatives are defined as
(20)
4;,, = [g(t)d(, + g’(t).G’@l. in which the subscript n denotes the submatrices (vectors) of one node in the element. Here the matrix [N’] holds the shape function derivatives with respect to the natural local coordinates of node n evaluated at the current integration point with coordinates {, r~ and t. The vector {u}, holds the displacement increments at node n (remember we have here nine degrees of freedom per node).
ell
(22)
(23)
The strain displacement matrix operator [mBL]nis obtained from premultiplication of matrix [“B& by [‘“HI and the strain increment of eqn (16a) is finally obtained as @,,I =
1
e22 en 2e12 2e13
2ez3
Fig. 6. StrainAisplacement derivatives relationship.
1086
Fig. 9. Displacement derivatives relationship for large displacements.
1087
Degenerated plate-shell elements
The final result is shown in Fig. 8. The following expressions were used for the sake of simplicity:
0% = (1 + “uz,,)&, + m4.*F;, +
mU3.2F;,
073
“%3G
=
(1
+
“u3,3)%
+
%,34,
+
D’i, = (1 + mu3,3)F;Z+ mu,.3F;z +
(25)
m~2,3F;Z.
The gradient of the total displacement, yet to be determined. The procedure is that used for the increment of displacement eqn (22), except that here, possible large rotations have to be considered
8. BOUNDARY
{“y:,) has similar to gradient, values for
The dimensions of matrix [“#‘I and vector {“zi} are slightly different from those used with their counterparts used in forming the gradient of the displacement increment, because of the large rotational characteristics of the total displacements. In such a case, due to complicated trigonometric relationship between the two rotational degrees of freedom $, and &, it is convenient to use three equations for the two degrees of freedom, instead of going for two separate complicated expressions. A detailed description of eqn (26) is shown in Fig. 9. 7. STRESS STATE MATRIX AND VECTOR AND (“9)
(“S]
There are only six components of stress state all together, but to comply with matrix multiplication requirements, they have to be manipulated into a square matrix form
CONDITIONS
Boundary conditions are constraints on certain force or displacement components. Simple constraints take the form of prescribed displacements or forces. The new degrees of freedom introduced herein cannot have nonzero prescribed displacements or forces. The associated nodal forces for the refined degrees of freedom are always zero, regardless of the loading condition on the structure. The major concern herein, is to determine for different cases whether to suppress these degrees of freedom at certain locations or not. It must be noted once again that the new degrees of freedom interact with ordinary degrees of freedom. The parabolic shear degrees of freedom & and fiZ not only contribute to the shear strain distributions, but also to rotations and, thus, bending moments. Therefore, an interaction with the rotational degrees of freedom tlr and c(~would be expected. Similarly, the unsymmetric shear degrees of freedom y, and yZ interact with the axial deformations. As a starting point it would be necessary to switch off the interaction phenomena at certain points to avoid undesirable strain fields. The parabolic shear strain degrees of freedom are rotational degrees of freedom. Thus, if a fixed rotational boundary condition is required, these degrees of freedom should be suppressed along with ordinary rotational degrees of freedom. The same reasoning, however, cannot be used for the unsymmetric shear strain degrees of freedom. Although the latter degrees of freedom can represent axial deformations as well as unsymmetric shear strain distributions, the middle surface of the element is not affected by these degrees of freedom. Thus, there is no need to suppress these degrees of freedom where the axial in-plane degrees of freedom are suppressed. Where a concentrated load is applied in the in-plane direction, the unsymmetric shear strain degrees of freedom, together with ordinary m-plane degrees of freedom, can represent the pinching effect on the middle surface of the cross section. Close attention I
m& 0
[“S] =
mS22
0
0
mS*9
mS,*
0
0
“Szz
0
mS*,
0
0
mS,,
mS,3
0
0
Y213
0
0
0
“S,,
0
0
0
“S,,
0
0
T,j
0
0
mS32 0
0
0
mS,,
0
0
mS23 0
The stress state vector (“j),
takes the traditional form
(“Q
“SlZ
= Y%
“SZZ
mS33
mS13 mS,,>.
(28)
(27) m&s33
mS33 0
v**
1
must be paid to this effect at point load locations. With no constraints on the y degrees of freedom, a concentrated load that is assumed to be applied at the
1088
M. Ziyaeifar and A. E. Elwi
middle surface is distributed nonuniformly on the cross-section. Therefore, the y degrees of freedom must be restricted where a concentrated load on the in-plane direction is expected on the structure. Support reactions with components in the in-plane direction are an example of this case. A node at which a concentrated load is applied with a component in the in-plane direction, is a potential candidate for restricting these degrees of freedom. On the other hand, suppressing this degree of freedom removes one of the most desirable aspects of refinement in the shell element. Therefore, it must be done only at the nodes at which a high magnitude of concentrated in-plane force is expected. The same phenomenon arises in the case of concentrated moments. In this situation there is strong interaction between the tl and B degrees of freedom. Therefore, the contribution of /3 to the final strain distribution across the thickness may not be proportional to the actual shear force. In this case the strain field is no longer acceptable. A remedy for this problem is to suppress the parabolic shear strain distribution degrees of freedom at the point of concentrated moment. Again, one should keep in mind that this solution may be advisable only in the presence of a high concentrated moment, because it eliminates the contribution of refinement on parabolic shear strain distribution. 9. EXAMPLES
In order to verify the performance of the refined element, a number of problems were solved and the solutions compared to those obtained from the original degenerated plate-shell element, hereinafter called the “ordinary element”, and other solutions. The examples include a shell structure, two square plates and two beams. In the following the examples are presented and the results discussed. The first example is a shallow spherical shell supported on four corners and subjected to a uniform normal pressure loading. The material property is assumed to be isotropic linear elastic. The consequences of using the refined element are manifested in the flexibility of the structure and the stress
0.1
1
“I
0
CenbtqOpoint, vertical
distribution. Both effects are thickness dependent. Therefore, this example has been solved with two thickness units (t = 200 and 2000 mm). The extreme thickness values have been used deliberately to bring out an exaggerated response that can be seen clearly. Figure 10 shows the load vs vertical deflection response of point B for both the ordinary element and the refined element in the case of the thin shell (t = 200 mm). The two solutions are almost identical. Figure 11 illustrates the normal and shear stress distributions at point A that is close to a support. Note that the difference in the normal stress distribution between the solutions of the ordinary element and the refined element is clear and, in contrast to the shear stress distribution, is highly thickness dependent as is shown in Fig. 12 for the case of the thickness shell structure (t = 2000 mm). This is expected, since the plane sections do not remain plane in the presence of high shear gradients in the refined element. The ratio of the shear stresses to the normal stresses in the thick shell at this location is much higher than that in the thin shell at the same location. The next example consists of a square plate supported at four corners and subjected to a uniform pressure. The material has a yield strength of 240 MPa and is modeled with an elastic perfectly plastic
shear stress MPa)
F
% 8 1 .z CI
50
4 0 1
1
-50
-100
(mm)
Fig. 10. Thin shell, vertical deflection of central point.
100
3 3 a ._
Al!
300
-50
.v F -100
Normal stress distribution Fig. 11.Thin shell, stress distribution at a point close to support.
1089
Degenerated plate-shell elements Onset of yield at extreme fibers
Onset of yield at middle of section
/ 0
30
10
40
50
Deflec2on (mm)
Normal
stress dist
@IF%)
Fig. 12. Thick shell, nortnal stress dist. close to support.
von Mises constitutive relation. The response of the refined element shows a different progression of yield in the structure
from that of the ordinary
element.
The refined element predicts yield in the structure earlier and also in a different mode compared to the oridnary element. The onset of yield takes place at the middle of the cross-section close to the support (point A in Fig. 13a) due to the high shear stresses at that point (Fig. 13b). This limits the load carrying capacity of the structures and constitutes a “punching” shear failure. On the other hand, when ordinary elements are used, the onset of yield is postponed to a higher level of load and it starts at the extreme fibers of the plate cross-section due to a combination of out-of-plane shear stresses and in-plane normal and shear stresses. The structure then exhibits a gradually softening response. Punching shear may not be detected if a different mode of failure takes place elsewhere such as a mid-span or at the free edge. The load-deflection curve is shown in Fig. 14. The next example shows the case of unsymmetric material behavior across the thickness. This is nrovided by applying a biaxial inplane normal force to a square plate supported at the central point. This represents the region of a flat slab around a column. The material model used here is the same as that of the previous model. In this case, the onset of yield in the top and bottom fibers of the plate, takes place at
Fig. 14. Plate supported at corners, vertical deflection at a point close to support.
different levels of transverse pressure. Figure 15 shows the normal and shear stress distributions at a point close to support that clearly depicts the effect of the unsymmetric component of shear strain across the thickness. As mentioned before, both the parabolic shear strain components, p, and the rotational degrees of freedom, u, contribute to the shear and flexure deformation modes simultaneously. At a discontinuity in the transverse shear forces a combination of both a and p could yield a non-zero transverse shear strain at the top and bottom of the plate surface. However, the integral of shear stresses across the thickness is always in total balance with the actual shear force on the cross section. Under such circumstances the Co continuity of the interpolation functions for refined degrees of freedom, /?, between elements (and higher orders of continuity for the nodes inside the element) exists. This does not allow the finite element model to imitate the discontinuity in the shear strain properly. In this situation a combination of the degrees of freedom a and /I gives a gradual, but steep transition of shear strains over the discontinuity region. If a three-dimensional model were to be used with the same constraints of plate or shell, the result of the strain field would be similar to that of the refined shell element. This phenomenon is illustrated in the last example.
250
Fig. 13. Plate supported at corners, stress dist. at a point close to support.
M. Ziyaeifar and A. E. Elwi
1090
Normal stress (MPa)
Fig. 15. Plate supported at centre, stress dist. at a pomt close to support
In this example, a beam supported at both ends and subjected to concentrated loads is analyzed. The beam is fairly thick to render the shear deformation mode dominant relative to the flexural mode. Two extreme cases in which a sharp discontinuity in the shear force exists in the beam are analyzed, one with two point loads applied at l/4 points of the span and one with a single mid-span load. These beams are shown in Fig. 16. Boundary conditions of plane stress are assumed in the width direction when refined
FEM model, 256 4 node Plane stress element I +
+
elements are used (12 bicubic elements). The same beams were constructed from 256 two-dimensional isoparametric eight-node plane stress elements also shown on the same figure. The large number of elements is necessary to examine the stress points as close as possible to the plane of force discontinuity. Figure 17a shows the results for the mid-span load case compared to the plane stress solution at different points in the vicinity of the discontinuity. At points close to the discontinuity (L/290), the shear stress distribution from both the plane stress and the refined element solutions are rather flat. A short distance away from the discontinuity (L/14) both solutions become much more parabolic. Figure 17b shows the two point load solutions for both types of analysis. The points at L/290 and L/14 are in the pure moment region. However, because of localized effects both solutions give a perturbation in the shear stress that is evident at L/290 and is negligible at L/14 from the load point in this region. The solutions at a distance of L/40 from the load point measured towards the support (nonzero shear force zone) are again almost identical and parabolic.
4
A7
10. CONCLUSIONS
FEM model, 12 bicubic “refined” elements Fig. 16. Beams loaded at mid-span and at l/4 of span.
A modified three-dimensional plate shell element that has a better transverse shear strain representation
-1
(a) Load case 1
(b) Load case 2
Fig. 17. Shear stress close to points of shear force discontinuity.
Degenerated
piate-she11 elements
presented. In this element two types of degrees of freedom are introduced, a symmetric set that governs a parabolic transverse shear strain distribution and a antisymmetric set that governs an unsymmetric contribution. The new degrees of freedom are treated as shape functions across the thickness and added in the manner of small rotations to the longitudinal displacement field representation. This effectively removes the assumption that plane sections remain plane and allows two new modes of deformation. The performance of the new element has been tested and proven to give a more realistic response in the case of elasto-plastic materials. The implications of removal of the plane section restriction are discussed at length and corresponding results are shown to be consistent, even in the presence of strong discontinuities m the transverse shear force distribution. has
been
Acknowledgments-This work has been carried out m the course of doctoral graduate research of the first author. The work has been funded in part by the Natural Sciences and Engineerma Research Council of Canada. The authors also wish to acknowledge the support of the Iranian Ministry of Culture and Higher Education. The first author wishes to thank Dr M Farshad for his motivating ideas.
51, 745-752
(1984).
1091
3. N. Buechter and E. Ramm, Shell theory versus Degeneration-a comparison in large rotation finite element analysis. Int. J. numer. Meth. Engng 34, 39-59 (1992). 4. J. N. Reddy, A refined nonlinear theory of plates with reansverse shear deformation. Inr. J. Solid Struct. 20, 881-896 (1984). rotation 5. J. N. Reddy, A small strain and moderate theory of elastic anisotropic plates. J. appi Mech. 54, 623-626 (1987). models of com6. J. N. Reddy, On refined computational posite laminates Int. J. num. Merh. Engng 27, 361-382 (1989). 7. Y. Basar, Y. Ding and R. Schultz, Refined sheardeformation models for composite lammates with finite rotations. Int. J. Solids Siruct. 30, 2611-2638 (1993). 8. D. R. J, Owen and 2. H. LI, A refined analysts of laminated plates by finite element displacement methods -I. Fundamentals and static analysis, 11. Vibration and Stability. Compur. Strucr. 26, 907-923 (1987). 9. S. Ahmad, B. M. Irons and 0. C Zienkiewtcz, Curved thick shell and membrane elements with particular reference to axi-symmetric problems, In Proc. ad Cony. Malrix Methods in Structural Mechanrcs, WrightPatterson Au Force Base. Ohio (1968) ‘nichtlineare ‘Elstostatik und 10. E. Ramm, Geometrisch finite Elemente. Habilitation, Berich Nr. 76-2, lnstitut fur Baustatik, Universitat Stuttgart, Germany (1976). 11. E. Ramm, A plate/hell element for large deflections and rotations. US-Germany Symp. on Formulations and Computational Algorithms in Finite Element Analysa, MIT 1976, MIT Press (1977). 12. H. Stegmuller, L. Hafner, E. Ramm and J. M. Satteie. Theoretische grundlagen rum FE-programm system NISA-80. Mitteilung Nr. I, Instttut fur Baustattk der Univetsitat Stuttgart, Germany (1983).
Pergamon 0045-7949(95)00449-l
DEGENERATED PLATE-SHELL ELEMENTS WITH REFINED TRANSVERSE SHEAR STRAINS M. Ziyaeifar and A. E. Elwi University
of Alberta,
Edmonton,
Alberta,
Canada
(Receined 18 November 1994) Abstract-Although traditional three-dimensIona plate-shell elements relax the constraint so that normal cross-sectlons remam normal to the neutral plane during transverse shear deformation, the section is still constrained to remain plane. The work reported here relaxes this constraint by introducing shape functions
across the thickness to approximate the transverse shear stram field in the thickness direction. These shape functions are treated in the manner of generalized angles undergoing small deformations. They are added as new degrees of freedom to the ordinary displacement field of the degenerated shell elements. The field is still able to simulate large deformation behavior of the element. Each shape function yields two independent variables, one in each direction. In this work, two types of shape functions are proposed allowing a parabolic transverse shear strain, as well as an unsymmetric transverse shear strain distribution in the tluckness directlon. These two modes of deformation are particularly important in the case of diffused material failure in shell structures. Displacement field representation and finite element formulation based on a total Lagrangian approach are given. Examples are presented demonstrating the applicability of this element in a variety of problems. Copyright 0 1996 Elsevier Science Ltd
1. INTRODUCTION
Early plate-shell elements were usually based on some basic assumptions, such as the Love-Kirchhoff hypothesis in which the cross-section remains plane and perpendicular to the neutral axis after deformation. Such assumptions reduce a three-dimensional continuum medium to a two-dimensional surface that could be a plate or shell. Although this approach reduces the required computational effort to analyze the structure by orders of magnitude, it costs too much accuracy when the hypothesis is no longer applicable. In conventional shell elements the stress picture is close to that of a plane stress field. Whereas, normal stresses in the direction of the thickness are excluded from the formulation, transverse shear forces are included. At best, however, the shear strain distribution is assumed uniform across the thickness. Strain fields based on prehminary assumptions could be very suitable for a large variety of problems dealing with general behavior phenomena (such as deflections, vibration frequencies, etc.). Also they are adequate if the structure is thin and material behavior is elastic. However, where the thickness of the shell increases or the material response depends on the strain level in the element, a more accurate representation for the strain field is required. The first stage of refinement that is currently in common practice, is referred to as first-order shear deformation shell elements. These elements that are formulated using independent interpolation functions for displacements and rotations, are commonly
addressed as Mindlin-Reissner plate elements. This type of formulation relaxes one of the main assumptions of the Love-Kirchhoff hypothesis in bending; the cross-section remains plane, but is not restricted to be perpendicular to the neutral axis. Therefore, a constant shear strain distribution across the thickness will result if a transverse shear force exists in the plate [ 11. In the next phase of refinement, other assumptions in classic bending theory can be relaxed. The crosssection need not remain plane after deformation. This means that, the cross-section can warp to comply with the shear strain distribution across the thickness. To have a parabolic shear strain distribution, a third-order variation of the in-plane displacement is required, which is named as third-order shear deformation theory [2]. Moreover, nonlinear material responses remove symmetry of the shear strain distribution across the thickness. This 1s pronounced if the differences in the material behavior in tension and compression are large. Analysis of diffused failure in shells and plates, particularly in concrete panels, is a case in point. To deal with this problem, the transverse shear strain treatment should allow unsymmetric distributions. This could be done for example, by adding an anti-symmetric function to the parabolic shear strain distribution. Traditionally there are two different approaches in dealing with plate-shell finite elements, shell theory based elements and continuum based elements (the so-called degenerated elements). The relationship between the two approaches was recently discussed by
1079
1080
M. Ziyaeifar and A. E. Elwi
Buechter and Ramm [3]. Several higher order (refined) shell or plate elements were formulated based on shell theory. Among those, the work by Reddy [2,4,5,6] is well known. Reddy’s work, which incorporates the third-order shear deformation theory, is particularly directed to laminated plates. In these problems in order to have a continuous shear stress distribution across the thickness, considering different material properties in each layer, a discontinuous shear strain distribution is required. Along the same lines higher order shell theories were expanded recently to large displacement applications by incorporation of finite-rotation shell theory instead of the usual moderate-rotation model [7]. A continuum based layered refined shell element for laminated shell elements in which a static condensation technique was used to reduce the computational efforts was introduced by Owen and Li [8]. Other than this particular effort, not much work has been done to refine continuum based shell elements to include out of plane shear strain distribution. This paper introduces a continuum based degenerated plate-shell element that allows parabolic and unsymmetric transverse shear strain distributions across the thickness. The displacement field and the finite element formulation are presented. Examples of the application of the element are given and compared with the performance of traditional degenerated plateshell elements. The accuracy of the resulted strain field using refined elements is demonstrated in comparison with a plane stress analysis. 2. FORMULATION
OF THE DISPLACEMENT
FIELD
The main concept of degenerated shell [9-l l] based on a three-dimensional continuum representation of the strain field, in which an isoparametric interpolation imposes kinematic constraints similar to those in the first-order shear deformation theories (i.e. Mindlin-Reissner theory). In other words, independent isoparametric interpolation is carried out for translational and rotational degrees of freedom of each node. Figure 1 shows the deformed elevation of a degenerated plate-shell element. The transverse shear strain at the middle surface, /lo, is defined as
Tii Bo
a
Fig. 1. Deformation of the cross-section considering constant shear strain distribution.
thickness, which is the ideal distribution of the linear elastic shear strain, satisfying the equilibrium conditions of shear stresses at the top and bottom fibers. The counterpart normal displacement field, ug, of such a strain distribution is a third-order function described as
in which fi,, is the shear strain at the middle surface, z is the coordinate along the thickness and h is the thickness of the element. The resulting strain field in the thickness direction is described as
Lx(Z)= eg_!=~o1_J!z2~
( >
(3)
Figure 2b illustrates eqns (2) and (3). The total longitudinal displacement is written as
u(z)=uo+$z -ufl(z).
(4)
Here, u, is the displacement of the middle surface. This is illustrated in Fig. 2a. Substituting for the derivative from eqn (1) eqn (4) takes the form u(z) =
Here, CIis the degree of freedom of rotation of the cross-section, while the first term on the right hand side is the rotation of the normal to the middle surface. Both terms are not functions of z, the coordinate along the normal to the middle surface. The direct consequence is the contribution of a constant shear strain distribution across the thickness. To have a more accurate strain field in the element, the assumption that the cross-section remains plane in bending theory must be relaxed. This should result in a parabolic shear strain distribution across the
tz
dwldx
110 +
ctz +
poz-
up(z).
Substituting for the parabolic shear adjustment term, uB, from eqn (2) yields u(z)=u,+az
(5) strain
+po $z3
( >
U(Z) = u,, + u,(z) + u, _P(z)
(symbolically).
(6b)
The first two components are a translational and a rotational degree of freedom that are already implemented into ordinary degenerated shell elements by
Degenerated plate-shell elements
1081
(4
@I
Fig. 2. Cross-section deformation due to parabolic shear strain distribution across the thickness.
imposing a constraint in which a cross-section normal to the middle surface before deformation remains plane, but not normal after deformation. The third component relaxes the requirement that the crosssection remains plane. Therefore, two additonal degrees of freedom are added, one in each direction. Since the degrees of freedom that generate the u, and u, contributions already exist, the additional degrees
yo=o ya=$+=130
Y,+
of freedom, associated with each out of plane shear strain components can be denoted &, and &,, such that BOX * YXZ and
BOY * r,, .
Figure 3 shows the different components displacement field in one dimension.
%z
Fig. 3. Normal displacement and shear strain distribution across the thickness.
of the
1082
M. Ziyaelfar
and A. E. Elwi
Conversely, the third-order shear deformation theory proposed by Reddy [2], possesses a displacement field of the form
h
in which the additional degree of freedom in each direction with respect to the classic shell elements is only a. In the above representation, the third term imposes a strict constraint on the parabolic shear strain distribution. Therefore, it is always consistent with the parabolic shear stram distribution and satisfies zero shear strain requirements on the top and bottom of the element. Involving the derivative of the normal displacement degree of freedom, W, is numerically inconvenient in finite element formulation. In the degeneration approach, the derivative is eliminated at the expense of introducing additional degrees of freedom. This may be viewed as a potential source of problems as will be discussed below. On the basis of the field described in eqn (6), the structure has two different channels to move into, in addition to axial deformation. Its share from each of those, is inherited from the basic concept of finite element. In other words, the structure would choose a combination from both of those channels of deformations in which the resulting strain distribution satisfies equilibrium and simulates the actual strain field. In the presence of transverse shear forces, the third displacement component in eqn (6) is directly proportional to the shear strain at the neutral axis, pi,. This component and its associated shear strain distribution are shown in Fig. 3. The second component of displacement in eqn (6) is representative of the rotational degree of freedom. Its derivative in the longitudinal direction provides the normal flexural stresses associated with bending moments. It also provides a constant shear strain contribution. In the following, the second and third terms of eqn (6) will be referred to as shape functions, albeit, across the thickness. The shape of the full shear strain distribution is a combination of the two contributions, and is depicted in Fig. 3. It is noted here that there is no strict constraint that imposes a perfect parabola shape on the shear strain distribution across the thickness, similar to that in the third-order shear deformation theory. Since each component of longitudinal displacement across the thickness having a contribution to the out of plane shear strain distribution is a totally independent degree of freedom, there is no assurance of a zero shear strain at the top and bottom fibers. Everything is relegated to the self adjustment capability of the finite element these new components formulation. Moreover, contribute to the longitudinal normal strains as well as the transverse shear strains. This may cause a discrepancy in the results in the case of discontinuity in moment and/or shear force in the structure. This issue will be discussed later.
,
uy (z> Fig. 4. Unsymmetric
r, (4
shear strain and its normal displacement distribution.
In the case of anisotropy, elasto-plastic behavior or material models where the behavior in tension is different from compression, unsymmetric shear strain distribution across the thickness may follow. These strain fields can be better represented by introducing more sophisticated shape functions in the same way as above. In this work a new shape function is added to the longitudinal displacement distribution across the thickness, eqn (6). This function must comply with zero transverse shear strain conditions on the top and bottom surface of the shell. Figure 4 is an example of the desired function and its contribution to the transverse shear strain distribution across the thickness that is anti-symmetric. This new shape function is constructed from spline functions to comply with the requirement of zero slope on top, bottom and middle surface:
iz2-$z’sign(z)
1 @aI (gb)
in which y0 represents the maximum transverse shear strain contribution at one fourth of the height of the cross-section. This shape function has no contribution to the flexural deformation of the cross-section. On the other hand, close attention must be paid to its contribution to the axial deformation of the element, which is normally represented only by the translational degree of freedom uO. This contribution is particularly pronounced at the discontinuity locations of axial forces such as point loads and support reactions. This problem will be addressed later. It should be noted that for each across-thethickness shape function, one degree of freedom is added in each direction at each node. Thus, for the parabolic and antisymmetric shear distributions, two additional degrees of freedom in each dtrection must be added to the ordinary five degrees of freedom at each node. The increase in computational effort is thus substantial. The use of the additional degrees of freedom should be based on the specific characteristics
1083
Degenerated plate-shell elements of each problem. If there is no possibility of having an unsymmetric behavior in the shell thickness, it is better to suppress or remove these degrees of freedom from the formulation. It is convenient to give the new displacement shape functions, the third component of eqn (6) and that of eqn @a), the same extreme value as that of the rotational degree of freedom contribution, the second term of eqn (6), at the top and bottom of the crosssection. This has no effect on our formulation because the additional components will be treated as shape functions. The intention is to provide the ability to compare the relative contributions of each component. Introducing t as the local coordinate in the thickness direction, the final equations for shape functions across thickness are as follows: u(t) = ug + U,(f) + u,_,(t)
+ u;.(t)
Fig. 5. Definition of rotational degrees of freedom. ence. The rotational degrees of freedom were defined as changes in the angles, $, the normal to the middle surface described with respect to the global frame of reference. Thus, “‘tl, and “‘~1~ are defined as
t = 2z/h u,(t) = ahtj2
“‘al = (“‘tj, - ‘$])
and
ma2= (“tiz - 0tj2), (10)
u1 -p(t) =&f(t) u,(t) = ?“g(t) f’(t) = ht’/2 g(t) = 3ht2/2 - ht3 sin(t).
(9)
Here g(t) and f(t) are shape functions. The concept of shape function that was used in the above representation can be easily expanded to the other type of required refinements in the shells, particularly in multi-layered shell elements. This can be achieved by introducing a piecewise smooth shape function across the thickness as a broken line. This function can accommodate the discontinuity of transverse shear strain across the thickness to provide the required continuity of shear force between layers. The substitution of such a function instead off(t) or g(t) does not need any special treatment in the formulation. 3. GEOMETRIC
TRANSFORMATIONS
The model developed above was implemented into Program NISA-80 that was developed at the University of Stuttgart by Stegmuller et al. [12]. This program was written in FORTRAN77 and contains the three-dimensional degenerated plate-shell elements developed by Ramm [ 10, 111. In expanding the formulation of the degenerated shell element, the notation used by Ramm [lo, 111,is used with minor changes. In the original element the degrees of freedom at each node were three translations and two rotations: %, , m~2,“‘u3, “‘a, and “‘u2. The translations take place at the middle surface and are defined with respect to the global axes of refer-
in which m denotes the current configuration and o denotes the original undeformed configuration. Figure 5 shows the definition of these angles. The four new degrees of freedom, “/3,, “f12, “y, and “‘yzintroduced in this work are naturally similar to, and are modeled after, the rotational degrees of freedom m~l and m~2. There are, however, two main differences. The new degrees of freedom have zero values in the original configuration. More importantly, since they have the magnitude of shear strain values they are not likely to have extremely large values. Therefore, they are treated as small displacement changes. This is in contrast to the contribution of formal changes in the angles of normal “$, and m$z (or m~l and Yx2), which existed in the original formulation, and can be as large as necessary. 4. FINITE ELEMENT DISPLACEMENT
FIELD
Ramm [I l] describes the mapping of a vector of length L along the normal to the middle surface in the current configuration that undergoes small changes in the angles “‘Ic/,and “‘tj2. The magnitudes of these small changes are CI,and a,. The quantities c(, and a2 are thus the incremental changes in angles “$I and “‘$*. The translations of the tip of the vector relative to the global frame of reference are thus A”x, A”x, 1 A”x, i 0 =
-sin(“$,) sin(“&) sin(“$,) ~os(~+~)
-sin(“ll/,) cos(?j,) cos(m+2) z: L cos(m$,) sin(“$,) 111 (1 la)
M. Ziyaeifar and A. E. Elwi
1084
{Amx} =[F(“$,,“‘$,]{cc}
(symbolically).
(llb)
If the changes in the original degrees of freedom of rotation are small from one step to the next, the form of eqn (11 b) can be used to obtain the corresponding contribution of the new degrees of freedom to the displacement field. Considering the above, the full displacement field of a point on the cross-section of an element at position t along the normal to the middle surface is described as
followed and material nonlinearities are allowed. Assuming that the response is linear in the increment between configurations m and m + 1, the principle is written as:
s
e,,Ct,kdekl dV +
‘“St,h, d V
Y
=
s s
m+ ‘f&,
m+ ‘fkduk dA +
A
-
“‘S,,e,, d V.
”
+ f
cwv, I- cos(“kh)
1
sin(?j,) cos(m+2) - sin(O$,) cos(Oi+$) sin(“tj,) sin(“&) - sin(“+,) sin(“&) 1
+f(t) [wk,
)
V2)1
1:;II 2
+ g(t)[wv,
7 “ti2)l
;I {
.
dV
sY
(14)
Here e is the linear part of the increment of Green-Lagrange strain tensor, rl is the nonlinear part, S is the second Piola-Kirchhoff stress tensor, C is the constitutive tensor, t is the vector of surface tractions and f is the vector of body forces per unit volume. The tensors e and 1 are defined in terms of the displacement derivatives and its incremental values as
(12)
I
To obtain the increment of displacement field, one notes that the corresponding increments in the rotational degrees of freedom are also small enough to be represented by a transformation based on eqn (11 b). If the same interpolation shape functions are used over the middle surface of the element, the increment of displacement field is written in terms of the nodal degrees of freedom as
In eqn (15) at a point defined by I, q and t in the element the strain tensor, e and increment of displacement gradient, u,,,, are written in matrix form in terms of the nodal displacement increments as
{e> =n&l[m~Lln{~n~
UW
{Vu1= ,;, [m4dn{~n~~
(16t.4
Substituting in eqn (14) and recognizing the arbitrary variation of the displacement increment, the equation of equilibrium is recovered as [m~L+mK~J(~} ={“‘+‘I? -mQ}
(13) Here, /IIt /12,y, and y2 are increments in the degrees of freedom “‘/l,, mjJ2, my, and “‘y2, respectively. The node number n varies from 1 to the total number of nodes in an element, I and [N], is the interpolation submatrix at node n. 5. FINITE
ELEMENT
(17)
in which the first bracket on the left hand side represents the linear (tangent) and geometric stiffness matrices, while the right hand side represents the current level of loads and the internal equilibrating forces. The stiffness matrices are recovered as
V’KJ =
[“&1T[C1[m4_1 dK sY
FORMULATION
In the following the element formulation is presented in the context of the incremental principle of virtual work. A total Lagrangian approach is
In the following, the forms of the matrices on the right hand side of eqns (16) and (18) are defined.
1085
Degenerated plate-shell elements 6. DERIVATION
OF THE STRAIN MATRIX
Next the transformation reference takes the form
DISPLACEMENT
The third and fourth terms of eqn (15) involve products of the increment of displacement gradient and the displacement gradient tensor in configuration”. Thus, the final representation of strain-displacement derivatives can be written as @I = rmfflPI
(19)
in which [“‘HI contains a template of the current displacement gradient, while {Vu} is the gradient of the increment of displacement. Equation (19) is shown in full in Fig. 6. Components of the gradient of the increment of displacement vector are obtained using the form of eqn (16b). Thus, the first step is the calculation of displacement derivatives matrix [I&]. Displacement derivatives must be first interpolated from nodal degrees of freedom in the natural local coordinate system (5, q, t), and then transformed to the global coordinates. For example, in order to evaluate the gradient of the increment of displacement vector, one must first evaluate the gradient with respect to the local coordinates:
i=l,2,3
and
{=l,q,t
to the global frame of
in which [J]-’ is the inverse of the Jacobean of transformation matrix. Combining eqns (20) and (21) and comparing to eqn (16b), the displacement gradient is obtained in the form
Pi
= 1 [Jl-'["N'(I3
9,
t>ln{U>,
“=I
=,i,
[m4-JL1n~~~“~
The vector on the right hand side represents the displacement degrees of freedom of node n, while [“‘&In stands for node n submatrix, and is shown in Fig. 7. In this figure R’; are the components of the transformation matrix defined in eqn (1 l), while the shape function derivatives are defined as
(20)
4;,, = [g(t)d(, + g’(t).G’@l. in which the subscript n denotes the submatrices (vectors) of one node in the element. Here the matrix [N’] holds the shape function derivatives with respect to the natural local coordinates of node n evaluated at the current integration point with coordinates {, r~ and t. The vector {u}, holds the displacement increments at node n (remember we have here nine degrees of freedom per node).
ell
(22)
(23)
The strain displacement matrix operator [mBL]nis obtained from premultiplication of matrix [“B& by [‘“HI and the strain increment of eqn (16a) is finally obtained as @,,I =
1
e22 en 2e12 2e13
2ez3
Fig. 6. StrainAisplacement derivatives relationship.
1086
Fig. 9. Displacement derivatives relationship for large displacements.
1087
Degenerated plate-shell elements
The final result is shown in Fig. 8. The following expressions were used for the sake of simplicity:
0% = (1 + “uz,,)&, + m4.*F;, +
mU3.2F;,
073
“%3G
=
(1
+
“u3,3)%
+
%,34,
+
D’i, = (1 + mu3,3)F;Z+ mu,.3F;z +
(25)
m~2,3F;Z.
The gradient of the total displacement, yet to be determined. The procedure is that used for the increment of displacement eqn (22), except that here, possible large rotations have to be considered
8. BOUNDARY
{“y:,) has similar to gradient, values for
The dimensions of matrix [“#‘I and vector {“zi} are slightly different from those used with their counterparts used in forming the gradient of the displacement increment, because of the large rotational characteristics of the total displacements. In such a case, due to complicated trigonometric relationship between the two rotational degrees of freedom $, and &, it is convenient to use three equations for the two degrees of freedom, instead of going for two separate complicated expressions. A detailed description of eqn (26) is shown in Fig. 9. 7. STRESS STATE MATRIX AND VECTOR AND (“9)
(“S]
There are only six components of stress state all together, but to comply with matrix multiplication requirements, they have to be manipulated into a square matrix form
CONDITIONS
Boundary conditions are constraints on certain force or displacement components. Simple constraints take the form of prescribed displacements or forces. The new degrees of freedom introduced herein cannot have nonzero prescribed displacements or forces. The associated nodal forces for the refined degrees of freedom are always zero, regardless of the loading condition on the structure. The major concern herein, is to determine for different cases whether to suppress these degrees of freedom at certain locations or not. It must be noted once again that the new degrees of freedom interact with ordinary degrees of freedom. The parabolic shear degrees of freedom & and fiZ not only contribute to the shear strain distributions, but also to rotations and, thus, bending moments. Therefore, an interaction with the rotational degrees of freedom tlr and c(~would be expected. Similarly, the unsymmetric shear degrees of freedom y, and yZ interact with the axial deformations. As a starting point it would be necessary to switch off the interaction phenomena at certain points to avoid undesirable strain fields. The parabolic shear strain degrees of freedom are rotational degrees of freedom. Thus, if a fixed rotational boundary condition is required, these degrees of freedom should be suppressed along with ordinary rotational degrees of freedom. The same reasoning, however, cannot be used for the unsymmetric shear strain degrees of freedom. Although the latter degrees of freedom can represent axial deformations as well as unsymmetric shear strain distributions, the middle surface of the element is not affected by these degrees of freedom. Thus, there is no need to suppress these degrees of freedom where the axial in-plane degrees of freedom are suppressed. Where a concentrated load is applied in the in-plane direction, the unsymmetric shear strain degrees of freedom, together with ordinary m-plane degrees of freedom, can represent the pinching effect on the middle surface of the cross section. Close attention I
m& 0
[“S] =
mS22
0
0
mS*9
mS,*
0
0
“Szz
0
mS*,
0
0
mS,,
mS,3
0
0
Y213
0
0
0
“S,,
0
0
0
“S,,
0
0
T,j
0
0
mS32 0
0
0
mS,,
0
0
mS23 0
The stress state vector (“j),
takes the traditional form
(“Q
“SlZ
= Y%
“SZZ
mS33
mS13 mS,,>.
(28)
(27) m&s33
mS33 0
v**
1
must be paid to this effect at point load locations. With no constraints on the y degrees of freedom, a concentrated load that is assumed to be applied at the
1088
M. Ziyaeifar and A. E. Elwi
middle surface is distributed nonuniformly on the cross-section. Therefore, the y degrees of freedom must be restricted where a concentrated load on the in-plane direction is expected on the structure. Support reactions with components in the in-plane direction are an example of this case. A node at which a concentrated load is applied with a component in the in-plane direction, is a potential candidate for restricting these degrees of freedom. On the other hand, suppressing this degree of freedom removes one of the most desirable aspects of refinement in the shell element. Therefore, it must be done only at the nodes at which a high magnitude of concentrated in-plane force is expected. The same phenomenon arises in the case of concentrated moments. In this situation there is strong interaction between the tl and B degrees of freedom. Therefore, the contribution of /3 to the final strain distribution across the thickness may not be proportional to the actual shear force. In this case the strain field is no longer acceptable. A remedy for this problem is to suppress the parabolic shear strain distribution degrees of freedom at the point of concentrated moment. Again, one should keep in mind that this solution may be advisable only in the presence of a high concentrated moment, because it eliminates the contribution of refinement on parabolic shear strain distribution. 9. EXAMPLES
In order to verify the performance of the refined element, a number of problems were solved and the solutions compared to those obtained from the original degenerated plate-shell element, hereinafter called the “ordinary element”, and other solutions. The examples include a shell structure, two square plates and two beams. In the following the examples are presented and the results discussed. The first example is a shallow spherical shell supported on four corners and subjected to a uniform normal pressure loading. The material property is assumed to be isotropic linear elastic. The consequences of using the refined element are manifested in the flexibility of the structure and the stress
0.1
1
“I
0
CenbtqOpoint, vertical
distribution. Both effects are thickness dependent. Therefore, this example has been solved with two thickness units (t = 200 and 2000 mm). The extreme thickness values have been used deliberately to bring out an exaggerated response that can be seen clearly. Figure 10 shows the load vs vertical deflection response of point B for both the ordinary element and the refined element in the case of the thin shell (t = 200 mm). The two solutions are almost identical. Figure 11 illustrates the normal and shear stress distributions at point A that is close to a support. Note that the difference in the normal stress distribution between the solutions of the ordinary element and the refined element is clear and, in contrast to the shear stress distribution, is highly thickness dependent as is shown in Fig. 12 for the case of the thickness shell structure (t = 2000 mm). This is expected, since the plane sections do not remain plane in the presence of high shear gradients in the refined element. The ratio of the shear stresses to the normal stresses in the thick shell at this location is much higher than that in the thin shell at the same location. The next example consists of a square plate supported at four corners and subjected to a uniform pressure. The material has a yield strength of 240 MPa and is modeled with an elastic perfectly plastic
shear stress MPa)
F
% 8 1 .z CI
50
4 0 1
1
-50
-100
(mm)
Fig. 10. Thin shell, vertical deflection of central point.
100
3 3 a ._
Al!
300
-50
.v F -100
Normal stress distribution Fig. 11.Thin shell, stress distribution at a point close to support.
1089
Degenerated plate-shell elements Onset of yield at extreme fibers
Onset of yield at middle of section
/ 0
30
10
40
50
Deflec2on (mm)
Normal
stress dist
@IF%)
Fig. 12. Thick shell, nortnal stress dist. close to support.
von Mises constitutive relation. The response of the refined element shows a different progression of yield in the structure
from that of the ordinary
element.
The refined element predicts yield in the structure earlier and also in a different mode compared to the oridnary element. The onset of yield takes place at the middle of the cross-section close to the support (point A in Fig. 13a) due to the high shear stresses at that point (Fig. 13b). This limits the load carrying capacity of the structures and constitutes a “punching” shear failure. On the other hand, when ordinary elements are used, the onset of yield is postponed to a higher level of load and it starts at the extreme fibers of the plate cross-section due to a combination of out-of-plane shear stresses and in-plane normal and shear stresses. The structure then exhibits a gradually softening response. Punching shear may not be detected if a different mode of failure takes place elsewhere such as a mid-span or at the free edge. The load-deflection curve is shown in Fig. 14. The next example shows the case of unsymmetric material behavior across the thickness. This is nrovided by applying a biaxial inplane normal force to a square plate supported at the central point. This represents the region of a flat slab around a column. The material model used here is the same as that of the previous model. In this case, the onset of yield in the top and bottom fibers of the plate, takes place at
Fig. 14. Plate supported at corners, vertical deflection at a point close to support.
different levels of transverse pressure. Figure 15 shows the normal and shear stress distributions at a point close to support that clearly depicts the effect of the unsymmetric component of shear strain across the thickness. As mentioned before, both the parabolic shear strain components, p, and the rotational degrees of freedom, u, contribute to the shear and flexure deformation modes simultaneously. At a discontinuity in the transverse shear forces a combination of both a and p could yield a non-zero transverse shear strain at the top and bottom of the plate surface. However, the integral of shear stresses across the thickness is always in total balance with the actual shear force on the cross section. Under such circumstances the Co continuity of the interpolation functions for refined degrees of freedom, /?, between elements (and higher orders of continuity for the nodes inside the element) exists. This does not allow the finite element model to imitate the discontinuity in the shear strain properly. In this situation a combination of the degrees of freedom a and /I gives a gradual, but steep transition of shear strains over the discontinuity region. If a three-dimensional model were to be used with the same constraints of plate or shell, the result of the strain field would be similar to that of the refined shell element. This phenomenon is illustrated in the last example.
250
Fig. 13. Plate supported at corners, stress dist. at a point close to support.
M. Ziyaeifar and A. E. Elwi
1090
Normal stress (MPa)
Fig. 15. Plate supported at centre, stress dist. at a pomt close to support
In this example, a beam supported at both ends and subjected to concentrated loads is analyzed. The beam is fairly thick to render the shear deformation mode dominant relative to the flexural mode. Two extreme cases in which a sharp discontinuity in the shear force exists in the beam are analyzed, one with two point loads applied at l/4 points of the span and one with a single mid-span load. These beams are shown in Fig. 16. Boundary conditions of plane stress are assumed in the width direction when refined
FEM model, 256 4 node Plane stress element I +
+
elements are used (12 bicubic elements). The same beams were constructed from 256 two-dimensional isoparametric eight-node plane stress elements also shown on the same figure. The large number of elements is necessary to examine the stress points as close as possible to the plane of force discontinuity. Figure 17a shows the results for the mid-span load case compared to the plane stress solution at different points in the vicinity of the discontinuity. At points close to the discontinuity (L/290), the shear stress distribution from both the plane stress and the refined element solutions are rather flat. A short distance away from the discontinuity (L/14) both solutions become much more parabolic. Figure 17b shows the two point load solutions for both types of analysis. The points at L/290 and L/14 are in the pure moment region. However, because of localized effects both solutions give a perturbation in the shear stress that is evident at L/290 and is negligible at L/14 from the load point in this region. The solutions at a distance of L/40 from the load point measured towards the support (nonzero shear force zone) are again almost identical and parabolic.
4
A7
10. CONCLUSIONS
FEM model, 12 bicubic “refined” elements Fig. 16. Beams loaded at mid-span and at l/4 of span.
A modified three-dimensional plate shell element that has a better transverse shear strain representation
-1
(a) Load case 1
(b) Load case 2
Fig. 17. Shear stress close to points of shear force discontinuity.
Degenerated
piate-she11 elements
presented. In this element two types of degrees of freedom are introduced, a symmetric set that governs a parabolic transverse shear strain distribution and a antisymmetric set that governs an unsymmetric contribution. The new degrees of freedom are treated as shape functions across the thickness and added in the manner of small rotations to the longitudinal displacement field representation. This effectively removes the assumption that plane sections remain plane and allows two new modes of deformation. The performance of the new element has been tested and proven to give a more realistic response in the case of elasto-plastic materials. The implications of removal of the plane section restriction are discussed at length and corresponding results are shown to be consistent, even in the presence of strong discontinuities m the transverse shear force distribution. has
been
Acknowledgments-This work has been carried out m the course of doctoral graduate research of the first author. The work has been funded in part by the Natural Sciences and Engineerma Research Council of Canada. The authors also wish to acknowledge the support of the Iranian Ministry of Culture and Higher Education. The first author wishes to thank Dr M Farshad for his motivating ideas.
51, 745-752
(1984).
1091
3. N. Buechter and E. Ramm, Shell theory versus Degeneration-a comparison in large rotation finite element analysis. Int. J. numer. Meth. Engng 34, 39-59 (1992). 4. J. N. Reddy, A refined nonlinear theory of plates with reansverse shear deformation. Inr. J. Solid Struct. 20, 881-896 (1984). rotation 5. J. N. Reddy, A small strain and moderate theory of elastic anisotropic plates. J. appi Mech. 54, 623-626 (1987). models of com6. J. N. Reddy, On refined computational posite laminates Int. J. num. Merh. Engng 27, 361-382 (1989). 7. Y. Basar, Y. Ding and R. Schultz, Refined sheardeformation models for composite lammates with finite rotations. Int. J. Solids Siruct. 30, 2611-2638 (1993). 8. D. R. J, Owen and 2. H. LI, A refined analysts of laminated plates by finite element displacement methods -I. Fundamentals and static analysis, 11. Vibration and Stability. Compur. Strucr. 26, 907-923 (1987). 9. S. Ahmad, B. M. Irons and 0. C Zienkiewtcz, Curved thick shell and membrane elements with particular reference to axi-symmetric problems, In Proc. ad Cony. Malrix Methods in Structural Mechanrcs, WrightPatterson Au Force Base. Ohio (1968) ‘nichtlineare ‘Elstostatik und 10. E. Ramm, Geometrisch finite Elemente. Habilitation, Berich Nr. 76-2, lnstitut fur Baustatik, Universitat Stuttgart, Germany (1976). 11. E. Ramm, A plate/hell element for large deflections and rotations. US-Germany Symp. on Formulations and Computational Algorithms in Finite Element Analysa, MIT 1976, MIT Press (1977). 12. H. Stegmuller, L. Hafner, E. Ramm and J. M. Satteie. Theoretische grundlagen rum FE-programm system NISA-80. Mitteilung Nr. I, Instttut fur Baustattk der Univetsitat Stuttgart, Germany (1983).