3.. .t. 4.!

‘1! u'l

n.||v pl.
[0)" it...
u Ill-"I

v I

VI. ”.4 . Krrll.

 

 

«2&8

RSILTY IBRAR

lHILHl

 

 

 

Hllzllll

 

Illlllllllllllllllllllll‘

31293 01694

This is to certify that the

dissertation entitled

NEW ZIG-ZAG COMPOSITE PLATE THEORIES AND 3-D
FINITE ELEMENTS FOR STATIC AND EXPLICIT DYNAMIC
ANALYSIS OF LAMINATED COMPOSITE AND SANDWICH PANELS

presented by

Yong Bae Cho

has been accepted towards fulﬁllment
of the requirements for

Ph.D Mechanics

degree in

 

 

W/W

Major professor

Date December 22, 1997

 

MSU is an Afﬂrmaliw Action/Equal Opportunity Institution 0‘ 12771

 

LIBRARY
Mlchlgan State
Unlverslty

 

 

 

PLACE IN RETURN BOX
to remove this checkout from your record.
TO AVOID FINES return on or before date due.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ma clone/Mam“

 

NEW ZIG-ZAG COMPOSITE PLATE THEORIES

AND 3-D FINITE ELEMENTS FOR STATIC AND EXPLICIT DYNAMIC

ANALYSIS OF LAMINATED COMPOSITE AND SANDWICH PANELS
By

Yong-Bae Cho

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Materials Science and Mechanics

1998

ABSTRACT
NEW ZIG-ZAG COMPOSITE PLATE THEORIES
AND 3-D FINITE ELEMENTS FOR STATIC AND EXPLICIT DYNAMIC
ANALYSIS OF LAMINATED COMPOSITE AND SANDWICH PANELS
By

Yong-Bae Cho

Reﬁned laminated plate theories based on linear or cubic plus zig-zag kinematics have
been developed, and new CO 3-D ﬁnite elements have been developed on the basis of those
theories for the purpose of static and dynamic analysis of sandwich and laminated
composite panels. In the ﬁrst order zig-zag plate theory, the in-plane displacement ﬁelds in
each sublaminate are assumed to be piecewise linear and vary in a zig-zag fashion through
the thickness of the sublaminate. In the third order zig-zag plate theory, the in-plane
displacement ﬁelds in the laminate are assumed to be piecewise cubic and vary in a zig-
zag fashion through the thickness of the laminate. In each case, the transverse
displacement ﬁeld is assumed to vary linearly through the thickness. The zig-zag functions
are evaluated by enforcing the continuity of shear stress at each interface. This in-plane
displacement ﬁeld assumption accounts for discrete layer effects without increasing the
number of degrees of freedom as the number of layers is increased. In order to maintain
C0 continuity of all variables, new rotation degrees of freedom are introduced, and
subsequently constrained via the penalty method. The transverse normal strain is
improved by assuming a constant transverse normal stress through a sublaminate or the
entire laminate, and hence Poisson’s ratio stiffening in associated with ﬁnite element

model is eliminated.

The proposed C0 ﬁnite elements have the topology of an eight-noded brick. Each node
has ﬁve engineering degrees of freedom - three translations and two rotations. In-plane
displacements and rotational degrees of freedom are approximated by the bilinear
Lagrange interpolation functions. For transverse deﬂection degrees of freedom, an
interdependent interpolation is utilized. The element stiffness coefﬁcients are integrated
exactly, yet the element exhibits no shear locking. This is achieved by using the consistent
transverse shear strain ﬁelds as well as the edge-consistency of the tangential shear strain
on any common inter-element edge. The elements are shown to be accurate, simple to use,
and compatible with the requirements of commercial ﬁnite element codes. Based on the
assumed displacement ﬁelds and the ﬁnite element approximations, consistent mass
matrices has been obtained. Making use of the consistent mass matrices, lumped mass
matrices have been derived according to HRZ lumping scheme. The developed elements
have been implemented into the explicit ﬁnite element dynamics code, NEPTUNE,
wherein the internal force vector is computed for the element.

Static, free vibration, and explicit structural dynamic analyses of sandwich panels and
laminated composites are performed using the elements. The numerical results show that

the elements are accurate, robust, and computationally efﬁcient.

To my late father, Hahk-Kyu Cho,
my mother, Eul-Soon Lee,

and my family

iv

ACKNOWLEDGEMENTS

I wish to express my sincere gratitude to my advisor, Dr. Ronald C. Averill for his
continuous encouragement, patience, suggestions and technical guidance throughout the
period of my Studies in Michigan State University and in the preparation of this
dissertation.

I would also like to place on record my appreciation to my doctoral committee
members, Dr. Nicholas J. Alterio, Dr. Dahshin Liu and Dr. Tien Yuen Li for their guidance
and suggestions during the course of this research.

I am grateful to Dr. Ronald Kulak and Dr. Ed Plaskacz of Argonne National
Laboratory who allowed me to use NEPTUNE, the in-house code for explicit dynamic
analysis, and for their valuable advice on its usage.

I gratefully acknowledge the Computational Mechanics Research Group (CMRG) for
their helpful discussions and contributions to my work. I specially thank Yuen Cheong
Yip, Venkateshwar Rao Aitharaju, Lorenzo Marco Smith and David Eby for their support
and encouragement.

Finally, I wish to thank my whole family for the emotional, moral and ﬁnancial

support, if not for which, this research would not have been possible

TABLE OF CONTENTS

TABLE OF CONTENTS ....................................................................................... vi
LIST OF TABLES ............................................................................................... viii
LIST OF FIGURES ............................................................................................... ix
CHAPTER I

INTRODUCTION ................................................................................................... 1
1 . 1 Introduction ................................................................................................. l
1.2 Goals and Objectives of Current Research ................................................. 4
1.3 Literature Review: Plate Theories ............................................................... 5
1.4 Literature Review: Three Dimensional Composite Shell/Plate Elements 11
1.5 Literature Review: Mass Matrices and Lumping Schemes ...................... 13
1.6 The Present Study ..................................................................................... 16
1.7 Organization of the Dissertation ............................................................... 18
CHAPTER 11

FIRST ORDER ZIG-ZAG BEAM THEORY BASED ON SUBLAMINATE
APPROXIMATIONS ............................................................................................. 21
2.1 Introduction ............................................................................................... 21
2.2 Kinematic Assumptions and Formulation of the Beam Theory ............... 22
2.3 Finite Element Formulation ...................................................................... 34
2.4 Numerical results ...................................................................................... 39
2.5 Summary .......................................................... _ ......................................... 45
CHAPTER III
PLATEELEMENTBASEDONFIRSTORDERZIG—ZAGSUBLAMINATEI'HEORY
................................................................................................................................ 57
3. 1 Introduction ............................................................................................... 57
3.2 Kinematic Assumptions ............................................................................ 58
3.3 Formulation of the Plate Theory ............................................................... 61
3.4 Governing Equations and Boundary Conditions ...................................... 69
3.5 Finite Element Formulation ...................................................................... 72
3.6 Poisson’s Ratio Stiffening Effect .............................................................. 97

vi

3.7 Numerical Results ..................................................................................... 99

3.8 Summary ................................................................................................. 104
CHAPTER IV

PLATE ELEMENT BASED ON CUBIC ZIG-ZAG THEORY .......................... 119
4.1 Introduction ............................................................................................. l 19
4.2 Formulation of Plate Theory ................................................................... 120
4.3 Finite Element Formulation .................................................................... 131
4.4 Numerical Results ................................................................................... 136
4.5 Summary ................................................................................................. 139
CHAPTER V

FREE VIBRATION ANALYSIS ......................................................................... 149
5. 1 Introduction ............................................................................................. 149
5.2 Mass Matrices ......................................................................................... 149
5.3 Free Vibration Analysis .......................................................................... 152
5.4 Numerical Results ................................................................................... 153
5.5 Summary ................................................................................................. 155
CHAPTER VI

EXPLICIT DYNAMIC ANALYSIS .................................................................... 167
6. 1 Introduction ............................................................................................. 167
6.2 Explicit Dynamic Analysis ..................................................................... 168
6.3 Critical Time Step ................................................................................... 170
6.4 Amplitude and Period Errors in Central Difference Approximation ...... 175
6.5 Geometrical Nonlinearities - Corrugation Approach .............................. 180
6.6 Numerical Results ................................................................................... 182
6.7 Summary ................................................................................................. 184
CHAPTER VII

CONCLUSIONS .................................................................................................. 192
7.1 Summary ...................................... 192
7.2 Conclusions ............................................................................................. 193
7.3 Recommendation .................................................................................... 195
APPENDICES

APPENDD( A ...................................................................................................... 196
APPENDD( B ...................................................................................................... 202
BIBLIOGRAPHY ..................................................... '. .......................................... 213

vii

Table 1.
Table 2.
Table 3.

Table 4.
Table 5.
Table 6.
Table 7.
Table 8.
Table 9.
Table 10.
Table 11.

LIST OF TABLES

Material properties used in the analyses. ................................................... 40
Ply stacking tables for the three laminates analyzed. ................................ 41
Material properties for TACOM’S composite armored vehicle (CAV) panel
.................................................................................................................. 102
Material properties for sandwich composite plate ................................... 105
Lamination Scheme for Sandwich Panel ................................................. 105
Material properties for damaged composite plate .................................... 140
Lamination Scheme for damaged composite plate .................................. 140
Material Properties for Sandwich Composite Plate ................................. 154
Lamination Scheme for Sandwich Panel ................................................. 157
The higher natural frequencies in the mesh ............................................. 157
Material properties ................................................................................... 185

viii

Figure 1.

Figure 2.
Figure 3.

Figure 4.

Figure 5.

Figure 6.

Figure 7.

Figure 8.

Figure 9.

Figure 10.

Figure 1 1.

Figure 12.

Figure 13.

Figure 14.

Figure 15.

LIST OF FIGURES

Schematic of sublaminate and layer divisions in the laminated beam theory.

.................................................................................................................... 23
Element topology and nodal degrees of freedom. ..................................... 36
Simply supported laminated composite beam subjected to sinusoidal loading.
.................................................................................................................... 39
Variation of normalized center deﬂection and maximum bending stress
versus number of elements in the mesh (Laminate 1, Uh=100) ................ 46
Variation of normalized center deﬂection versus beam aspect ratio
(Laminate 1). .............................................................................................. 47
Inplane displacement versus normalized thickness coordinate (Laminate 1,
L/hﬂ). ....................................................................................................... 48
Bending stress versus normalized thickness coordinate (Laminate 1, Llh=4).
.................................................................................................................... 49
Variation of normalized center deﬂection versus beam aspect ratio
(Laminate 2). .............................................................................................. 50
Inplane displacement versus normalized thickness coordinate (Laminate
2, Llh=4). ................................................................................................... 51
Bending stress versus normalized thickness coordinate (Laminate 2, Llh=4).
.................................................................................................................... 52
Variation of normalized center deﬂection versus sandwich beam aspect
ratio (Laminate 3). ..................................................................................... 53
Inplane displacement versus normalized thickness coordinate (Laminate 3,
Llh=4). ....................................................................................................... 54
Bending Stress versus normalized thickness coordinate (Laminate 3, Llh=4).
.................................................................................................................... 55
Bending stress in face sheets versus normalized thickness coordinate
(Laminate 3, Llh=4). .................................................................................. 56

Schematic of sublaminate and layer divisions in the laminated plate theory.
.................................................................................................................... 59

ix

Figure 16.

Figure 17.
Figure 18.
Figure 19.
Figure 20.
Figure 21.

Figure 22.

Figure 23.

Figure 24.

Figure 25.

Figure 26.

Figure 27.

Figure 28.

Figure 29.

Figure 30.

Figure 31.

Figure 32.

Figure 33.

Figure 34.

Figure 35.

A composite under transverse tension load (a) and the corresponding

transverse strain distribution through the thickness (b) ............................. 67
Element topology and nodal degrees of freedom. ..................................... 73
Quadrilateral element coordinates and deﬁnition of rotation .................... 76
Coordinate transformation between edge and global coordinates ............. 80
Schematic and deﬁnition of rotation and deformed shape ......................... 92
Simply supported laminated composite plate subjected to sinusoidal loading.
.................................................................................................................. 100
Variation of normalized center deﬂection versus plate aspect ratio at the
center of the plate (CAV Panel, alh=4). ................................................... 106
In-plane displacement, u, versus normalized thickness coordinate (CAV
Panel, alh=4). ........................................................................................... 107
Bending stress, on, versus normalized thickness coordinate at the center
of the plate (CAV Panel, alh=4). .............................................................. 108
In-plane Shear stress, 6x» versus normalized thickness coordinate at x, y =
0, (CAV Panel, a/h=4). ............................................................................. 109
Transverse shear stress, 012, versus normalized thickness coordinate at x =
0, y = 2/a (CAV Panel, alh=4). ................................................................. 110
Transverse normal strain, ea, versus normalized thickness coordinate at
the center of the plate (CAV Panel, alh=4). ............................................. 111
Transverse normal stress, ow versus normalized thickness coordinate at
the center of the plate (CAV Panel, alh=4) .............................................. 112
In-plane displacement, u, versus normalized thickness coordinate
(Sandwich Panel, alh=4). ......................................................................... 113
Bending stress, on, in face sheets versus normalized thickness coordinate
(Sandwich plate, alh=4). .......................................................................... 114
In-plane shear stress in the face sheets, 0,), versus normalized thickness
coordinate at x, y = 0, (CAV Panel, alh=4). ............................................. 115
Transverse shear stress, on, versus normalized thickness coordinate at x =
0, y = 2/a (Sandwich Panel, alh=4). ......................................................... 116
Transverse normal strain, 82,, versus normalized thickness coordinate at
the center of the plate (Sandwich Panel, alh=4). ..................................... 117

Transverse normal strain, ca, versus normalized thickness coordinate at
the center of the plate (Sandwich Panel, alh=4). ..................................... 118

Schematic of a laminate and layer divisions in the cubic zig-zag laminated
plate theory. .............................................................................................. 121

Figure 36.
Figure 37.

Figure 38.

Figure 39.

Figure 40.

Figure 41.

Figure 42.

Figure 43.

Figure 44.

Figure 45.

Figure 46.

Figure 47.

Figure 48.

Figure 49.

Figure 50.

Figure 51.

Figure 52.

Figure 53.

Figure 54.

Element topology and nodal degrees of freedom. ................................... 132

Simply supported laminated composite plate subjected to sinusoidal loading.
.................................................................................................................. 137

Variation of normalized center deﬂection versus plate aspect ratio at the
center of the plate. .................................................................................... 141

In-plane displacement, u, versus normalized thickness coordinate (alh=20)
.................................................................................................................. 142

Bending stress, on, versus normalized thickness coordinate at the center
of the plate (alh=20). ................................................................................ 143

In-plane Shear stress, 6,), versus normalized thickness coordinate at x, y =
0, (alh=20). ............................................................................................... 144

Transverse shear stress, on, versus normalized thickness coordinate at x =
0, y = 2/a ................................................................................................... 145

Transverse shear strain, ‘sz. versus normalized thickness coordinate at x =
0, y = 2/a ................................................................................................... 146

Transverse shear stress, 0'”, versus normalized thickness coordinate at x =
2/a, y = 0 ................................................................................................... 147

Transverse normal strain, ea, versus normalized thickness coordinate at
the center of the plate (alh=20). ............................................................... 148

Variation of normalized lowest natural frequency versus plate aspect ratio
by consistent mass matrix and lumped mass matrix (FZZ3D). ............... 158

Variation of normalized lowest natural frequency versus plate aspect ratio
by consistent mass matrix and lumped mass matrix (CZZ3D). ............... 159

The ﬁrst eigen mode of transverse deﬂection d.o.f of simply supported
sandwich panel (alh=4) ............................................................................ 160

The second eigen mode of transverse deﬂection d.o.f of simply supported
sandwich panel (alh=4) ............................................................................ 161

The third eigen mode of transverse deﬂection d.o.f of simply supported
sandwich panel (alh=4) ............................................................................ 162

The fourth eigen mode of transverse deﬂection d.o.f of simply supported
sandwich panel (alh=4) ............................................................................ 163

The ﬁfth eigen mode of transverse deﬂection d.o.f of simply supported
sandwich panel (alh=4) ............................................................................ 164

The second eigen mode of transverse deﬂection d.o.f of simply supported
isotropic plate (alh=4) .............................................................................. 1.65

The third eigen mode of transverse deﬂection d.o.f of simply supported
isotropic plate (alh=4) .............................................................................. 166

xi

Figure 55.
Figure 56.
Figure 57.
Figure 58.

Figure 59.
Figure 60.

Figure 61.
Figure 62.

Amplitude and period errors in central difference solution ..................... 176

The corotational coordinate system ......................................................... 181
The loading and boundary conditions and the geometry of plate ............ 186
The time history of Center deﬂection of orthotropic square plate under

suddenly applied pressure loading. .......................................................... 187
The transient response at the center of two-layered angle-ply (-45°/45°.)

plate under suddenly applied pressure loading. ....................................... 188
The effect of material orthotropy on the transient center deﬂection analysis
.................................................................................................................. 189
An Irnpulsively Loaded Cylindrical Sandwich Panel .............................. 190

The time history of displacement of the sandwich cylindrical panel at the
center of outer surface and at node 13. .................................................... 191

xii

CHAPTER I
INTRODUCTION

1.1 Introduction

Fiber reinforced laminated composites have been used in structural members due to
better stiffness-to-weight and strength-to-weight ratios than those of conventional mono-
lithic materials. The primary characteristics of laminated composites are that the material
properties are orthotropic in the laminate plane, the ratio of transverse shear modulus to
in-plane modulus is low and layers are laminated through the thickness. Therefore, unlike
homogeneous isotropic materials, the behavior of laminated composites is unique and
complicated at both the global and the local levels. Moreover, the early laminated compos-
ites were applied to thin structures and were designed to withstand mainly in-plane tension
loads. But recently they are also employed for primary load-carrying structural members
in automobile front aprons, bridge frames, submarine hulls and aircraft because of
advances in their manufacturing techniques. Such recent applications of composite materi-
als involve the use of thick-section laminates containing perhaps over one hundred layers
and/or sandwich construction that contains a thick core between face sheets. Furthermore,
such structures are typically designed to withstand complex mechanical loading states,

e.g., static loadings and dynamic impact loading, in the presence of harsh environmental

conditions.

Because the importance of laminated composites in aerospace, automotive, civil and
marine structures continues to increase, the development of accurate structural theories
and efﬁcient computational models for analyzing these structures is currently a critical
area of research. Especially important is to develop accurate, efﬁcient, and convenient
models that can be used for static, dynamic, and progressive failure analysis of thick and
thin laminated composite and sandwich structures. Thus, for the purpose of modeling the
behavior of laminated composites efﬁciently under the various loading conditions in ser-
vice, many reﬁned plate theories have been developed and pTOposed on the basis of vari-
ous kinematic assumptions.

The analytical solutions of such proposed theories are conﬁned by the conﬁgurations
of plate, boundary conditions and loading, and are usually not available for most real prob-
lems. The way of overcoming such obstacles is to make use of the ﬁnite element method, a
prevalent tool used to obtain approximate solutions in complicated engineering problems.
The quality of a ﬁnite element can be signiﬁcantly affected by the kinematic assumptions
made throughout its formulation. Thus, every aspect of the problems to be solved should
be reﬂected in the kinematic assumptions of ﬁnite element formulation.

In the analysis of moderately thick laminated composites, local effects are so impor-
tant that the three dimensional kinematics should be used on the development of ﬁnite ele-
ment models. Moreover, when the loading is excessive, the composites undergo local
damages, e.g., ﬁber breakage and/or ﬁber kinking, and local delarnination. Thus, the ele-
ment should have capabilities to capture more accurate distribution of local strains and

stresses and predict damage in composite structures associated with proper failure criteria.

However, conventional three-dimensional ﬁnite elements are computationally too expen-
sive to be used in the analysis of laminated composites to capture local effects of a com-
posite structure, if one element per layer is used. Thus, a lot of research has been carried
out to develop a robust ﬁnite element model that is more computationally efﬁcient than the
conventional displacement based continuum ﬁnite element and describes adequate local
kinematics.

When external loading is varying with time or suddenly applied to a composite struc-
ture, e.g., seismic loading and impact loading, the response of the composite structure is
time dependent and can be obtained using structural dynamics techniques. Structural
dynamics is one of the rapidly expanding areas of application of the ﬁnite element method.
Many problems in structural dynamics cannot be solved effectively by analytical methods,
and the advances in both the digital computer and various associated numerical techniques
of analysis have signiﬁcantly enhanced solution capability of the ﬁnite element method.
Therefore, the ﬁnite element method is as valuable a tool in dynamics as it is in the static
problems.

In general structural dynamics problems, the governing equations are second order or
fourth order partial differential equations in the space domain, and hyperbolic in the time
domain. In this case, the solution is a function deﬁned both on the space domain and on
the time domain. Therefore, as a numerical solution procedure in explicit dynamics analy-
sis, the partial differential equations in space and time are ﬁrst discretized in space, yield-
ing a system of decoupled ordinary differential equation in time. The disretization in space
is usually achieved by employing the ﬁnite element method, and called a ﬁnite element

semidiscretization, because though displacements are discrete functions of space, they are

still continuous functions in the time domain. To solve a system of decoupled ordinary dif-

ferential equations the ﬁnite difference method is adopted in the time domain.

1.2 Goals and Objectives of Current Research

To analyze static and dynamic responses of composite structures more accurately, ﬁrst
of all, it is required to develop a composite plate theory where the local distributions of
stresses, Strains and displacement ﬁelds within a laminate are taken into account. On the
basis of the developed composite plate theory a robust ﬁnite element model that is more
computationally efﬁcient and accurate needs to be developed. The main objective of the
current research is to develop a composite plate theory and the associated elements to pro-
vide static and dynamic analysis of composite laminate and sandwich structures. Thus the
plate theory and element must possess the following features.

. Using 3-D kinematics, a zig-zag discrete layer-wise plate theory wherein local effects
can be accounted for very well is to be developed.

. The developed plate theory should be C0 continuous in displacement ﬁelds. Thus, a C0
ﬁnite element should be developed.

- A shear correction factor Should not be necessary in the developed plate theory.

. The number of degrees of freedom of the ﬁnite element model associated with the
plate theory should be independent of the number of layers in the composites.

- The developed element should have only engineering degrees of freedom (translations
and rotations) to facilitate implementation into commercial ﬁnite element code.

- The developed element should be free of shear locking and Poisson’s ratio stiffening

effects.

- To be used in structural dynamics problems, the consistent mass matrix should be
developed. On the basis of the consistent mass matrix, the lumped mass matrix should be
developed.

- In structural dynamics problems, the internal force vector is calculated for the devel-
oped element.

- The developed element should have a critical time step (related to natural frequencies)

associated with explicit dynamic analysis that is as large as possible.

1.3 Literature Review: Plate Theories

Many plate theories have been developed on the basis of various kinematic assump-
tions to predict the response of laminates more accurately. Comprehensive reviews of lit-
erature on the development of modern plate theories have recently been carried out by
several authors [1-6]. They usually classify the plate theories by the assumptions of their
displacement ﬁelds. Such models can be roughly divided into three categories: equivalent
single layer theories, layerwiSe theories, and zig-zag in-plane displacement theories.

In the equivalent Single layer approach, the displacement components are assumed to
be C1 continuous through the thickness of the plate and functions of in-plane coordinates.
Thus, the material properties of all the layers are “smeared,” and the laminate is modeled
as an equivalent single anisotropic layer.

The classical laminated plate theory (CLPT) based on Kirchhoff - Love classical plate
theory (CPT) was the ﬁrst laminated plate theory employed in the early days of study. In
this theory, it was assumed that the in-plane displacements vary linearly through the thick-

ness, while the transverse displacement is constant through the thickness. The shortcom-

ing of this theory is that the predictions of displacements and stress become inaccurate
when the aspect ratio of the plate becomes small (less than 20) or the transverse shear
modulus is low. The shortcoming mainly comes from the fact that CPT is based on the
Kirchhoff hypothesis that lines originally straight and normal to the reference surface
remain straight and normal to the reference surface during deformation so that the traverse
shear deformation is not allowed.

Timoshenko ﬁrst recognized that rotary inertia and transverse shear deformation are of
importance in the vibration analysis of elastic bodies, and considered the effects of rotary
inertia and transverse shear deformation on the vibration. of elastic bars [7]. Then Reissner
[8] and Mindlin [9] expanded this concept to isotropic plate bending problems. On the
basis of these theories, an improved composite plate theory was developed by Mindlin,
called the ﬁrst order shear deformation theory(FSDT) [10]. This theory assumes that a line
originally straight and normal to the reference surface remains straight during deformation
but not necessarily perpendicular to the reference surface. This theory yields good predic-
tions of overall laminate behavior (e.g., deﬂections, natural frequencies, and buckling
loads) and inplane stresses provided the plate is thin and the material properties of adja-
cent layers do not differ signiﬁcantly. However, FSDT does not account for warpage of the
cross section, which may be signiﬁcant in laminated composites, because FSDT accounts
for only average kinematics. Furthermore, FSDT yields constant transverse shear strains
through the thickness of the plate, therefore a shear correction factor(s) is required to com-
pute the shear strain energy accurately [11,12].

In order to reduce the inaccuracies of the FSDT, higher order shear deformation theo-

ries (HSDT) were proposed (e.g., [13-17]). In these models, it is assumed that the dis-

placements are of higher order polynomial form and are C1 continuous through the
thickness. Using the shear traction boundary conditions on the top and bottom of a lami-
nate, the higher-order undetermined functions in the displacement ﬁeld are evaluated and
are eliminated. The assumptions of displacement ﬁelds permit nonlinear variations of dis-
placements, strains and stresses through the thickness, and thus warpage of the cross sec-
tion is allowed. However, even though they often can predict well the gross behavior of
thin and some moderately thick laminates, all equivalent single layer theories have a com-
mon shortcoming. Since they are unable to account for the sometimes severe discontinuity
in transverse shear strains that occur at interfaces between two adjacent layers with drasti-
cally different stiffness properties, they preclude the satisfaction of continuity of trans-
verse stresses at interfaces between adjacent layers with different material properties, and
do not reﬂect the proper kinematics in laminates. Therefore, the local deformations and
stresses, and sometimes even the overall laminate response, are not well predicted.

If the transverse shear stresses are needed, they can be recovered indirectly by integrat-
ing the equilibrium equations [18-22]. However, this procedure leaves controversies in its
accuracy for the analysis of highly irregularly shaped plates since it involves higher order
in-plane differentiation of the displacement functions. Moreover, the higher order in-plane
differentiation requires to use at least quadratic ﬁnite element interpolation functions to
obtain transverse shear stresses and cubic interpolation functions to capture transverse
normal stresses.

In an effort to overcome the shortcomings of the equivalent single layer approach, dis-
crete-layer (or layerwise) theories have been proposed [23-32]. These theories are based

on a unique displacement ﬁeld for each layer, and enforce interlaminar continuity of dis-

placements and sometimes of transverse stresses, as well. In addition, the through-layer
variation of displacement is interpolated by one dimensional Lagrangian linear [23,24,27]
or quadratic shape functions [30] or a high order hierarchical element [31]. Reddy [25]
assumes that the displacement ﬁelds for each layer are expanded through the thickness
direction using Lagrangian shape functions employed in the conventional ﬁnite element
method. Therefore, displacement components are C0 continuous, so that transverse strains
are piecewise continuous through the thickness. Toledano and Murakami assume indepen-
dent ﬁelds for displacements and stresSes in each layer, and employ Reissner’s mixed vari-
ational principle in the formulation of model. They Show that the in-plane displacement
can be improved by this technique in the presence of transverse shear effects. In these the-
ories, some of the dependent variables are eliminated in the development of the model by
enforcing displacement continuity across the layer interfaces. In most of these theories,
the continuity of transverse shear stresses is met in an integral sense through the potential
energy formulation. When these theories are compared with the conventional 3-D dis-
placement ﬁnite element model, the most signiﬁcant feature of these theories is that it has
a computational time - saving data structure while yielding exactly the same results for
comparable meshes. These theories predict excellent global and local distributions of in-
plane and out-of-plane displacements and stresses. The major shortcoming in these theo-
ries is their large computational expense, for as the number of layers increases, the number
of degrees of freedom increases proportionally.

A new class of laminate theory, called here the ﬁrst order zig-zag theory (FZZT), was
developed by DiSciuva in the mid 1980’s [33,34]. In this theory, in-plane displacements in

a laminate are assumed to be piecewise (layerwise) linear and continuous through the

thickness, yet the total number of degrees of freedom is only ﬁve (not dependent on the
number of layers). This is accomplished by analytically satisfying the transverse shear
stress continuity conditions at each interface in the laminate. This theory was shown to be
very accurate for many cases, especially symmetric laminates. On the basis of the concept
introduced in [33,34], DiSciuva as well as other researchers made signiﬁcant improve-
ments to the FZZT [35-42]. The primary improvement was achieved by superimposing a
piecewise linear variation of in-plane displacements on a continuous cubic function of the
transverse coordinate [35,36,40-42], creating a displacement ﬁeld that can better account
for the warping that occurs during bending of unsymmetric laminates. These latter theo-
ries, denoted here as higher-order zig-zag theories (HZZT), also satisfy the homogeneous
shear traction boundary conditions at the top and bottom surfaces of the laminate in order
to maintain the number of degrees of freedom at ﬁve. This class of theories appears to
have an ideal combination of accuracy and efﬁciency, making them well-suited for use in
computational simulations. However, these theories have the unfortunate shortcoming that
the transverse deﬂection degree of freedom wo is required to be C1 continuous, so that
Hermitian-type interpolation of wo must be used within the ﬁnite element models [33,35].
Thus additional rotational degrees of freedom (gradients of wo) are present in the ﬁnite
element models, making it inconvenient, if not impossible, to implement the ﬁnite element
models based on these theories into commercial ﬁnite element software packages that
allow only six degrees of freedom -- three translations and three rotations.

Averill proposed a generalized form of the ﬁrst-order zig-zag theory [38] and a gener-
alized form of the high-order zig-zag theory [40] for beams to alleviate the requirement of

Cl continuity on the transverse deﬂection. A new variable was introduced along with an

10

associated constraint condition which was enforced by employing an interdependent
(anisoparametric) element interpolation scheme and the penalty method. On the basis of
the theories, C0 two-noded elements were developed. These elements were shown to be
simple, robust and accurate for application to thick and thin laminated beams. However,
these models still contained one additional rotational degree of freedom, and it was found
that for very thick and complex laminates, more reﬁnement through—the-thickness was
rmeded.

Yip and Averill developed a model for laminated beams [41] and plates [42] that com-
bined the beneﬁts of the discrete-layerwise and high-order zig-zag theories. This model
allows the laminate to be subdivided into a number of sublarninates, with each sublami-
nate containing several, even many, physical layers. \Vithin each sublaminate, very accu-
rate high-order zig-zag kinematics are assumed in which degrees of freedom describe
displacements, rotations, and transverse shear stresses (tractions or interlaminar stresses)
at the top and bottom surfaces of the sublaminate. Each ﬁnite element represents one sub-
laminate, and, if cast in the form of a four-noded quadrilateral (for beams) or an eight-
noded brick (for plates), reﬁnement of a model by through-thickness discretization can be
achieved without the use of any special constraints. When only one sublaminate is used
through the entire thickness of the laminate, nodal degrees of freedom are three transla-
tions and two rotations. However, when multiple elements (sublaminates) are needed to
model a laminate, additional shear stress degrees of freedom are present, making the ele-
ment unsuitable for implementation into many of the current commercial ﬁnite element

codes.

11

1.4 Literature Review: Three Dimensional Composite Shell/Plate
Elements

Laminated composites are usually modeled by using plate and shell elements which
are based on plate and shell theories. There have traditionally been three approaches to
model general shell structures in ﬁnite element representation [43,44]:

1) The ‘faceted” model of ﬂat plate element

2) The “curved” shell element formulated on the basis of curved shell theory

3) The “degenerated” isoparametric element

Among the above elements, the Ahmad type [45] degenerated isoparametric elements
where the interpolation functions for rotational and translational d.o.f are independent are
most popular. The degenerated elements, since the thickness of plate or Shell is very small
when compared with other dimensions, are obtained by constraining adjacent thickness-
direction nodes on solid elements to have the same thickness-direction displacement and
merging the top and bottom nodes to a mid-node.

Recently, a number of three dimensional composite plate/shell elements that have a
thickness dimension have been proposed for the analysis of composite structures [47-53].
Such three dimensional elements are classiﬁed as “regenerated,” in contrast to the popular
two dimensional “degenerated” shell elements. This regenerated approach has been
employed successfully, and renders the very important advantage of allowing discretiza-
tion in both inplane and through-thickness directions without imposing any special multi-
point constraints. The regenerated elements must have just 3 translational degrees of free—
dom (u, v and w) to describe the deformations of a plate, whereas the degenerated element

must have at least 5 degrees of freedom (3 translational d.o.f + 2 rotational d.o.f). Thus, as

l2

far as the discretization through-thickness direction, the description of geometry and kine-
matics of deformation are concerned, the three dimensional regenerated element is pre-
ferred over the two dimensional degenerated element. The capability of discretizing in the
through-thickness direction enables one to simulate the separation of the laminated com-
posites due to delamination into two parts in the thickness direction.

Most proposed regenerated 3-D composite elements are similar to continuum elements
that have only 3 translational degrees of freedom per node. Moreover, contrary to the
Structural element cases, in-plane displacement ﬁelds are approximated using the shape
functions, while special functions in some cases are used for transverse displacement to
improve the continuity of transverse stresses at the interface between two adjacent layers.

Adding composite analysis capability to Ausserer and Lee’s isotropic shell solid ele—
ment [46], Kim and Lee [47] presented an 18-node solid isoparametric element using
assumed independent strains based on Hellinger-Reissner principle in an attempt to allevi-
ate locking. However, the transverse strains on this element in are continuous through
thickness direction, so that the continuity of transverse stress is questionable. Kong and
Cheong [52] developed an element incorporating the advantages of the layerwise plate
theory and the equivalent single layer plate theory. Layerwise local Shape functions in the
regions were used where transverse shear stresses are of interest and importance, while a
special through-thickness global-local interpolation is used in the zone where only global
responses are of interest. The compatibility of the displacement between these two region
is satisﬁed by introducing a transition region. Wanthal and Yang [48] proposed a family of
three quadrilateral isoparametric elements based on the layerwise modeling approach to

include the effects of transverse shear and normal stress for the analysis of delamination,

13

holes and joints. The crucial shortcoming of these elements is that the total number of
d.o.f depends on the number of layers in the laminate. Han and Hoa [50] developed a
hybrid formulated plate element by taking three in-plane stresses and three transverse
stresses as the basic variables. The continuity of the transverse stresses at the interface
between two adjacent layers are satisﬁed by introducing a partial stress ﬁeld parameter.
The element is a super element composed of some three dimensional eight noded ﬁnite
elements (layer element), so that the total number of d.o.f does not depend on the number

of layers.

1.5 Literature Review: Mass Matrices and Lumping Schemes

If a structure vibrates freely or excitation frequencies are higher than 1/3 of the lowest
natural frequency of the structure, inertial effects must be included when investigating the
response of the structure. When we try to analyze the dynamic response of a structure
using the ﬁnite element method, the mass matrix accounts for inertia and represents dis-
cretely the continuous distribution of mass in a structure. Thus, in general structural
dynamic response problems, the stiffness matrix is used to describe the elastic characteris-
tics in terms of force - displacement relationships, and the mass matrix is exploited to
deﬁne the inertial characteristics of the structure in terms of inertial force - acceleration
relationships at a ﬁnite number of nodal points chosen to represent the behavior of the
structure. In this case, the accuracy of the analysis is Signiﬁcantly affected by the accura-
cies of the stiffness matrix as well as the mass matrix.

Up to the early 19605, the mass matrix was usually obtained by employing a lumped

parameter method [54], in which the structural mass was physically grouped as a series of

l4

concentrated mass at the nodes. The derived mass matrix is diagonal, and a simple tech-
nique is required to formulate and solve the equations. However, the computed natural fre-
quencies and mode shapes might be greatly different from the exact solution, if a large
number of elements are not used.

A consistent mass matrix was proposed by Archer in an attempt to improve the accu-
racy of the dynamic analysis [55]. The mass matrix was derived employing equations cor-
responding to the Rayleigh-Ritz approach, and the resulting mass matrix is “consistent”
with the actual distribution of mass throughout the structure. Later on, Archer also a con-
sistent mass matrix for Timoshenko beam element. The mass matrix was derived by super-
posing linear inertias on rotary inertias which were determined assuming that the mass
angular acceleration is consistent with the angular motion of the transverse ﬁbers of the
beam [56].

Even though the consistent mass matrix represents a truer distribution of mass
throughout the structure compared to a lumped mass matrix, it results in a coupled system
of equations for explicit structural dynamic analysis as well as a need to construct a global
mass matrix. Thus the solution procedure to solve the system of equations becomes com-
plicated and computationally expensive. In addition, lumping lowers the highest wave
speed in the ﬁnite element mesh, and increases the critical time step for explicit integration
schemes. These facts motivate analysts to ﬁnd rational ways to develop diagonal mass
matrix by lumping appropriately. Several lumping schemes based on a consistent mass
matrix have been proposed.

The easy-to-use and well-known lumping techniques are the row-sum lumping, the

optimal lumping and proportional lumping (so called HRZ lumping) schemes. In row sum

15

lumping, the components in each row are added and lumped on the diagonal using the
property of the interpolation functions that the sum of the interpolation functions is 1. This
lumping scheme sometimes causes negative masses, especially for comer nodes of the
eight-node serendipity element. On the other hand, in the optimal lumping scheme pro-
posed by Fried and Malkus [57], the lumping of mass was accomplished by placing the
location of nodes of an element at the integration point of a quadrature rule. This lumping
scheme is practically useful only for low order Lagrange elements with only translational
degrees of freedom, because this scheme produces a block diagonal lumped matrix for
elements with rotational d.o.f, and zero or negative nodal masses for quadratic and higher
order elements.

A special lumping technique was presented by Hinton, Rock and Zienkiewicz [58].
This lumping scheme is sometimes called as the proportional lumping scheme or the HRZ
lumping scheme. In the HRZ (proportional) lumping scheme, lumping is achieved by cal-
culating the diagonal elements of the lumped mass matrix to be proportional to the diago-
nal elements of the consistent mass matrix While preserving the total mass of the element.
The most important feature of this lumping scheme is in that it always produces positive
lumped masses without regard to the element type, because the diagonal terms in the con-
sistent mass matrix must be positive by virtue of the positive-deﬁniteness property. For
low order elements, the accuracy of this lumping scheme in the prediction of natural fre-
quencies often surpasses that of the consistent mass matrix. For Lagrange linear and qua-
dratic elements with constant material density, the lumped mass matrix computed by the
proportional lumping is the same as that obtained by the row-sum lumping.

The lumping scheme for two-dimensional degenerated general shell elements and axi-

16

symmetric shell elements with rotational d.o.f was developed by Surana [59]. In this
lumping scheme, the rotary inertia terms are obtained by assuming that the masses in
translational d.o.f at the top and bottom surfaces of the plate are rotating with respect to
the 'mid surface, while the mass terms in lumped mass {matrix for translational d.o.f. are
derived by the same procedure as in the HRZ lumping scheme.

It was shown in [60,61] that the reduced integration rule to evaluate the element mass
matrix may result in a singular mass matrix. The existence of a singular mass matrix
implies that there are certain unacceptable modes of deformation associated with the ele-
ment which produce zero kinetic energy. Thus, full Gauss quadrature rule is recommended
to be used to evaluate the element mass matrix.

The the rates of convergence of eigenvalues and eigenfuctions between consistent
mass and lumped mass formulations was studied in [62-64]. Tong, Pian and Bucciarelli
[62] found that the lumped mass formulation suffers no loss in order of convergence for
the membrane and rod cases employing Simple interpolation functions. However, for prob-
lems with higher order equations such as the beam or plate, or for the second order equa-
tions with more accurate interpolation functions, the consistent mass formulation yields a
better rate of convergence.

The effects of negative and zero masses and stability conditions reproduced by optimal
lumping scheme in transient ﬁnite element analyses were investigated by Malkus and Ple-

sha [65], and Malkus, Plesha and Liu [66].

1.6 The Present Study

Reﬁned plate theories based on linear or cubic plus zig-zag kinematics have been

17

developed, and new C0 3-D ﬁnite elements have been developed on the basis of those the-
ories. The in-plane displacement ﬁelds in each sublaminate (for linear zig-zag) or in the
laminate (for cubic zig-zag) are assumed to be piecewise linear or cubic functions and
vary in a zig-zag fashion through the thickness of the sublaminate or the laminate. The
transverse displacement ﬁeld is assumed to vary linearly through the thickness direction.
The zig-zag functions are evaluated by enforcing the continuity of shear stress at each
interface. This in-plane displacement ﬁeld assumption accounts for discrete layer effects
without increasing the number of degrees of freedom as the number of layers is increased.
The rotational variables in the in-plane displacement ﬁeld are eliminated in favor of the
translational variables, the inplane translations at the top and bottom surfaces of the sub-
laminate or laminate, in order to facilitate the development of convenient ﬁnite element
models. For the cubic zig-zag model, high order rotational variables are removed by
employing zero shear traction conditions at the top and bottom of the laminate. Moreover,
since the derivatives of transverse deﬂections appear in the displacement ﬁeld, their sec-
ond derivatives will be present in the Strain energy functional, so Cl continuity of trans-
verse deﬂection is required. Such a requirement is alleviated by introducing new rotational
degrees of freedom. The constraints between the derivatives of transverse deﬂections and
new rotational degrees of freedom are imposed via the penalty method. To remove Pois-
son’s ratio stiffening effect, the transverse normal strain is improved by assuming a con-
stant transverse normal stress through a sublaminate or the entire laminate.

The equations of motion, the essential and natural boundary conditions, as well as the
displacement-based ﬁnite element model are developed from Hamilton’s principle. The

new C0 ﬁnite element is developed with the topology of an eight-noded brick. Each node

18

has ﬁve engineering degrees of freedom - three translations and two rotations.

Inplane displacements and rotational degrees of freedom are approximated by the
bilinear Lagrange interpolation functions. For transverse deﬂection degrees of freedom, an
interdependent interpolation concept developed by Tessler and Hughes [75] is utilized.

-The above interdependent interpolation scheme alleviates the shear locking prOblem
but does not eliminate it totally. The element developed using this scheme still locks in the
very thin regime. Thus, the consistent transverse shear strain ﬁelds as well as the edge-
consistency of the tangential shear strain on any common inter-element edge were
achieved utilizing Prathap and Somashekar’s approach to prevent the element from lock-
ing [76,77]. The element stiffness coefﬁcients are integrated exactly, yet the element
exhibits no shear locking. The element is shown to be accurate, simple to use, and compat-
ible with the requirements of commercial ﬁnite element codes.

Based on the displacement ﬁelds and the ﬁnite element model, a consistent mass
matrix has been obtained. Making use of the consistent mass matrix, lumped mass matri-
ces for explicit structural dynamic analyses have been derived by several lumping
schemes.

The developed elements have been implemented into the explicit ﬁnite element
dynamics code, NEPTUNE [78]. In the code the internal force vector is computed for the

element. Using NEPTUNE, the dynamic analysis of composite plate was performed.

1.7 Organization of the Dissertation

The dissertation will be focused on the development of robust and accurate composite

beam and plate elements for the static and dynamic analyses of laminated composites and

19

sandwich structures. The dissertation is composed of 7 chapters, and the features and con-
tents of each chapter are as follows.

In Chapter 2, ﬁrst order layerwise beam theory and the two dimensional beam element
associated with the beam theory is presented. Numerical results are obtained for cylindri-
cal bending cases and compared with the exact solutions by Pagano [72].

A laminated plate theory with ﬁrst-order zig-zag sublaminate approximations and a
new 3-D ﬁnite element based on that theory are developed in Chapter 3. To alleviate lock-
ing, interdependent interpolation, consistent transverse Shear strain ﬁelds, and edge-con-
sistency of the tangential shear strain will be introduced. To remove Poisson’s ratio
stiffening effect, the transverse normal strain will be improved by assuming a constant
transverse normal stress through each sublaminate. Numerical performance of the current
element is investigated for a composite armored vehicle panel and a sandwich panel with
low and high aspect ratios.

Chapter 4 will introduce a reﬁned laminated plate theory based on cubic zig-zag kine-
matics and the associated element. The discretization in the through-thickness direction is
not allowed in this element, because homogenous shear traction conditions at the top and
bottom surfaces are assumed. The numerical experiments will be carried out by simulating
a damaged composite plate under a double sinusoidal loading.‘

In Chapter 5, on the basis of the developed elements, the consistent mass matrices and
the lumped mass matrices will be derived. Utilizing the derived stiffness and mass matri-
ces, free vibration of a laminated plate will be analyzed.

The structural dynamic analyses of composite structures will be performed using the

proposed elements in Chapter 6. Since the proposed elements have been implemented into

20

NEPTUNE, which utilizes the explicit time integration, the stability condition and critical
time step for the explicit time integration will be also discussed.
The conclusions and recommendations of the current research will be present in Chap-

ter 7.

CHAPTER H
FIRST ORDER ZIG-ZAG BEAM THEORY
BASED ON SUBLAMINATE APPROXINIATION S

2.1 Introduction

In this chapter, a new beam ﬁnite element based on a new discrete-layer laminated
beam theory with sublaminate ﬁrst-order zig-zag kinematic assumptions is presented and
assessed for thick and thin laminated beams. The model allows a laminate to be repre-
sented as an assemblage of sublaminates in order to increase the model reﬁnement through
the thickness, when needed. Within each sublaminate, discrete-layer effects are accounted
for via a modiﬁed form of DiSciuva’s linear zig-zag laminate kinematics [33], in which
continuity of interfacial transverse shear stresses is satisﬁed identically. In the computa-
tional model, each ﬁnite element represents one sublaminate. The ﬁnite element is devel-
oped with the topology of a four-noded rectangle, allowing the thickness of the beam to be
discretized into several elements, or sublaminates, if necessary, to improve accuracy. Each
node has three engineering degrees of freedom, two translations and one rotation. Thus,
this element can be conveniently implemented into general purpose ﬁnite element codes.
The element stiffness coefﬁcients are integrated exactly, yet the element exhibits no shear

locking due to the use of a consistent interdependent interpolation scheme. Numerical per-

21

22

formance of the current element is investigated for an arbitrarily layered beam, a symmet-

rically layered beam and a sandwich beam with low and high aspect ratios.

2.2 Kinematic Assumptions and Formulation of the Beam Theory

It is assumed that a laminate is composed of N perfectly bonded layers, with each layer
being of independent thickness and having independent Stiffness properties. The laminate
will be modeled as M sublaminates (1 S M), with each sublaminate containing Nm layers

such that:

N = 2 Nm (2.2.1)

In order to facilitate the development of the ﬁnite element model, all expressions in the
following derivation pertain to the mth sublaminate, and the sublaminate number designa-
tion is omitted for brevity hereafter.

The constitutive equation for a plane stress state in the kth layer of the mth sublaminate

can be written as:

 

 

 

 

 

 

' (If ' I i (k)‘

911 CY? CE? 0 811
(k) _ (k)
(k) (k) (k)

913, _ 0 0 C55_ 1713 J

where the principal material axes are assumed to coincide with the X,Z axes, and the sub-

scripts 1,3 correspond to the X,Z or x,z directions, respectively. The global X- and local x-

23

directions coincide and are taken along the length of the beam. The global Z-axis is taken
perpendicular to X and has its origin at the bottom of the laminate, while 2 is a local sub-

laminate coordinate in the direction of Z with its origin at the bottom of the sublaminate as

shown in Figure 1.

 

 

 

 

 

 

 

 

 

 

ll 2, 2
2M
i b ‘
2’2 ZN2 (2)
_ _ EuﬁarﬁnLTe? — ‘
z2 (2)
z1 \ Z1(2)
V/J/J/j/J/Jz/M/Jx z, <2>
Sublaminate l
a) xxxxxxxxxxxxxxxxxXxxxxxxxxxxxxxxxxxxxx

 

 

X, x
Figure 1. Schematic of sublaminate and layer divisions
in the laminated beam theory.

Coefﬁcients of the material stiffness matrix are given by:

 

(k) (k)
C(k) = Ei C(k) = E3
11 _1 (k) (k) ’ 33 1 (k) 00
"V13 V31 ‘V13V31 223
(“E“) (10E ( ' ' )
C(k) = V13 1 = V31 3 COO _ 00¢)
13 1 (k) (k) 1 (k) (k) ’ 55 ‘ 13
‘V13 V31 ‘V13V31

where E,- are the Young’s Moduli, and VU are the Poisson’s ratios in kth layer, and Gij are

shear moduli.

The strain - displacement relationship for a plane stress state can be written as:

24

(k) (k) (k) (k)
(k) 3“ (k)_auz (k)_auz +3“;

811 = 3;)“ , 33 —a—Z- , 13 —-3; a—z (2.2.4)
The assumed displacement ﬁelds are initially assumed in the following form:
k k- 1
u; )(x, z) = ub(X) + 2W + 2 (z - z,)§.-
i= 1 (2.2.5)

ugk)(x, z) = wb(x)(l — IE1) + w,(x)(’§l)

where h is the thickness of the sublaminate, ub is the axial displacement at z = 0, u! is the

rotation at z = 0, and wb and w, are the transverse deﬂections of the bottom and top sur-

(’0

x varies In a piece-

faces, respectively, of the mth sublaminate. Thus, it is assumed that u

(’0

Z varies linearly through

wise linear fashion through the thickness of a sublaminate and u
the thickness.
The parameters Q in (2.2.5) are eliminated by enforcing the continuity of shear stress

at each interface. The shear stress continuity condition between the kth layer and the

k+1th layer can be expressed as follows:

"" - 0““) at z: z, A (2.2.6)

21:" 2x

Substitution of the constitutive equation (2.2.2), the strain - displacement relationship

(2.2.4) and the assumed displacement ﬁelds (2.2.5) into Eq. (2.2.6) gives:

25

 

 

 

 

 

(k) (k) (kt-1) (k+1)
= G
013 sz 2:th 13 72x z=z,
k
dw z dw z
_ (k+l) ' __b ___5 __i_k
‘ G” {W 2§‘+dx (1 h)+dx (12)} (2.2.7)
1': l _
Z - 2*
U” 1) (k)
= 013 (72x z=zi+§k)
From the above (2.2.7), ﬁk is:
(k)
g G 13 (k)
k G(k+1) u z = z,

13

(k) k 1 (2'2'8)
Gl3 — (1W1, z), dw, zk

- [gm—r)” 1%“ 251’? (“FIVE (T)
13 i: 1
Solving the recursive relationship for 1;, ﬁnally 7;,- is found to be:
,. A Wb A dwr
g,- = aiw + b9; + C'E (22.9)

where

 

 

6111;) Z. 1-1
. _ , ..
b,- -( (M)— l][l—z+ 2 13,] (2.2.10)

From Eqs. (2.2.9) and (2.2.10) it can be seen that ﬁt. depends on the ratios of shear proper-

26

ties between two adjacent layers and the shear deformation in each sublaminate. If either
of these quantities is small, then discrete-layer effects will also be small.
The modiﬁed displacement ﬁeld is given by:
k ' 1 dw

.. dw
i=1 (2.2.11)

“(2k)”, 2) = wb(x)(l - i) + MOE)

In order to facilitate the development of convenient ﬁnite element models, the rota-
tional variable 1|! is now eliminated in favor of the translational variable u,, the inplane
translation at the top surface of the sublaminate. Thus, rather than describing the inplane
displacement ﬁeld by a translation and rotation at one point, it can more conveniently be
described here by the translation at two points. Introducing the new variable up u! can be

obtained as follows:

1 13de _dw,
W - 5|:“t-“b— 3; _CE] (2.2.12)

where ab and u, are the axial displacement at the top and bottom of each sublaminate, and

can be described mathematically as:

= u (I) u = u (N) (2.2.13)

u
b 2: 2:0 I x z=h

The coefﬁcients a, I.) and E are deﬁned as:

27

N—l
a = h+ 2 (h—zﬁdi
i=1
N-l
B = 2 (h—z,)13,. (2.2.14)
i=1
N-l A
Z‘ = 201—200,-

i=1

(k)

x can be written as follows:

Therefore, u

x ’ dx dx

(2.2.15)
k- 1 .

* dw dw dw dw

a - b - r - b .. t

+.21(Z-Zi){5(ut“ub-b-d—x- — 2;)4'b-d—x- +025}
1 =
. . . . dwb th . .

Since the derivatives of transverse deﬂections, a and E , appear in the displacement

ﬁeld, their second derivatives will be present in the strain energy functional, so Cl continu-
ity of w), and w, is required. It is desirable to alleviate such a requirement by introducing

new rotational degrees of freedom as follows:

_ = 91»— : 6, (2.2.16)

The above relationship still must be enforced within the model, as discussed later. The

ﬁnal form of the axial displacement ﬁeld can now be written as:

28

k—l k—l
“(x10 = ub{l — aoz- 2 (Z “ 290,-} ‘I' ut{aoz + Z (Z —' 290,-}
"1 (2.2.17)

i=1

k-l k-l
+ 6,,{boz + 2 (z — 20171} + 9,{coz + 2 (z - mg}

i=1 i=1

where z,- is the coordinate of the ith interface in the sublaminate. The coefﬁcients

ao, b0, CO, at, bi and ci are given by:

1 ..
b E _ x 01-. .. _ ..
0 — “'5 = ‘aob, b. = bi—Eb = bi—aoba' (2.2.18)
- Z- — _ - A ai_ _ A _ A
CO — -;_'1 — -aoc, Ci — i-EC - Ci—aocai
Thus, the strains are deﬁned as:
du k-l du k-l
k b
3(11) = d; {1- aoz- 2 (z — 2,.)ai} +E'{aoz + 2 (z - 2001}
1‘: 1 i=1
+3; {boz + 2 (z — zi)bi} +E'{coz + 2 (z — 2091
i=1 i=1
(k) Wb ‘ “’1 (2.2.19)

 

{ZHXHXHX}

i=1 i=1 i=1 i=1

The governing equations, the essential and natural boundary conditions, as well as the

displacement-based ﬁnite element model can be derived from the principle of total poten-

29

tial energy. Again, attention will be limited to a single sublaminate. The potential energy
. functional is modiﬁed here to include the imposition of the constraints (2.2.16) via the
penalty method [68-70]. The modiﬁed total potential energy for the mth sublaminate can

be deﬁned as follows:

XX exx

k=1v. (2.2.20)
dw 2 dw 2
I __b I __‘
+1329!) dx ) dx+j:2(e, dx)dx

where N", is the number of layers within the mth sublaminate, Vk represents the volume of

“(m)_ (108(k) (k) (k) (k) (1‘)
—N2 IZIG + 0'2: 822 + 021sz lde-j:(qrwr + qbwb)dx

the kth layer in the sublaminate, 'y is the penalty parameter, and qb(x) and q,(x) are the dis-
tributed transverse loads applied on the surfaces of the sublaminate (if these surfaces also
represent surfaces of the beam). As 7 is assigned successively larger values, the con-
straints of (2.2.16) are satisﬁed more exactly in the least squares sense. Substituting Eqs.
(2.2.2), (2.2.19) into (2.2.20), pre-integrating the internal strain energy terms through the
thickness and through the width of the beam, taking the variation of It?) , and integrating

by parts, (2.2.20) can be rewritten as:

30

L 2
0 dx
dr, d9 d w de dwb
+8WIN ”q, dx +7 2;— d 2 +66b Sb-Z—x'x b+Y(eb—E )
x
dM dw i
+69,{S, 7; " + 7(9, -2—x-')}]dx (2.2.21)

d

+ {a — 7(1), — g‘)}8w,]

where N xb, N n, M xb, M n, N zb, N U, Qb, Q” Sb, 5,, Tb, T, are the beam sublaminate

L

=0
0

 

stress resultants deﬁned as:

m k-l k—l
k k
Nx(b,t) = Z I: six)b( ){[1—aoz— 2 (Z-Zi)ai], [002+ 2 (Z—Zi)ai]}dx
k=1 "“ i=1 i=1
Nil k k k-l k-]
k
wa.” = r Gig)“ ){[boz+ 2 (Z-Zﬂbi], (COz-F 2 (z—zi)ci]}dx
k =1 z,_, i=1 i=1
N" 1 1
* (k) (k)
Nz(b,i) = £._.Ozz b {(‘er (2)}dx
[:1 k 1 k 1 (2.2.22)
" t (k) (k)
Q(b,r) = 2i: 02: b {[ 0' 2 01'], [40+ 2 01]}dx
k=1 "“ i=1 i=1
N... (k) (k) k-l k-l
k
S(b,l) = 2g 62X b {[bo'I’ bl]’ [Co'I’ 2 Cl]}dx
k=l *" 1:] i=1
NM
._ " (k) (’0 Z Z
- Emu" {(1 7.1 (71)}...

31

Substituting Eqs. (2.2.2), (2.2.4) into (2.2.22), the sublaminate stress resultants are

described in terms of the generalized displacement degrees of freedom by:

N A A0 dub Aa B du‘ Ab Ba dOb Ac Bb d9!
x(b,r)—( ’ )2;+( 9 )2;+( , )2; +( , )2;
+ (C, C“)w,,—(C, C”)w,
du du b d9 b (19
Mx(b.:) = (Ab: Ac)?§b+(D, 3b13';t+(Da,D )EbHD .185!
+ (Cb: Cc)Wb-(Cb, CC)W’
dub dur b b deb c (19,
N (Cl-C)??? +(Ca,-Ca)-d-; +(C ,—C ).d_ +(C ’-CC)d—x

z(b, r) " x

+ (G, -G)Wb + (-G, G)w,
a b b (2.2.23)
QW, = (H, -H)ub+(-H, H)u,+(H”, —H )0b+(H ,-H )6,

dw dw
c c b d d t
+(H,—H )a +(H,—H )d—x
S“, 0 = (Ha, Hb)ub + (—H", —H”)u, + (I, I“)e,, + (1", J)e,
dw dw
b {a b c b t
”1’ )3 ”1’15?
T — (H‘ H")u +(—H" —H")u +(I” I‘)9 +(J“ J”)e
(b, r) " ’ b 3 r ’ b 3 r

dw dw
+ (K. K537” + (K“. mg;

where the coefﬁcients in (2.2.23) are deﬁned as follows:

32

k-l \ - k-I
(A, A“, Ab, AC) = 2N1]: C(lk)b(k)[l —a Oz— 2 (z—zﬁai ( l—aoz— Elk-2.1611]
' i=1 )

b i=1

k-l k-l - P k—l
,[aoz + Z (2 — 200,], [boz + 2 (z — z,)b,- coz + 2 (z - 2,)ciDdx

i=1 1:] J '- i=1

9

 

 

N»l k—l k—l
(B, B“, B”) = 2 f‘ C‘,"’b“"[a a0z+ 2 (z-z,)a,.]-[[aoz+ Z (z—z,)a,],

Z
"‘ i=1 i=1

k=1

k— 1 k-l
[boz + 2 (z — 2,)bi], [coz + z (z — z,)ci]]dx

i=1 i=1

_ N, 1-1
(C, Ca, Cb, Cc) = 2 ft C(lg)b(k)(— It). [[1 —aoz— 2 (Z-Zi)at]’

Z
k=l "' i=1

k-l k-l k-l
[aoz + 2 (z — 21M], [boz + 2 (z - 2,)bi], [coz + 2 (z - Iraq-Dru

i=1 i=1 i=1

i=1 i=1

k- l k-l
(D, D, Db ) = [‘21]: C(lk)b(k)[b0z+ 2 (z- z )b ][[aoz+ 2 (z—z,)a,-:|,
‘ (2.2.24)
k-l k-l
[boz + Z (z — zi)b,], [coz + 2 (z - z,)c,-Ddx

i=1 i=1

k- 1
(F): Nil]: C‘,"’b""[coz+z(z- z,)c,-]-[coz+2(z z)c,)dx

i=1 i=1

(G) = N2 J: |C§k)b(k)(%)o (ﬂax

k=l

N, k-l k-l
(H, Ha, Hb, He, Hd) = 2 J41 C(Ic)b(k)[_aO _k2 11“] 6-00 -2 ai], [b0+ 2 bi],

k=l i=1 i=1

[so+"i'c.-]11-ziiziidx

i-l

k-l k-l
(1,1“,1”,f)- = 2 j: Cg’;’b""[b0+ 21,17]. ﬂbw 2 b,], [c0+ 2 0.],

i=1 i=1 i=1

[13112])“

33

Mg cgk>i<k>[.,.“ 2 q.) [[219 -5], [2]}.

i=1 i=1

:r W1) (II-ziiziidx

Z
1 k- 1

Off) -

I
211N3M2

R-
2"

(L) j: lcgk’b‘k’m (2)213:

k l

The governing equations and boundary conditions for the mth sublaminate are obtained
from (2.2.21). If the entire laminated beam is modeled using a single sublaminate, then the

beam equilibrium and boundary conditions are:

dN b
Dub: —E;x +Qb = 0
(IN
5a, —‘—,;“+Q, = 0
dT d6 d w
8wb: sz—ab" qb(x x)+y[d—b—d_:.b] =
2x (2.2.25)
dT, d9, dw,
5W,: N21- d—x- -q,(X) +7 d—x —-d—x7 '-

deb dwb

dM dw
66,: -$ "+S,+y(9,—a') = O

34

Essential B.C Natural B.C.
“b Nxb
U, er
9b M xb
9, M x, (2.2.26)
dw,,
Wb Tb - Y(9b ‘ 3; )
dw,
WI T: " 7(91 - 2'; )

2.3 Finite Element Formulation

From (2.2.23), it can be seen that the primary variables are functions of x only, thus we
need to consider only one dimensional shape functions for the ﬁnite element model. The
most obvious path of development would be toward a simple two-noded beam element
with six degrees of freedom per node. However, such an element would not fall within the
constraints of having three translational and three rotational degrees of freedom. Even
more importantly, a major advantage of the present formulation would be lost if a tradi-
tional two-noded element were developed. Because the degrees of freedom represent
quantities at the surfaces of a sublaminate, it is most convenient to cast the ﬁnite element
model in the form of a four-noded element with three degrees of freedom per node as
shown in Figure 2. This element is classiﬁed here as a “regenerated,” element. This
approach has been taken by other investigators as well [41,42,67], and has the very impor-
tant advantage of allowing discretization in both inplane and through-thickness directions
without the need for any special multipoint constraints. Therefore, the element can be used

to predict the global response of a laminate using only one (or a small number) of ele-

35

ments through the thickness, and local effects such as interlaminar stresses can be cap-
tured by reﬁning the mesh in the thickness direction near the interface(s) of interest.
Inplane displacement and rotation degrees of freedom are approximated by the linear

Lagrange interpolation functions:

2 2
“17,: = Z ubi,riPi 9b.: = Z ebmipi (2-3-1)

i=1 i=1

where
P1 = %(1—§) P2 = %(1+§) (2.3.2)

and g is the element natural (local) coordinate.

The constraints between the derivatives of transverse deﬂection and the rotation degree
of freedom is imposed in the strain energy functional via the penalty method. As men-
tioned, the rotation degrees of freedom are approximated by the linear Lagrange interpola-
tion functions. Thus, if the linear Lagrange interpolation functions are employed for the
transverse deﬂection degrees of freedom, the constraints in the Strain energy functional
enforce the relationships between the linear variations of rotation degrees of freedom and
the constant part of the deﬂection approximation, because the ﬁrst derivative of a linear
function is constant. This causes locking of the element. In oeder to alleviate locking,
transverse deﬂection degrees of freedom are approximated by an interdependent interpola-
tion concept developed by Tessler and Dong [71]. Initially, the transverse deﬂections, wb

and w,, are approximated by quadratic Lagrange interpolation functions.

36

3
“’12.: = X wbi,riNi (2.3.3)

(—

where
N1 = %(§2-§), N2 = %(§2+§), N3 = i—gz, —1 3:51 (2.3.4)

The virgin element is shown in Figure 2, and has a mid-node at the top and the bottom sur-

faces.

(a) virgin element
at, Wt, 6! wt “19 Wt! 9!

 

 

 

 

 

“in WI» 9b “’1: ub, wt, 9,,

(b) constrained element
u,, Wt, 9: at. W. 9:

“b, va eb ub’ Wb. 6b

 

 

 

Figure 2. Element topology and nodal degrees of freedom.

The resulting mid-side deﬂection degrees of freedom can be eliminated by imposing the

37

high-order part of the constraints in (2.2.16):

9

d dwot
dx “—2; J

= 0, or = 1(bottom), 2(top) (2.3.5)

Substituting Eqs. (23.1-23.4) into (2.3.5):

d dwa _ 0
dx “-5 ) -
2
1d0a 1 d we,
- “5 Jidt’
2 2 2 (2.3.6)
z(dPl sz ) 4 le sz dN3
=——0 +—6 -——w +—w +—w
he d1"; “1 d5 “2 h: d§2 “1 d§2 012 (1:2 a3
"9011 +9112) 4
= —— -—(w +w —2w )
( he kg a1 012 013
where J,is line Jacobian and he is the element length.
If we solve for we,3 in (2.3.6), we,3 is:
1 he
W013 = -2-(wa] +wa2) +§(9a1 -9a2) (2.3.7)

The substitution of we,3 into (2.3.3) results in the new interpolation functions for trans-

verse deﬂection degrees of freedom:

2 2
“’12,: = Z Wbr, n‘P 1+ 2 abruNr (2.3.8)

i=1 i=1

.. h x h
WhereN1= -85(l-§2), N2 = -§e(1-§2)

38

This approach improves the element behavior very effectively without introducing
additional degrees of freedom. Since the interdependent interpolation scheme identically
satisﬁes the linear part of the constraint in (2.2.25), the penalty parameter need only
enforce the constant part of the constraint. It has been found that a penalty parameter on
the order of the largest material shear stiffness is all that is needed to impose the constraint
sufﬁciently, so the possibility of numerical ill-conditioning due to large penalty parame-
ters is eliminated. Further, shear locking is obviated so that an exact order of integration
can be used.

An interesting property of the present model is uncovered when it is used to model iso—
tropic beams. Note that for the case of an isotropic beam the displacement ﬁeld in (2.2.15)
does not contain any derivatives of the transverse deﬂection degrees of freedom, since the
coefﬁcients of those terms vanish. The theory thus requires C1 transverse deﬂection
degrees of freedom for laminated beams in which any two adjacent layers have different
transverse shear stiffnesses, but requires only C0 transverse deﬂection degrees of freedom
for isotropic beams. As formulated, the present ﬁnite element is CO, with the constraints in
(2.2.16) imposed via the interpolation scheme and the penalty method. When used to
model isotropic beams, the current model yields predictions almost identical to those of
the one-point-integrated four-noded isoparametric plane stress element. The present ele-
ment may also be viewed as a planar element with drilling degrees of freedom, though the
rotations deﬁned in (2.2.16) are in terms of only deﬂection degrees of freedom. The
present element is designed speciﬁcally for the analysis of beam bending, not for general

plane elasticity problems.

39
2.4 Numerical results

The accuracy of the proposed laminated beam theory and the appropriateness of the
ﬁnite element model were assessed by simulating the response of three different simply-
supported laminated beams subjected to a transverse sinusoidal loading as shown in Fig-
ure 3. The elasticity solution of Pagano [72] was suitably modiﬁed to provide comparison
with the current beam results. The material properties used in the analyses are listed in

Table 1, and the lamination schemes studied are deﬁned in Table 2.

p = sin(1l:x/L)

 

 

 

 

ll
1

Figure 3. Simply supported laminated composite beam
subjected to sinusoidal loading.

40

Table 1. Material properties used in the analyses.

 

 

 

 

 

 

 

1. E11 E33 G13 V13 J
Material 1 1.0066 1.0066 0.2066 0.25
Material 2 25.066 1.0066 0.5066 0.25
Material 3 33.066 21.066 8.0066 0.25
Material 4 32.666 10.666 8.2166 0.10
Core 50.063 50.063 21 .763 0.15

 

 

 

 

 

 

 

 

For a given location along the length of the beam, explicit expressions for the through-
thickness variation of stresses are provided by the present model. Along the beam axis,
stresses were evaluated at the centroid of each element, where stress is super-convergent
[73]. The bending stresses presented in the forthcoming discussion were evaluated at the
center of the element that is nearest the center of the beam. The transverse deﬂections are
obtained at the top surface of the center of the beam where load is applied, while the in-
plane displacement is obtained at the end of the beam where x = 0. Due to the symmetry of
the problem, only half of the beam was modeled. Unless noted otherwise, ten elements
were used along half the length of the beam, and the aspect ratio of the beam was set to
four.

An aspect ratio of four was chosen in the present examples in order to highlight the
ability of the present model to obtain accurate results for very thick laminates via through-
thickness reﬁnement of the mesh. Analyses of beams with the same lamination schemes
but with an aspect ratio of ten were also performed. In most of these cases, the predictions
using one sublaminate were nearly indistinguishable from the exact elasticity solution.

These results are omitted for brevity.

Table 2. Ply stacking tables for the three laminates analyzed.

41

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Laminate Layer Material Thickness
4 j

1 1 2 0.125
2 1 0.125
3 2 0.125
4 1 0.125
5 1 0. 125
6 2 0. 125
7 1 0.125
8 2 0.125

2 1 2 0.300
2 1 0.200
3 2 0.150
4 3 0.250
5 1 0.100

3 1 1 0.010
2 4 0.025
3 2 0.015
4 1 0.020
5 2 0.030
6 Core 0.800
7 2 0.030
8 1 0.020
9 2 0.015
10 4 0.025

 

 

 

 

H
#

 

 

0.010

 

 

42

Laminate 1:

The ﬁrst example is denoted as Laminate 1 and has the lamination sequence (0/90/0/
90),. In Figure 4, normalized predicted center deﬂection and normalized predicted in-
plane stress is plotted as a function of the number of elements used in the mesh. In this set
of results, Llh=100, and only one sublaminate was used. The error was less than 1% when
seven elements were used in half the length of the beam. The inﬂuence of aspect ratios on
transverse deﬂection is shown in Figure 5. Nine different aspect ratios were considered to
check a wide range of beam bending behavior (Uh: 4, 8, 10, 20, 30, 50, 100, 500, 1000,
and 10000). When one sublaminate was used, the error in predicted deﬂection was less
than 5 percent even for the case Llh=4. The use of two sublaminates of equal thickness
resulted in more accurate predictions, but the trends were almost the same as those when
one sublaminate was used. When eight sublaminates (one element per layer) were used,
the error in the predicted deﬂection was less than one percent for all aspect ratios consid-
ered.

In Figure 6, the variation of in-plane displacement through the thickness of the lami-
nate is plotted for a beam with aspect ratio of four. The exact distribution of in-plane dis-
placement is approximately linear within each layer, and globally piece-wise linear and of
zig-zag pattern due to the difference in the shear moduli of adjacent layers. For the model
containing only one sublaminate, the difference between the present model predictions
and the elasticity solution at the top surface was about 13 percent. However, when eight
sublaminates were used, the numerical predictions were in excellent agreement with the
elasticity solution, with error of about 2.75 percent. The stress distribution through the

thickness is illustrated in Figure 7. It can be observed that the prediction of in-plane stress

43

()3.Jr converges rapidly to the exact solution as more sublaminates are used.

Laminate 2:

Laminate 2 is a random ﬁve layer laminate in which the thickness and material proper-
ties vary considerably from layer to layer, a challenging problem for any model to solve
accurately. In this case, most predicted results are locally relatively poorer than those for
Laminate 1 which has a symmetric lamination sequence. Predictions of normalized trans-
verse deﬂection versus aspect ratio are represented in Figure 8. It can be observed that
when the aspect ratio is greater than 20, a one sublaminate approximation gives quite
accurate estimates of center deﬂection. However, ﬁve sublaminates (one element per
layer) are needed to obtain accurate predictions of deﬂection for an aspect ratio of four.
The in-plane displacements predicted by the present model and by elasticity are plotted
versus the thickness of the beam in Figure 9. Even though the overall prediction of in—
plane displacement by the current model agreed well with the elasticity solution, a rather
large error in the prediction of in-plane displacement occurred in the ﬁrst layer of the
beam when either one or two sublaminates were used. This is an example of the shortcom-
ing of ﬁrst-order zig-zag theory in modeling very unsymmetric laminates. However, the
predictionusing ﬁve sublaminates was very close to that of the elasticity solution. The pre-
dicted distribution of in-plane stress as a function of thickness is plotted in Figure 10.
Overall, the predictions of the current model are very good, with the maximum error

occurring at the top of the laminate.

Laminate 3:

The last example is a sandwich panel, denoted here as Laminate 3, which consists of a

44

thick core between two sets of ﬁve face sheets. As can be seen in Table 2, the core occu-
pies eighty percent of the thickness of the beam, while each set of face sheets contains ﬁve
layers (with a lamination sequence very similar to that of Laminate 2) and occupies ten
percent of the total thickness. Three types of through-thickness discretization were used in
the analyses: (1) a single sublaminate, (2) three sublaminates -- one for the core and one
each for the face sheets, and (3) ﬁve sublaminates -- three for the core and one each for the
face sheets.

The predicted normalized deﬂections for various aspect ratios are plotted in Figure 11.
When the aspect ratio was greater than 4, only one sublaminate was needed to yield accu-
rate results. In Figure 12 the predicted through-thickness distributions of in-plane dis-
placement are illustrated. All predictions were very accurate in the upper and lower face
sheet regions, but there were signiﬁcant errors in the core region when only one or three
sublaminates were used. These errors were due to the very large transverse shear strains
and transverse squashing in the core of this thick sandwich beam. The nonlinear variation
of inplane displacement in this region could not be captured accurately by a single sub-
laminate in the core region. For thinner beams (i.e., L/h greater than about 10), the inplane
displacements throughout the thickness are predicted satisfactorily using a single sublami-
nate through the entire thickness of the laminate.

The bending stress predictions shown in Figure 13 are so accurate that it is difﬁcult to
distinguish between analytical and ﬁnite element results. To assess these more clearly, the
stress distribution in the upper face sheet region is enlarged and shown in Figure 14. It can
be observed that the ﬁnite element model using only one sublaminate predicts the stress

distribution in the face sheets very accurately.

45

2.5 Summary

An improved ﬁrst order zig-zag theory and an associated ﬁnite element model are pre-
sented. The model allows the representation of a laminate as an assemblage of sublami-
nates in order to increase the mode] reﬁnement through the thickness, when needed.
Within each sublaminate, an accurate ﬁrst-order zig-zag sublaminate kinematic approxi-
mation is made which minimizes the need for multiple sublaminates for many problems of
practical interest yet easily accommodates through-thickness reﬁnement for predicting
local variations of stress and deformation. The ﬁnite element model is “regenerat ” in the
form of a four-noded planar element with three degrees of freedom (two translations and
one rotation) per element. Thus, it suitable for implementation into commercial ﬁnite ele-
ment codes. Numerical results demonstrate that the current model is accurate, efﬁcient,

and robust for analysis of a wide variety of thick or thin laminated beams.

46

 

 

 

 

 

 

1.005 . T . . . . T
1 L /—— on] 6x1: elasticity
—————————— *—-—————————____
—ax
0.995 - ‘ _
W/W elasticity
0.99 - -
l
I
0.985 - , _
l
1
at:
0.98 l 1 1 1 , I 1 l
0 10 20 30 40 50 60 70 80
Number of Elements

Figure 4. Variation of normalized center deﬂection and maximum
bending stress versus number of elements in the mesh
(Laminate 1, Llh=100).

W / Wexact

47

 

 

 

 

 

 

 

 

 

 

 

1,02 ....., . 2 ..... .j . . .....f,
1.01 - 2
1e _ _ _
,,. 3K ---- "* ‘ﬁjZ'Z'_-_-_-_‘:* _-.-_-*-.': ::'-".-_'_ _ _ -i
#3:: 4s—
0. ..
0.98 _
0.97 -
— 1 Sublaminate
0.96 - - 2 Sublaminates 1
-- - - 8 Sublaminates
0.95 -
094 .....i . . W; . . ......i 4 J
101 102 103 10‘

ASPECT RATIO L / h

Figure 5. Variation of normalized center deﬂection
versus beam aspect ratio (Laminate 1).

§“IN

48

 

 

   

 

 

 

 

 

 

0.9 - F
0.8 ~ ..
0.7 - _
0.6 - _
0.5 - -
0.4 - -
0.3 _ — EIastncuty q
----- 1 sublaminate
02 _ - - 4 sublaminates _
-- - - 8 sublaminates
0.1 - s
-1.5 -1 -0 5 0 0.5 1 1.5
x 10'6

Figure 6. Inplane displacement versus normalized thickness
coordinate (Laminate 1, Llh=4).

§“IN

49

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 v I ‘
0.9 -
0.8 r
0.7 ~
0.6 ~
0.5 -
0.4 -
— Elasticity
0.3 - ‘ ----- 1 sublaminate
0 2 _ - - 4 sublaminates
- -- - - 8 sublaminates
0 5 10 15 20

 

Figure 7. Bending stress versus normalized thickness coordinate
(Laminate l, Llh=4).

W / W exact

50

 

0.95 -

0.85 r

 

 

 

 

—- 1 Sublaminate
- - 2 Sublaminates
-- - - 5 Sublaminates

 

 

 

A I. L I 1 4L j A I I I l I A I A l I I A l J

 

 

0.8

10‘ 1o2 10
ASPECT RATIO Uh

Figure 8. Variation of normalized center deﬂection
versus beam aspect ratio (Laminate 2). ‘

1O

51

 

0.9

l

0.8

0.7

I

0.4 *-

0.3 r

0.2 r-

0.1 r

 

 

  
     

— Elasticity

----- 1 sublaminate
— - 2 sublaminates
-- - - 5 sublaminates

 

 

 

 

 

Figure 9. Inplane displacement versus normalized thickness
coordinate (Laminate 2, Llh=4),

3'IN

52

 

0.6 -

0.5 r

0.4 -

0.3 %

0.2 '-

0.1-

 

 

—— Elasticity

----- 1 sublaminate
- - 2 sublaminates
-- - - 5 sublaminates

 

 

 

 

 

 

 

 

 

15

Figure 10. Bending stress versus normalized thickness coordinate
(Laminate 2, Llh=4).

20

W / Wexact

53

 

1.02L

’.—.-_o—-_o_‘_c-c-g_.--_-.-n-.—-~.—-o-_--.-.-_...—-o-—O—°—'-
.

‘

 

.0
(D
on

0.96

 

— 1 Sublaminate
- - 3 Sublaminates
-- - - 5 Sublaminates

 

0.94

 

 

 

 

 

 

0.92 A l A AA] A A A A A A AA] A ##A A A AA]
10 10 10
L/h

Figure 1 1. Variation of normalized center deﬂection
versus sandwich beam aspect ratio (Laminate 3).

1O

3‘IN

 

54

 

 

 

 

 

 

 

1 " T 1 l T r I
0.9 '- :1 _
.3/ -
.‘ / !
0.8 '- C. / ! —
‘_ Elasticity '/ / l
0.7 _ '''' 1 sublaminate // _I q
- - 3 sublaminates ~/‘ /
0-5 ” '- - - 5 sublaminates ‘3‘, 1‘ _
,- /
.- / I.
0.5 ” ’ .l‘ 4
’ I
I I.
0.4 1- .3., / -
.‘ /
."/
0.3 *- -‘ / . q
/ I
. / ‘1'
0.2 - ./ /’.’ _
:7 '

0.1 ~ _

o l 1 ¥ 1 l J ‘
'4 ‘3 '2 '1 0 1 2 3 4
x 10'6

u

Figure 12. Inplane displacement versus normalized thickness
coordinate (Laminate 3, Llh=4).

55

 

 

 

 

 

 

 

 

 

 

 

 

Figure 13. Bending stress versus normalized thickness coordinate

(Laminate 3, Llh=4).

1 i I l I I . ,4 :x _n .
l
0.8 - i J
l
0.7 ~ _
0.6 - _
5 0.5 - -
h
0.4 - -
0,3 - -— Elasticity _
.- - - 1 sublaminate
0.2 ~ - - 3 sublaminates -
0.1 "' -/J/’ _.
0 I: 1 1 —A L l L L l
-100 -80 -60 —40 —20 O 20 4O 60 80 100

3‘IN

0.95

0.9

56

 

 

 

 

l

 

-- - - 1 sublaminate

!
1 -— Elasticity
l
: - - 3 sublaminates

 

 

 

 

 

 

 

 

Figure 14. Bending stress in face sheets versus normalized thickness
coordinate (Laminate 3, Llh=4).

100

CHAPTER III
PLATE ELEMENT BASED ON
FIRST ORDER ZIG-ZAG SUBLAMINATE THEORY

3.1 Introduction

In this chapter, after estimating the shortcomings and merits of the ﬁrst order zig-zag
beam theory and the beam element, a reﬁned plate theory based on sublaminate linear zig-
zag kinematics and a new 3-D ﬁnite element based on that theory are developed. The in-
plane displacement ﬁelds in each sublaminate are assumed to be piecewise linear func-
tions and vary in a zig-zag fashion through the thickness of the sublaminate. The zig-zag
functions are obtained by satisfying the continuity of transverse shear stresses at layer
interfaces. This in-plane displacement ﬁeld assumption accounts for discrete layer effects
without increasing the number of degrees of freedom as the number of layers is increased.
The transverse normal strain predictions are improved by assuming a constant variation of
transverse normal stress in each sublaminate. In the computational model, each ﬁnite ele-
ment represents one sublaminate. The ﬁnite element is developed with the topology of an
eight-noded brick, allowing the thickness of the plate to be discretized into several ele-
ments, or sublaminates, where each sublaminate can contain more than one physical layer.

Each node has ﬁve engineering degrees of freedom, three translations and two rotations.

57

58

Thus, this element can be conveniently implemented into general purpose ﬁnite element
codes. The element stiffness coefﬁcients are integrated exactly, yet the element exhibits no
shear locking due to the use of an interdependent interpolation scheme and consistent
shear strain ﬁelds. Numerical performance of the current element is investigated for a
composite armored vehicle panel and a sandwich panel with low and high aspect ratios.
Comparison of numerical results with elasticity solutions shows that the element is very

accurate and robust.

3.2 Kinematic Assumptions

In this plate theory, a laminate will be modeled as M sublaminates, and each sublami-
nate is assumed to consist of N,,, perfectly bonded layers which have independent thick-
ness and independent stiffness properties. Mathematically, the total number of layers, N, in

the laminate is written as:

N = 2 N... (3.2.1)

The global Z-axis is taken perpendicular to the inplane X, Ycoordinate axes and has its ori-
gin at the bottom of the laminate, while 2 is a local sublaminate coordinate in the direction
of Z with its origin at the bottom of the sublaminate (see Figure 15). In the following deri-
vation all expressions are related to the mth sublaminate in order to facilitate the develop-
ment a computationally convenient ﬁnite element model. All indices indicating the

sublaminate number are omitted for brevity.

59

Sublaminate 2

 

Figure 15. Schematic of sublaminate and layer divisions
in the laminated plate theory.

The constitutive relations for a three dimensional stress state in the kth layer of the mth

sublaminate can be written as:

 

 

'05:" CW C(k) C(12) 0 0 Cone eff
“i? C“) C“) cg? 0 0 eggs 83;)
021%: CE’S’C‘Z'Q’C‘!) o o cg“, 82"} (3.2.2)
0;? o 0 0 C34“ cg? o 7;?
cg? o o 0 Cf,’;’C§§) 0 7""
~ 61(3)} C(k) C(zk)C(k) 0 0 C212 H733

 

 

 

 

For monoclinic materials, 13 coefﬁcients of the material stiffness matrix are independent,

but for most practical laminated composites, nine coefﬁcients of the material stiffness

60

matrix are independent and the others are obtained through a transformation law [79].

These nine coefﬁcients are given as follows:

 

 

 

(k) (k) (k) (k) (k) (k) (k) (k)
—(k) (k)1‘V23V32 —(k) __ (k)"21 +V31V23 —(k) _ (k)V31 +V21V32
C11 = 1 T’ 12 'El (k) ’ C13 - 1 (k)
A A A
(k) (k) (k) (k) (k) (k) (k)
50:) _ (k)1’V13V31 €213) _ (k)V32 +V12 V31 1‘13) _ E(k)1-V12V21
22 - 2 ——' - 2 ’ - 3 ——
A“) AU!) A“)
—(k) (k) -(k) (k) —(k) (k) (3-2-3)
C44 = 023’ C55 = G13’ C66 = 012
(k) (k) (k) (k) (k) (k) (k) (k) (k) (k)
A = 1“’12V21 ‘V23 V32 "’31 V13 ‘2V21V32V13
(k) (’0
vi. v..
J = j] . A =
(k) (k) "1 1’2’3

, . . . . k . , .
where E?“ are the Young s moduli in x, y, 2 directions and vgj) are the Porsson s ratios
associated with transverse strain in the j direction due to axial strain in the i direction in
the kth layer. From the transformation law of a fourth order tensor, 13 coefﬁcients are

derived using the transformation law of a fourth order tensor:

(’0 —(k)
”kl = a- ajnakpalqcmnpq (3.2.4)

1} rm

where am stands for the direction cosine associated with the material and the global
axes.

The inﬁnitesimal strain - displacement relationships are deﬁned as:

(k) (k) (k)
80:) _ Ex 80:) _ 8“y 80:) _ 5“:
xx " ax ’ yy '3; ’ zz -a—Z
(3.2.5)
(’0 (k) (k) (k) (k) (k)
(k) _ 8“2 at? (k) _ auz +§x (k) _ 8“y + 3“x

”2'23 +82’ “'3'; 82’ Icy-3; ﬁ

61

3.3 Formulation of the Plate Theory

The displacement ﬁelds in a sublaminate are initially assumed in the following form:

k-l
k
u; )(x, y, z, t) = ub(x,y,t)+zw,(x.y.t)+ 2 (24.9%.-

i=1
k-I
k
u; )(X, y, Z, I) = Vb(x, y: t) +ZWy(x, y! t) + 2 (Z—zi)ni

i=1

(’0 Z 2
uz (x, y, z, t) = wb(x, y, t)(1 — h) + w,(x, y, 0(5)

(3.3.1)

where h is the thickness of the sublaminate, ub and vb are the axial displacements in x and
y directions, respectively, at z = 0, w, and my are rotations of the normal at z = 0, and Wb

and w, are the transverse deﬂections of the bottom and top surfaces, respectively, of the

(k)

x

(k)
and uy

mth sublaminate. Thus, it is assumed that u vary in a piecewise linear fashion

)

through the thickness of a sublaminate and uik varies linearly through the thickness.

The parameters ﬁt. and I], in Eq. (3.3.1) are eliminated by enforcing the continuity of
shear stress at each interface [33, 34]. The transverse shear stress continuity conditions at

the kth interface can be expressed as:

(k) _ 0“”) 600 - 00””, (3.3.2)

yz-yz’ zx’zx

Assuming that inﬁnitesimal strain theory holds, the form of the assumed displacement

0‘) and “(1‘)

ﬁelds in (3.3.1) is such that the inplane displacements u x y

contribute only con-

(’0

(k)
yz and o",z .For con-

stant (i.e., no variation in 2) terms to the transverse shear stresses, 0'
(k)

sistency, then, the transverse displacement uz

should also contribute only constant

terms to these stress components. In order to achieve this consistency, terms in the trans-

62

verse shear strain expressions that involve the transverse displacement variables are evalu-
ated at the mid-plane of the sublaminate, thereby taking a thickness-averaged value of
these terms. This approach is equivalent to ignoring the effect of transverse normal strain
on the transverse shear strains, a very good assumption for most laminates. Thus, the
transverse shear strains in the kth layer are:

k 113wb+BW)

72:) = + 2% + gab +331:
'=‘ (3.3.3)

3w
(’0
sz = y+k21ni+ %(+ayt)

i=1

From (3.3.3), the relationships between transverse shear strains in the (k+1)th layer and

the kth layer are:
k
(k+ 1) 1 3W1, 3W
73" "+ Zgi+§£8x +3xt)
i= 1
k—l
1 awb 3w (3.3.4)
=Wx+2§i+§(ax a I) a]:
1': 1
k
= 'Yix) 1' £1
By the same token,
k 1 k
7;: ) =7? + nk (3.3.5)

From the constitutive equation in (3.2.2), the transverse shear stress-strain relationships

are:

63

(k) (k) (k) (k) (k)
cyyz = C44 sz + C45 yxz

(k) _ (k) (k) (k) (k)
xz " C45 yz + C55 yxz

(3.3.6)
0

Making use of (3.3.4), (3.3.5) and (3.3.6), the transverse shear stress continuity conditions

at the kth interface can be recast as:

(k

k k k k k1 k k
1CL’1§,’+C£5’1§,’>M =16” ’11”+n,.)+Cf.51v’,.z’+§,.)}|_

 

1* (3.3.7)

k k k k kl k k k
1C25’1§,’+C§5’v§,’) = {Ci * ’1v‘zs’+n.)+C§’1 111"+§,.)}|z=z

 

Solving (3.3.7) for 11k and 5,, , and using (3.3.3), we can develop the following recursive

equations.

11 +21- «oz—:1 :—:1211;212:1213—1141}

 

 

k 1 (3.3.8)
.. ’13wb+aw13wb+aw,
5" = ck1‘"y+.zl"‘+§(ay ay )1 W + 2:: + 5(5 311)}
1: 1:
where
ii" = A1(C(Sl;+1)Cg4k)_C£§-l-1)C(§))_‘ 1
[1+1
131 = —(C(5k+1)C(k)—C,(,5+1)C(k))l
A[11-1-1
a, = (Cf,4"*”cf,’;’-Cf,’;*”cf,4"’) 13.3.9)
Ak:-1-l
ak_ A _(C(k+l)C(5k)_ £§+1)C(k))—l
k+l

(1+1) (k) (k+l)
A1+1 = C44 C55 ~C45

64

Solving the recursive relation, (3.3.8), i,- and 11,- are then found to be:

a — +.1_(a_wb+§f‘) +d +1(3Yb+.a;.vt)
i-c"wy23y 3y iw‘ZBx 3x

(3.3.10)
11'. = ai{\|ly+§(a—y +a_y.’)}+bi{\|lx+-2-(§-x— +§;)
where
1-1 1-1
a' = a, 1+ XajJ+EiZcJ
J'= i=1
i: 1: (3.3.11)

.53 .9"
II II
.5." P"
+ -l-
Ml. MT '
JP .5“

\——/ \_a/
+ +
.50 .5»
MI. MT
'3 e-
N x.

From Eqs. (3.3.9-3.3.1 1) it can be seen that £1 and 11,. depend on the ratios of shear prop-
erties between adjacent layers and the shear deformation in each sublaminate. If either of
these quantities is small, then discrete-layer effects will also be small.

The rotational variables W): and my in the displacement ﬁelds are now eliminated by
introducing translational variables u, and v,, the inplane. translations at the top surface of
the sublaminate [41,42] in order to expedite the development of versatile ﬁnite element
models. Thus, rather than describing the inplane displacement ﬁeld by a translation and
rotation at one point, it can more conveniently be described here by the translation at two

points.

65

Moreover, the displacement ﬁelds still need to be manipulated to develop a C0 ﬁnite

element, because the derivatives of transverse deﬂections,

aw, aw, aw, 3w,

appear in the displacement ﬁeld. The existence of the derivatives of transverse deﬂections
in the displacement ﬁelds results in their second derivatives appearing in the strain energy
functional, so Cl continuity of wb and w, is required. It is desirable to alleviate such a

requirement by introducing new rotational degrees of freedom as follows:

=—9 9

3 W3; = - ”,5; =9xb’ay n (3.3.13)

These relationships will be constrained in the strain energy via the penalty method. There-

fore, the displacement ﬁelds can be written as follows:

u(k) = (1)1111)“ “b + (Dyan + (1,1210% + (bygv, + (by? 9x!) + (bade + «1,119,, + e‘k’ey,
= 111+” +111... 211+“ 11.2 ”31111,. 1121».
k) (10 (3.3.14)
+ ‘1‘,“ 9y}; + ‘1’42ﬂ6

(k)

z = wbﬂl “”1292

Employing the indicial notation, the modiﬁed displacement ﬁelds can be expressed in

more concise forms:
11) - 11:) (k)_ k 1
u = u Be 3, u M113, u§’=anB, a = 1...4;13 = 1,2 (3.3.15)

where the index [3 is used to denote the top and bottom of the sublaminate with 1: bottom

66

and 2: top, and (PL? , W3; and RB are shape functions in the thickness direction, and will

be deﬁned in Appendix A. Unless noted otherwise, summation on repeated indices is

implied. The variables Bali are:

‘7113 = “B , (.423 = VB , 17313 = 6x5 , 1743 = Gyﬁ (3.116)
Since through-the-thickness shape functions, (1):}; , ‘l’g‘g and 913’ which are only func-

tions of z, are included in displacement ﬁeld, nodal variables in the ﬁnite element can be
approximated by two dimensional shape functions in x, y. Thus, the strains are obtained

from Eq. (3.3.15) as:

(k) _ 1k) 1k) - 1k) 1k)
exx = “aﬂx‘baﬂ eyy = “0113wa 822 = WBQBJ
1k) - 1k) 1k) _ - (k)
7,. = “aa‘l’ap,z+wa,y913 1. 'sz -“a3¢aﬁ.z+wﬁ.xaﬁ ,, (3.3.17)
Z = i Z = i

k - k - (k
yky) = uaB,y¢ExB)+uaB,qua[i a = 1-"43ﬁ = 1:2

where tensor notations are employed and a coma in the subscript implies partial differ-
entiation. In (3.3.17), the contribution of the transverse deﬂection to the transverse shear
strains 7;]? and 7;? are evaluated at the middle of the sublaminate to ensure constant
shear stresses through-the-thickness of a sublaminate (as discussed previously).

From (3.3.17), it can be seen that the transverse normal strain,

k W -W .
2:2) = was)“ = ' h b (3.3.18)

 

is constant in a sublaminate. Such a uniform normal strain may give rise to substantial

errors in composites which have a soft core or layer. Let’s consider a composite which is

67

composed of 3 layers, one of which is very compliant. When the composite is under load-
ing as shown in Figure 16 (a), the uniform transverse normal strain distribution through
the thickness is predicted from (3.3.18) as shown in Figure 16 (b).

11111111111111

actual

/ r uniform

N
»

Compliant

 

 

 

(a)

15;

Figure 16.A composite under transverse tension load (a) and the corresponding
transverse strain distribution through the thickness (b)
However, the actual distribution is not uniform, so that we can see that there is a big devi-
ation between the predicted strain and the actual one. Thus, it is desirable to improve the
distribution of the transverse normal strain through the thickness. If the transverse strain
ﬁeld is improved, the Poisson’s ratio locking can be prevented as well (to be discussed
later). The improvement can be achieved by instead assuming a constant transverse nor-

mal stress, (T-

u , through-the-thickness of a sublaminate. Tessler and Seather attempted to

improve the transverse strain, but they used a ﬁeld consistent polynomial [80]. (Tu can be
determined using Reissner’s mixed variational principle [81]. From (3.2.2), the transverse

normal strain is obtained as:

68

(k) (k) (k)
-(k) _C_1_3 (k)_ C__2_3 (k) C36 (k)+1 —

ﬁzz = -C—(—];) xx _Cg-T) yy _C—g§)yxy+ CT]? 22 (3.3.19)
o—zz is then determined by analytically solving the following relation:
k
_k k
0 = z I (822 - Ezz)dz
k = lzk_ l
(3.3.20)

 

(k) (k) (k)
2 }{W:; Wb_ C13 80:) C23 80:) C__367(k)+ 1 _;]}dz
cm H cm W cm U can zz

k=lzkl

Substituting the strain - displacement relations in (3.3.17) into (3.3.20), the constant trans—

verse normal stress o—zz is found to be:

_ Bub Bu, av, av, aeyb

on = 37 (P1+56)+Tx(P2‘56)+$(P3+S7)+a—x(‘P3+SB)+a—x (P4+Sg)
86 39,, 86 ea, Bu

+a—xy '37(P4+S9)+ Xb-é—(P5+Sw)+ x’"(P5+S]O)+a—y b1(Q,+S “a? (-Q1+Sz)

(3.321)

av, av aeyb 86y, 66x1,

+—b(Q2+S3)+a—y (Q3— 53)”??? yb(Q4+S4)+—y 'M'Sﬂ'i' Jrb(Q5'*'Ss)
89x!

+a—y- (Q5+SS)+wb(—T])+wt(Tl)

where the coefﬁcients are deﬁned in the Appendix A. The newly deﬁned transverse nor-

mal strain can now be written as:

k
.g = uaB xxggﬂ+uaﬂ yx§ﬁﬂ+wﬂxﬁ a = 1...4; B = 1,2 (3.3.22)

where the functions I.“ B) and x3 are deﬁned in the Appendix A.

The important properties of the present laminate theory are:

69

1) The laminate thickness may be subdivided into several sublaminates for reﬁned analy-
sis, if needed.

2) Discrete-layer effects are captured within each sublaminate by allowing a piecewise
(layerwise) continuous through-the-thickness variation of the in-plane displacements.

3) Continuity of transverse shear stresses through the thickness of each sublaminate is sat-
isﬁed.

4) Transverse normal strain distribution through the thickness of sublaminate was
improved by assuming constant transverse normal stress.

5) Accuracy and efﬁciency of the theory are adaptable, depending on the number of sub-
laminates used.

6) In the displacement ﬁelds, only ﬁve engineering d.o.f. (three translations + two rota-

tions) are present.

3.4 Governing Equations and Boundary Conditions

The equations of motion, the essential and natural boundary conditions, as well as the
displacement-based ﬁnite element model can be developed from. Hamilton’s principle.
Again, attention will be limited to a single sublaminate. The constraints between the deriv-
atives of transverse deﬂection and the rotational d.o.f in (3.3.13) should be imposed in the
energy functional. There are two ways to impose the constraints - either Lagrange multi-
plier method or penalty method. Since the Lagrange multiplier method introduces addi-
tional d.o.f., the penalty method has been used to maintain computational efﬁciency and
convenience. Therefore, the energy functional is modiﬁed here to include the imposition

of the constraints via the penalty method. Hamilton’s principle for the mth sublaminate

70

can be deﬁned as follows:
5 '{U—K—W+P}dt=0 (3.4.1)
‘0

where U is the internal strain energy:

xx 8xx+oyy€yy+ozzezz +0”sz +6 u'sz +0

N
"ll/ck kk(k)k kk(k)k kk

=2J‘_2_[O. O.£_:()() ()() () ()() () ()Y(y)]de (3.42)
=1v

Vk is the volume of the k-th layer. K is the kinetic energy:

NM
2
= 2 I; Mkll -2+(k) “(’02 +1221!) lde (3.43)
=1V,‘

pm is the density of the k-th layer and a superposed dot indicates differentiation with
respect to time. W is the work of external forces:
w = IwaPadQ+ Itiuidl"

(3.4.4)

I =oijnj;oc=1, 2; i=1, 2

P is the penalty constraint:

P = $11211(u3B-gﬂf+(u4ﬂ+::p)2}]dn (3.4.5)

where y is the penalty parameter. As 7 is assigned successively larger values, the con-

straints of Eq. (3.3.13) are satisﬁed more exactly in the least squares sense. In the internal

71

strain energy expression, (3.2.4), we need to pay attention to the transverse normal stress
terms. The constant 0:; leads to an asymmetric stiffness matrix. Thus, Cu from the con-
stitutive equation, (3.2.2), instead of the constant transverse normal stress should be
employed.

Substituting Eqs. (3.2.2, 3.3.15, 3.3.17, 3.3.22, 3.4.2, 3.4.3, 3.4.4, and 3.4.5) into (3.4.1),

taking the variation of (3.4.1), and integrating by parts, (3.4.1) can be rewritten as:

t
' _ _aNmB _aNyaB _a_NwB_ _a_QwB aRmB_ aRymB
J J[{5u“ 8x 3y 8x 8y +anﬂ+QmB 3y 3x

'1

2
8W5 32 ur a a? V
(—)( — .M Jim}

2 2
BR BR 817 a w 317 8 w 82 w
+ 5WB{RZB ay ax PB Y[ay ayz J + 'Y[ + x J+— a—'t2 215p}]dﬂ

(3.4.6)
+ §[5§“B{Nman + NymBny + Nzaan + QzaBny + nyaBny + Ryxaan}

_ awB _ BWB
+8wB Rany+Rxan+yny(u3B—-é—y- j—ynx(u4B+§; ) ds

1!:
N" (k) - (k) (k) -
+12 [(511515 «95015435215 5 aa‘l’aa’)dzds]dr = 0

(a,r=1...4;[3,p=1,2)

Where NxaB’ NyuB’ N zaB’ NzaB’ QzaB’ anB’ Qxaﬂ’ nyaﬁ’ RyxaB’ RzB’ RyB’ RxB are the

sublaminate stress resultants shown in the Appendix A, nx, n are the directional cosines,

)’

Ill

aBrp’ [LBW 13p are the plate inertias deﬁned in the Appendix A, and SaB is the Kro-

necker delta. Considering that variations of the generalized displacement degrees of free-

dom are independent for to < t < t1 and zero for t = to, II , the equations of motion are

72

derived as:

.. . aNxthaNyaBaNzaBa_Qz(1B anyaB+ aRyxaB
5““13' 3_x +8y +3x +a—y - anﬁ - Qmﬂ + 3— +3x

2_
awB BWB Bu am 32 umv
_ _ ._ I
+ 7(u uaB ay B)831 + 7(1“ GB + ax “)5“ 3‘2 IaBrp+ at2 aBrp (3-4-7)
312 an a- 32 a a a2
u w u w w
5 _yB+_xB_R +1; + _3B__B _4B+ _B =_p1 W
WB . 3y 3 x zB B Y[By ay2 7 a x +3x2 a: IPB
The corresponding boundary conditions are:
Essential B.C Natural B.C.
aaB (NxZan + NyaBny + Nullan 4' Qany + nyaBny + Ryxaan)
(k)n (Ion (k)
+ X 2! (Gun jcp ffB + 0‘2} nJ‘PaB)dz (3.4.8)
=12»:
wB Rany + Rxan + yny(u3B -a—y )- -'yn (“413 +— a? )

3.5 Finite Element Formulation

Since the primary variables are functions of x, y and I only in (3.4.6), we need to con-
sider only two dimensional shape functions for the ﬁnite element model. Thus we can con-
sider that the most obvious element to be developed would be a simple four-noded
Ahmad-type degenerated plate element with 10 degrees of freedom per node (5 bottom
surface d.o.f + 5 top surface d.o.f). However, such an element would not fall within the
constraints of having three translational and three rotational degrees of freedom. (Most

available commercial codes do not allow more than six d.o.f per node). Even more impor-

73

tantly, a major advantage of the present formulation would be lost if a traditional four-
noded element were developed. Because the degrees of freedom represent quantities at the
surfaces of a sublaminate (this was not by chance), it is most reasonable to develop the
ﬁnite element model in the form of an eight-noded element with ﬁve degrees of freedom
per node (see Figure 17). Such elements are classiﬁed here as “regenerated,” compared to
the popular “degenerated” shell elements obtained by imposing constraints on solid ele-

ments.

u,, v,, w,, 9,5,, GB, 14,, v,, w,, 9,“, 0y,

 

 

   
 

 

 

r. c-‘u 4.4,: ',' ; ;.‘ .«'.'“.
_._.... .-..
"43': 4121’.- .-'Tr“.3-.r4m.
-‘-.~.-.-'.~. .453

_ 0'
.. , . . .
W ‘r—x‘p ~... 1-
J— " "' - 1
q
.,

rl-l'

. . . '. .i ;‘_

. ‘. 3'
.|_‘.._ Bf .'-"‘~"' Bz.__ “ , B, ..I:.- “-7. . ,
.B...',Q.-Vf*‘.:".‘ﬁ;.:k,h>. . 17-34453‘ . B . . 1 . _. -- ~., .. ; ..:4‘.‘.'-‘ 1

  

.'-' .. *".4; ':..‘ a; 22‘." m'“,
41. :32..:a-‘ru.3 eats-32‘ {3.4—.41 ~

I If"? t r': \ ._).’,.‘..m’—1r.'¥‘.‘tf'.‘t‘:‘1§;__§f
‘4: - "‘2 , 4‘ B”. -. 5.7.34

   

“1» V1» W1» 9:2» eyb “1» Vi» W1» 9x1» eyb

Figure 17. Element topology and nodal degrees of freedom.

This approach has been used with great success [41,42], and has the very important advan-
tage of allowing discretization in both inplane and through-thickness directions without
the need for any special multi-point constraints. Elements are assembled in the standard
way. Therefore, the element can be used to predict the global response of a laminate using
only one (or a small number) of elements through the thickness, and local effects such as

interlaminar stresses can be captured by reﬁning the mesh in the thickness direction near

74

the interface(s) of interest. The capability of discretization in the through-thickness direc-
tion enables one to simulate the separation of the laminated composites due to delamina-
tion into two parts in the thickness direction.

Inplane displacement and rotation degrees of freedom are approximated by the bilin-

ear Lagrange interpolation functions:

4 4
‘716 = 2145,10,. 7‘26 = Zvﬂipi
i=1 i=1
4 4 (3.5.1)

1733 = 2 exﬁipi 1743 = 2 eyBiPi

i=1 i=1

where

P.- = ,l,(1+§,-§)(1+nm) (3.5.2)

and g, n are the element natural (local) coordinates. For transverse deﬂection degrees of
freedom, an interdependent interpolation concept similar to that developed by Tessler and

Hughes [75] is utilized.

Interdependent interpolation

One of objectives in the current research is to develop a C0 element. However, the ﬁrst
derivatives of the transverse deﬂection are present in the displacement ﬁelds due to evalu-
ating the zig-zag functions. In consequence a Cl element needs to be developed unless we
relax such a requirement. The requirement of CI continuity was eliminated by introducing

new rotational d.o.f. as discussed previously. Then the constraints between the derivatives

75

of transverse deﬂection and the rotation degrees of freedom are imposed in the strain
energy functional using the penalty method. Recall that the rotational degrees of freedom
were approximated using bilinear shape functions. If we employ bilinear shape functions
to approximate the transverse deﬂection d.o.f. the constraints in the strain energy func-
tional must enforce the relationships between the bilinear variations of rotation degrees of
freedom and linear variations of transverse degrees of freedom. As a result, the element
locks. Thus, the order of shape functions to be used for transverse deﬂection d.o.f. is
biquadratic. The simplest shape functions to be used is the quadratic Serendipity shape
functions which have mid-side nodes. However, such a ﬁnite element is not of convenient
form. Thus, the interdependent interpolation scheme introduced by Tessler and Hughes
[75] is utilized in the current research. This approach improves the element behavior very
effectively without introducing additional degrees of freedom.

The virgin element is assumed to have mid-side nodes with only deﬂection degrees of
freedom at the top and bottom surfaces (see Figure 18), and then the mid-side nodes are
removed by writing the deﬂection degrees of freedom in terms of rotational degrees of
freedom according to interdependent interpolation scheme. The interdependent interpola-
tion scheme enforces the linear part of the constraints of the form:

i‘rb, t aWb t

ax +9yb,,=0, _’ -9 , =0 (3.5.3)

where the coma in the subscript does not imply any differentiation. Generally, in quadri-

lateral element coordinates, the constraint (3.5.3) is recast into:

(.31: _ 9,.) = 0 (3.5.4)

76

Lets designate the coordinate running along an element edge and denote n the coordinate

running normal to an element edge as shown in Figure 18.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

xij = xi—xja yij=yi-yj‘

Figure 18.Quadrilateral element coordinates and deﬁnition of rotation

The xy coordinate and the ﬁn coordinate systems are referred to as the “global” and
“local” coordinate systems, respectively. Initially, a Serendipity (eight-node) transverse

deﬂection distribution is assumed:

8
wb,t = 2 Ni(&n n)wbi,ti

i=1

4 4 (3.5.5)
x = 2 Pi(§,n)x,- y = 2,10,03,71)»

i=1 i=1

77

where Pi(§, n) , and N ,(i, n) are bi-linear and Serendipity shape functions respectively,
and wbi, n. denote transverse nodal degrees-of-freedom at the bottom and top surfaces. The

biquadratic serendipity functions are [82]:

N, = gamma+n,-n>-(1-§2)(1+n,—n)-<1-n2)<1+§.§)1§%n3
(3.5.6)

+%(1-§2)(1+n,~n)(1-§.-2)n.-2+%(1"12)(1+§i§)(1“li2)§i2

where —1 .<_ {3, n S l , and 111,};- are the location of node in the local coordinate. Rotation
degrees of freedom are approximated by the bilinear Lagrange interpolation functions as
in (3.5.1).

For side 1-2, the constraint along this side (1] = -l) is:

aw“
(a: -e..)

where a = l, 2 imply bottom( 1) and top(2) of the sublaminate. Since the different shape

= 0 a = 1,2 (3.5.7)
n=-l

 

functions are used, the number of d.o.f at the comer nodes is not same as those at mid-side
nodes. In order to have a uniform four-node topology at the top and bottom surfaces, the
four mid-side node transverse deﬂection degrees-of-freedom at edges are condensed out

by using four differential constraints at the edges:

= o (o s z s 1) (3.5.8)

aw
i(_ a - Gnu)
§n=il

as as

 

where hk is the edge length. For the edge 1-2, the differential constraint is given by:

78

a
a W“- em) = 0 (3.5.9)

333

 

Due to the property of the shape functions, the evaluation of the derivatives of transverse

deﬂection at the edge leads to:

I

 

 

 

8
Wa|n=—1 = ZNIWO‘I (3.5.10)
i=1 11:4
' = = = N = N = O,
smce N3|n=—1 N4In=—1 N6In=-1 7ln=-1 8|"="
1 l
wal _ = -§(§-1)w,+-§(§+1)w2+(1-§2>w5
11-- 2 2
.a_wa — (_l+§\w +(l+§)w -2§W l
as _-1 2 J 1 2 2 5 1: (3.5.11)
82
Wu 1
_ = W +W -2W —
n=—] 71:"

 

 

where ds = J [I 1d§’ and J 1 is the line Jacobian. The normal rotation at the edge is

n = -
approximated using bilinear shape functions. Then the derivatives of the normal edge rota-

tion along side 1-2 due to the property of the shape functions is given by:

= Z 12.9,")n = -1 (3.5.12)

79

4
9n n =_ = 2 Premln = -1

1:1

1

= i[(1—§)e,,1+(1+§)e,,2] (3.5.13)
39,, 1 1
g - - 5(6n2—enl)7l
‘n—-l n=—l

Substituting (3.5.11) and (3.5.13) into (3.5.9), solving for the mid-side W5 degree-of-free-

dom, then w5 is found to be:

11 1

Since the geometry of the element can be approximated by bilinear shape functions, the

line jacobian, J 1 , is obtained as follows:

 

 

 

 

 

 

3x
Jc|n=_1 = .zlPixi , Tan=4 — 2(x2—x1)
l— "=4
4
yl = 21’) ' a—y -(y y)
TI=-l . 1 It 9 .. 3§n__1 2 2 1
I: 71:4 -
J" = [ax ]2+[3_y )2 (3.5.15)
"=‘1 §Eﬂ=-l agn-_1

 

 

 

II
%
H
N
NI
.3‘
\_/
N
+
/—\
‘<
N
NI
.3
\_/
N

The other mid-side W6, W7, W3 degrees-of-freedom can easily be obtained by the same

procedure, and are found to be:

8O

1 l h2
“’6 = §(W2+W3)-Z(9,,3-9,,2)—2-
1 1 h
W7 = §(w3+w4)—Z(6n4-9n3)§3 (3.5.16)

_ 1 l ’14
W8 - 5(W4+W1)‘Z(9n1‘9n4)'2'

Now we need to transform the rotations at the edges into ones in global coordinates. Tak-

ing the edge 1 - 2, we can derive the transformation law.

 

 

 

Figure l9.Coordinate transformation between edge and global coordinates

From Figure 19, the transformation law is derived as:

[n = cos(n, x) cos(n, x) x (35.17)
s cos(s, x) cos(s, x)

where (n, x) is the angle between the two axes. Then on the 1 - 2 edge,

81

 

9:11 = cos(n, x) cos(n, x) 9x1 (3.5.18)
631 cos(s,x) cos(s,x) El),l
where
cos(n, x) = cosG = cos(-B+90) = sinB = )’_2’_:_ﬂ
’ (3.5.19)
. x1 "x2
cos(n, y) = cos(9+90) = —sm9 = -cosB = h
1
Thus, (3.5.18) becomes:
yZ-yl Jr1‘3‘2 yZ-yl x1'x2
6n] = -—h—1—9x1+T9y1, an = TQM—ray, (3.5.20)

By the same procedures, the coordinate transformations of the other sides are obtained:

on the 2 - 3 edge,

y3-y2 xz‘xs Y3‘Y2 x2-x3

 

 

 

 

9,12 = —h—2—9x2 + TOW 9n3 = ’12 6x3 T9y3 (3.5.21)
on the 3 - 4 edge,
y4‘y3 x3 ‘x4 y4—y3 x3 "x4
on the 4 - 1 edge,
_ yl-Y4 x4-x] YITY4 x4-x1

82

Substitution of Eqs. (35.20-35.23) into (3.5.14) and (3.5.16) yields the following expres-

sions for the mid-side node transverse deﬂection degrees-of-freedom in global coordi-

nates:

w5 = gov.+w2)-§t(y2—y,)(e.2—e.l)+(x,—x2)(e,2—9,,>1

W5 = %(W2 + W3) " $03 " Y2)(9x3 " 9x2) + (x2 ‘ x3)(9y3 — 9,2)1
(3.5.24)

W7 = %(W3 + W4) ‘ %[()’4 - Y3)(ex4 — 9x3) + (x3 — x4)(9y4 " 9y3)]

W3 = %(W4 + W1) - $101 'Y4)(9x1 — 9x4)+(x4—x1)(9y1—6y4)]

Finally, interdependent interpolation functions are derived by substituting the mid-side

node transverse deﬂection degrees-of-freedom expressions, (3.5.24), back into (3.5.4):

4
w = 20>,ng .9 .+Ny,e

I! X!

y.) (3.5.25)

i=1

where N xi and N yi depend on N i + 4 and the projections of each edge onto x and y axes

and are given by:

1 1
le = —§(>’12N5-)’41N8)’ Nyl = §(x12N5—x41N8)
1 l
NxZ = ‘§(>’23N6‘>’12N5)4 Ny2 = §(x23N6-x12N5)
1 1 (3.5.26)
Nx3 = _§(y34N7_y23N6)’ Ny3 = §(x34N7-x23N6)
1 l
Nx4 = —§(>’41N8‘>’34N7)’ Ny4 = §(x41N3‘x34N7)

where

83

N5 = %(1-§2)(1-n)5 N6 = %(1-n2)(1+§)
N7 = $0 —€2)(1 +11), N8 = %(1-Tl2)(1-§) (3.5.27)
xij = xi—xj’ yij = yi—yj

Furthermore, the mid-side transverse deﬂection d.o.f. in terms of covariant rotations is

readily obtained from Eqs. (3.5.14, 3.5.16) as:

1 1
ws 2 5(w1+w2)—Z(9n1-9n2)
l 1
(3.5.28)

1 1
W7 = 5(W3+W4)’Z(9n4’9n3)

1 l

The interdependent interpolation functions in terms of covariant rotations can be written
as in the following expression by substituting (3.5.28) into (3.5.5) (these interdependent
interpolation functions will used to derive ﬁeld consistent shear strains):

4
w = Z (Pl-wt. + Ngieg. + Nniem) (3.5.29)

i=1

where N g1 and ”hi are given by:

1 _ 1
”:2 ZN6’ an - 3N5
1 1 (3.5.30)
N 1N N 1N
§4 - -2 8’ 114 _ -— 7

Since the interdependent interpolation scheme is employed, the transverse deﬂection
d.o.f. is coupled with the rotational d.o.f. Therefore, the transverse distributed load acting

on the top and bottom surfaces should be applied to the appropriate nodal d.o.f.

Field Consistent Shear Strain

The locking phenomenon has been a big issue in the development and application of
the ﬁnite element method. In an attempt to alleviate shear locking, the interdependent
interpolation scheme was adopted. However, it turned out that it did not remove the shear
locking completely, and the developed element still locks in the very thin regime.

To avoid shear locking, a lot of works have been accomplished. One popular method
employed is reduced/selective integration [83-88]. However, reduced/selective integration
has had limited success. In some elements the locking phenomenon is not completely
avoided. Moreover, when the order of integration is reduced enough to alleviate locking,
spurious kinematic or zero strain energy modes arise in the element stiffness matrix. These
kinematic modes can be suppressed by adding a stabilization matrix to the stiffness matrix
[89,90]. Thus, some researchers try to alleviate the locking problem by special treatment
of the transverse shear strain components [91-94], while using full integration to preserve

the correct rank of the element stiffness matrix.

85

Prathap and Somashekar found that the source of the shear locking phenomenon
comes from inconsistency of the transverse shear strain ﬁelds with respect to the inplane
coordinates, x and y [72,73]. There are also inconsistencies in the thickness coordinate, 2,
but terms associated with the thickness coordinate are small when the aspect ratio is large,
so that the transverse interpolation inconsistencies are not so crucial in locking problems
compared to the in-plane interpolation inconsistencies [95]. Thus, in order to prevent the
shear locking phenomenon, we need to make only the transverse shear strain ﬁelds be ﬁeld
consistent with respect to only the in-plane coordinates.

Since a general quadrilateral element is non-uniformly mapped from a Cartesian (glo-
bal) system into a natural (local) coordinate system, it is difﬁcult to see the consistent form
of the transverse shear strains clearly. We need to transform them to the natural coordinate
system. The transverse shear strains in the natural coordinate system are called covariant
shear strains. Any ﬁelds can be transformed from the covariant to Cartesian coordinate

systems by the following expression:

1952 13x2
{ } = [J]{ } (3.5.31)
1371?. ﬁyz

where [J] is the J acobian matrix and is deﬁned as:

3x a

[J] = [1” J ‘2] = (3.5.32)
121 122

‘1

 

 

a— a
iay
_an a

The inverse Jacobian is obtained as:

 

86

-1 7 i 1 j -J°
[J] = 111112 =m e2 .12 =
121 122 “121 111

where I] I is the determinant of the Jacobian matrix.

(3.5.33)

1°) wla'

lg” S’lg’ .

 

 

'cv
cu
E

The transverse shear strain in the Cartesian coordinates, 7” is mapped into the natural

coordinates as:

aux
sz(§21’l) = a—z' +3}!
(3.5.34)
aux Buzag Eula"

= a7 *(‘aE5§*ﬁa—x)

Considering that the inverse J acobian terms are included in Eq. (3.5.34), some new covari-
ant in-plane displacement ﬁelds, u: and “n are introduced in order to derive the covariant

shear strains, so that the in-plane displacement ﬁeld, u x can be rewritten as:

BE 31’] (3.5.35)

ux = “Ea—x-‘i' “"3;

As the modiﬁed displacement ﬁelds can be expressed in more concise forms in Eq.

(3.3.15), covariant in-plane displacement ﬁelds, u: and “n can be expressed by:
_ (k) _ (k) - .
“g = u§a8¢§aﬁ, U1.I = unaBCDnaB a. = l, 2, B = 1,2 (3.5.36)

where the new set of degrees-of-freedom will be represented as:

87

a :17 ,a = u ,9
fa". _ng _m :5 g” B = l(bottom),2(top) (3.5.37)
"naﬂ' “nlﬁ’unzﬂ = “M3: 9116 .

Thus, u: and “n can be deﬁned in terms of their new set of degrees-of-freedom, (igaﬂ ,

“naB) 1n the covariant frame as:

ux = amBqaggg—E” Mofﬂﬂ-gg a : 1, 2; [3: 1, 2 (3.5.38)

where (1)223 and (Pm nab are still shape functions of the thickness coordinate (2) only and

are given in detail in the Appendix A. Substituting Eq. (3.5.35) into Eq. (3.5.34), the Car-

tesian based shear strain, 'sz’ can be obtained in terms of the covariant shear strains,

7&2’7112 as
+23“ a: 8“z 311
"ME n)= (.3.? +a§ —)3—x+ (g: +331): (3.5.39)
=7§3111+Yn2112

where

_ _8_Zu§+ a_uz _ dun duz

Now, the covariant shear strain is derived in terms of covariant degrees of freedom by sub-

stituting Eq. (3.3.15) and Eq. (3.5.38) into Eq. (3.5.40):

7:2 = “élpq’gla, z 4’ 9gzpq’g23, z + quﬂg = aéab‘pgap, z + WB’énﬂ a = 1, 2 (3 5 41)
7n: = “nlﬁ‘pnlb.z+9n213‘pnzb.z+WB’nQB = ar1<=:13“’1r1m13.z+“’B’nna B = 1, 2

88

The shear strain ﬁeld, 'sz , can be expressed in terms of covariant degrees of freedom

by the same procedure as 7“:

Bu au
sz(§, Tl) = a—y +a—'z
Z y
(3.5.42)
= Buy + cal-33+ 8153-“)
BE; ay+ 811 By
u y can be written as:
_ . 3i . 3n
uy — 14:53:4' (Ln-a; (3.5.43)
Covariant in-plane displacement ﬁelds, it: and an can be expressed as:
uy = agaﬁwgaﬁ 1’2, «Landﬁwnaﬂ 132 a = 1, 2; [3: 1, 2 (3.5.44)

where through-the-thickness shape functions, was and \Pnaﬁ are given in the Appendix
A. Then the shear strain, yyz, can be expressed in terms of covariant shear strains,

“7&2 and 7m. Mathematically, this is as shown below:
(g )_ (Bug +114“) +(a:n+ iu)-
7” .n _ 3% ’2‘ an ’22 (3.5.45)
= 7:2121 +7nz122

where

_ (k)
‘ “zaa‘l’wz

~2
W
N

I

             

II
NI}

(3.5.46)

.. (k)
?nz = unaﬁwnaﬁ,z+wﬂsn “- = 1, 2, B = 1,2

     

89

ﬁéaB and anus of Eq. (3.5.37) can be approximated by bilinear Lagrangian interpolation

functions, Pi(§, n) whiCh are used to interpolate 17,43 degrees-of—freedom in Eq. (3.5.1):

“eons = P “éaﬁi’ W " Pia “nal” (3.5.47)
i=1, ..,4;a=1,2,;[3=1,2

For “’6 , the interdependent interpolation functions in the covariant frame as in Eq.
(3.5.29) should be used. Using Eqs. (3.5.25), (3.5.41), (3.5.47), the covariant shear strain

component, 7&2 in 7“ of Eq. (3.5.39) can be recast as:

we s-z)

           

_ (Pi. éwﬁi + ”4449415” Nni. 4911139

I
NI}

_C__-8712)
—{ enlB+9n23‘en3B+ 9714

     

NII'}

+%[(1— n){¢§1a 204.5213 ' “4213) + (9:23 ‘ eélbwém 2 “61115 " e

       

+(1+"){‘1’§m.z(“§3a‘ “2:46) + (936" Gunman: + (91133-9

       

 

                                             

=(WZB‘ WIB)}

531':-

 

{q’m z(“§BB + “4415) + @4211: _ (W313 — W419}

                   

NI}

Since the coefﬁcients associated with the quadratic term (1 - 112) in Eq. (3.5.48) include
only 9111' rotational degrees-of-freedom, they apparently yields an inconsistent shear strain
ﬁeld, 75,. Those terms do not vanish in the thin regime where “7&2 tends to zero, leading to

spurious constraints and shear locking. It was identiﬁed that the inconsistent terms came

[_‘1

90

from the interdependent interpolation scheme for degrees-of-freedom WB' Eliminating the

inconsistent terms, a consistent 752 can then be obtained in a concise form:

71:2: Piaigaa‘l’éanzHP 4W5.- +N§i§9wima

 

2:

(3.5.49)

NI}

i=1,...,4; 0L=1,2; B=1,2

Thus, we can obtain the consistent transverse shear strain, Yﬁz of Eq. (3.5.49) by approxi-

mating WB by the modiﬁed interpolation function as follows:

WB = Pini+N§ie§Bi i = 1, ..., 4, B = 1,2 (3.5.50)

The other consistent covariant strain, 7m can be derived from the inconsistent shear strain,

711: in Eq. (3.5.39), in a similar way as

.(bu‘) + (P

=P'u naﬂz

mal: (3.5.51)

l-{l
.3

naBi i’Tl WiB+ NTIi’TI 6115i

hill}-

with interpolation function for the WB degrees-of-freedom now given by:
EB = Pini+Nn OMB! i = 1, ...,4,B = 1,2 (3.5.52)

The other transverse strain sz of Eq. (3.5.46) can be written in terms of the other two
covariant consistent transverse shear strains, 3:2 and $112 using the inconsistent covariant
shear strains, 3%, and '7,” of Eq. (3.5.45). The other two covariant consistent transverse

shear strains are:

91

1’42 = Pi ”"igaa‘l’mazﬁpvawm +N§v§9§mma|z= ,,

N"

(3.5.53)

Ynz = Piunaﬂi‘yﬂaﬁ 2+(Pi’ﬂw Bi +N1‘Ii’ﬂ6 713i

       

NII'}

i=l,...,4;0t=1,2;|3=1,2

The derived consistent transverse strain ﬁeld is based on the concept that all transverse
shear strain ﬁelds should vanish in the thin regime of plate due to Kirchhoff constraints,
and the inconsistent ﬁelds which don’t vanish should be eliminated. This method can be
viewed as one of the assumed shear strain ﬁelds approaches [91-94], and it is analogous to
the selective/reduced integration concept except that it only removes the inconsistent
terms. However, selective/reduced integration removes even necessary terms in the strain

ﬁelds and causes spurious zero energy modes.

Consistent Penalty Constraints

To develop a C0 laminate theory and the ﬁnite element associated with the theory,
four constraints between rotation d.o.f and the ﬁrst derivatives of transverse deﬂection
were imposed in the energy functional by employing the penalty method. The constraints

are given as follows:

Cw aw, e Cyﬂ aw“ e -12 3554
=W =(a7+5) P” (")

These constraints are very similar to the form of the shear strains in ﬁrst order shear defor-
mation plate theory. This can be simply veriﬁed. If we consider the following rotation con-

vention and deformed shape,

92
/ 6y

/deformed shape

U/ —-

\original shape

 

 

 

 

Figure 20.Schematic and deﬁnition of rotation and deformed shape

then displacement ﬁeld becomes:

u = uo—za—x = u0+z6)

aw _ (3.5.55)

Thus, the transverse shear strains are given by:

_aw 8u_awo
”farm-5; “’1

_aw av _ aW0
'sz =8y +3z 3y 6"

(3.5.56)

We can see clearly that the constraints of Eq. (3.5.54) are exactly the same as the trans-
verse shear strains in ﬁrst order shear deformation plate theory. Therefore, just like the

above transverse shear strain ﬁelds, we need to check if any of these constraints are ﬁeld

93

inconsistent, and hence causes any spurious constraints. Using the relation from Eq.

(3.5.31), the constraints in Eq. (3.5.54) can be written in terms of their covariant form as:

aw
3y ‘5
aw aw
_ 5 B 5 5 ’ 7
_ “Va—g +1125 ~1129§p+1119nﬁ (3-5-57)

7 3W3 7 aWB
111(5E + 91.3) + 112(ﬁ — 953)

i11C§p+j12CnB

After substitution of shape functions into the covariant penalty constraints, C§B can be

l'CCﬂSt as:

aWB
C413 = 53 +9113)

_(l-nz)
- 3 (9452‘9§BI+9:B4‘9§133)

(1- 1+
4n)(W52 _ WI“ + 6an + 91132) + LIT-12(WB3 - “’84 + 91:53 + Ogﬂ4)

(3.5.58)

 

 

+

If we examine Eq. (3.5.58), the coefﬁcients of the quadratic term ( l — n2) are found to be
inconsistent because they are composed of only the rotational degrees-of-freedom 955,-.
Thus, they must be removed to make all the ﬁelds consistently vanish without inducing
spurious or locking constraints. The consistent covariant penalty component, C55 can be
obtained by interpolating the WB d.o.f. in a similar form as in Eq. (3.5.50) above instead of
Eq. (3.5.25) while keeping the interpolation scheme for the rotational degrees-of-freedom

unchanged. 6&5 can then be expressed in a succinct form,

94

5:6 = Pvgwm + (Nun: + Poegpi (3.5.59)
Similarly, the other covariant penalty constraint Cm; can be derived as:

with the W5 degrees-of-freedom interpolated as in Eq. (3.5.52).

Edge Consistency

Now we have developed ﬁeld consistent shear strains by examining transverse shear
strains in the covariant frame and eliminating inconsistent terms. However, even though
transverse shear strain ﬁelds are ﬁeld-consistent, the element still locks if edge-consis-
tency of transverse shear strain on any common inter-element edge is not met, especially
for general quadrilateral elements. It is crucial to ensure matching of the tangential shear
strain on any common inter-element edge. Mismatch of tangential shear strain will give
rise to spurious constraints on the edges when the shape of elements is a general quadrilat-
eral. This stems from the fact that when the ﬁeld consistent covariant transverse shear
strain is mapped back into global cartesian coordinates from the covariant frame, jacobian
transformations of two adjoining elements at their own integration points will redeﬁne this
consistency, and cause a mismatch of tangential shear strain at the common edge. Thus,
the choice of sampling points of the jacobian transformations plays a very important role
in edge consistency of tangential shear strain at the common edge. Prathap and
Somashekar [72] showed that the consistency of tangential shear strain along a common
edge between two elements is maintained by transforming ﬁeld consistent shear strains in

covariant frame into global cartesian ﬁelds at each node of the element instead of at Gauss

95

integration points. In this way, edge consistency in the ﬁeld consistent covariant penalty
constraint ﬁelds can be treated as well.

Replacement of Eqs. (3.5.49) and (3.5.51) with Eq. (3.5.39) gives the ﬁeld consistent

shear strain in the covariant frame as:

7xz(§.11) = 2:21:11 “an1712
_ k 7' 7
= (Piuéaﬂi¢(§ozﬂv z + (Pirawn: + N §,.§0§B.-)Qp)1 11

(3.5.61)

_ k " ’

nab: 1105.2 t
where 1‘: 1,...,4; a=1,2; B: LEE-25:95

2 =1'
2
To maintain the edge consistency on the common inter-element edge, the covariant nodal

d.o.f, 17:3,, Ting,- of Eq. (3.5.37), should be mapped into their Cartesian coordinate deﬁni-

tions, ﬂag,- of Eq. (3.3.16) at each node [42,76], i.e.:

§=-1,n=—l; ﬁ=l,n=—1; §=l,n=1; §=-1,n=1 (3.5.62)

Then, using Eqs. (3.5.32), (3.5.61), after the nodal transformation is accomplished, the

consistent Cartesian based shear strain, 7“ is given by:

4
- . (k) 7 . (k) 7
73:45:") = Z [uﬂi{(111)ipi¢§lﬂ, 2111 + (1203091113, 2112}

i=1

‘ k 7 . 7 7 7 _.
+ vBi{(112)iPi(D(§l)B, 2111 + (122)ipi¢f1k1)ﬁ, 2112}+ WBi(Pi’§./11 + Pi’n112)nﬂ
(3.5.63)

. (k) 7 " 7 . k 7 " 7
"' exai{(122).-(Pi¢§2p, 2111 + NgvnQBJn) - (112),(P,-¢$,2)3, Z112 + Nnv§93111)}

. (k) 7 7' 7 . k 7 '— 7
+ eyﬂi{(111)i(Pi¢nzﬂ, 2112 + angﬂﬁln) _(121)i(Pi¢(§2)B.2111 + Ngi’n08112)}:l

96

iii are used to map Cartesian d.o.f into the covariant frame d.o.f. Therefore, it must be
evaluated at the integration points. j i]. are employed to map covariant frame d.o.f back to

Cartesian d.o.f, and must be assessed at the nodal points. Similarly,

4
- . (k) 7 . (k) 7
'sz(51 11) = 2 [uﬂi{(111)ipiql§lﬂ, z121 "' (121)1Pi‘11n1m122}

i=1

. k 7 . k 2 7 _
+ ”Bi{(1 12)1PI‘P(§1)B, 2121 + (122),P.:‘P$11)3, 2122} + W55(Pi’§b21 4’ Pi’n122)aﬂ
(3.5.64)

. k 7 — 7 . k 7 — 2
4’ 9x5i{(122),-(Pi‘yéz)p, 2121 + N gngBJzz) 7 (112),(P 31412)), z122 4' angﬂplzﬁ}

. k 7 — 7 . k 7 — 7
+ Gypi{(111),-(Pi‘1’£12)3, z122 + angﬂpjzl) - (121)i(Pi‘P(§2)B, z121 + Ngi,nQB]22)}]

Field consistent penalty constraint 6x13 is obtained by expressing Eq. (3.5.57) in terms of
the ﬁeld consistent covariant penalty constraints, 6&3 and C115 from Eqs. (3.5.59),

(3.5.60):

(15$. T1)= 5:131:11 + (71151712
i= 1,...,4; [3: 1,2

By nodal transformation, the nodal covariant rotations, 6g,- and 9111' are transformed back
to the Cartesian coordinate deﬁnitions of the rotations, 9n- and Gyi. Thus, the consistent

covariant penalty constraint, ExB can be written in terms of 9x: and 0,,- explicitly as:

97

4
i=1 . . , , (3.5.66)
+ eyBi{(Jll)i(Nni9§ + 19,91117(121),-(N§,-,,1 — P9112}
+WBi(Pi’§jll+Pivnj12)] [3 = 1.2

By the similar procedure, the consistent covariant penalty constraint, 5343 are:

4
ayﬁ<§9 TI) = Z [eth{—(j12)i(NT|i’§+Pi)j21 +(j22)i(N§i91]—Pi)j22}

i: l . 7 , ... (3.5.67)
+ eysi{(11 1 ).-(Nn.w§ 4' P0121 - (121)i(N§i’n — P9122}

+W5i{Pi’§j2l +Pi’nj22}] B = 1, 2

The consistent shear strain ﬁelds, '7“ and '7” in Cartesian coordinates as well as the con-
sistent penalty constraints, CxB and Cy]; derived in this method ensure to alleviate locking
phenomenon, therefore an exact numerical integration rule can be used to evaluate the

shear strain energy.

3.6 Poisson’s Ratio Stiffening Effect

The element with uniform transverse normal strain predicts deﬂections that are about
10 % too stiff when one sublaminate (one element through thickness) is used in the thin
regime, depending on Poisson’s ratio, even though the element is developed using ﬁeld
consistent shear ﬁelds and edge consistency to prevent shear locking. This phenomenon
was identiﬁed as Poisson’s ratio stiffening effect [77,96,97].

From the assumption of transverse deﬂection displacement, Eq. (3.3.1), the original

form of transverse normal strain is constant through-the-thickness as:

 

e = ’ (3.6.1)

This transverse normal strain can be also obtained from the constitutive equation Eq.
(3.2.2) and is given by:
(k) (k) (k)
(k) _C__13 (k)_ C_2__3 (k)__ C___3_6 (k)+__ 1 —o (3 62)
ea = C—'(k) xx C(k) W Can 0 Can 22 ° '
When the laminate is thin, the laminate is in a plane stress state, and the transverse normal

stress approaches to zero. Eq. (3.6.2) becomes:

(’6) (k)
(k) C 13 (k)_ C23 —e(k)- C36 (’0 (363)
82.2: _C(k) xx _Cg") yy €735)ny . °

Due to the assumption of displacement ﬁeld, Eq. (3.3.1), eggs” "" and 7"" in Eq. (3.6.3)
vary linearly in the thickness direction, while 82:) in Eq. (3.6.1) is constant through the
thickness by the deﬁnition of transverse strain. Thus, the element can not capture the vari-
ation of egg) in the thickness direction. The difference in the through-thickness variation
between transverse normal strain by Eq. (3.6.1) and by Eq. (3.6.3) produces a stiffening
effect of the order of 10%. This effects goes away rapidly when more than one element is

used through the thickness. In the current element, this stiffening effect is removed by

using the assumed transverse normal stress approach as discussed previously.

99

3.7 Numerical Results

The accuracy of the proposed laminated plate theory and the appropriateness of the
ﬁnite element model were assessed by simulating the response of two different simply-
supported laminated plates subjected to a transverse double sinusoidal loading as shown in
Figure 21. The elasticity solution of Burton and Noor [98] and the predictions based on
ﬁrst-order shear deformation theory (FSDT) [10] were used to provide comparison with
the current plate results (FZZBD). For a given location in the plane of the plate, explicit
expressions for the through-thickness variation of stresses are provided by the present
model. In the plane of the plate, stresses were evaluated at the centroid of each element,
where stress is super-convergent [73]. The bending stresses presented in the forthcoming
discussion were evaluated at the center of the element that is nearest the center of the
plate. The transverse deﬂections are obtained at the top surface of the center of the plate
where load is applied. The in-plane displacements are obtained at the edges of the plate
where x = O and y = a/Z (for u x displacement) or y = O and x = a/2 (for uy displacement),
and a is the length of each side of the plate. Due to the symmetry of the panels, only a
quarter of the plate was modeled in each case. A uniform mesh of 8x8 elements was uti-
lized in the plane of a quarter of the plate, and the through-thickness disoretization was
varied to assess the utility of the present models. A length-to-thickness aspect ratio of four
was chosen in the present examples in order to highlight the ability of the present model to
obtain accurate results for very thick laminates via through-thickness reﬁnement of the

mesh.

 

Figure 21. Simply supported laminated composite plate
subjected to sinusoidal loading.

Example 1: CAV Panel

The ﬁrst example is a panel from the US Army TACOM’s Composite Armored Vehi—
cle (CAV). Material properties used in the analysis are given in Table 3. The laminte con-
sists of 55 layers and is nearly two inches thick. The laminate can be divided into three
sections, as described below:
- Inner Shell: Four plies of S-2 Glass/Phenolic fabric with stacking sequence
{[902/02] } (ply thickness = 0.0] in.). Thirty-seven plies of S-2 Glass/855340 Epoxy
tows with stacking sequence {[0/90] , [45/—45/0/90]4 , [45/0/45] ,
[90/ 0/—45/45]4 } (ply thickness = 0.021 in.).
0 Armor Core. One layer of EPDM rubber (ply thickness = 0.0625 in.) and one layer of
ceramic tile with inserts (ply thickness = 0.7 in.).

0 Outer Shell. TWelve plies of S-2 Glass/855340 Epoxy fabric with stacking sequence

101

[0/ 90/ 45/ —-45/ 0/ 90] 5 (ply thickness = 0.01 in.).

To allow comparison of predictions of the current model with the elasticity solution, all
layers are oriented at either 90° or 0° from the reference axis by replacing all 45° plies with
90° plies and -45° plies with 0° plies.

Two types of through-thickness reﬁnement were used in the analyses: (1) a single sublam-
inate, and (2) three sublaminates with a sublaminate for the inner shell, a sublaminate for
the rubber layer, and a sublaminate for the ceramic tile layer and the outer shell.

The inﬂuence of aspect ratios on transverse deﬂection is shown in Figure 22. Eight dif-
ferent aspect ratios were considered to check a wide range of plate bending behavior (Uh:
4, 8, 10, 20, 30, 50, 100, 500, and 1000). When one sublaminate was used, the error in pre-
dicted deﬂection was about 15 percent for the case alh=4. The use of three sublaminates
resulted in more accurate predictions, with an error of less than one percent. It was
observed that the prediction of transverse deﬂection converges rapidly to the exact solu-
tion as more sublaminates are used. In Figure 23, the variation of in-plane displacement,
ux, through the thickness of the laminate is plotted at x = 0 and y = a/2 for the case of
aspect ratio of four. The exact distribution of in-plane displacement is approximately lin-
ear within each layer, and globally piece-wise linear and of zig-zag pattern due to the dif-
ference in- the shear moduli of adjacent layers. It is shown that the current model
containing only one sublaminate captures the trends of the distribution of displacement
very well, while FSDT is unable to simulate such behavior. The bending stress distribution
through the thickness is illustrated in Figure 24. The stress was evaluated at the Gauss
point nearest to the center of the plate. It can be observed that the prediction of in-plane

normal stress by the current model agrees well with the elasticity solution.

102

Table 3. Material properties for TACOM’s composite armored vehicle (CAV) panel

 

 

 

 

 

 

 

 

 

 

 

S-2 Glass] S-2 Glass] S-2 Glass/ EPDM
Phenolic 8553-40 8553-40 Rubber Ceramic
=L Fabric Tow Fabric

E1 1 3.0e6 6.2e6 3.0e6 3.0e3 5.0e6
E22 3.0e6 1 .0e6 3.0e6 3.0e3 5.0e6
E33 1.2e6 1 .0e6 1.1e6 3.0e3 1 .25e5
v12 0.13 0.29 0.13 0.45 0.15
v23 0.18 0.37 0.18 0.45 0.15
v13 0.18 0.29 0.18 0.45 0.15
(312 1.0e6 0.3e6 1.0e6 1.0e3 2.5e6
623 4.6e5 0.3e6 3.9e5 1.0e3 8.5e3
613 4.6e5 0.3e6 3.9e5 1.0e3 8.5e3

 

 

 

 

 

 

 

 

 

The prediction of in~plane shear stress is shown in Figure 25. There is signiﬁcant error in
the armor core region, but predictions of the current model are very good in other regions
of the laminate. The distribution of transverse shear stress, (13,z through thickness is shown
in Figure 26. As assumed, this component of stress is constant in each sublaminate. It is
shown that the use of multiple sublaminates allows the current model to capture the trends
of the through-thickness distribution of shear stress. The transverse normal strain, 8w is
predicted in Figure 27. When three sublaminates are used, the distribution of transverse
strain is nearly indistinguishable from the elasticity solution. In Figure 28, transverse nor-
mal stress, ow versus normalized thickness coordinate at the center of the plate is illus-
trated. As assumed, the model captures about the average value of the exact distribution in

a sublaminate.

103

Example 2: Sandwich Panel

The second example is a sandwich panel which consists of a thick core between two
sets of ﬁve face sheets. The material properties used in the analysis are listed in Table 4.,
and the lamination scheme is deﬁned in Table 5. The core occupies eighty percent of the
thickness of the plate, while each set of face sheets contains ﬁve layers and occupies ten
percent of the total thickness. Two types of through-thickness discretization were used in
the analyses: (1) a single sublaminate, (2) three sublaminates with equal thickness. In Fig-
ure 29, the predicted through-thickness distributions of in-plane displacement are illus-
trated. All predictions were very accurate in the upper and lower face sheet regions, but
there were signiﬁcant errors in the core region when only one sublaminate was used.
These errors were due to the very large transverse shear strains and transverse squashing
in the core of this thick sandwich plate. The nonlinear variation of inplane displacement in
this region could not be captured accurately by a single sublaminate in the core region. For
thinner plates (i.e., a/h greater than about 10), the inplane displacements throughout the
thickness are predicted satisfactorily using a single sublaminate through the entire thick-
ness of the laminate. To assess the bending stress, on, more clearly, the stress distribution
in the upper face sheet region is enlarged and shown in Figure 30. It can be observed that
the current ﬁnite element model using only one sublaminate predicts the stress distribution
in the face sheets fairly accurately, and even greater accuracy is achieved with three sub-
laminates. In Figure 31, in-plane shear stress in the face sheets are predicted. It is shown in
Figure 32 that the average transverse shear stress is captured well by the current model.
Figure 33 depicts the distribution of transverse normal strain, 82,, through the thickness of

the laminate. When three sublaminates were employed, the strain distribution through the

104
thickness is predicted relatively well, capturing the trend very well in the core region.

Transverse normal stress, ow versus normalized thickness coordinate at the center of the
plate is shown in Figure 34. It can be observed that it takes about the average value of the

exact distribution in a sublaminate, and captures the distribution trends as more elements

are used through thickness.

3.8 Summary

An improved ﬁrst order zig-zag theory and an associated ﬁnite element model are pre-
sented. The model allows the representation of a laminate as an assemblage of sublami-
nates in order to increase the model reﬁnement through the thickness, when needed.
Within each sublaminate, an accurate ﬁrst-order zig-zag sublaminate kinematic approxi-
mation is made which minimizes the need for multiple sublaminates for many problems of
practical interest yet easily accommodates through-thickness reﬁnement for predicting
local variations of stress and deformation. The ﬁnite element model is “regenerat ” in the
form of an eight-noded three dimensional element with ﬁve degrees of freedom (three
translations and two rotations) per node. Thus, it is suitable for implementation into com-
mercial ﬁnite element codes. Numerical results demonstrate that the current model is
accurate, efﬁcient, and robust for analysis of a wide variety of thick or thin laminated

plates.

105

Table 4. Material properties for sandwich composite plate

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

L 4 material 1 material 2 material 3 core
E11 1.0e6 r 32.6e6 : 25.e6 50.e3
E22 25.0e6 1 .0e7 1 .0e6 150.e3
E33 1.0e6 l .0e7 1 .0e6 50.e3
V12' 0.01 0.10 0.25 0.01
v23 0.25 0.25 0.25 0. 15
v13 0.25 ' 0.10 0.25 0.15
612 0.50e6 8.21e6 0.5e6 21 .7e3
623 0.50e6 0.5e6 0.2e6 42.0e3
613 0.20e6 8.21e6 0.5e6 21 .7e3

 

Table 5. Lamination Scheme for Sandwich Panel

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. Layer No. dilateﬁal Thickness
1 1 0.010
2 2 0.025
3 3 0.015
4 1 0.020
5 3 0.030
6 Core 0.800
7 3 0.030
8 1 0.020
9 3 0.015
10 2 0.025
11 1 0.010

 

 

 

106

 

1.5

1.4r

1.3

I

 

 

\ — FZZSD: 1 sublaminate
‘. - - F223D: 3 sublaminates

 

 

 

 

 

 

 

-§ 1.2” \ "
-§ -\ '---FSDT
'4
.3. t
N .
E x
.\

1— _______________ \_'_‘_ :..-—— ............
o.9~ _
0.8- .

10‘ 102 1c

a/h

Figure 22. Variation of normalized center deﬂection

versus plate aspect ratio at the center of the plate (CAV Panel, alh=4).

3‘IN

107

 

0.9

l

0.8

0.7 r

0.3 r

0.2 -

0.1 r

 

 

 

    
    

— Elasticity
- - F223D: 1 sublaminate

-— - - FZZSD: 3 sublaminates
----- FSDT

 

 

     
 

l L l L

 

—0.8 -O.6 -0.4 -0.2 O 0.2 0.4 0.6 0.8
U Displacement x 10

Figure 23. In-plane displacement, u,
versus normalized thickness coordinate (CAV (Panel, alh=4).

 

3‘IN

108

 

 

 

 

  
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.9 - -
0'3 ’ — Elasticity a/h 4 -
0 7 _ - - FZZSD: 1 sublaminate - r
-- . - FZZSD: 3 sublaminates
0.6 _ ----- FSDT ‘
0.5 - ’ - -
0.4 ~ 1
0.3 - _
0.2 - -
0.1 - ﬂ
/
o ‘ .
-30 -20 '10 0 10 20

SIGMA xx

Figure 24. Bending stress, om versus normalized thickness
coordinate at the center of the plate (CAV Panel, alh=4).

30

§“IN

109

 

 

 

 

 

 

 

 

 

 

 

0.9 r
0.8 -
0.7 -
0.6 -
0.5 -
0.4 - — Elasticity
— - F230: 1 sublaminate
0.3 r - — . - F230: 3 sublaminates
..... FSDT
0.2 '-
0.1 -
o ' '
_4 4 6

Figure 25. In-plane Shear stress, on, versus normalized
thickness coordinate at x, y = 0, (CAV Panel, alhﬂ).

 

B‘IN

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

110

 

I

I

I

 

 

— Elasticity
- - FZZSD: 1 sublaminate
-- - - FZZSD: 3 sublaminates

 

 

 

 

 

l l l l l J 4

 

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
SIGMA xz

Figure 26. Transverse shear stress, on, versus normalized
thickness coordinate at x = 0, y = 21a (CAV Panel, alh=4).

 

:‘IN

111

 

 

 

 

 

 

 

 

 

 

 

l_
0.9 - ' I -
I I
I i
0.8 - : ‘ _
I
0.7 ~ I .
I I a/h = 4
I i
0.6 r I l r
I i
0.5 " , ’ v / ..
0.4 - -
0-3 " — Elasticity ‘
- - F2230: 1 sublaminate
0'2 b - - - - FZZSD: 3 sublaminates -
I
0.1 - _
o l l l l I l
..1 2 3 4 5 6 7
Epsilon 22 x 10-5

Figure 27. Transverse normal strain, 82,, versus normalized thickness

coordinate at the center of the plate (CAV Panel, alh=4).

3‘IN

112

 

 

 

 

 

 

 

1 I I I I II I r I VI
I .
. E
0.91. I .
l I
I I
0-8’ — elasticity : '
- - FZZSD: 1 sublaminate . i
07" -—-- Fzzeo: 3sublaminates : i
I i
06* I i
I .
I _!
0.5*' '._._.|__._____._._
I l
0.4- i I
l I
0.3- I '
: I
. I
0.2 ' | I
I l
0.1 - I '
i I
- I
0 J L 1 pl 1L l l J l
O 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
SIGMAZZ

Figure 28.'Iransverse normal stress, ow versus normalized thickness

coordinate at the center of the plate (CAV Panel, alh=4)

 

3‘IN

113

 

 

 

 

 

 

 

 

1 _ I I I
0.9 " ‘ :l \. ..
/
0.8 _ . , ’ _
— Elasticity 5 / ’
.' / _,
0-7 ’ — - F230: 1 sublaminate / /
._ - - FZZ3D: 3 sublaminates ,5 / .
0.6 - : ’ - ,I r
----- FSDT I / _,
I -’
0.5 - [7 x'/ -
a”! = 4 / I: ./.
0.4 ' / I . /ll -
’ .' /'
/ 3 '
0.3 "' / / .5. ‘I ’ ..
/ .. I . /
/ g"
0.2 e I / ., , .- _
/ / ’ /
/ _ I
0.1 - .. J- g -
-1.5 -1 -0 5 0 0 5 1.5
U Displacement —6

Figure 29. In-plane displacement, u, versus normalized
thickness coordinate (Sandwich Panel, a/hd).

3‘16!

114

 

0.95 -

 

 

 

— Elasticity

- - FZZ3D: 1 sublaminate
-- - - F223D: isxblaminates
..... FS

 

 

 

 

 
  
  

 

 

0.85
-30

—20 -1 O O 10 20
SIGMA xx

Figure 30. Bending stress, on, in face sheets versus
normalized thickness coordinate (Sandwich plate, alh=4).

115

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 I I I I I 3 I I I
0.98 r d
0.96 r ' ‘
g aﬂl = 4 g
h 0.94 ~ 5 _
El
— Elasticity 3'
0-92 * — — FZZ3D: 1 sublaminate '1 «
-- - - FZZ3D: 3 sublaminates '1
..... FSDT El
0.9 e " d
0.88 ‘ ‘ J 1 l I
-25 -2o _15 -10 -5 0 5 ‘0

SIGMA xy

Figure 31. In-plane shear stress in the face sheets, 6,), versus normalized
thickness coordinate at x, y = 0, (CAV Panel, a/hﬂ).

15

116

 

 

 

 

 

 

 

 

 

 

 

 

 

1 I M I I I I II T I I

l'

o 9 - .............. I?) _
I
II

0.8 - i ‘
i

0.7 - I '
l

0.6 - I ‘
— Elasticity I

0.5 ~ — - F230: 1 sublaminate "
-— - - FZZSD: 3 sublaminates l

0-4 ” ----- FSDT : d

0.3 ~ I ‘
l
I.

0.2 ‘ II C
I II

0.1 .. ............... - II -
I!

O l l/—-— I l l 1 ll 1 l l
-0.2 -0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7

SIGMA xz

Figure 32. Transverse shear stress, on, versus normalized thickness

coordinate at x = 0, y = 2Ia (Sandwich Panel, alh=4).

0.8

3‘IN

117

 

0.9 '-

O.8 '-

0.7 r

0.6 L

0.5 -

I

0.4

I

0.3

I

0.2

0.1 -

 

   
    

 

I
|
a/h=4 I
I
I

- - F230: 1 sublaminate I!
-- - - FZZ3D: 3 sublaminates

 

 

 

 

 

Epsilon 22 x 10-6

Figure 33. Transverse normal strain, aw versus normalized thickness

coordinate at the center of the plate (Sandwich Panel, alh=4).

3‘IN

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

118

 

 

 

    

— Elasticity
- - FZZ$Dz 1 sublaminate
-- - - F230: 3 sublaminates

 

 

 

I

 

b - — — - — — — - —.---

 

 

 

.5 0 0.5 1

SIGMA 22

Figure 34. Transverse normal strain, ow versus normalized thickness
coordinate at the center of the plate (Sandwich Panel, alh=4).

1.5

 

CHAPTER IV
PLATE ELEMENT
BASED ON CUBIC ZIG-ZAG THEORY

4.1 Introduction

Composite materials are being considered for use in the front end structures of vehi-
cles to help reduce overall vehicle mass and, thus, improve fuel efﬁciency. Acceptance of
composite materials in structural members will depend on their ability to do crash energy
management. Numerical simulations can greatly aid in the design of these critical struc-
tures and reduce the number of crash tests. When the laminated composites undergo crash,
the fracture of composite materials results from four possible major reasons: (i) ﬁber frac-
ture, (ii) matrix cracking, (iii) ﬁber/matrix interface debonding, and (iv) ﬁber buckling or
kinking. These modes of failure don’t occur simultaneously in all constituents but initiate
as micro-cracks at the constituent level, and unite or propagate as external loading (for
example, impact force, heat and moisture) is increased. As the micro-crack propagates
parallel to the ﬁber direction through resin rich region, delamination occurs. The simplest
way to predict delamination is to evaluate the interlaminar shear stress/strain states in each
layer, neglecting the stress/strains in the constituent level. The discretization in the thick-

ness direction gives a good prediction of the interlaminar shear strains, but makes the com-

119

120

putational cost expensive (especially for dynamic analysis). For this purpose, the
development of more accurate theories and elements to predict interlaminar shear stress/
strain is indispensable.

A new ﬁnite element, which is based on laminated plate theory with cubic zig-zag
approximations, has been developed to model the relevant mechanics that occur in com-
posite materials during crash events. In the plate theory, the in-plane displacement ﬁelds in
a laminate are assumed to be piecewise cubic functions and vary in a zig-zag fashion
through the thickness of the laminate. The zig-zag functions are obtained by satisfying the
continuity of transverse shear stresses at layer interfaces. This in-plane displacement ﬁeld
assumption accounts for discrete layer effects without increasing the number of degrees of
freedom as the number of layers is increased. The transverse normal strain predictions are
improved by assuming a constant variation of transverse normal stress through the thick-
ness in a laminate. The ﬁnite element is developed with the topology of an eight-noded
brick. Each node has ﬁve engineering degrees of freedom, three translations and two rota-
tions. Thus, this element can be conveniently implemented into general purpose ﬁnite ele- ,
ment codes. Numerical performance of the current element is investigated by simulating a

plate with a damaged layer.

4.2 Formulation of Plate Theory

A laminate is assumed to be composed of N perfectly bonded layers. Each layer in the
laminate may have independent thickness and material properties. The global Z-axis is
taken perpendicular to the inplane X, Y coordinate axes and has its origin at the bottom of

the laminate (see Figure 35).

121

 

     

 
 

‘0

  

.OYO_0.0.0.0.0.0...0.0.0.0.0.0.0.0...0.0.0.0.-

   
    

 

     
  

    
      

 

 

sansmmx
Z] M
99.3.0°3’9’3’0’3’9’3’0’3’0’3’0’3’0’3’.’3’e’£ '3?"
20 39303030?of.to203030fofofe?03030302030302.30365‘ _
X

Figure 35. Schematic of a laminate and layer divisions
in the cubic zig-zag laminated plate theory.

To develop the element, three dimensional kinematic assumption were utilized, and the

three dimensional constitutive equation in the kth layer can be written as:

.3: "cl: Ci? c2? 0 0 CE? 8553
61? egg) cg; egg) 0 o cg? e33
(’0 k k (k)
<o::)> = C(13) C(23) Cg]? (3k) (3k) Calgtéi)» (4.2.1)
Cyz 0 0 0 C44 C45 0 sz
of? o o 0 C5,? cg? o yﬁﬁ’
as pit) c9: cl: 0 o c221 Liz:

 

 

 

 

 

 

For most practical laminated composites, the 13 coefﬁcients of the material stiffness

matrix are obtained from 9 independent material constants and a transformation law. The

122

stress - strain relationship is deﬁned as:

(k) (k) (k)
8(k) _ 93x 8(k) _ 33y a“) _ 3‘1:
xx“ x ’ w‘a ' zz’az
(k) (k) y(k) (k) (k) (k) (4'22)
Y(k) = ﬁe + ﬂy (’0 = 3_“z + is , (k) = fl“), +911.
yz y 82 ’ ‘2 8x 82 ’ ’0’ 8x By
The displacement ﬁelds are initially assumed in the following form:
2 3 k-l
k
“(x )(X, y: Z, t) = ub+sz+Z ax+z Bx+ 2(Z_Zi)§i
i= 1
2 3 It—l
k
"(y )(x, y, z, t) = v,, + zwy + 2 ply + 2 By + z (z — zimi (4-2-3)

i=1

(’0 Z Z
uZ (x, y, z, t) = wb(l — II) + w,(,—1)

where h is the thickness of the laminate, ab and vb are the axial displacements in x and y
directions, respectively, at z = 0, w, and w, are rotations of the normal at z = 0, and wb and
w, are the transverse deﬂections of the bottom and top surfaces, respectively, of the lami-
nate. ax, cry, Bx and By represent the high order rotations of the normal at z = 0, and
account for warping of the normal, which occurs when the material properties in, each
layer are signiﬁcantly different. Comparing Eq. (4.2.3) with Eq. (3.3.1), it can be observed
that this plate theory is reduced to the ﬁrst order zig-zag plate theory by setting these high
order rotations to zero. It is assumed in Eq. (4.2.3) that a piecewise linear zig-zag function

(’0 (k)
x and u),

k)

is superimposed on cubic displacement function, so that u vary in a piecewise

cubic fashion through the thickness of a laminate, and it: varies linearly through the

thickness. The resulting transverse shear strains vary a quadratic fashion through thick-

123

ness, so interlaminar shear stresses can be predicted accurately. The parameters 5, and 11,.
of Eq. (4.2.3) are removed by enforcing the continuity of shear stress at each interface [33,
34]. At the kth interface in the laminate, the shear stress continuity condition can be given

by:

(k)_ (k+1) (k)_ (k+1)
yz - Gyz _ 62x " O.er _ (4-2-4)
Z - Z} Z " zt

0'

The transverse shear stress - strain relationship can be obtained from the constitutive equa-

tion (4.2.1) and is deﬁned as:

(k) _ (k) (k) (k) (k) '
0'yz ' C44 sz +C45 7x2

(k) _ (k) (k) (k) (’0
0x2 - C45 sz + C55 sz

(4.2.5)

By differentiating displacement ﬁelds with respect to z, x or y properly, the transverse
shear strains are obtained, and the relationship between the transverse shear strains in the

kth layer and (k+1)th layer are derived as:

(k+l) 2 k aw” Z aw‘z
'sz = Wx+220Ix+3z Bin“ 2§i+a—x (1';)+§} (5)
i=1

0‘)
= sz +§k

k
(“U __ 2 3W1; 2 8w, 2
yyz — wy+22ay+3z By-I- 211,43 (l-Z)+$ (’7)

i=1

(4.2.6)

(k)
= sz + TI];

The transverse shear stress continuity condition is rewritten as:

124

 

 

 

 

k k k k k l k It It
(C C()7;z)+ (32.-2:37;)” _ {C( + )(Y( z)+TIk)+C(5('Y) )+§k)}I(z=z
(k) (k) C(k) (k) (k l) (k) (k)( ,YUc) (4'2'7)
4.
(Caryn szsssz)| _ = {C45 (7,, +n,.)+C55 z..+§,.)}|_
- zk
The following equation is obtained by solving (4.2.7) for 11k and Q:
k k
7115 = ay;2)+by() 428
(k) (k) ( ' ' )
gk = 6"sz +d‘sz
where
. l k k 1 k
a, = At l(C5’5+ ’C‘ 44’-Cf,5*c ’ 5’)—l‘
+
.. 1 k l k
b5 = Ak11(cg’5+ )C( 5’—Cf,’;* ’C5 5’)
+
Ck = A 1(C(k+1)C(k)_ £§+1)Cg:)) (4.29)
k-i-l
3k = A 1(C(k+l)C§k)_C£§+l) C(k))_1
k+1
Ah] _ C(k+l)C(5k)__ C(k+l)

Substitution of Eq. (4.2.6) into Eq. (4.2.7) gives recursive equations about E,- and 11,- , and

solving the recursive equations for £1- and 11,. renders:

(4.2.10)

125

where

aik
J'=I j-l

bik =

a,[f ,,+za]+b,.zc
125+; IN:

MI
. .,( 5.

(4.2.11)
]+dkkzlc:
j=1 j-l
dik=dk[f fik+kzljdj+ckkzlbij9 (i=192’394)
j-1 =1
fik = "(Z/Ji-

From Eqs. (4.2.9)-(4.2.11) it can be observed that Q and n,- depend on the ratios of shear
properties between adjacent layers and the shear deformation in the laminate. If the differ-
ence in shear stiffness quantities between two adjacent layers or the shear deformation in
the laminate is small, then discrete-layer effects will also be small.

Convenient ﬁnite element models (i.e. regenerated elements) can be developed by
writing the rotational variables W5 and 115 in the displacement ﬁeld in terms of the trans-
lational variables a, and v,, the inplane translations at the top surface of the laminate
[41,42]. Thus, in-plane displacement is described more conveniently by the translation at
two points using surface quantities. Moreover, since the existence of the ﬁrst derivatives of

transverse deﬂections,

!
37 , $6 "a? and 8_y (4.2.12)

in the displacement ﬁelds result in the presence of their second derivatives in the strain

126

energy functional, Cl continuity of w,, and w, is required. In order to facilitate the develop-
ment of more convenient ﬁnite element models, it is desirable to alleviate such a require-

ment by introducing new rotational degrees of freedom as follows:

(4.2.13)

These relationships will be imposed as constraints in the energy functional. Furthermore,
ax, (1y , Bx and By can be eliminated by using the free shear traction conditions at the

top and bottom of the laminate as follows [15]:

(1)

 

 

 

 

 

= (ataxia-c2355.?) - o

55:.) -—- <01§’v§’:’+ c1251?) — 0

1:2,) =(Cﬁﬁ’y§§’+c§§’y§’;’)l =0 (4.2.14)
z=h z=h

55:.“ = (castucliiitblzz, = 0

Substituting the modiﬁed displacement ﬁelds into Eq. (4.2.14), solving for ax, 01y , x
and By , we can attain ﬁnal modiﬁed displacement ﬁelds. Here, we need to pay attention to
the fact that since the free shear traction conditions were utilized, the discretization

through thickness is impossible in ﬁnite element model. Therefore, the displacement ﬁelds

can be written as follows:

(k) - (k) (k) - (k) (k)
= u B‘Daﬂ , My =uaB‘I’aB, “Z =WBQB, (1 = 1...4,B = 1,2 (4.2.15)

where the index [3 is used to denote the top and bottom of the laminate with 1: bottom and

2: top, and (1)23, if; and “B are shape functions in the thickness direction, which are

127

explained in Appendix B. Unless noted otherwise, summation on repeated indices is

implied. The variables Cap are:
l-llB = “3, I225: VB, E3B=GIB,T¢4B=9yB (4.216)

The transverse normal strain, which is constant in the laminate due to the assumption
of displacement for thickness direction, can be improved by instead assuming a constant
transverse normal stress <1—ZZ through-the-thickness of the laminate. 6'; can be deter-
mined using Reissner’s mixed variational principle [81]. From Eq. (4.2.1), the transverse
normal strain is obtained as:

e“) - figs“ - £8“) — 92¢“ + io— (4.2.17)

22 ‘ (15) xx (k) W (k) I? (k) 21
C33 C33 C33 C33

o—ZZ is then determined by analytically solving the following relation:

N 1k
..k k
0 = 2 j (ea—cum (4.2.18)
k=lzk—l

Substituting the strain - displacement relations and Eq. (4.2.17) into Eq. (4.2.18), the con-

stant transverse normal stress (I, is found to be:

128

_ aub

Bu, avb 3v, 39y];
ozz = 8—x (P1+56)+$(P2—56)+3;(P3+S7)+a—x(—P3+SB)+§; (1044.59)

39),, 39x), 39x, aub aut
+5; (P4+59)+-a; (P5+S‘°)+$ (P5+Slo)+§;(Q1+51)+§;(’Q1+52)
(4.2.19)
+-a—y(Q2+S3)+$(Q3—S3)+$ (Q4+S4)+a—y (Q4+S4)+a_y (Q5+Ss)
39x:
+9.; (Q5+S5)+Wb(-T1)+W:(Tl)

where the coefﬁcients are deﬁned in the Appendix B. The newly deﬁned transverse normal

strain can now be written as:

(k) _ (k) — (k) . _
822 = “GB: xlxaﬂ + uaB’ ylyaﬁ ‘l' WBXB a = 1. . .4, B — 1, 2 (4.2.20)
where the functions 2.225 , 14,23 and x5 are deﬁned in the Appendix B.

The current plate theory has the following features:
1) A piecewise linear zig-zag function is superimposed on a cubic in-plane displacement

(k) (k)
x and u),

k)

function, thus, u vary in a piecewise cubic fashion through the thickness of

a laminate, while u: varies linearly through the thickness.

2) Discrete-layer effects are captured within a laminate by allowing a piecewise (layer-
wise) continuous through-the-thickness variation of the in-plane displacements.

3) The transverse shear strains vary in a piecewise quadratic fashion through thickness,
and interlaminar shear stresses can be predicted accurately.

4) Continuity of transverse shear stresses through the thickness of a laminate is satisﬁed.

5) By introducing rotational degrees of freedom, the C0 plate theory was developed.

6) Because of utilizing homogeneous shear traction boundary conditions, reﬁnement

through the thickness can not be accomplished.

129

6) Transverse normal strain distribution through the thickness of a laminate was improved
by assuming constant transverse normal stress.
7) In the displacement ﬁelds, only ﬁve engineering d.o.f. (three translations + two rota-
tions) are present.
The equations of motion, the essential and natural boundary conditions, as well as the
displacement-based ﬁnite element model can be developed from Hamilton’s principle.
The energy functional is modiﬁed here to include the imposition of the constraints Eq.

(4.2.13) via the penalty method. Hamilton’s principle can be deﬁned as follows:
8 l{U—K—W+P}dt= 0 (4.2.21)
‘0

U is the internal strain energy:

N...
lkk kk(k)k kk(k)k kk
U = 2 J §[o( )e( ) +0( )8( ) +6 8( ) +0373 +0er ygx)+o§y)7§y)]de (4.2.22)
k

xx xx yy yy 22 zz

where Vk is the volume of the k-th layer. K is the kinetic energy:

N

1 k .k2 .k2 .k2 .

K: 2 I§p()[ui) +14) +u§)]dvk (4.2.23)
k=1V,, '

where pm is the density of the k-th layer and a superposed dot indicates differentiation

with respect to time. W is the work of external forces:

w = JwaPadQ+ ItiuidI‘
Q r (4.2.24)

130

P is the penalty constraint:

P = all zl{(u3p_ggﬂ)2+ (u4B+g_:B)2}:|d.Q (4.2.25)

where y is the penalty parameter. As 7 is assigned successively larger values, the con-
straints of Eq. (4.2.13) are satisﬁed more exactly in the least squares sense.

Substituting Eqs. (4.2.2), (4.2.15), (4.2.20), (4.2.22), (4.2.23), (4.2.24), and (4.2.25) into
Eq. (4.2.21), taking the variation of Eq. (4.2.21), and integrating by parts, Eq. (4.2.21) can

be rewritten as:

_yxaﬁ

aNmB_ aNyaB_ anB {:wa a_nyaB 3R
:1 1““ 6"— ay ax 197 ”WWW 37 a7

de 32 u, 32 u,
““7( (NB a_'y MM: N( “013+?— B)84tz'*'— 2pIaBrp+ _2pIaBrp}w

at2 3:2
3R 8R ya—u _322 w 3_u +3_2 w +8—22 wp

(4.2.26)
+ §[5aaB{Nxaan + N yaB”y + N :0,an + szﬂny + nyapny + Rymbnx}p
S

__ aWB _ 3W5 ‘
+ 5wB Ryany + Rxan + Yny(u3B —$ )- ynx(u45 +$ ) ds

N 21
k) - k
+J‘ 2 I (0(11- njau (15(1) (43+ C(21)"). -6uaﬂwgcﬂ))dzds:|dt = 0

3k: lzk-l

(a, r =1...4;[3,p =1,2)

where NxaB’ NyaB’ NZGB’ NZGB’ QZGB’ anp, Qxap, nyaB’ Ryxaﬁ’ RZB’ RyB, Rxﬂ are the

stress resultants shown in the Appendix B, nx, ny are the directional cosines,

131

1:5”), 1:15”, 13;, are the plate inertias deﬁned in the Appendix B, and 8013 is the Kro-
necker delta. Considering that variations of the generalized displacement degrees of free-

dom are independent for r0 < t < t1 and zero for t = t0, t1 , the equations of motion are

derived as:

anyaB 3R

_ny13
NB Qde + 3)? +37

31v “+3 aN 3352,1539

_ . zaB_
5““B' ax ay ax ay Q

8w 8w 3 it, u 32 u,
+ 7(‘7045 ‘ 533)531+ 707043 "' 3; BJSM = _2' I” «W +— at2 2pIaBrp (4-2-27)

The corresponding boundary conditions are:

Essential B.C Natural BC.
3013 (Nxaan '1' Nyaﬁny + Nzaﬂnx + QZaBny + nyaﬂny + Ryxaﬂnx)
N z (k) (k) (on (k)
+ z I (011-711(1) B + oz}. nJ‘PaB)dz (4.2.28)
k ' 1231— 1
W13 Ryﬁny + Rxan + yny(u3ﬂ -a—y )— —7n (“43 +— a—x )

4.3 Finite Element Formulation

Because the degrees of freedom represent quantities at the surfaces of the laminate, to
take a major advantage of the present formulation, the ﬁnite element model was cast in the
form of an eight-noded element with ﬁve degrees of freedom per node (see Figure 36).

From Eq. (4.2.3), it can be seen that the primary variables are functions of x, y and t only,

132

thus we need to consider only two dimensional shape functions for the ﬁnite element

model.

14,, v,, w" 9”, By, 14,, v,, w,, 9,“, 9y,

 

       
     

m ﬁr—wqw~\- "... -_¢.- __.-. 1.... r
‘ mm 5 3""
r . '

ActéeudLW—i‘nau. “34.52.355-41 53.13.... .--

 

  

 

   

"3.4.

. f
. ,. «v.2. . :. -41.“? gr. —-,_ ..ta'
”crab-1%)); '{ try-31011. VIZ: ”M::1Lﬁt% 221* ﬁ’ "5'... ‘4... _... ._.. .
.53.” ‘3; W!“ :F “.2 x $15; In: ‘5.

-_'r"":-"::‘Tr .31; as". dwarf...

l
6

     

uba vb, Wbr exbr eyb ub, Vb, Wb, 9x1” eyb

Figure 36. Element topology and nodal degrees of freedom.

Inplane displacement and rotation degrees of freedom are approximated by the bilin-

ear Lagrange interpolation functions:

4 4
7‘15- 2%”: '7213 - 2%th
I: i:
4 4 (4.3.1)

where

= in +§.§)(1 +n.~n) (4.3.2)

133

and 2;, n are the element natural (local) coordinates. For transverse deﬂection degrees of
freedom, an interdependent interpolation concept similar to that developed by Tessler and

Hughes [75] is utilized. The interdependent interpolations functions are:

4
WB = 2 (Piwﬂi + inexBi + Nyieyﬂi) (4.3.3)

i=1

where N“. and Nyi are:

l 1

in = ‘§()’12N5 ‘3’41N8)’ Nyl = §(x12N5 ‘x41N8)
1 1

NxZ = "§()’23N6 'y12N5)’ Ny2 = §(x23N6 ‘x12N5) (4 3 4)
1 1 ° '

Nx3 = -§()’34N7 ‘Y23N6)’ Ny3 = §(x34N7 _x23N6)

l 1
Nx4=-§(y4lN8_y34N7)’ Ny4= §(x4lN8-x34N7) xij=xi-xj; yij=Yi—yj

where N i are the quadratic serendipity shape functions.

The developed element locks in the thin regime even though the interdependent inter-
polation scheme was utilized. This shear locking phenomenon is mainly attributed to
inconsistent terms in the transverse shear strain ﬁelds with respect to the in-plane coordi-
nates, x and y [73,74]. The inconsistent terms do not vanish in the thin regime of the ele-
ment, which causes locking. The inconsistencies in the thickness coordinate, 2, also exist
but can be neglected relative to the inconsistencies in in-plane interpolation. Therefore, we
need to ensure only the transverse shear strain ﬁelds to be ﬁeld consistent with respect to
the inplane coordinates in order to prevent shear locking. This is performed using the

approach presented in [73,74]. Through the similar procedure as in the previous chapter,

134

the consistent shear strain ﬁelds in covariant frame are obtained as follows:

- (k) 1'
YXZ(§’ 1‘) = [Piuéaﬁfpgaﬁ’ Z + (Pi’gwﬂi + N§i9§6§Bi)QB]111
.. (k) -.
+ [Piunaai‘bnaﬂz + (Pi’ﬂwﬂi + Nni'nenﬂi)ﬂl31112

where i=1,...,4; 0t=1,2; B=1,2

- k ,
7,.(5. n) = [Piugaﬁpréga , + (Pi’ﬁwﬁi + N§,.§e§,,,-)n,31121
+ [Pianamwnaﬂ z + (mem + Nnmenﬂimp] in

where i=1,...,4; a: 1,2; [3: 1,2

The consistent constraint ﬁelds are recast into the following form:

+ [Pi’nwﬂi + (Néivn - Pi)enﬂi1j12

Cyﬁ(§a Tl) = [Pi’éwﬂi + (Nnigg + Pi)9§B,]j21
+ [Pi’nwﬁi + (”gm - Pi)9nﬁi]j21

(4.3.5)

(4.3.6)

(4.3.7)

(4.3.8)

To make sure that an element is free of locking, still the derived consistent shear strain

ﬁeld and consistent constraints are required to satisfy the edge-consistency on the inter-

element edges. The edge consistency is fulﬁlled by mapping back the covariant-based

shear strain ﬁelds and the covariant-based penalty constraint ﬁelds to Cartesian coordinate

at element nodal coordinates. The edge consistent shear Strains are:

135

4
_ . k 7 . (k) 7
'sz(§1 11) = 2 [uﬂi{(111)ipid)(§l)ﬂ, 2111 '*' (121),?1‘1’1115, 2112}

i=1
° P111”) 7 +(' )P-<1>(k) ‘ + .P- f11+P-2j120

+VBi (112); 1 {3113,2111 122,- 1 7113,2112 W51 1’: In 13

(4.3.9)

. k 7 7 . k 7 7
+ exﬂi{(122)j(Pi¢(§2)B, 2111 + ngnﬂpllz) 7(112),-(P5¢f12)p,2112 + Nni2§93111)}

. k '. 7 . k a. 2
+ ey5i{(111)i(Piq>$12)B,z-’12 + Nni’éﬂﬂjll) ' (121 ).-(Pid’(§2)5, 2111 + N§i9n03112)}]

4
— . (k) 7 . (k) 7
yyz(§’ '1) = 2 [“ﬁi{(l11)2-Pi‘¥§1|3, 2121 +(121)2-P1‘Pn13.2122}

i= 1
. k 7 . k 7 7
+ vﬂi{(112)ipi‘y(§1)a, 2121 + (122),-P2"P1(q1)3,2122} + W5;(Pi’§b21 "’ P 1221122)!)‘3
(4.3.10)

. k 7 7 . k 7 7
+ 9.2.{1122),.(P9P§2’2, zJ21 + N2...,a.2122) - (1,2),(P2a'gz’p, .122 + N,,-.2522122)}

. k ‘7 7 . k 7 ..-
+ 6y3i{(111),-(P.-‘I’$,2)5, 2122 + angﬂﬂm) - (1201032313, 2121 + Ngimﬂﬂ 122)”

where jij and ii}- are deﬁned in Eqs. (3.5.32) and (3.5.33).

The edge consistent constraints are:

4
aid“; 11) = 2 [exﬂi{—(j11)i(an§+Pi)jll + (j22)i(N§i2n -P,-)_;12}

i= 1
+ OyBi{ (111)2(Nnvg + P291711 7 (1'21),-(N:;,-,11 - P9112}
+WBi(Pi’§jll +Pi2nj12)] B = 1, 2

(4.3.11)

136

4
a)“; 71) = 2 [9x3i{-(j12)i(Nm-,§ + P0172: + (j22)i(N§i’T| 7 P0132}
i: 1 . 7 , , (4.3.12)
+ eyﬂi{(]ll)i(Nni’§ + Pi)1217(121),-(N§,-.n - P9122}
4' WBi{Pi,§j21 + P;.nJ-'22}] B = 1. 2

Since the adopted method ensures to alleviate locking phenomenon, therefore an exact
number of integration points can be used to evaluate the shear strain energy without caus-

ing locking.

4.4 Numerical Results

The response of simply-supported laminated plates ‘with delamination or ply damage
under a transverse double sinusoidal loading as shown in Figure 37 is presented. In order
to investigate the ability of such a laminate model to capture the local kinematic variations
associated with damage, a very thin layer with compliant material stiffnesses is placed
between two adjacent layers. The material properties used in the analysis are listed in
Table 6., and the lamination scheme is deﬁned in Table 7.. The current plate results
(CZZBD) were compared with the elasticity solution of Burton and Noor [98]. For a given
location in the plane of the plate, explicit expressions for the through-thickness variation
of stresses are provided by the present model. In the plane of the plate, stresses were eval-
uated at the centroid of each element, where stress is super-convergent [73]. The bending
stresses presented in the forthcoming discussion were evaluated at the center of the ele-
ment that is nearest the center of the plate. The transverse deﬂections are obtained at the

top surface of the center of the plate where load is applied. The in-plane displacements are

137

obtained at the edges of the plate where x = O and y = 312 (for ux displacement) or y = O
and x = alZ (for uy displacement), and a is the length of each side of the plate. Due to both
the geometric and boundary condition symmetries of the panel, only a quarter of the plate

was modeled.

 

Figure 37 . Simply supported laminated composite plate
subjected to sinusoidal loading.

A uniform mesh of 8x8 elements was utilized in the plane of a quarter of the plate. A
length-to-thickness aspect ratio of twenty was chosen in the present example in order to
simulate the automotive-type ﬁber reinforced composite thin wall structures. The inﬂu-
ence of aspect ratios on transverse deﬂection is shown in Figure 38. Eight different aspect
ratios were considered to check a wide range of plate bending behavior (Uh: 4, 8, 10, 20,
30, 50, 100, 500, and 1000). The error in predicted deﬂection was about 3 percent for the

case a/h = 10. In Figure 39, the variation of in-plane displacement, ux, through the thick-

138

ness of the laminate is plotted at x = 0 and y = a/2. The exact distribution of in-plane dis-
placement is approximately linear within each layer, and globally piece-wise linear and of
zig-zag pattern due to the difference in the shear moduli of adjacent layers. Moreover it
can be observed that there is a big change in the distribution of in-plane displacement in
the compliant layer. The bending stress, on, distribution through the thickness is illus-
trated in Figure 40. The stress was evaluated at the Gauss point nearest to the center of the
plate. It can be observed that there is a sudden change of distribution of in-plane normal
stress in each interface region, and the prediction of in-plane normal stress by the current
model agrees well with the elasticity solution. The prediction of in-plane shear stress, on,
is shown in Figure 41. It is shown that current model captures the in-plane stress distribu-
tion in interface region very accurately. The distribution of transverse shear stress, (1’Jrz
through thickness is shown in Figure 42. This component of stress is quadratic through the
thickness direction in each layer, and is under predicted compared with the elasticity solu-
tion. In Figure 43 the distribution of the transverse shear strain, 7,2, is shown. It can be
observed that the shear stress is large in the compliant layer, and the prediction of the cur-
rent model is in good agreement with the elasticity solution. Thus, the developed element
can capture delamination with a proper composite interface failure criterion when the
delamination occurs at an interface. Figure 44 shows the distribution of transverse shear
stress, 6)., through the thickness. It can be observed that 0),, varies in quadratic fashion
like on, however it is over-predicted compared with the elasticity solution. The transverse
normal strain, 8w is predicted in Figure 45. The trend of distribution of transverse strain is
predicted very well, but the strain in the compliant layer is overpredicted by approximately

a factor of two.

139

4.5 Summary

An improved cubic zig-zag theory and an associated ﬁnite element model are pre-
sented. The in-plane displacement ﬁelds in a laminate are assumed to be piecewise cubic
functions and vary in a zig-zag fashion through the thickness of the laminate. This
assumption accounts for discrete layer effects without increasing the number of degrees of
freedom as the number of layers is increased. The ﬁnite element model is “regenerated” in
the form of an eight-noded three dimensional element with ﬁve degrees of freedom (three
translations and two rotations) per node. The ﬁnite element is free of shear locking, robust,

and accurate. Thus, it is suitable for implementation into commercial ﬁnite element codes.

140

Table 6. . Material properties for damaged composite plate

 

 

 

 

 

 

 

 

 

 

 

material 1 material 2 material 3]
E111 139.e9 4.3—0e9ﬁ 4.30e6
1322 9.86e9 4.30e9 4.30e6
13.33 9.86e9 4.30e9 4.30e6
v23 0.37 0.34 0.34
v13 0.30 0.34 0.34
v12 11 0.30 0.34 0.34
(31;, 5.24229 1.60e9 1.60e6
(3,2 5.24e9 1.60e9 1.60e6

 

 

 

 

 

 

 

Table 7. . Lamination Scheme for damaged composite plate

 

 

 

 

 

 

 

 

 

 

 

 

 

Layer N 14Eaterial Thickness Angle
1 l 0. 16583 0.00
2 2 0.00100 0.00
3 1 0. 16583 90.0
4 3 0.00100 0.00
5 1 0.16583 0.00
6 2 0.00100 0.00
7 1 0.16583 0.00
8 2 0.00100 0.00
9 1 0.16583 90.0
10 2 0.00100 0.00
1 l 1 0. l 65 83 0.00

 

 

 

 

 

 

 

141

 

1.1 V V V V I I V I j Y I I U T T V I

1.05 7

 

0.95 7

WM elasticity

P
to
I

0.85 7

 

L A L A l

 

 

0.75 2 2 2 2 2 2 2 l1 2 2 2 2
10 10

a/h

Figure 38. Variation of normalized center deﬂection versus
plate aspect ratio at the center of the plate.

10"

ﬁ‘IN

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

142

 

I

I

I

I

 

    

 

— Elasticity
~- - - czzao

  
  

 

 

 

l l

 

 

-O.5 0 0.5 1
U Displacement

Figure 39. In-plane displacement, u,
versus normalized thickness coordinate (alh=20)

143

 

0.9

 

 

0.8 7 7

0.7 +- _

 

 

.0
00
I
1

— Elasticity
-- - - CZZSD _

I

 

 

0.2

 

 

I

0.1

 

 

l l l

0 l
—2.5 -2 —1.5 -1 -0.5 O 0.5 1 1.5 2
SIGMA xx x 106

 

 

Figure 40. Bending stress, om versus normalized thickness
coordinate at the center of the plate (alh=20).

144

 

 

 

 

1 I I I I
0.9 7 _
— Elasticity C}
0-8 ” - - - - czzao ‘

 

 

 

 

 

 

 

x10

Figure 41. In-plane Shear stress, ox}, versus normalized
thickness coordinate at x, y = 0, (alh=r-20).

R‘IN

145

 

0.9 7

0.8 7

0.7 7

0.5

0.3 7

0.2 7

0.1

 

 

 

— Elasticity
~ — - - czzso

 

 

 

 

 

 

-6 -5 -4 -3 —2 -1 0
SIGMA xz

Figure 42. Transverse shear stress, om versus normalized
thickness coordinate at x = 0, y = 21a

146

 

0.9 7

0.8 7

I

0.7

I

0.6

3"
.°
0'!

I

0.3 7

0.2 7

0.1 7

 

 

— Elasticity
-- - - CZZSD

 

 

 

 

 

 

 

0
—0.03

-0.025 -0.02 -0.015 -0.01 -0.005
GAMMA xz

Figure 43. Transverse shear strain, 7w versus normalized
thickness coordinate at x = 0, y = 71a

0.005

147

 

0.9

I

0.8

0.7 7

0.6 7

0.5 7

 

37IN

I

0.4

 

 

 

I

0.2

I

0.1

 

 

 

LO

 

Figure 44.Transverse shear stress, 0,2, versus normalized
thickness coordinate at x = 21a, y = 0

0.9

0.8

0.7

0.6

7 0.5

0.4

0.3

0.2

0.1

148

 

I

I

 

 

— Elasticity
'- ° - CZZSD -

 

 

 

 

 

 

l l l L l I l l

 

—7 -6 —5 -4 -3 -2 -1 O 1
Epsnlon 22 x 10—4

Figure 45. 'ﬁ'ansverse normal strain, ea, versus normalized
thickness coordinate at the center of the plate (alh=20).

CHAPTER V
FREE VIBRATION ANALYSIS

5.1 Introduction

In free vibration analysis, the mass matrix accounts for inertia and represents dis-
cretely the continuous distribution of mass in a structure. In this chapter, consistent mass
matrices for the previously developed elements are derived, and the lumped mass matrices
associated with the consistent mass matrices are obtained employing the HRZ lumping
scheme. On the basis of the derived mass matrices, free vibration analysis of simply sup-
ported laminated composites will be accomplished. Since the derived lumped mass matri-
ces render the decoupled dynamics governing equations in explicit dynamic analysis, they

will be utilized for explicit dynamic analysis in the next chapter.

5.2 Mass Matrices

Inertial effects must be included when analyzing the response of a structure if it
vibrates freely or if excitation frequencies are higher than 1/3 of the lowest natural fre-
quency of the structure. In the ﬁnite element analysis of the dynamic response of a struc-

ture, the mass matrix accounts for inertia of a structure, and represents discretely the

149

150

continuous distribution of mass in a structure.

A consistent mass matrix can be obtained from the kinetic energy Eq. (3.4.3) using the
displacement ﬁelds Eq. (3.3.14) (for cubic zig-zag plate element, Eqs. (4.2.23) and
(4.2. 15)) and the same shape functions are employed as those used to develop the stiffness

matrix [55,56]. The consistent mass matrix can be written as:

N
[M] = 2 H: P(k)[N(k)lT[5/(k) 14de (5.2.1)
k = In "“
where N m is the number of layers in the mth sublaminate or the laminate, pm is the den-
sity of the kth layer, and [II/(1‘)] is the shape function matrix in x, y and z coordinates.

The lumped mass matrix is derived on the basis of the consistent mass matrix. Among
the lumping schemes that produce a diagonal mass matrix, the row-sum lumping, the opti-
mal lumping [57] and proportional lumping (so called HRZ lumping) schemes [58] are
well-known. However, since the current ﬁnite element model includes rotational d.o.f., the
HRZ lumping scheme which can treat rotational d.o.f. reasonably is employed in this
research. The basic idea of the HRZ lumping scheme is that only the diagonal terms of the
consistent mass matrix are scaled in such a way that the total mass of the element is main-
tained. The procedure of the HRZ lumping scheme is as folloWs. First, we need to obtain
the total mass of the element:

N».

Z:
m = j ] p‘k’dadz (5.2.2)
k = Izk_‘ﬂ

Second, calculate only diagonal components ( M ii) of the consistent mass matrix.

151

[M,,]= N; H: p""’[N,"" ’T1[N,-""’1dzdn (5.2.3)

k=ln '

Third, compute total masses, (3", 3,, sw) in each translational direction (u, v, w) respec-

tively:

40 N...
- " (k) (k)
s“- 2 2 [ﬁt—1p [N,""’1[N, ldzdn
1:51-4 1:19
N
2,: E)" .. if" p‘k’[1'i(,""’1[N,-""’1dzdn 524
j=5i—3 k=ln 2"" ""’
40 ..
_ 1.2" (“W-(“ITW-(“l d d!)
sw ' Z 2,1 0...") J 1 Z
j=5i—2 k=lQ
where z: 1...8

Fourth, scale all the diagonal components by multiplying them by the appropriate ratios.

Forud.o.f:
p(k) 7 (k) T 7 (k)
114.. .—. N- N- dd!)
1 11] "nuNEiEp " , 1" I 1" (5.2.5)
where j = 5i-4,i = 1...8
Forvd.o.f:
k) 7 (k) T 7 (k)
[M..]= " [N ][N- 121de
”31.1%"me 1 1 (5.2.6)
where = 5i-3,i = 1...8

For w d.o.f:

152

N,
— ﬂ
[My] 7 SW 2 ‘11 2“ (5.2.7)

For rotational d.o.f, attention should be paid to which .9 will be used. In this dissertation,
the factor, 3, for rotations was obtained by averaging three factors as follows:

_ (su+sv+sw)

Sex = Se), — 3 (5.2.8)

 

Then, by similar procedures, the rotational d.o.f. components in the lumped mass matrix
are obtained. The most important feature of this lumping scheme is that it always produces
positive lumped masses without regard to the element type, because the diagonal terms in
the consistent mass matrix must be positive by virtue of the positive-deﬁniteness property.
For low order elements, the accuracy of this lumping scheme in the prediction of natural

frequencies often surpasses that of the consistent mass matrix.

5.3 Free Vibration Analysis

The natural frequency of a structure can be obtained by solving the dynamic equations

without damping and external force as follows:

[M]{l'i(x, y, z, t)}+[K]{u(x, y, z, 1)} = {0} (5-3-1)

where [M] and [K] are global mass and stiffness matrices respectively. A solution to the

above equation may be obtained by assuming a displacement vector as

153
{u(x, y, z. t)} = {17(x, y, z)}e""°" (5.32)

where {17(x, y, 2)} is the modal vector and a) is the natural frequency of the plate. The
acceleration vector can be derived by differentiating Eq. (5.3.2) twice with respect to time

and can be written as:

{12(x. y, z. a} = -w2{a(x.y.z)}e“°’ (5.3.3)

Substitution of Eq. (5.3.2) and Eq. (5.3.3) into Eq. (5.3.1) leads to the eigenvalue problem:

([Kl-coleMﬁ} = {0} (5.3.4)

Eq. (5.3.4) is solved by using the standard eigenvalue "solution schemes after boundary
conditions are imposed. Thus, the modal vectors ({ﬁ}) and the natural frequencies ((0)

are obtained.

5.4 Numerical Results

The accuracies of the mass matrix and lumping scheme were evaluated by simulating
the response of a simply-supported laminated plate. The example is a sandwich panel
which consists of a thick core between two sets of ﬁve face sheets. In this example, a uni-
form 16x16 mesh was used for a full plate, which is equivalent to an 8x8 mesh for a quar-
ter plate. The reason why a full plate was used for free vibration analysis is that both the
geometric and boundary condition symmetries of the panel do not imply symmetry of all
vibration modes, so that all antisymmetric modes would be excluded by imposing symme-

try. The material properties used in the analysis are listed in Table 8, and the lamination

154

scheme is deﬁned in Table 9. The core occupies eighty percent of the thickness of the

plate, while each set of face sheets contains ﬁve layers and occupies ten percent of the

total thickness.

Table 8. Material Properties for Sandwich Composite Plate

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Material 1 Material 2 Material 3 Core
E1 1 1.0e6 3.26e7 25.e6 50.e3
E22 25.0e6 1 .0e7 l .0e6 150.e3
E33 1 .0e6 l .0e7 1 .0e6 50.e3
623 0.50e6 0.5e6 0.2e6 42.0e3
Gl3 0.20e6 8.21e6 0.5e6 21 .7e3
G12 0.50e6 8.21e6 0.5e6 21.7e3
v23 0.25 0.25 0.25 0.15
v13 0.25 0.10 0.25 0.15
v12 0.01 0.10 0.25 0.01
p 0.065 0.0542 0.0578 0.0047

 

 

The predictions of the lowest natural frequency by FZZBD and CZZBD are shown in
Figures 46 and 47. The lowest natural frequencies were obtained at the various plate
aspect ratios using both a consistent mass matrix and a lumped mass matrix, and normal-
ized with the solution of Burton and Noor [98]. In Figure 46, it can be observed that the
prediction of the lowest frequency by ad hoc lumping scheme is very poor when the aspect
ratio is small. This comes from the fact that the rotary inertia terms are excluded in ad hoc

lumping. Thus we can recognize that keeping rotary inertia terms in the lumped mass

155

matrix is very important when the plate is thick. The biggest error in the prediction of the
lowest natural frequency is less than 1.6% for both the consistent mass matrix and the
lumped mass matrix as seen in Figures 46 and 47, and it can be concluded that treatment
of the total mass factor for rotational d.o.f was correct and that the HRZ lumping scheme
is effective and accurate. In Table 10, the higher natural frequencies in the 16x16 mesh
based on the consistent and the lumped mass matrices are listed. It can be observed that
the lumped mass matrix yields a lower the highest natural frequency than the consistent
mass matrix. Therefore, in explicit dynamic analysis, the lumped mass matrix not only
enables the dynamic equations to be decoupled but also gives rise to larger time step (this
fact will be discussed in the next chapter).

Figures 48, 49, 50, 51 and 52 show the ﬁrst, second, third, fourth and ﬁfth eigen mode
shapes of transverse deﬂection d.o.f. of the simply supported sandwich plate at the aspect
ratio of a/h = 4. From the ﬁgures we can see that all the eigen modes are linearly indepen-
dent. As the plate considered above is an orthotropic plate, we observe different eigen val-
ues for the second and thrid eigen modes given in Figure 49 and Figure 50. But in the case
of an isotropic plate, the mode shapes shown in Figure 53 and Figure 54 degenerate into a
single mode with the same eigen value. This results from the orthotropic material proper-

ties of the sandwich plate.

5.5 Summary

The consistent mass matrices have been derived from the previously developed ele-
ments and the lumped mass matrices associated with the consistent mass matrices have

been also derived. Free vibration analysis of simply supported sandwich panel has been

156

carried out. From the free vibration analysis, it has been shown that the lowest frequency
predicted by the consistent and the lumped mass matrices are in good agreement with elas-
ticity solution. It was observed that the predicted highest natural frequency by lumped
mass matrix was lower than that by the consistent mass matrix. Also from the results it is
seen that the eigen modes for orthotropic materials are distinctly different from conven-

tional isotropic materials.

157

Table 9. Lamination Scheme for Sandwich Panel

 

 

 

 

 

 

 

 

 

 

 

 

 

Layer No. Material Thickness
l 1 0.010
2 2 0.025
3 3 0.015
4 1 0.020
5 3 0.030
6 Core 0.800
7 3 0.030
8 1 0.020
9 3 0.015
10 2 0.025
1 1 1 0.010

 

 

 

 

Table 10. The higher natural frequencies in the mesh

 

 

 

 

 

 

 

 

 

Mode No. consistent HRZ lumping
2478 0. 1375E+08 0.3574E+07
2479 0.1379E+08 0.3581E+07
2480 0.1380E+08 0.3581E+07
2481 0.1380E+08 0.3583E+07
2482 0.1381E+08 0.3583E+07

 

 

158

 

 

 

 

 

 

 

 

 

I * r ' I
1.2— ..
1~ _ ........ I ___—___ _
l
g I
x I
u I
8 0.8- , ..
E 1
I
l
I
0.6- I q
I I
, — Consustent
I ------ HRZ lumping
I - - - Ad Hoc
I
0.41' I -
I
I
I
I
I
0.2 . ...I.I 4 . . .11.41 4 4 14L14_L
10‘ 102 103
a/h

Figure 46.Variation of normalized lowest natural frequency versus plate aspect ratio
by consistent mass matrix and lumped mass matrix (FZZBD).

159

 

 

 

 

 

 

 

 

 

 

I I l
1.2 ~ 3
Consistent
- - - HRZ lumping
1.1 - a
8
><
o
8
a 1 - M —————————— ..
0.9 ~ _
0.8 )- _
O7 . . i A; I 1 r . L I . I I 14 A , A , L , , ,
101 102 103
a / h

Figure 47 .Variation of normalized lowest natural frequency versus plate aspect ratio
by consistent mass matrix and lumped mass matrix (CZZ3D).

 

 

Figure 48.The ﬁrst eigen mode of transverse deﬂection d.o.f of
simply supported sandwich panel (alh=4)

161

 

 

Figure 49.The second eigen mode of transverse deﬂection d.o.f of
simply supported sandwich panel (alh=4)

162

 

 

Figure 50.The third eigen mode of transverse deﬂection d.o.f of
simply supported sandwich panel (alh=4)

163

 

 

Figure 51.The fourth eigen mode of transverse deﬂection d.o.f of
simply supported sandwich panel (alh=4)

164

 

 

Figure 52.The ﬁfth eigen mode of transverse deﬂection d.o.f of
simply supported sandwich panel (alh=4)

165

 

 

Figure 53.The second eigen mode of transverse deﬂection d.o.f
of simply supported isotropic plate (alh=4)

166

 

 

Figure 54.The third eigen mode of transverse deﬂection d.o.f
of simply supported isotropic plate (alh=4)

CHAPTER VI
EXPLICIT DYNAMIC ANALYSIS

6.1 Introduction

The governing equations in general structural dynamics analysis are second order or
fourth order partial differential equations in the space domain, and hyperbolic in the time
domain. In this case, the solution is a function deﬁned both on the space domain and on
the time domain. Methods of dynamic analysis have been developed to solve the dynamic
equations on the time domain, and the popular methods are modal methods and direct
integration methods. In modal method, the dynamic equations are uncoupled, and then
each of the equations can be solved independently of the others [99,100]. In direct integra-
tion method the dynamic equations are discretized in time, and a sequence of simulta-
neous algebraic equations is obtained [99-101]. The direct integration methods can be
divided into two types - explicit direct integration and implicit direct integration methods.

Since the developed elements are implemented into Argonne National Laboratory’s in-
house code, NEPTUNE [102-104], which utilizes explicit direct integration method, only
the explicit direct integration method is studied in this chapter. In NEPTUNE the internal

force vector is calculated from the developed elements at each time step.

167

168
6.2 Explicit Dynamic Analysis

In structural dynamic analysis, the major goal is to investigate how a structure moves
with time under external excitation, that is the time-history of movement of the structure.

The dynamic equations for the undamped forced vibration of an elastic system under-
going small displacements can be derived by the principle of virtual work and rewritten as
follows by describing displacement ﬁelds in terms of nodal displacement and shape func-

tions.

[M1{a(r)}+[K1{u(t)} = {R"“(r)} (6.2.1)

where [M] is the global mass matrix, [K] is the global stiffness matrix, {Rw(t)} is the

external time varying force vector due to body forces, traction, concentrated forces:

{R"“} = [V [N1T{F}dV+ IS [N1’{‘11}dS+ 2 {p}.-

i=1

(6.2.2)

{F}: body force vector, {‘I’}: traction vector, {p}: concentrated load vector

{u(t)} and {ii(t)} are global displacement and acceleration vectors respectively. The
above dynamic equation is a system of coupled, second order and ordinary differential
equations in time. Although displacement and acceleration vectors are discrete functions
of space, they are still continuous functions of time. Methods of dynamic analysis "have
been developed to solve this equation, and direct integration methods are the popular
methods to deal with a large problem. In the direct integration method, the ﬁnite element
method is used in the space domain, and the ﬁnite difference method is applied in the time

domain. Thus, the dynamic equations are discretized in time, and a sequence of simulta-

169

neous algebraic equations is obtained [99-101]. At a speciﬁc instant of time, the dynamic

equation is:

[M1{ﬁ(t)},,+[K]{u(t)},, = {Rmunn (6.2.3)

where n represents time MM and At is the size of time step.

The explicit direct integration method utilizes the central difference method to solve a
sequence of simultaneous algebraic equations. In the central difference method, velocity
and acceleration at time (MM) is approximated as follows by expanding {u(t)}n +1 and
{ u(t) }n _1 in Taylor series about the time nAt.

{120)}, = flauuunMI—{uunhn

1 (6.2.4)
—,({u(r)},..1-2{u(t)},,+{u(r)},,_,)
At

{‘70)}n
Substituting Eq. (63-4) into 139- (6.2.3), the dynamic equation can be described as:
gl—ZIM]{u(I)}n +1 = {Rext(t) }n - [K]{u(t)}n + Zl—ZlMKZh‘U) }n — {u(t)}n_1X6.2.5)
" t

In Eq. (6.2.5), the term, [K ]{u(t)}n is called the internal force vector and represents
loads at nodes caused by straining of the material. The internal force vector can be rede-

ﬁned as:
{R(t)""'}. = [Kl{u(t)},, = JV [BlrlDllBldV'{u(t)},, = jv [B]T{6},,dV (6.2.6)

where [B] is the strain - displacement matrix, {0'} is the stress vector, and [D] is the

170

material property matrix. It can be seen in Eq. (6.2.6) that the internal force vector at each
time step can be calculated by using the strain - displacement matrix and the stress vector,
and the element stiffness matrix does not need to be calculated and assembled into a glo-
bal stiffness matrix. Thus, the computing time can be decreased by reducing the number of
matrix multiplications, and computer storage can be saved because element and global
stiffness matrices need not be stored. In NEPTUNE the internal force vector is calculated
for the developed element at each time step. Since the internal force vector needs to be cal-
culated at each time step, the use of reduced integration quadrature can decrease computa-
tional cost, but in this case, hour glass control or a stability matrix is required to add
artiﬁcial stiffness to the element to suppress the zero-energy modes [106,107]. Thus, the
reduced integration quadrature will not be used.

If [M] is a diagonal matrix, then the Eq. (6.2.5) is a system of decoupled linear alge-
braic equations. By solving Eq. (6.2.5), the displacements at the end of the time step, n
are obtained.

Since {0'} in Eq. (6.2.6) could be a nonlinear function of strain, Eq. (6.2.6) is valid for
both linear and nonlinear material behavior. Therefore, the explicit method is well suited

to treating material nonlinearity.

6.3 Critical Time Step

The central difference algorithm is relatively simple; however, it is only conditionally
stable. To assure a stable solution, the computational time step must be equal to or less
than the critical time step. Thus, the determination of the critical time step is a key to

reducing computational time and maintaining a stable solution. Overly conservative esti-

171

mates will increase computing times and estimates that are too large will produce numeri-

cal instabilities.

Stability Analysis of the Equations of Motion

The method of assessing stability of the equations of motion based on direct integra-
tion methods can be roughly divided into two categories: spectral (Fourier) and energy sta-
bility analysis. In spectral stability analysis, a single equation of motion is derived by
modally uncoupling the original global equations of motion, and then the effects of a time
integration method on a single equation of motion are examined. While energy stability
method utilizes the original global stiffness and mass matrices, and leads to the conditions
that a norm of the solution at time, t = nAt can be bounded a norm of the solution at
time, t = 0. Since spectral stability analysis usually provides insight to interpret the
results (stability conditions) and sometimes predicts more precise results, the spectral sta-
bility analysis will be discussed.

If the equations of motion based on direct integration methods are stable without any
external loading, they will also be stable with external loading. Thus, we need to consider
only the homogeneous form of the equations of motion in stability analysis. The modally

uncoupled, undamped, and homogeneous equations of motion are given by:
Z’, + (0,22,. = 0 (6.3.1)

where Z, and 0), are the modal amplitude and natural frequency of mode i, respectively.
The subscript i will be omitted for brevity hereafter. Applying the central difference

method Eq. (6.2.4) to Eq. (6.3.1), a central difference approximation of Eq. (6.3.1) at time

172

t = nAt can be written as:
2,,+1 + ((62/512— 2)z,, + z,,_1 = 0 (6.3.2)

The exact solution of a central difference approximation, Eq. (6.3.2), consists of a linear

combination of two parts as follows:

2 = C,7\.’,"+C27L; when 11.0.2 (6.3.3)

2 = C,7t’,'+nAtczxg when 3., = 2., (6.3.4)

where Cl and C2 are constants and determined using initial conditions, and

2.. = . 2““ i = 1, 2 (6.3.5)

I

in which 11, and 112 are generally complex. Z n of Eq. (6.3.3) will decay or remain steady
only if 2., :1: 3.2, Pull S 1 and [22] S l . Therefore, the criterion of time step size, At, pro-
viding a stable computation can be derived by assuming that A, at 12 and selecting initial

conditions such that C 1 or C2 is zero. Zn of Eq. (6.3.4) becomes:
Z = C71." where A" = 81mm 1 (6.3.6)

Substituting Eq. (6.3.6) into Eq. (6.3.2) and dividing by C2." - 1 , the following characteris-

tic equation can be obtained:

2.2 + (m2A12 — 2)>.+ 1 = 0 (6.3.7)

173

The solutions of Eq. (6.3.7) are:

2 2 (2 2
7*. _ 2—(1) At imAt a) At —4 (6.3.8)

1,2 _ 2

Since the characteristic equation of Eq. (6.3.7) is quadratic, we can recognize that the mul-
tiplication of the solutions of the quadratic equation, 2.12.2 is 1. If the radicand of Eq.
(6.3.8) is positive, that is, 2,1(132At2 —4 > 0, then both 3., and 2.2 are real and 2.13.2 = 1
implies that the absolute value of one I. is less than unity while that of the other is greater
than 1. Thus instability results because one 9. does not satisfy Ill 5 l . If the radicand is
negative, that is, A/(DZAIZ — 4 < 0 , then 3., and 2.2 are complex conjugates and each is of
unit modulus by virtue of 2.13.2 = 1. Thus the dynamic equation is always stable. If
MAR—4 = o, 2., = 2.2 = —1 is obtained. This violates 2., $2.2, which is required
for stability. The stablity condition of the equations of motion associated with central dif-

ference method is such that:
At < — (6.3.9)

However the equation of motion of Eq. (6.3.1) is only one of many uncoupled equations.
Thus we need to evaluate Eq. (6.3.1) for every possible modeand select the most conser-

vative critical time step as follows:

At < (6.3.10)

 

max

where (om,Jr is the maximum natural frequency. From Eq. (6.3.10) it can be observed that

174

the size of the time step depends on the highest frequency of the mesh. Moreover, it is
known and mathematically proven that the highest frequency of the mesh is bounded by
the maximum frequency of a single, unsupported element. For structural and continuum
elements a bound on the critical time step can normally be obtained from the highest
eigenvalue of the unconstrained elements. Simple formulas based on element eigenvalues
have been derived for some elements [111,112], however, isotropic elastic material prop-
erties were assumed. The upper bound of the highest frequency of an element can be

obtained by Gerschgorin’s theorem [113].

Gerschgorin ’s Theorem
The approach followed here to obtain an estimate for the critical time step is to apply
Gerschgorin’s theorem to the unassembled element stiffness matrix of each element in the
mesh. The time step, Are”, for the ith degree-of-freedom of node I of element e can be
expressed as [114]:

0.5
z(milil)

mi, 3 (6.3.11)

 

0.5
Kilil+ Z |Kiljll
it j
I at]
where m,,,-, is the element nodal mass of node I in the ith direction, Kim is the diagonal

stiffness, and KiljJ is an off-diagonal stiffness. The maximum time step for the mesh is

determined from:

At = min(Ate.-1) (6.3.12)

crit

175

In order to decrease computational cost, it is desirable to decrease the highest frequency of
a single element. It is known that a lumped mass matrix provides for larger time steps.
Moreover, a lumped mass matrix provides uncoupled dynamic equations and generally

more accurate results than a consistent mass matrix.

6.4 Amplitude and Period Errors in Central Difference Approximation

In order to analyze the errors in using the central-difference method, the modally
uncoupled, undamped, and homogeneous equation of motion, Eq. (6.3.1), is utilized. The

exact solution of Eq. (6.3.1) is given by:

exact
Zn

= C1(coswt + isinmt) + 632(cosmt — isinwt) ' (6.4.1)

where C, and C2 are determined from initial conditions, and i = J3. However, the
approximate solution of Eq. (6.3.4) obtained by the central difference method may have
amplitude error and period error. Amplitude error can be either amplitude increase or
amplitude decay. The amplitude increase can be interpreted as instability, and amplitude
decay can be viewed as artiﬁcial damping or viscosity. Period error can be either period
elongation or period contraction [99]. These errors are illustrated in Figure 55. Since 2.,
and 2.2 of Eq. (6.3.5) are complex conjugates, It, and 1.12 are also conjugate and can be

written as:
It, 2 a+ib 1.12 = a—ib (6.4.2)

Thus,

176

 

 

 

     
      
   

 

 

    

 

 

 

1.5 l r I I I
1 Amplitude 1

*7' 0.5

t:

"D Possible

E central 1’ —>

o ' .

o o - gigflfleézréce I Period

«1 ’ contraction ‘
I—t \ I \

a. ’ \

m ’ \
...4 \

-o 5 -
Q Exact ‘
harmonic
response
..1 -
_1 . 5 l l l I 1
o 0.5 1 1.5 2 2.5
Ti m c

Figure 55.Amplitude and period errors in central difference solution

177

unAt
3.","=el -

A 'b A At . .
ea" "e" n " = ea" (cosbnAt+tsmbnAt)

n uznAt anAt -ibnAt anAt . . (6.4.3)
2.2 = e = e e = e (cosbnAt-tsrnbnAt)

Substitution of Eq. (6.4.3) into the exact solution of a central difference approximation,

Eq. (6.3.3) gives:

N
ll

,, C 2."+C 2."
" 1 2 2 (6.4.4)

C, eanA"( cosbnAt + i sin bnAt) + CzeanA"( cosbnAt — i sinbnAt)

When the exact solution of the central difference approximation, Eq. (6.4.4) is compared
with the exact analytical solution of Eq. (6.4.1), the possible errors can be distinguished. If

a at 0 , amplitude error occurs, and if b ¢ (1) , period error result.

Period Error

Period error, P, is deﬁned as:

_ 21t/(o _

P- 21t/b '-

(6.4.5)

WIS

where 21t/o) and 21t/ b are the period of the actual system and the period of the system
created by the central difference algorithm. In order to characterize the period error, we
need to determine b of Eq. (6.4.5). b is obtained by using Eq. (6.4.3) with n = 1 in the fol-
lowing procedure:

ImO.) _ sin(bAt)

Re(2.) 7 _cos(bAt) = tan(bAt) (6.4.6)

Using Eq. (6.3.8) (remember that the radicand in Eq. (6.3.8) is complex), the following

178

equation is derived:

Im(7L) _ toAtA/4- (02At2 (6.4.7)

R60") — 2-(112At2

Equating Eq. (6.4.6) and Eq. (6.4.7), b is found to be:

2 2
b = itam[""A"“4'"° A" ] (6.4.8)

A 2 — (112m2

Substitution of Eq. (6.4.8) into Eq. (6.4.5) yields the period error of the central difference

method as:

2 2
P = (wAt)/atan[mAw4-m A" ] (6.4.9)

2 — (oz/5:2

The possible period error can be one of the followings:

P > 1 period elongation
P = 1 no period error (6.4.10)
P < 1 period contraction

From Eq. (6.4.9), we can see that the period error is subject to the value of At. Further-
more, the frequency, (1) is the natural frequency for mode i among many possible modes
that a multi-d.o.f structure can have. Therefore the appropriate selection of the time step,
At , can prevent period error only for one mode. In addition, even though At is chosen to
be small enough to maintain the stable condition of the equations of motion, the period

error occurs. Hence it can be concluded that the period error always happens in the central

179

difference method. However it is observed that period error is partially compensated when

a lumped mass matrix is used with the central difference method.

Amplitude Error

The following amplitude errors are possible:

a > 1 amplitude growth (instability)
a = 1 no amplitude error (6.4.11)
(1 < 1 amplitude decay (artificial damping)

Amplitude error in the central difference method can be checked as in the following proce-
dure by obtaining the expression for a. The sum of it, and 3.2 using Eq. (6.4.3) with n = 1

is expressed as:

2. +2. = eaA"(cosbAt+isinbAt) +eaA"(cosbAt—isinbAt)
1 2

 

 

(6.4.12)
= 2eaAtcosbAt
While the sum of 2., and 2.2 using Eq. (6.3.8) is derived as:
2. 2. _ 2 - (112m2 + (11At2,/c02At2 - 4 2 — cozAt2 - (DAIJCDZAIZ — 4
, + 2 — +
2 2 (6.4.13)
= 2 — (oz/5:2
Equating Eq. (6.4.12) to Eq. (6.4.13) provides:

ZeamcosbAt = 2 — 612A? (6.4.14)

Substitution of the expression for b of Eq. (6.4.8) into Eq. (6.4.14) leads to the following

expression:

180

1 — (mzAt2)/2

i
f 2 2
C08["" [mAt 4-6) At ]]

2 — (1)2At2

(6.4.15)

 

aAt = ln<

 

 

If Eq. (6.4.15) is evaluated for 0 S toAt < 2 , it is shown that aAt = 0. Hence we can con—
clude that the central difference method has no amplitude error. However, we need to pay
attention to the fact that this does not imply that the amplitude predicted by the central dif-
ference method agrees with the exact one. When the time step is very close to the critical
time step, the radicand of Eq. (6.3.8) approaches to zero. Thus, 3., and 2.2 are almost
equal and the expressions in the parentheses of Eq. (6.4.4) are almost linear combinations
of one another. In this case, C, and C2 of Eq. (6.4.4) may be very different from C1 and
C2 of Eq. (6.4.1), and hence the beating phenomenon occurs. Thus, zero amplitude error
in central difference method implies that the numerical solution provides a mean value of

amplitude that does not grow or decay in comparison with the exact solution.

6.5 Geometrical Nonlinearities - Corotation Approach

The use of the present element in crashworthiness simulations requires the capability to
handle geometrical nonlinearities. Geometrical nonlinearities can be accounted for using a
co-rotation approach [108-111], wherein a local coordinate system is attached to the ele-
ment and moves in space with the element. For this eight-noded hex-plate element, a coro-
tational coordinate system is embedded in the mid-plane of the element and translates and

rotates with the element. The coordinate system is developed as follows. A set of mid-

181

plane pseudo-node coordinates, 27,-, I7,- and Z,- , is deﬁned between each pair of comer

nodes as follows:

— 1
._ 1 .
Y1: i(y,-+y,-+4) l= 1...4 (6.5.1)
- l
1'= 5(Zi+zi+4)

The orthonormal basis vectors e,, e2, and e3 are computed from the pseudo-nodes. The e,
basis vector is taken to be along the line connecting pseudo-nodes l and 2. The e3 basis is
normal to the cross product of the diagonal vectors between pseudo-nodes 1 and 3 and
psuedo-nodes 2 and 4. The remaining basis vector, e2, is obtained from the cross product

of e3 and e,.

 

pseudo nodes

 

 

global coordinates

Figure 56.The corotational coordinate system

182

6.6 Numerical Results

The numerical convergence and accuracy of the transient behavior of the previously
developed element are investigated through several cases. In order to prove the validity of
the present formulation, the results obtained are compared with those available in the liter-

ature. In the examples, it is assumed that the plates remain elastic.

Example 1. Simply Supported 0" Orthotropic Plate

A simply supported 0° orthotropic plate subjected to a suddenly applied uniformly dis-
tributed load of 0.1 MPa is considered. The dimensions. of the plate are 0.25 m x 0.25 m
with aspect ratio 5. A 4 x 4 mesh is used in a quarter plate model due to the biaxial sym-
metry of the problem. The material properties of material 1 in Table 11 are used. The load-
ing and boundary conditions and the geometry are shown in Figure 57. From the
Gershgorin’s theorem, Eq. (6.3.11) the critical time step, At = 1.2337e-7 sec. was
obtained. However, the larger time step At = 8.542e-7 sec. derived by the trial and error
method was used in the analysis. The deﬂection at the center of the plate is obtained and
compared with Kant et al.’s solutions using nonlinear ﬁrst and high order shear deforma-
tion theories [115], and Swaddiwudhipong et al.’s nonlinear results [116]. It is shown in
Figure 58 that the amplitude and the period response at the center of the plate are in good
agreement with Kant et al. and Swaddiwudhipong et al.’s predictions. The small deviation
in the amplitude and period of the response may come from the fact that the plate is very
thick (a/h = 5) so that the transverse normal effect, which only the current element has the
ability to predict, is signiﬁcant. The top and bottom deﬂections are different, and a

through-thickness (barrelling) mode of vibraiton can be observed in Figure 58.

183

Example 2. Simply Supported Two-Layered Angle-Ply Plate

The behavior of a simply supported two-layered angle-ply (45°/45°) plate is investi-
gated in this example. The thickness of each layer is 0.025 m. The material properties, the
loading and the geometry are the same as in the example 1. However, the biaxial symme-
try can not be utilized even when the plate geometry and boundary condition are symmet-
ric, because the bending is coupled with stretching in an angle ply. Thus, the full plate is
modeled with 8 x 8 mesh. The purpose of this example is to demonstrate the effect of lam-
ination scheme on the transient response of a composite plate. The transient response of
center deﬂection of plate is shown in Figure 59, and it indicates good agreement with Kant
et al. and Swaddiwudhipong et al’s prediction. If Figure 59 is compared with Figure 58, it
is observed that the difference in the lamination scheme causes a big difference in both the

amplitude and period distribution in the dynamic response of a laminated plate.

Example 3. The Eﬂ‘ect of Material Orthotropy on the Transient Center Deﬂection
Analysis
In order to examine the effect of material orthotropy on the transient center deﬂection
analysis, the following case was analyzed. The loading, boundary conditions and geome-
try used are as shown in Figure 57. Due to the biaxial symmetry of the problem, a 4 x 4
mesh is used in the quarter plate. The Materials 2 and 3 in Table 11 is used for orthotropic
plate and isotropic plates, respectively. It is shown in Figure 60 that even though only E , ,’
s are different among the material properties, the effect of orthotropy results in decreasing

both amplitude and period of center deﬂection [117].

184

Example 4. An Impulsively Loaded Cylindrical Sandwich Panel

The panel is clamped at both ends along its length and simply supported at both ends
along the circumferential direction as shown in Figure 58. The dynamic response of the
sandwich panel was obtained in the elastic range when the impact load is applied on the
panel at the initial velocity of 5650 in/sec over a 3.08 in x 10.205 in area. The material
properties and the lamination scheme of the panel are shown in Table 8 and Table 9,
respectively. A mesh of 16x8 elements was utilized in the panel. Using Eq. (6.3.11) At =
2.92459E-07 was obtained as the critical time step. However, the larger time step At =
2.60E-06 was derived by the trial and error method.

In Figure 61, the time history of displacement in the y direction of the panel subjected
to impulsive loading is represented at the middle of the outer surfaces of the panel and at
node 13 (x = -1.31 134E-07, y = 3, z = -9.42). By comparison of examples 1, 2 and 3, it is
shown that the response of the sandwich panel is irregular in the periodicity of amplitude

and period, primarily due to geometric nonlinearity.

6.7 Summary

For explicit structural dynamics analysis, the elements are, implemented into Argonne
National Laboratory’s in-house code, NEPTUNE wherein the internal force vector at each
time step is calculated for the developed element. Gerschgorin’s theorem predicts the
lower bound of critical time step of an element, and the larger time step can be obtained by
try and error method based on the prediction of Gerschgorin’s theorem. Numerical results
demonstrate that the current models are accurate, efﬁcient, and robust for dynamic analy-

sis of laminated plates.

185

Table 11. Material properties

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Material 1 Material 2 Material 3
J
13,, 5 25e11=1 5.25e11 21e9
1322 21e9 21e9 21e9
E33 2169 2169 2169
v12 0.25 0.25 0.25
v23 0.25 0.25 A 0.25
v13 0.25 0.25 0.25
(3,, 1.05e10 8.4e9 8.4e9
c323 1.05e10 8.4e9 8.4e9
G13 1.05e10 8.4e9 8.4e9
p 0.8e3 0.8e3 0.8e3

 

 

 

 

 

 

186

 

Region Modeled

9%"
"1111
coo

 

 

 

 

 

 

 

 

 

0.25 m
l
|
l
| simply supported I
l I thickness: 0.05 m
l |
_L L
4
0.25 m >1
A P(t)
P = 0.1 MPa
>
0 Time, I

 

Figure 57.The loading and boundary conditions
and the geometry of plate

Displacement (m)

187

 

 

x 10
3 h — CZZSD: top
— - - 0Z23D: bottom
2 7 o o Kant et aleSDT '
D U Kant et aleSDT
1 ~ A A Swaddiwuhipong at al -

 

 

 

 

 

 

 

Lto- 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

TIME (second) " ‘°

Figure 58. The time history of Center deﬂection of orthotropic square plate
under suddenly applied pressure loading.

188

 

 

 

 

x 10"5
3 " — czzao: top I
- - - C2230: bottom
2 r U 0 Kant et aleSDT ‘
A A Swaddiwuhipong et al

 

Displacement (m)

 

 

 

l 1 l J I 1 l
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
x 10

 

LML- 1

TIME (second)

Figure 59.The transient response at the center of two-layered angle-ply (45°/45°)
plate under suddenly applied pressure loading.

x10

189

I I

 

I
0|

Displacement (m)

-15 7

 

 

 

orthotropic: top
orthotropic: bottom
isotrOpic: top
isotropic: bottom

 

 

 

...Vf. '/ .
\ k" .,," ""5.
-20..
l I l g I I 1 I
o 0.5 1 15 2 2.5 3 as 4 45
1110‘4
TIME (second)

Figure 60.The effect of material orthotropy on the transient

center deﬂection analysis

 

 

190

simply supported

 

v0 =5650 in/sec

impact area

 

R = 2.9375 in

Figure 61.An Impulsively Loaded Cylindrical Sandwich Panel

191

 

Displacement (in)

 

 

— CZZSD: center
- - - CZZSD: node 13

 

 

 

-5—

-5_ _

 

 

 

TIME (seconds)

Figure 62. The time history of displacement of the sandwich cylindrical panel
at the center of outer surface and at node 13.

CHAPTER VH
CONCLUSIONS

7.1 Summary

The primary goals of the current research to develop the ﬁrst order and the third order
zig-zag C0 laminated plate theories and associated ﬁnite element models are
accomplished. On the basis of the developed elements, static, free vibration, and explicit
dynamic analysis for sandwich and laminated composite panels were performed. The ﬁrst
order zig-zag model allows the representation of a laminate as an assemblage of
sublaminates in order to reﬁne the model through the thickness, when needed. In the third
order zig-zag model, discretization through the thickness of a laminate is not allowed.
However, since in-plane displacement ﬁelds are assumed to vary in a piecewise cubic
fashion through the thickness of a laminate, accurate analysis. of composite behavior can
be achieved without any discretization through the thickness. In particular, the associated
element is very efﬁcient to analyze and predict the crack (delamination) propagation as
time marches in explicit dynamic analysis of composite plates under impact loads. In these
plate theories, all variables are required to be CO continuous, and the total number of
degrees of freedom is independent of the number layer in a laminate. The ﬁrst order zig-

zag theory has some common features of other approaches of laminated plate theories,

192

193

namely the discrete layerwise and the zig-zag plate theories. The cubic zig—zag plate
theory has some common characteristics of conventional cubic zig-zag theories. However,
both the ﬁrst order and the cubic order zig—zag theories can be distinguished from the
already existing theories in that the variables are expressed in terms of surface quantities
in order to facilitate the development of convenient ﬁnite element models. In each of these

models, the transverse normal strain is improved by assuming a constant normal stress.

7.2 Conclusions

With regard to the use of the proposed laminated plate theories and associated ﬁnite
element approximations for static and explicit dynamic structural analysis of sandwich
and laminated composite panels, the following conclusions can be made:

. Zig-zag discrete layer-wise plate theories wherein local effects are accounted for very
well are developed on the basis of 3-D kinematics.

. The continuity conditions of transverse shear stresses are explicitly satisﬁed in the der-
ivation of displacement ﬁelds, and hence transverse shear stresses are continuous
through the thickness.

- The transverse shear stresses in the ﬁrst order zig-zag theory are constant within a sub-
laminate.

- The number of degrees of freedom of the ﬁnite element model associated with the
plate theory are independent of the number of layers in the composites.

- The developed ﬁnite elements are “regenerated” elements in the form of an eight-
noded three dimensional element, and they have only engineering degrees of freedom

(3 translations and 2 rotations) to facilitate implementation into commercial ﬁnite ele-

194

ment codes.

The constant transverse normal strain is improved by assuming constant transverse
normal stress instead, and the Poisson stiffening effect in the associated ﬁnite element
is thus alleviated.

The element stiffness coefﬁcients are integrated exactly, yet the element exhibits no
shear locking. This is achieved by using the consistent transverse shear strain ﬁelds as
well as the edge-consistency of the tangential shear strain on any common inter-ele-
ment edge.

To be used in structural dynamics problems, the consistent mass matrices are devel-
oped. On the basis of the consistent mass matrices, the lumped mass matrices accord-
ing to the HRZ lumping scheme are developed.

The predicted highest natural frequency for the given mesh by employing a lumped
mass matrix is lower than that by employing a consistent mass matrix.

When the plate is thick, rotary inertia terms are found to play a very important role in
the prediction of natural frequencies. Thus, the rotary inertia terms should not be
neglected in the lumping scheme.

Gershgorin’s theorem predicts a conservative critical time step for explicit structural
dynamic analysis.

In structural dynamics problems, the internal force vector is calculated for the devel-
oped element.

Numerical results demonstrate that the current models are accurate, efﬁcient, and
robust for analysis of a wide variety of thick or thin laminated plates in static and

explicit dynamic analysis.

 

195

7.3 Recommendation

There are many needs to improve the proposed plate theories and associated ﬁnite

elements in order to obtain better prediction of static and dynamic responses of laminated

plates and sandwich panels.

The developed elements have the capability to simulate shell analysis by allowing the
in-plane integration points of Gauss quadrature to be varied in each layer. The shell
geometry can not be described exactly by this method, however, reasonable results can
be obtained when the curvature radius is not too small. Thus, by introducing transfor-
mation matrix to transform stiffness matrix and force vector in local coordinate to glo-
bal coordinate, shell analysis can be performed.

In the deveIOped plate theories, it is assumed that the thickness of a laminate is con-
stant. To account for the change of thickness of a laminate due to ply drop-off etc., it is
desirable to add the capability to model non-uniform thickness through in-plane direc-
tions.

For computationally efﬁcient explicit dynamic analysis, an element with a larger criti-
cal time step needs to be developed. This can be accomplished by adding rigid motion
in the transverse direction between the top and bottom nodes.

In order to predict the location and time of delamination commencement due to
dynamic loading, it is desirable to obtain the time history of stress in the composite

plate.

 

APPENDICES

APPENDIX A

l 96
APPENDIX A

- Shape Functions:

The rotational variables 111,, and my are eliminated in favor of the translational variables u,
and v,, the inplane translations at the top surface of the sublaminate. In that process the fol-
lowing variables are obtained.

N,,, N...
A = h+ 201721)",- B = 201—2001

i=1 i=1
N. N... A.(l)
C: 201—296,. D=h+2(h—z,)a, A=A-D—B-C

i=1 i=1

where h is the thickness of mth sublaminate, and a,, b,, c,, d, are deﬁned in Eq. (3.3.11).

On the basis of the variable of Eq. A.(1), new variables are introduced as:

 

 

 

or, = I; on, ___ _g, 03 = {B.C—D_-(A—h)}’ a4=_B'h
A A 2A
A.(2)
B =_g ,_4 B=-C_4 B={B-C-A-<D-h>}
‘ Z’ 2 [1’ 3 23’ 4 271
k-l k-l k—l k-l
i=1 i=1 i=1 i=1
k—l k-l k—l . 1-1
Bl =a2+1322€i+a22di BZ=BZZciZi+a22dizi
i=1 i=1 i=1 i=1
k-l lk-l k-l lk-l
C, = 013+B32c,+(a3+§)2d, C2 = B32c,z,+(a3+-2-)2d,z, A.(3)
i=1 i=1 i=1 i=1
lk-l k-l lk-l k-l

i=1 i=1 i=1 i=1
k-l k-l k—l

k-l
Al=Bl+Blzai+alei A2=Blzaizi+a12bizi
i=1

i=1 i=1 i=1

197

k-l k-l

k- —1 '
B,=B2+B22a,+a22b, 32=B220,z,+a22b,z,

i=1 i=1 i=1 i=1

k-l lk—l k—l 11-1
6, = B3+B3Za,+(a3+§)2b, C2 B32a,zi+(a3+§)2b,z,

i=1 i=1 1:1 i=1

lit-1 k-l lk—l k-l
D1=B4+(B4+2)zai+a42bi D2 = (B4+§)Xaizi+a42bizi

i=1 i=1 i=1 i=1

The shape functions for ail") are:

(bgk) = -AIZ+1+A2, (13(12)=A12-A2, (by: 1)=-§1Z+§2, ¢(2I2)=§12—§2
A.(4)

og’” = o"’2’- = 1212-132, <1>f,"’=<1>"4’2’=—C,z+Z‘2
(k)

The shape functions for uy are:

‘1’",k’=—A z+A2, ‘P",2’=Alz— A2, ‘P"2,’ =—B,z+82+1, ‘1’"22)=Blz- 82

(k) ,,,__(k) (k) (k)_ MS)
The shape functions for ugk’ are:
2 z

The shape functions associated with $2z based on a constant transverse normal stress

assumption are:

2.",’,- = E12+E2, 2.",’2=F,z+F2, 2.22’,=Glz+(—}2, k"22=H,z+H2
11,31 = "132 = 111+12ik7~i§i= A1’32 =71:+72

2% = Elz+Ez, 2.y,’2 = F12+F2, 2.y2’,— -Giz+G2, 2.;"2’2 = le+Hz A-(7)
2.",’;’, = 2.322 = A“), 2.32,: 2.-‘,’f,’2 _ 212+);

Xi = “K1, X2=K1

198

where coefﬁcients associated with the shape functions are:

(k)
.13-

C(3k)‘l

(k)
36 —

(k)C"’

(1k) (1‘)

 

_ C36 +P, + 56
E" = ‘ (210421" C702” C-—"—"’
C33 33
(k1 (k1
— C13 - C36 ~ P 2 S6
F2 C(k1A2 + ‘(71’42 +— (_1.)
33 C33 C33
0:) (k1
Gz—--V)Bz-Tk)(82+1)+ (k)
33 33 33
— C§§’- Cé'fs’ - P3 + 58
H2: (1:132 "(k—132+ (k1
C33 C33 C33
- C","§’_ C¥2C P4 + $9
"2 " C""’ ”Cm 2'" C""’
33 33 33
(k1 (k1
- C13— C36 Ps‘l'Sro
"FWDZ + (7102+ (k1
C33 C33 C33
(k1 (k1
C2 Q, + S,
E; - - (z—k-SAZ— (3-—k)(A2 + 1)+ (k) A-(3)
33 33 33
,.. c121. c_g'2- - o .5,
2 371""2 C""’ 2 + (k1
33 33 C33
(k1 (k1
. 23 36 — Q2 “l” 53
62 ‘ (“(32 + 1’ (1:132 C""’
33 33 33
H- C_(2_3_k)3+ i32- Q3 -
2: C""’ 2 C7132 +_ (k)
33 C33
22 _C_(2_,3)C 13,2- Q4 + S4
= (k1 2 ("7102 + (k1
]_ 1:12:04. Ci;- Qs "’ 55
2: C""’ 2 C—(Ic1D2 (k)
33 33 C33

199

The coefﬁcients deﬁned Eq. (3.3.21) are given as follows:
P, = P3={N 2 R‘3" ,’—[ A x2+x, (1 +Az)]}/A P2 ={N3 2 R§"’[A,x2—x, A2]}/A
k =

N
k

l
R""[ A "
l

IIMz

R-

R§§’[—B,x2+x,(32+ 1)]}/A ng)[31K2-K,Bz]}/A
1

:M2 FMz

R§"’[D,x2 — 1c ,D2]}/A
l

N
2 ng)[-A1K2+K,(1+A2)]}/A ng)[A1K2—K,A2]}/A

R'-

IIMz "Mz

{ .
{ R§"’[— Clx2+x C2]}/AA (9)
k 1

a.
ll
u—o

3M2

R(k)[-131K2+ K 1A2]}/A
l

_

={
=1
={

P4
Q1
Q3
Rg;)[— C1K2+K,C2]}/A Q5
52
S4
S6
53

M2

57: 2R“):31K2+K,(§2+1)]}/ARgi)[31x2—K,BZ]}/A

a-
II

1

N
(k ‘ ‘ k
59 = {ZR Rg4)[- C1K2+ C2]}/A SlO— _{ 21R R()[b1K2-K,D2]}/A
k =1 k =
T1= {1: 2 ng’x ,/h}/A
(k) (k) (k)
where Ry?) - —23 (3’3) — .93 (3’3) — _1. R“) - C36
- (k) " (k) ' (k) 34 - ‘77,
C33 C33 C33 C33
2 2 N
(Zk —Zk- l) (k)
K2 = 2 K] = (Zk-Zk-l) A=-2R33K1

k=l

200

- Plate Inertias:

’aBrp" N2 I p‘”¢‘"’<l>‘f,’dz, 13M_ 2" I p(k1‘,,(k1.,,£:1 dz
= 11:: 1 = 17-1—1
N,— ZR
1‘”,3 = 2 j pmﬂﬁﬂpdz a,r = 1...4; [3,11 = 1,2
k = 11,“,

- Shape Functions in the Covariant Frame.

(7.3.1)

The shape functions deﬁned in Eqs. A.(4), A.(5) are expressed in covariant frame through

the transformation matrix Eqs.(3.5.38) and (3.5.44) as follows:

(10+ 111 (k) _j22 (k) 112 0:)
(DUB = (1’14'1351—2‘122131 ‘I’azB - j—‘I’BB ‘j—‘I’4B
11 11
(11+ 122 (k) 121 (k) 111 (k)
%B=‘I’1+B —¢2l3’ ¢n2B=‘j—¢3B‘j— ‘1’4B
12
111 (k1+ (k) 122 (k) 112 (k)
TUB: —‘*’1+B ‘2’2B’ ‘P:2B= 713‘? 3'75‘2’4B

j
_1_:‘ng1 (Pm

121 (k)+111 (k)
1"B 213’ ‘Pn2B=‘j—2"'3+JB ‘2‘24‘1’

wnlﬂ = 4B

A.( 10)

A.(“)

° Stress Resultants:

NxaB = i IOUQCDUde

1:17-21

N=zaﬂ N2 } 6(k);\'(ka)de,

5:12-74:

R215 = N2 Iaif’xﬁdz,

k = 12:,
RyB ‘1‘? }T(k)§gdz,
= 13:1
RxB = kNZ 3‘“ Tx JQBdZ:
= lzk- 1

201

N=yaB E I 01") ‘1’“)

k=l:k- 1
N" 0:) (k1
QzaB= 2 Ion A de
k=:zk- 1
N" 0:) (k)
anB= 2 3:1” ‘I’aﬁzdz A.(12)
k=lzkz— 1
N" (k1 (k1
QxaB= z Zfrx z¢aBZdz
k=lzk- 1

nyaﬂ= 2'” I 10011300212, RyxaB' E JT(k)‘i’(k)(Iz

k=lZ‘l “1211

APPENDIX B

202
APPENDIX B

- Shape Functions:

The rotational variables 11],, and w), are eliminated in favor of the translational variables u,
and v,, the inplane translations at the top surface of the laminate. In that process the fol-
lowing variables are obtained.

=2(h— z)d,,,B,= z(h- z)c,,,C,=2(h—~ z)b,,,D,= 2(h- 2111,,

i=1 i=1 i=1 i=1

N N N N
A2 = 2 (h—z,)d2,,B2= 201—2,112,, C2: 2 (h—z,)b2,,D2 = 201-2,1112,

i=1 i=1 i=1 i=1

N N N A.(13)
243:2(h 21")d31133=2(h Zi—)c3i’C3=Z(h Z1)b31’D3=2(h 29031

i=1 i=1 i=1 i=1

N N N N

A4: 20' Z1“)d41' ’34: z(h'ziﬁ‘u’ C4: z(h—Zi)b4i’D4= 202-Z1304“

[:1 i=1 '=l i=1

 

 

 

 

 

 

 

 

 

 

A=A,+h, D=D,+h, A=A-13—B,.C,
where hrs the thickness of the laminate, and au, bij’ C1}, (1,1-i , = 1...4 are deﬁned in Eq.
(4.2.11).
13 B, B,C2-f)(h2-A2) B,(h2+D2)-f)B2
(1,= :, (1.2 = —T, (13 = _ , (X4 = _
A A A A
B,C —D(h3 +A3) B,(h3+D3)-bB3 b(A4-A,1+B,(C,—C41
(15- — (I6 = _ , (17 = _
A A A
bus —B 1+3 (D -D B C -13A B D _133
as: 4 1 _ 1 1 4), a9 = 1 4_ 4’ 0‘10 = 1 4_ 4 A.(14)
A A A
C, A C,(h2+A2)-AC2 C,B2-A(h2+D2)
B, = —-:_—, [32 = :, B3 = __ , B4: _
A A A A

3 ~ ~ :-
=C,(h +A3)-AC3 B =B3C,—A(h3+D3) [5 =A(C4—C,)+C,(A,-A4)
A" ’ 6 71 ’ 7 A

 

 

 

 

B3:

3|

WI

(‘11!

:‘ﬂ

Ql

SI

:1
I

203

 

 

A 9-

“1312613011,:2411.’ 24.2

i=1 i=1]

zﬁﬁzicuwzkzdw 32

i=1 i=1

k-l k-l k-

 

, 10=

DI

A

k-l k-l
612 c,,-z,+a, 2 d121,-

i=l i=1

k-l k-l

= B22cuzi+a22duzi
i=1 i=1

1

018+(138+ 1)2c,,+a82d,,— 2C4,

1:1 [:1 i=1

k-l k—l k-

l

-(33 +1); 01121—“3 2 (111214” 2 64121

i=1 i=1 i=1

k-l 1-1 1-1
—B72c,,—(1+a7)2d,,+ 2"“

k-

B72 c,,-Z,+(l+(17)2 (1],-Z,-

i=l i=1
1 k—l

i=1

1—1
2 d4iz1'

i=1 i=1 i=1

1W

1+2d3,

3+kB3Zlcli+a3kzldli+kzld21’E

i=1 i=1 i=1

4+k54261ﬁa4kzdli+ [(2621,

i=1 i=1 i=1

5+kBS2161i+a5kzldli+kzld3v

i=1 i=l i=1

6"”1‘362611'1'0‘61‘2‘1114' kzcm»

i=1 i=1 i=1

i=1 ] i=1 i=1

i=1 i=1 i=1

k— 1 k- 1A(15)
[i3kzlicuz +a32duz + 2d2iz,
i=1 1:] i=1
1-1 k-l k-l

7’2 = 542 €1121+a4 z dlizi+ Z 32121

i=1 i=1 i=1

1-1 k-l 1—1
02 = 352cuzi+a52duzi+zd3izi

1:1 (=1 [:1

k—l . k-l 1-1
= 136 Z C1121+ 0‘6 2 dlizi+ 2 “3111'

i=1 i=1 i=1

5|

11-1 11-1
+3103: 611+ k2 C41412 =‘Blokz 6112102102‘111'7-1“ 2 04121

1:] i=1 i=1

1-1 1—1 1—1 _ 11-1 k-I k-l
1+ 2‘2”]“592911' 2d4,,12 = B92C1121+a92d1121+ 2‘141'7-1'

i=1 i=1 i=1

 

204

1-1 k-l 1—1 1-1
31 = 51+B1 20114411 21211 A2 = Bl zalizi+al Zhuz,
[:1 i=1 t=l i=1
1-1 1-1 k-l 1.1
B1 = 52+522011+a22b11 32 = 52201121+a22b1121
. i=l 1:] 1:] i=1
k-l k-l k-l
6] = 53+(Bs+1)2011+a82b11‘ 2041'
1:1 i=l i=l
1-1 k-l k-l
C2 = -(I33+1)201121‘a82 b11212" 204111
1:] i=l i=1
1—1 k-l k—l

bl = 437437201141 +0‘7)2b11+ Eb“

i=l i=1 i=l

k-I k-l 1-1
D2 = B72a,,z,+(1+a7)2b,,z,- 212412.-
1'=1 i=1 i=l
k-l 11-1 11-1 k—l 11—1 11-1
£1 = B3+B32a,,+a32b,,+ szv 52 = B32alizi+a3zblizi+ 2122121
i=1 1:1 i=1 i=l i=l i=1
11-1 1—1 k-l 11-1 11-1 11-1

F1= B4+B42a,,+a42b,,+ 2021', F2=B4201121+a42b1121+ 202121

i=l i=1 i=l i=l i=l i=l
k—l k-l k-l

k-l 1—1 1—1
G] = B5 + 55 2 “11+ 0‘5 2 1211+ 2 1231" G; = 35 2 “1121'" 0‘5 2 121121" Z 123121

i=l i=l i=l i=l i=l i=l
k-l k-l k-l

1-1 1-1 k—l
H1=B6+B6Zali+a62bli+ 2031', HZ: 562011z1+a62b11z1+ 203121

i=1 i=l i=l i=l i=l i=l

k—l k—I 11-1 k-l 11-1 k—I
21=Bro[1+ z “11]+ [3102 1211+ 2 “41172 = “5102 “1121'““102 blizi_ 2 “4121'

i=l i=l i=l i=l ' i=l i=l

11-1 k-I 11-1 11-1 11-1 11-1
~71=-B9{1+ Zari]’°‘92b11‘ 2b41~72 = B9Zalizi+a9zblizi+ 2124121

i=l i=l i=l i=l i=1 i=l

The following variables are obtained when Eq. (4.2.14) is evaluated.

ab 1, utl, vb], vtl, bel, 01:11, Bybl, Gytl, ub2, utZ, vb2, vt2, 0xb2, 0x12, 0yb2, 0yt2,
ub3, ut3, vb3, vt3, 0xb3, 0x13, 0yb3, 0yt3, ub4, ut4, vb2, v14, 0xb4, 0x14, 0yb4, 9yt4

 

205

The shape functions for “(‘11) are:

(DY? = ub2-23+ub1~z2+(—Zl+ubl~El+ub3-F,+ub2-(_}1+ub4~H1)°z
+1+Az—ubloE‘z-ub3oF2-ub2-62-ub4oi—12

(DY? = ut2-z3+ut1~22+(A,+utl~El+ut3°F1+ut2-§l+ut4oI—il)-Z
—Az—utl-EZ—ut3-F2-ut2-(-;2—ut4-I72

1119;) = vb2-z3-1-vbl~22+(-§,+vbl-E,+vb3~F,+vb2-6,+vb4-I7,)-z
+§2—vb1-E2-vb3-F2—vb2-32—vb4-172

(D212) = v12-z3+vt1-22+(§,+vtl-E,+vt3-F,+v12--C-}1+vt4-F11)-z
—Bz—vtl ~E2—vt3-F2-vt2 - Gz—vt4-H2 A.(16)

111$): 9xb2-z3-1-0xb1~22+(Cl+0xbl~l—i'1-1-0xb3-F1+9xb2-C—}1+0xb4-I—11)z
+E‘z—0xbl~E2-0xb3-F2-0xb2o(_?2—0xb4-I72

111;}; = 91112-23 +0xtl -zz+(71+0xtl ~l-i‘1-1-0xt3-F, +ex12-C,+ex14.17,1z
+72-0xt1~52-0xt3~F2-0x12-32—0xt4ol72

off? = 9yb2-z3+0ybl -zz+(I-),+9yb1 .E,+eyb3-F,+eyb2-C,+eyb4.17,1z
+132—eyb1-E2—9yb3-F2—eyb2-(’;2—eyb4-172

«115,? = 0y12-23+(9yt1)-22+(.71+0ytl T5,+9y13-F‘,+6y12-(—}1+9yt4-I71)z
+72-0yt1-F2—0yt3-F2—9yt2-52-0yt4-ﬁ2

The shape functions for 11;“ are:

 

206

11“,? = ub4-z3+ub3-zz+(—A1+ub1-E1+ub3-F1+ub2-Gl+ub4oﬂl)-z
+A2-ubl-Ez—ub3-F2-ub2-Gz—ub4-Hz

‘I’(,’2) = ut4~23+u13~22+(A1+u11-El+ut3-F1+utZ-Gl+ut4-f11)-z
-A2—utl'32—u13-F2-u12-Gz—ut4-Hz

T2,?) = vb4-23+vb3'zz+(—Bl+vbl-El+vb3-F1+vb2-Gl+vb4-H1)~z
+1+Bz—vb1-E2—vb3-F2-vb2-Gz-vb4-Hz

T2? = vt4°z3+vt3-zz+(31+vtl~51+vt3-F1+vt2-Gl+vt4-Hl)~zz
—32—vtl -£’2-vt3 -F2—v12°02—vt4-H2 A.(17)

1115'? = 0xb4 . 23 + (0xb3) - 12 + (C, + bel -E, +0xb3-F1-1- 0xb2- 61+ 0x174 - 161,1
z+Cz—0xb1~82-0xb3-F2-6xb2-Gz-9xb4-Hz

W212) = 0xt4-z3+(0xt3)-z2+(71+0xt1~El+0x13-F1+0x12-Gl+0xt4-Hl)z
+72—9xtl-Ez—6x13~F2—6x12-Gz—9xt4-Hz

1143’? = 6yb4 - 23 + 0yb3 12 + (D1 + Gybl -E, + Obe - F, + Gbe - G, + 0yb4 . 19112
+D2-6yb1~52-9yb3-F2-9yb2-Gz—0yb4-Hz

‘11:? = 6yt4-z3+6y13~22+(31+6y11 ~31+6yt3-F1+6y12-Gl+6yt4-Hl)z
+32—0ytl-Ez—0y13~F2—0yt2'Gz-0yt4-Hz

The shape functions for “£111 are:

{22 = A.(18)

B'IN

The shape functions associated with 8,2 based on a constant transverse normal stress

assumption are:

207

2$11: R(k)d>(,k)+R(k)‘P(,k)+ R(k)(P, +199)
21112 _ Rgli1¢1lk1+ R13k1,,,1lk1 11:1,,D2 +510)
11921- = R‘,"’<1>‘2*2+R§4"’\P‘2"’+R§" 3(P3+S,,)
1:21— = 4141:) 41111;) saw.»
2“1‘1311__,,1k1(,,1k1 R(k)‘l’(k)+ 11:1, P5 + 513)
11(32_ R0041“) +R(k)‘P(k)-I- R(k)(P6 +314)
lit), = R(k)¢(k)+ R(3k)\P(k)+ R(k)(P7+Sls)
1133’2— = R§’§’¢§*’+R§’§"P§"’+R§§’ (P8+S,6)
253—1 _ R3211) R‘k’o‘“ Know, +5 1)
23:22 = ng’wf’” + R‘k’of") + R§"’(Q2 + 521
133)] = Rg§)\P(k)+ Rg§)¢(k)+ R(k)(Q3+S3)
,1212 _ R131511P1k1 R134k1¢1k1 R("2(Q4+S41
14?, = R""~P""+R§"’<I>§,"’+R§" 3(Q5+551
1:5; _ R(k)\P(k)+ R(k)¢(k)+ R(k)(Q6+S6)
71% _ R‘k’w‘“ R‘“<I>""+ R(k)(Q7 + 57,
1:522-R(k)‘P(k)+R(k)<b(k)+R(3k)(Q8 +S8)

K, = R§"’(T,1

where coefﬁcients associated with the shape functions are:

P1

_T,kleg"’{u12 1hk4+1111 K3+(A,+utl 131411113 F1+ut2 G,

+ur4-H,)1c2+(—Az-utl1722—1113-F2-utz-Cz—ut4ol72)x,}

A.(19)

_T,k21R(3k){ub2 thk4+ub1 K3+(—A1+ub1 1314111113 F1+ub2 G,

+ub4-H1)K2+(1+K2—ubl~52—ub3-I_~'2—ub2--C_?2—ub4-712)1c,}

A.(20)

Q:

Qz

208

- —T,kle§k){vb2 K4+vb1 K3+(—B,+vbl E,+vb3 F,+vb2 G,

+vb4-H,)x2+(E2—vb1-E2—vb3-E2—vb2-62—vb4-H2)x,}
_T,21R§"’{v:2 K4+vt1 x3+(E,+vzl E,+vr3 E,+v:2 G,
+vt4-H,)x2+(—E2—v:1-E2—vt3-E2—vt2c';2—vt4-Hz)x,}

-1, 2R§"’{exb2 K4+6xbl K3+(C,+6xbl E,+exb3 F,+exb2 G,
+ 6xb4-H1)1(2 +162 -6xb1-E‘2—6xb3-f‘2— 9xb2 - 62 - 9xb4 . 172119}
_1, 2 R§’§’{9x:2 - x4 + Oxtl .113 + (i, + Oxtl ,' E, + 6xt3 - E, + 9xt2 - E,
+9xt4 ~1'L'1{,)1<2 + (iz—thl -E2—6xt3 - 755-9th - (72-91::4 - I72)1c,}
_T,2R""{eyb2 K4+6ybl K3+(D,+6ybl E,+eyb3 E,+eyb2 G,
+ eyb4 - H, 1112 + (1‘)2 — eybl -E2 — 6yb3 - EZ— 11yb262 — eyb4 - mm}
.1, 2Rg"){ey:2 K4+9ytl K3+(J,+9yt1 E,+ey:3 E,+ey12 G,
+9yt4-H,)K2+(.72—9yt1~52—9yt3-F2—6yb2~62—9yt4-l72)xl}
_T,2R§"’{ub4 K4+ub3 x3+(—A,+ub1 El+ub3~F1+ub2 G,
+ub4-H,)x2+(A2—ub1-E2—ub3-E2—ub2-Gz-ub4.H2)x,}
_le2R(k){ut4 K4+ut3 K3+(Bl+utl EH-utS- F1+ut2 G1
+ut4-HI)K2+(32—utl~52-ut3-F2—ut2-GZ-ut4-Hz)xl}
—T,21R‘3’;’ {vb4 K4+vb3 K3+(-B,+vb1 E,+vb3 E,+vb2 G,
+vb4-H,)K2+(1+Bz—vbl-Ez—vb3-F2—vb2-Gz-vb4~ﬂz)xl}

—T1k21ng){vt4 K4+vt3 K3+(81+vt1 El+vt3 F1+vt2 G]

Q6

Q7-

Q3-

209

+vt4oﬂl ).x2+(—Bz—v11-Ez—vt3.F2—v12.Gz—vt4-Hz)-x,}
_Tllegk’wxm x4+6xb3 K3+(C,+6xbl E,+exb3 F,+exb2 G,

+ be4 - H,)1<2 + (C, — (3be F, - 6xb3 - F2 — 9xb2 . G, - 9xb4-Hz)1c, }

- —T121R(k){9xt4-K4+9xt3 K3+(7,+9xt1- E,+6xt3 F,+6xt2 G,

+Oxt4-Hl)-K2+(72—9xt1~52-9xt3-F2-9xt2-32—9xt4-92)K1}
_leng){6yb4 K4+Oyb3 K3+(D,+6yb1 E,+eyb3 F,+eyb2 G,

+ 6yb4 . H, )x2 + (b, — Oybl .52 —9yb3 - F2 — eybz - Gz- 6yb4 - £12m}
_leng){9yt4 K4+6yt3 x3+(3,+9y:1 E,+eyt3 F1+9yt2 G,

+6;t4-H,)-K2+(.72—0yt1oEZ—Gyt3-F2-Gyt2-Gz—9yt4-HZ)K1}
—T, 2R§§’{ub2 K4+ubl K3+(—A,+ub1-E‘,+ub3-F1+ub2-(_3,

+Ub4'H1)°K2+(1+Zz-ub1'Ez-ub3'F2-ub2°az-ub4°I—12)K1}

- —T12ng){ut2 K4+utl K3+(A,+utl E1+ut3 F1+ut2 G,

+u:4'H1)'K2+(-Zz-utl -Ez—ut3~72—ut2-52—ut4-I72hc1}

- —T121R(k){vb2 K4+vb1- K3+(—B,+vb1 E,+vb3 F,+vb2 G,

:Vb4'H1)'K2+(§2-Vbl 'Ez—Vb3'Fz-Vb2‘62-Vb4'ﬁzhcl}

- -T1k21Rg:){vt2 K4+vtl K3+(B,+vtl E,+vt3 F,+vt2 G,

+vt4-H1)-K2+(-§2—vtl .Ez—vt3-F2—v12-(‘;2—vt4-ﬁ2)x,}

- _1, zng’wxbz K4+9xb1 K3+(C,+9xb1 E,+exb3 F,+exb2 G,

+ 9xb4 411,)1:2 + (‘62 - exbl - 1'32 — exb3 . F2 — 9xb2 - (‘22 — 9xb4-I-12)K,}

210

N
—T1 2 R§§){6xt2 K4+9xtl K3+(i, +9xt1 -E‘, +6xt3 - F, +9th - Z7,

+6xt4-I—11) - K2+(72-6xt1 -E2—6xt3 - 172—6th - Gz—exm -I72)1c,}

- —T,2R‘3’f,’ {6yb2 K4+6ybl K3+(D,+6ybl E,+9yb3 F,+6yb2 G,

:fyb4-H,)K2+(l_)2-6ybl~Fz—6yb3-72—9yb2-62-0yb4-712)K1}

- —T,21R§4) {9yt2 K4+9ytl-K3+(.71+9ytl E,+6yt3- F,+6yt2 G,

+0yt4-I7,)-K2+(.72-6yt1~E2—9yt3-F2—6yb2-C—;2—9yt4-T12)K,}

— —T,2R§"’{ub4 K4+ub3 K3+(-A,+ub1 E,+ub3 F,+ub2 G,

+ub4.H, )-x2+(Az—ub1 -Ez-ub3-F2-ub2-Gz-ub4-F12)K,}

- —T,21R(k){ut4 - K4+ut3 - x3 +(1‘3, +ut1 .H, +ut3 ~F1+ut2° G,

+kut4-ﬂ,).K2+(§2—utl-Ez-ut3~F2—ut2-Gg-ut4~ﬂz)xl}
-T,k21g‘,R)-{vb4 K4+vb3 K3+(— B,+vb1 E,+vb3 F,+vb2 G,
+vb4-H,)-K2+(1+Bz—vb1~Eg—vb3-F2-vb2-Gz-vb4-Hz)x,}
_leleg:){vt4 K4+vt3 K3+(B,+vt1 E,+vt3 F,+vt2 G,
+vt4-FI,)-K2+(—Bz—vtl~Ez-vt3~F2—vt2-G2-vt4-thcl}
—T, k2 R(34k){6xb4-K4+9xb3 - K3+(C’, +9xb1 ~E, +9xb3oF, +9xb2-C’,
+ 0xb4 . H,)x2 +(G2—exb1~Hz—exb3-F2—exb2-G2—exb4-H2)x,}
_T, Zng’wm K4+6xt3 x3+(i,+extl H,+ex:3 F,+ex:2 G,
+9xt4-F11) - x2+(22—ext1 -E'2-9xt3 -F2—9x12- G2-9x14-H2)1c,}

-1, 2 R(k){9yb4-1c4 + 9yb3 . x3 + (b, + eyb1 .H, + Obe - F, + 6yb2 - G,

 

211

+6yb4-H,) - K2+(D2—9ybl -Ez—0yb3 - FZ—eybz - Gz—Gyb4-Hz)x,}

N
5,6 = —TI 2 R§§){9yt4 K4+9yt3 x3+(], +6yt1 -E, +6yt3 -F, +6yt2 - G,

 

k=l
+0yt4-H,)eK2+(.72-9ytl-E’2—9yt3-F2-6y12-Gz—Byt4-thcl}
1 4 1 3 3 1 2 2
T, = N , K4 = z(b,-21-1). K3 = 3(Zk-zk-l)’ K2:3,(Zk‘zk-1)
2 Kg“ thkl
k=1
(k) (k) (k)
_ R1k)_ C31 R1k)_ C32 Rm 1 R00 C36
K1“(Zk‘zk—1)’ 31 ' ‘71,”), 32 ‘ "W, 33 "m, 34 ‘75
C33 C33 C33 C33
Plathnertias:
N Q
(k) (k) (k) V _ (k) (k) (k)
IaBrp= Z 3“) ¢a13r¢pdz’ IaBrp ' 2 JP ‘1’ aB‘Prp dz
k=l k=l
Q 1 Q-1 .AJZI)
N Q
k
1;, = Z Ip( ’52,!)de a,r = 1...4; B,p = 1,2
k=lzk-—l
Stress Resultants:
(k) (k)dz (k) (k)d
N=xaB ZJGHOIBM N=yaB 2}°yyabdz‘l’
=11: 1kN=l:k-1
(k) (k) (k) (Dd
N=Zaﬂ k2}OMXBdZ’ QzaB=kz}ouWABd
N=lz—Q 1 k=le- 1
_ (k) (1‘)
R2,, - £3] a: kxﬁdz, an5= 2 j r wapzdz A.(22)
:IQ 1k=le-1
_ (k) (k) (k)
R”- 2 It” dez, (2,0,3: 2 j: (Dapzdz
k;l7-:- 1 k=lQ— 1

R4; = 2 ]' 1(k)Qde,nyaB_ 2 ]‘ munch“) dz, Ryxap- 2 3' 1:(Ic),l,(k)dz

k=lQ1k=lQ1k=lQ|

212

. Shape Functions in the Covariant Frame.
The shape functions deﬁned in Eqs. A.(16), A.(17) are expressed in covariant frame

through the transformation matrix Eqs.(3.5.38) and (3.5.44) as follows:

_ (1k) 111 (k) _ (141112 (k)
A.(23)
(1k) 122 (k) _ 121 (k) 111 (k)
111 (k) (k) _ 122 (k) 112 (k)
‘Pm = 3'11? 15 “'23, T5213 - E‘I’asf "7‘11, 4B
A.(24)

J
jl:\P(lk)+ ‘1’“)

‘I’nw =1+p 213, ‘I’nza = -—

 

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[81

[9]

[10]

[11]

[12]

BIBLIOGRAPHY

Kapania, R.K. and Raciti, 8., “Recent Advances in analysis of laminated beams and
plates Part I: Shear effects and buckling”, AIAA Journal, Vol. 27, 1989, pp. 923-934

Noor, AK. and Burton, W.S., “Assessment of Shear Deformable Theories for Multi-
layered Composite Plates”, Appl. Mech. Rev., Vol.42, 1989, pp. 1-12

Noor, A.K. and Burton, W.S., “Assessment of Computational Models for Multilay-
ered Composite Shells”, Appl. Mech. Rev., Vol.43, 1990, pp. 67-97

Noor, A.K., “Mechanics of Anisotropic Plates and Shells - a New Look at an Old
Subject”, Computers and Structures, Vol. 44, 1992, pp. 499-514

Reddy, J.N., “A Review of Reﬁned Theories of Laminated Composite Plate”, Shock
V1br. Dig., Vol. 22, 1990, PP. 3-17

Jemielita, G., “On Kinematical Assumptions of Reﬁned Theories of Plates: A Sur-
vey”, J. Appl. mech., Vol. 57, 1990, pp. 1088-1091

Timoshenko, S., Vibration Problems in Engineering, 2nd ed, D. Van Nostrand Com-
pany, Inc., New York, 1937

Reissner, E., “The Effect of Transverse Shear Deformation on the Bending of Elastic
Plates”, J. Appl. mech., Vol. 12, 1945, pp. A69-A77

Mindlin, R.D., “Inﬂuence of Rotary Inertia and Shear on ﬂexural motions of isotro-
pic, elastic plates”, J. Appl. Mech., Vol. 18, 1951, pp. 31-38

Yang, RC., Norris, CH, and Stavsky, Y., “Elastic Wave Propagation in Heteroge-
neous Plates,” Int. J. Sol. Struct., Vol. 2, 1966, pp. 665-684.

Whitney, J.M., “Shear Correction Factors for Orthotropic Laminates Under Static
Loa ”, J. Appl. Mech., Vol. 1973, pp. 302-304

Vlachoutsis, S., “Shear Correction Factors for Plates and Shells”, Int. J. Num. Meth.

213

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

214

Engng., Vol. 33, 1992 pp.1537-1552

Whitney, J.M. and Sun, C.T., “A Higher Order Theory for Extensional Motion of
Laminated Composites, Journal of Sound and Vibration, Vol. 30, 1973, pp.85-97

Lo, K.H., Christensen, RM. and Wu, E.M., “A High-Order Theory of Plate Defor-
mation - Part 2: Laminated Plates,” J. Appl. Mech., Vol. 44, 1977, 669-676.

Reddy, J.N., “A Simple Higher-Order Theory for Laminated Composite Plates”, J.
Appl. Mech., Vol. 51, 1984, pp. 745-752.

Barbero E]. and Reddy J.N., “An Accurate determination of stresses in thick lami-
nates using a Generalized plate theory”, Int. J. Numer. Meths. Engrg., Vol. 29, 1990,
pp. 1-14

Tessler, A., “A Higher-Order Plate Theory with Ideal Finite Element Suitability”,
Comp. Meth. Appl. Mech. Eng, Vol. 85, 1991, pp. 183-205.

Lo, K.H., Christensen, RM. and Wu, E.M., “Stress Solution Determination for High
Order Plate Theory”, Comput. Struct., Vol. 14, 1978, pp. 655-662

Chaudhuri, R.A., “An Equilibrium Method for Prediction of Transverse Shear
Stresses in a Thick Laminated Plate”, Comput. Struct., 1986, Vol. 23, pp. 139-146

Lajczok, M.R., “New Approach in the Determination of Interlaminar Shear Stresses
for the Results of MSC/Nastran”, Comput. Struct., 1986, Vol. 24, pp. 651-656

Chaudhuri, R.A. and Seide, B, “An Approximate Semianalytical Method for the Pre-
diction of Interlaminar Shear Stresses in an Arbitrarily Laminated Thick Plate”,
Comput. Struct., 1986, Vol. 25, pp. 627-636

Byun, C. and Kapania, R.K., “prediction of interlaminar stresses in a laminated plate
using global orthogonal interpolation polynomials”, AIAA Journal, Vol. 30, 1992,
pp. 2740-2749

Srinavas, S, “A Reﬁned Analysis of Composite Laminates”, J. Sound Vibr. Vol. 30,
1973. PP. 495-507

Owen, DR]. and Li, 2.11., “A Reﬁned Analysis of Laminated Plates by Finite Ele-
ment Displacement Methods - 1. Fundamentals and Static Analysis”, Comput.
Struct., Vol. 26, 1987, pp. 907-914

Reddy, J.N., “Generalization of Two-Dimensional Theories of Laminated Compos-
ite Plates”, Comm. Appl. Num. Meths. Vol. 3, 1987, pp.173-180.

Toledano, A. and Murakami, H., “A Composite Plate Theory for Arbitrary Laminate

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

[40]

215

Conﬁguration”, J. Appl. Mech., Vol. 54, 1987, pp.18l-l89.

Reddy, J.N. and Barbero, E.J., “A Plate Bending Element Based on a Generalized
Laminate Plate Theory”, Int. J. Num. Meth. Engng. Vol. 28, 1989, pp. 2275-2292

Lu, X. and Liu, D., “An Interlaminar Shear Stress Continuity Theory for Both Thin
and Thick Composite Laminates”, J. Appl. Mech., Vol.59, 1992, pp.502-509.

Lee, CY. and Liu, D., “An Interlaminar Stress Continuity Theory for Laminated
composite analysis”, Comput. Struct., Vol. 42, No.1, 1992, pp. 69-78.

Robbins, DH. and Reddy J .N., “Modeling of Thick Composites Using a Layerwise
Laminate Theory”, Int. J. Num. Meth. Engng., Vol. 36, 1993, pp. 655-677

Orth, E]. and Surana, K.S., “p-Version TWO-Dimensional Curved Beam Element
Using a Piecewise Hierarchical Approximation for Laminated Composites, Comput.
Mech., Vol. 11, 1993, pp. 141-161

Lee, J., Gurdal, Z. and OH. Grifﬁn Jr, “Layerwise Approach for the Bifurcation
Problem in Laminated Composites with Delaminations”, AIAA Journal, Vol. 31,
1993,1313. 331-338

DiSciuva, M., “Development of an Anisotropic Multilayered Shear deformable
Rectangular Plate Element”, Comput. Struct., Vol. 21, 1985, pp.789-796.

DiSciuva, M., “An Improved Shear-Deformation Theory for Moderately Thick Mul-
tilayered Anisotropic Shells and Plates,” J. Appl. Mech., Vol. 54, 1987, pp. 589-596.

DiSciuva, M., “A General Quadrilateral Multilayered Plate Element with Continu-
ous Interlaminar Stresses,” Comput. Struct., Vol. 47, 1993, pp. 91-105.

Cho, M. And Parmerter, R.R., “Efﬁcient Higher Order Composite Plate Theory for
General Lamination Conﬁgurations”, AIAA Journal, Vol. 31, 1993, pp. 1299-1306.

Murakami H., “Laminated Composite Theory with Improved In-plane Response”, J.
Appl. Mech., Vol. 53, 1986. PP- 661-666.

Averill, R.C., “Static and Dynamic Response of Moderately Thick Laminated
Beams with Damage”, Composite Engineering, Vol. 4, 1994, pp. 381-395.

Cho, Y.B. and Averill, R.C., “An Improved Theory and Finite Element Model for
Laminated Composite and Sandwich Beams Using First Order Zig-Zag Sublaminate
Approximations”, will be appeared in Composite structures in 1997

Averill, RC. and Yip, Y.C., “Development of Simple, Robust Finite Elements Based
on Reﬁned Theories for Thick Laminated Beams”, Comput. Struct., Vol. 59, l996pp.

 

[41]

[42] ‘

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[50]

[51]

[52]

216

529-546.

Averill RC. and Yip, Y.C., “An Efﬁcient Thick Beam Theory and Finite Element
Model with Zig-Zag Sublaminate Approximations”, AIAA Journal, Vol. 34, 1996,
pp. 1627-1632.

Ylp, Y.C., and Averill, R.C., “A 3D Laminated Plate Finite Element with Zig-Zag
Sublaminate Approximations”, submitted for publication, 1996.

Gallagher, R.H., “Problems and Progress in Thin Shell Finite Element Analysis”, in
Finite Elements for Thin Shells and Curved Members, Eds. Ashell, D.G. and Galla-
pher, R.H., John Wiley & Sons, New York, 1976, pp. 1 -l 4

Huang, H.C., Static and Dynamic Analyses of Plates and Shells: Theory, Software,
and Applications, Springer-Verlag, London, 1988, pp. 1 - 43

Ahmad S., Irons, E.M., and Zienkiewicz, O.C., “Analysis of Thick and Thin Shell
Structures by Curved Finite Elements”, Int. J. Num. Meth. Engng., Vol. 2, 1970, pp.
419—451

Ausserer, M.F. and Lee, S.W., “An Eighteen-Node Solid Element for Thin Shell
Analysis”, Int. J. Num. Meth. Engng., Vol. 26, 1988, pp. 1345-1364

Kim, Y.H. and Lee, S.W., “A Solid Element Formulation for Large Deﬂection Anal-
ysis of Composite Shell Structures”, Comput. Struct., Vol. 30, No. 1/2, 1988, pp. 269
- 274

Wanthal, SP. and Yang, H.T.Y., “Three-Dimensional Finite Element Formulation for
Laminated Plates”, Journal of Reinforced Plastics and Composites, Vol. 10, 1991,
pp. 330 - 355

Barbero, E.J., “3-D Finite Element Modeling of Laminated Composites by Layer-
wise Constant Shear Elements”, NDE-Vol. 10, Enhacing Analysis Techniques for
Composite Materials ASME, 1991, pp. 229 - 237

Han, J and Hoa, S.V., “A Three-Dimensional Multilayer Composite Finite Element
for Stress Analysis of Composite Laminates”, Int. J. Num. Meth. Engng., Vol. 36,
1993, Pp. 3903 - 3914

Zinno, R. and Barbero, E.J., “A Three-Dimensional Layer-wise Constant Shear Ele-
ment for General Anisotropic Shell-Type Structures, Int. J. Num. Meth. Engng., Vol.
37, 1994, pp. 2445 - 2470

Kong, J. and Cheung, Y.K., “Three-Dimensional Finite Element Analysis of Thick
Laminated Plates”, Computers and Structures, Vol. 57, No. 6, 1995, pp. 1051 - 1062

 

[53]

[54]

[55]

[56]

[57]

[58]

[59]

160]

[61]

[62]

[63]

[64]

217

Dorfmann, A. and Nelson, R.B., “Three-Dimensional Finite Element for Analyzing
Thin Plate/Shell Structures”, Int. J. Num. Meth. Engng., Vol. 38, 1995, pp. 3453 -
3482

Bisplinghoff, R.L., Ashley, H., Halfrnan, R.L., Aeroelasticity, Adddison-Wesley
Publishing Co., MA., 1955, pp.125 - 131

Archer, J .S., “Consistent Mass Matrix for Distributed Mass Systems”, Journal of the
Structural Division, Proceedings of the American Society of Civil Engineers, Vol.
89, No. ST4, Aug. 1963, Part 1, pp. 161-178

Archer, J .S., “Consistent Matrix Formulations for Structural Analysis Using Finite
Element Techniques”, AIAA Journal, Vol. 3, No. 10, 1965, pp. 1910-1918

Fried, I. and Malkus, D.S., “Finite Element Mass Matrix Lumping by Numerical
Integration with No Convergence Rate Loss”, Int. J. Solids Structures, Vol.11, 1975,
pp. 461-466

Hinton E., Rock, T. and Zienkiewicz, DC, “A Note on Mass Lumping and Related
Processes in the Finite Element Method”, Earthquake Engng. Struct. Dynamics, Vol.
4, No. 3, 1976, pp. 245-249

Surana, K.S., “Lumped Mass Matrices with Non-zero Inertia for General Shell and
Axisymmetric Shell Elements”, Int. J. Num. Meth. Engng., Vol. 12, No. 11, 1978,
pp. 1635-1650

Bossak, M.A.J., and Zienkiewicz, O.C., “Free vibration of Initially Stressed Solids
with Particular Reference to Centrifugal-Force Effects in Rotating Machinery”,
Journal of Strain Analysis, Vol.8, No. 4, 1973, pp. 245-252

Hinton, E., and Rock, T., “Free Vibration and Transient Response of Thick and Thin
Plates Using the Finite Element Method”, Earthquake Engineering and Structural
Dynamics, 1974, Vol. 3 pp.51-63

Tong, P., Pian, T.H.H. and Bucciarelli, L.L., “Mode Shapes and Frequencies by
Finite Element Method Using consistent and Lumped Masses”, Comput. Struct.,
Vol. 1, 1971, pp. 623-638

Li, Y.N., Liang, R.Y. and Wang, D], “On Convergence Rate of Finite Element
Eigenvalue Analysis with Mass Lumping by Nodal Quadrature”, Computation
Mechanics, Vol. 8, 1991, pp. 249-256

Chan, H.C. Cai, CW. and Cheung, Y.K., “Convergence Studies of Dynamic Analy-
sis by Using the Finite Element Method with Lumped Mass Matrix”, Journal of
Sound and Vibration, Vol. 165, No. 2, 1993, pp. 193-207

 

[65]

[66]

[67]

[63]

[69]

[70]

[71]

[72]

[73]

[74]

[75]

[761

[77]

218

Malkus, D.S. and Plesha, M.E., “Zero and Negative Masses in Finite Element Vibra-
tion and Transient Analysis”, Comp. Meth. Appl. Mech. Engng., Vol. 59, No. 3,
1986, PP. 281-306

Malkus, D.S., Plesha, ME. and Liu, M.R., “Reversed Stability Conditions in Tran-
sient Finite Element Analysis”, Comp. Meth. Appl. Mech. Engng., Vol. 68, No. l,
1988, pp. 97-114

Flanagan, G., “A Sublaminate Analysis Method for Predicting Disbond and Delami-
nation Loads in Composite Structures,” J. Reinforced Plastics and Composites, Vol.
12, 1993, Pp. 876-887.

Lynn, PP. and Arya, S.K., “Finite Elements Formulated by the Weighted Discrete
Least Squares Method”, Int. J. Num. Meth. Engng., Vol. 8, 1974, pp. 71-90

Zienkiewicz, O.C., Owen, D.J.R., and Lee, K.N., “Least Square Finite Element for
Elasto-Static Problems. Use of Reduced Integration”, Int. J. Num. Meth. Engng. Vol.
8, 1974, pp. 341-358

Reddy, J. N., “A Penalty Plate-Bending Element for the Analysis of Laminated
Anisotropic Composite Plates”, Int. J. Num. Meth. Engng. Vol. 15, 1980, pp. 1187-
1206

Tessler, A. and Dong, S.B., “On a Hierarchy of Conforming Timoshenko Beam Ele-
ments”, Comput. Struct, Vol. 14, 1981, pp. 335-344.

Pagano, N .J ., “Exact Solutions for Composite Laminates in Cylindrical Bending”, J.
Campos. Maten. Vol. 3, 1969, pp. 398-411.

Barlow, J. “Optimal Stress Locations in Finite Element Models”, Int. J. Numer:
Meths. Engrg., Vol. 10, 1976, pp. 243-251.

Cho, Y.B and Averill, R. C, “Zig-Zag Sublaminate Finite Element Model for Lami-
nated Composite and Sandwich Structures”, The 12th Annual Advanced Composites
Conference & Exposition (ACCE ‘97)

Tessler, A and Hughes, T.J.R., “An Improved Treatment of Transverse Shear in the
Mindlin-Type Four-Node Quadrilateral Element”, Comp. Meth. Appl. Mech. Engng.,
Vol. 39, 1983. Pp. 311-335.

Prathap, G. and Somashekar, B.R., “Field- and Edge-Consistency Synthesis of a 4-
Noded Quadrilateral Plate Bending Element”, International Journal for Numerical
Methods in Engineering, Vol. 26, 1988, pp. 1693-1708.

Prathap G., “The Displacement-type Finite Element Approach - From Art to Sci-
ence”, Progress in Aerospace Science, Vol. 10, 1994, pp. 295-405

 

[78]

[79]

[80]

[31]

[82]

[83]

[84]

[85]

[86]

[87]

[88]

[89]

[90]

219

Kulak, RF. and Fiala, C., “Neptune: A System of Finite Element Programs for
Three-Dimensional Nonlinear Analysis”, Nuclear Engineering and Design, Vol.
106, 1988, pp. 47-68

Reddy, J.N., Energy and Variational Methods in Applied Mechanics, John Wiley &
Sons, New York, 1984, pp. 51-55

Tessler, A. and Saether, E., “A Computationally Viable Higher-Order Theory for
Laminated Composite Plates”, Int. J. Num. Meth. Eng., Vol. 31, 1991, pp. 1069-1086

Reissner, E., “On a Mixed Variational Theorem and on Shear Deformation Plate
Theory,” Int. J. Num. Meth. Eng., Vol. 23, 1986, pp. 193-198

Ergatoudis, I., Irons, B.M., and Zienkiewicz, O.C., “Curved, Isoparametric, “Quadri-
lateral” Elements for Finite Element Analysis”, Int. J. Solids Structures, Vol. 4,
1968, pp. 31-42

Zienkiewicz, O.C., Taylor, R.L., and T00, J .M., “Reduced Integration Technique in
General Analysis of Plates and Shells”, Int. J. Num. Meth. Eng., Vol. 3, 1971, pp.
275-290

Hughes, T.J.R., Cohen, M., and Haroun, M., “Reduced and Selective Integration
Techniques in the Finite Element Analysis of Plates”, Nuclear Engineering and
Design, Vol. 46, 1978, pp. 203-222

Pugh, E.D.L., Hinton, E., Zienkiewicz, O.C., “A Study of Quadrilateral Plate Bend-
ing Elements with ‘Reduced’ Integration”, Int. J. Num. Meth. Eng., Vol. 12, 1978,
pp. 1059-1978

Stolarski, H. and Belytschko, T., “Membrane Locking and Reduced Integration for
Curved Elements”, J. Appl. Mech., Vol. 49, 1982, pp. 172-176

Hughes, T.J.R., Taylor, R.L., and Kanoknukulchai, W., “A Simple and Efﬁcient
Finite Element for Plate Bending”, Int. J. Num. Meth. Eng., Vol. 11, 1977, pp. 1529-
1543 '

Parisch, H., “A Critical Survey of the 9-Node Degenerated Shell Element with Spe-
cial Emphasis on Thin Shell Application and Reduced Integration”, Comp. Meth.
Appl. Mech. Engng., Vol. 20, 1979, pp-323-350

Belytschko, T. and Ong, J S] ., “A Consistent Control of Spurious Singular Modes in
the 9-Node Lagrange Element for the Laplace and Mindlin Plate Equation”, Comp.
Meth. Appl. Mech. Engng. Vol. 44, 1984, pp. 269-295

Belytschko, T., Liu, W.K., Ong, J S] ., and Lam, D., “Implementation and Applica-
tion of A 9-node Lagrange Shell Element with Spurious Mode Control”, Comput.

[91]

[92]

[93]

[94]

[95]

[96]

197]

[98]

[99]

220

Struct, Vol. 20, No. 1-3, 1985, pp. 121-128

Crisﬁeld, M.A., “A Four-Noded Thin-Plate Bending Element Using Shear Con-
straints-a Modiﬁed Version of Lyons’ Element”, Comp. Meth. Appl. Mech. Engng.,
Vol. 38, 1983, 93-120

Bathe, K.J. and Dvorkin, E.N., “A Formulation of General Shell Elements - the Use
of mixed Interpolation of Tensorial components”, Int. J. Num. Meth. Eng., Vol. 22,
1986, pp. 697-722

Hinton, E. and Huang, H.C., “A Family of Quadrilateral Mindlin Plate Elements
with Substitute Shear Strain ﬁelds”, Comput. Struct. Vol. 23, pp.409-431, 1986

Huang, H.C., Static and Dynamic Analyses of Plates and Shells: Theory, Software,
and Applications, Springer-Verlag, London, 1988, pp. 46 - 50

Robbins, D.H., Reddy, J .N., and Krishna-Murty, A.V., “On the Modeling of Delami-
nation in Thick Composites”, in Enhancing Analysis Techniques For Composite
Material, NDE-Vol. 10, Winter ASME Conference 1991, Eds. Schwer L. Reddy J .N .
and Mal A., pp. 133-149.

Prathap, G., “The Curved Beam/Deep Arch/Finite Ring Element Revisited”, Int. J.
Num. Meth. Eng., Vol. 21, 1985, pp. 389-407.

Prathap, G., “The poor bending response of the four-node plane stress quadrilateral,
Int. J. Num. Meth. Eng., Vol. 21, pp. 825-835.

Burton, W.S. and Noor, A.K., “Three-Dimensional Solutions for Thermomechanical
Stresses in Sandwich panels and Shells”, Journal of Engineering Mechanics, ASCE,
Vol. 120, No. 10, 1994, pp. 2044-2071.

Cook, R.D., Malkus, D.S. and Plesha, M.E., Concepts and Applications of Finite
Element Analysis, 3rd Ed., John Wiley & Sons, New York, 1989, pp. 367-418

[100]Bathe, K.J., Finite Element Procedures, Prentice Hall, New Jersey, 1996, pp. 768-

837

[101]Belytschko,'T. and Hughes, T.J.R., eds., Computational Methods for Transient Anal-

ysis, North Holland, Amsterdam, 1983, pp. 2-65.

[102]Kulak, RF, and Fiala, C., “NEPTUNE: A System of Finite Element Programs for

Three-Dimensional Nonlinear Analysis,” Nuclear Engineering and Design, Vol. 106,
1988, PP. 47-68.

[103] .Kulak, R. F., Plaskacz, E. J. and Pfeiffer, P. A., “Structural Mechanics Computations

on Parallel Computing Platforms,” Fluid-Structure Interaction and Structural

221

Mechanics 1995, C.Y. Wang, ed., ASME Publication PVP-Vol. 310, 1995, pp.129-
133

[104].Kulak, RF. and Pfeiffer, P.A., “A Computational Engine for the Nonlinear, Tran-
sient Response of Preloaded Structures”, Development, Validation, and Application
of Inelastic Methods for Structural Analysis and Design, Sammataro, RF. and
Ammerman, D.J., eds., ASME Publication PVP-Vol. 343, 1996, pp.231-238

[105]Belytschko, T., Lin, 1.1. and Tsay, C.S., “Explicit Algorithms for the Nonlinear
Dynamics of Shells”, Comp. Meth. Appl. Mech. Engng., Vol. 42, 1984, pp. 225-251

[106]Belytschko, T., Ong, J.S.J., Liu, W.K. and Kennedy, J.M., “Hour Glass Control in
Linear and Nonlinear Problems”, Comp. Meth. Appl. Mech. Engng., Vol. 43, No. 3,

[107]Flanagan, DD. and Belytschko, T. “A Uniform Strain Hexahedron and Quadrilateral
with Orthogonal Hour Glass Control”, Int. J. Num. Meth. Enggng., Vol. 17, No. 5,
1981, pp. 679-706

[108]Hsiao, KM. and Hung, H.C., “Large-Deﬂection Analysis of Shell Structure by
Using Corotational Total Lagrangian Formulation”, Comp. Meth. Appl. Mech.
Engng., Vol. 73, 1989, pp. 209-225

[109]Hsiao, K.M., “Corotational Total Lagrangian Formulation for Three-Dimensional
Beam Element”, AIAA Journal, Vol. 30, No. 3, 1992, pp. 797-804

[110]Jiang, L. and Chemuka, M.W., “A Simple Four-Noded Corotational Shell Element
for Arbitrarily Large Rotations”, Comput. Struct, Vol. 53, No. 5, 1994, pp. 1123-
1132

[111]Key, S.W., “Transient Response by Time Integration: Review of Implicit and
Explicit Operators,” Advanced structural dynamics, J. Donea, ed., Applied Science
Publishers LTD, 1979

[1 12] Flanagan, DP. and Belytschko, T. “Eigenvalues and Stable Time Steps for the Uni-
form Strain Hexaheron and Quadrilater ”, Journal of Applied Mechanics, Vol. 51,
1984, pp.35-40

[113]Isaacson, E. and Keller, H.B., Analysis of Numerical Methods, John Wiley and Sons,
New York, 1966

[114]Kulak, R. F., “Critical Time Step Estimation for Three-Dimensional Explicit Impact
Analysis,” Structures under shock and impact: proceedings of the ﬁrst international

conference, Cambridge, Massachusetts, USA, July 1989, RS. Bulson, ed., Computa-
tional Mechanics Publications and Elsevier, 1989, pp. 155-63.

222

[1 15] Kant, T, Ravichandran, R.V., Pandya, B.N., and Mallikarjuna, “Finite Element Tran-
sient Dynamic Analysis of Isotropic and Fibre Reinforced Composite Plates Using a
Higher Order Theory”, Composite Structures, Vol. 9, 1988, pp. 319-342

[116]Swaddiwudhipong, S. and Liu, Z.S., “Response of Laminated Composite plates and
shells”, Composites Structures, Vol. 37 , No. 1, 1997, pp. 21-32

[117]Reddy, J .N., “Geometrically Nonlinear Transient Analysis of Laminated composite
Plates”, AIAA Journal, Vol. 21, No.4, 1983, pp. 621-629

"7'1111111111'111111