.I . 4s..T.I
. ..1..." .y;

4..
«a...

1....
a. flaw:

e . e an”,
weee.
«wear.
a

.. e we
. w e, .
vdlir .Jnufirhflh

. v . A
3%.

i
s

e :5
.Enu
. era; .
, . . «u. a... ..
. . : fly. a???
:. Hun—haul #. a»
_ . awe .e .
4 . . , H. a}
e :L , L. , 5'
3.53.35. . w . . .
gnaw: .
«a: . 3%
e a?! .we, s
.. .1 ex.

th-quhu 3
«in «mm?

. a .e A;
r s 4”“. .
. ”(I

it
. e.”

3”

 

, Been; A; ,
u V1.1.
A!) w u .
e .r'i‘ - e
. . e
rkll
. .
133. Pat» ..
7.3V. .n. .954.
b. .
.af
.-
e 5.0.“ 4

. . 3.1L
.. .. v2 ; .
.

, E3325»?

51 ”mafia...

 

THESIS

 

 

 

 

2.
:3 C r: !
LIBRARY
Mlchigan State
University I
This is to certify that the
dissertation entitled

COUPLED MULTIOBJECTIVE OPTIMIZATION OF
CRASHWORTHINESS AND MODAL FREQUENCIES IN STRUCTURES
USING GENETIC ALGORITHMS

presented by

David J. Eby

has been accepted towards fulfillment
of the requirements for

Ph.D. degree in Mechanics

Mam

Major professor

 

Date mm

MS U is an Affirmative Action/Equal Opportunity Institution 0-12771

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINE return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

“0v 0 4 2003

 

villi? ‘

 

 

 

 

 

 

 

 

 

 

 

 

 

11/00 chIRCJDateDuopfiS-p.“

COUPLED MULTIOBJECTIVE OPTIMIZATION OF CRASHWORTHINESS AND
MODAL FREQUENCIES IN STRUCTURES USING GENETIC ALGORITHMS

By

David J. Eby

AN ABSTRACT OF A DISSERTATION

Submitted to
Michigan State University
in partial fiilfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Materials Science and Mechanics Department

2000

Professor Ron C. Averill

ABSTRACT

COUPLED MULTIOBJECTIVE OPTIMIZATION OF CRASHWORTHINESS AND
MODAL FREQUENCIES IN STRUCTURES USING GENETIC ALGORITHMS

By

David J. Eby

This dissertation presents an approach to concurrent optimization of
stnictures by using an island injection Genetic Algorithm for: 1) crash energy
management, 2) modal vibration frequency response, 3) peak crush force, 4)
weight reduction, and 5) manufacturability, while considering stochastic variability
of design variables and loading conditions. Existing automated design technology
has never simultaneously addressed all of these competing mechanisms. In
addition, this dissertation presents the development of a geometrically nonlinear
first-order zig-zag sublaminate theory and finite element model which accounts
for moderately large displacements and moderate rotations using a Total
Lagrangian formulation. The accuracy of the model is demonstrated by
comparing its structural response predictions with results from previous
experimental investigations and with numerical tests. Finally, the optimal design
of sandwich panels was approached with a simple Genetic Algorithm and the
developed geometrically nonlinear first-order zig-zag sublaminate finite element

model.

Copyright by
David J. Eby
2000

DEDICATION

I would like to dedicate this body of work to all of my family and friends.
This includes my mother, father, three brothers, two sisters and many friends. In

particular, Karen Kneibes for her love, patience, and understanding.

iv

ACKNOWLEDGMENTS

I would like to express my gratitude to Dr. Ron Averill, Dr. Erik Goodman,
Dr. William Punch Ill, Dr. Alejandro Diaz, and Dr. Baisheng Yan for serving on
my defense committee. Also, John Madakacherry of General Motors provided
insurmountable insight and leadership during the duration of this work. I would
also like to acknowledge the support of the Computational Mechanics Research
Group and the Genetic Algorithms Research and Application Group at Michigan
State University for all of their patience and advise. I would like to thank my
advisor, Professor Ron Averill, for helping me develop the tools needed to
overcome the many challenges of life and research - without his guidance,

understanding, and foresight none of this work would have been completed.

TABLE OF CONTENTS

List of Tables ................................................................................. vii
List of Figures ................................................................................. viii

Chapter 1: Introduction

1.1 Introduction .................................................................... 1
1.2 Literature Review ............................................................. 3
1.3 Present Study .................................................................. 14
Chapter 2: The Mechanics of Genetic Algorithms
2.1 Introduction ..................................................................... 16
2.2 Simple Genetic Algorithms .................................................. 16
2.3 Pros and Cons of Genetic Algorithms .................................... 18
2.4 Island Injection Genetic Algorithms ....................................... 25

Chapter 3: Optimization of Crashworthiness and Modal Frequencies
of Automotive Structures

3.1 Introduction ..................................................................... 28
3.2 Optimization Problem Statement .......................................... 28
3.3 Optimization Objective Functions ......................................... 32
3.4 Optimization Constraints .................................................... 33
3.5 Application of Island Injection Genetic Algorithm ...................... 36
3.6 Results ........................................................................... 38
3.7 Conclusions ..................................................................... 41
Chapter 4: Development of a Geometrically Nonlinear First Order Zig-Zag
Finite Element Model
4.1 Introduction ..................................................................... 65
4.2 Basic Ingredients Needed For The First-Order
Zig-Zag Plate Theory ......................................................... 68
4.3 Virtual Work ..................................................................... 73
4.4 Finite Element Model ......................................................... 75
4.5 Newton-Raphson Method ................................................... 78
4.6 Finite Element Approximations ............................................ 80
4.7 Numerical Results ............................................................. 82
4.8 Conclusions ..................................................................... 89
Chapter 5: Optimization of Sandwich Panels using Genetic Algorithms
5.1 Introduction ...................................................................... 109
5.2 Optimization Problem Statement .......................................... 109
5.3 Results ............................................................................ 112
5.4 Conclusions ..................................................................... 1 13
Chapter 6: Summary and Conclusions ................................................ 121
Appendix A .................................................................................... 122
References .................................................................................... 134

vi

I;

LIST OF TABLES

Table 3.1. Target values of inequality constraints for fitness function ......... 62
Table 3.2. Weight values used for fitness function ....................................... 62
Table 3.3. Genetic operator rate for each island in the island injection
Genetic Algorithm. .................................................................. 62
Table 3.4. Primitive and thickness design variable information .................... 63
Table 3.5. Crush evaluation times per island and load case ........................ 63
Table 3.6. Number of full and partial design evaluations .............................. 63
Table 3.7. Thickness design variables for each best design ......................... 63

Table 3.8. Objective and constraint data for the best design found in three
different independent runs ....................................................... 64

Table 3.9. Stochastic variability of internal energy for the best designs found
(20 independent evaluations) with respect to thickness .......... 64

Table 4.1. Material properties for square single-layer orthotropic plate ....... 107

Table 4.2. Material properties for square cross-ply laminated plate

(0/90/90/0) ............................................................................... 107
Table 4.3. Material properties for sandwich panel ........................................ 107
Table 4.4. Lamination scheme for sandwich panel ....................................... 107

Table 4.5. The effect of transverse refinement considering small

displacements ......................................................................... 108
Table 5.1. Ply material properties ................................................................. 119
Table 5.2. Core material properties .............................................................. 120

Table 5.3. Baseline and optimized sandwich panel design variables.
All material names are in reference to Tables 5.1 and 5.2. ....120

Table 5.4. Baseline and GA optimized sandwich panel objectives and
constraints ............................................................................... 120

vii

LIST OF FIGURES

Figure 2.1. Simple Genetic Algorithm Structure ............................................ 27

Figure 2.2. A simple iiGA topology that simultaneously searches various
levels of resolution of the design space on three separate
computational nodes (islands) ................................................. 27

Figure 3.1 Biasing of cross sectional primitives along the length of a rail
structure ................................................................................... 43

Figure 3.2 Surface of rail structure is developed with a linear blending of
the cross sectional primitives. The surface was smoothed at
each of the cross sectional primitive locations. ...................... 44

Figure 3.3. Potential thickness changes occurred at the placement of
each cross sectional primitive .................................................. 45

Figure 3.4. Cross sectional primitive library that contains a set of
knowledge-based manufacturable cross sections ................... 46

Figure 3.5. Finite element model and of rail structure running into a
stationary rigid wall for the first load case. A lumped mass
was placed directly behind the rail structure ........................... 47

Figure 3.6. Finite element model and of rail structure running into a
stationary rigid wall for the second load case. An eccentric
lumped mass was placed behind the rail sthcture ................. 48

Figure 3.7. An iiGA topology that breaks down the design space into
three levels of design variable and evaluation time based
resolutions with multiple load cases ......................................... 49

Figure 3.8. An iiGA topology that monitors islands the performance of
the best designs for the two objectives defined by two
structural load cases ............................................................... 50

Figure 3.9. Crush zone internal energy as a function of generation
number for load case 1, island 9 ............................................. 51

Figure 3.10. Normalized constraints as a function of generation number
for island 9. A feasible design has normalized mass, force,
and manufacturability constraints which are less than or
equal to unity. A feasible design also has normalized first
and second modal frequency constraints which are greater

viii

 

 

Fl.

than or equal to unity ............................................................... 52

Figure 3.11. Crush zone internal energy as a function of generation
number for load case 2, island 18 ........................................... 53

Figure 3.12. Normalized constraints as a function of generation number
for island 18. A feasible design has normalized mass, force,
and manufacturability constraints which are less than or
equal to unity. A feasible design also has normalized first
and second modal frequency constraints which are greater
than or equal to unity ............................................................... 54

Figure 3.13. Crush zone internal energy as a function of evaluation
number for load cases 1 and 2 from monitoring island ............ 55

Figure 3.14. Normalized constraints as a function of evaluation
number for monitor islands. A feasible design has
normalized mass, force, and manufacturability constraints
which are less than or equal to unity. A feasible design
also has normalized first and second modal frequency
constraints which are greater than or equal to unity ................ 56.

Figure 3.15. Animation of deformation for the best design found

(load case 1) ............................................................................ 57
Figure 3.16. Animation of deformation for the best design found

(load case 2) ............................................................................ 58
Figure 3.17. Surface models for best design found, run 1 ........................... 59
Figure 3.18. Surface models for best design found, run 2 ............................ 60
Figure 3.19. Surface models for best design found, run 3. ......................... 61

Figure 4.1. Schematic of ply stacking and coordinate system used in
first-order sublaminate theory ................................................. 90

Figure 4.2. Element topology for the eight-noded brick and the

six-noded wedge finite elements based on first order

zig-zag sublaminate theory ...................................................... 90
Figure 4.3. Load versus deflection for a square simply-supported

single-layer orthotropic plate under uniform pressure at

various levels of in-plane and through-thickness refinement. 91

Figure 4.4 Load-deflection curves for the clamped square cross ply

ix

 

Fl

Fig

Fig

Flgu

Flglll

Flgui

“gun

plate (0/90/90/0) under uniform pressure at various levels of
in-plane and through-thickness refinement ............................. 92

Figure 4.5. Linear and nonlinear load versus deflection at the top of a
simply-supported square sandwich panel under uniform
pressure using first-order zigzag sublaminate theory .............. 93

Figure 4.6. Linear and nonlinear load versus deflection at the bottom
of a simply-supported square sandwich panel under uniform
pressure using first-order zigzag sublaminate theory .............. 94

Figure 4.7. Through-thickness variation of total X1 displacement at the
plate midspan and X1 =0 .......................................................... 95

Figure 4.8. Through-thickness variation of total X2-displacement at the
plate midspan and X2=0 ......................................................... 96

Figure 4.9. Through-thickness variation of the normal Cauchy stress
(in the X1-direction) at the Gauss points nearest to the
center of the plate .................................................................... 97

Figure 4.10. Through-thickness variation of the normal Cauchy stress
(in the X2-direction) at the Gauss points nearest to the
center of the plate .................................................................... 98

Figure 4.11. Through-thickness variation of the normal Cauchy stress
(in the X1 -direction) for the top five layers of the sandwich
panel at the Gauss points nearest to the center of the plate...99

Figure 4.12. Through-thickness variation of the normal Cauchy stress
(in the X2-direction) for the top five layers of the sandwich
panel at the Gauss points nearest to the center of the plate...100

Figure 4.13. Through-thickness variation of the normal Cauchy stress
(in the X3-direction) at the Gauss points nearest to the
center of the plate .................................................................... 101

Figure 4.14. Through-thickness variation of transverse normal strain (833)
at the Gauss points nearest to the center of the plate ............ 102

Figure 4.15. Through-thickness variation of the transverse Cauchy stress
(023) at the Gauss points nearest to the center of the plate ..... 103

Figure 4.16. Through-thickness variation of the transverse Cauchy stress
(013) at the Gauss points nearest to the center of the plate ..... 104

Figure 4.17. Through-thickness variation of the in-plane Cauchy stress
(0'12) at the Gauss points nearest to the center of the plate... 105

Figure 4.18. Measurement points for the displacements, strains, and
stresses for the sandwich panel .............................................. 106

Figure 5.1. Finite element mesh including applied boundary conditions.
Location of maximum transverse displacement and failure
index. The sandwich panel is 24 inches wide and 84 inches
long ......................................................................................... 115

Figure 5.2. The mass of the best design found as a function of
generation number .................................................................. 116

Figure 5.3. Normalized constraints of the best design found as a function
of generation number .............................................................. 117

Figure 5.4. Linear and geometrically nonlinear load deflection response

curves for the baseline and simple GA optimized sandwich
panels ...................................................................................... 118

xi

va
or;
slr
ba

Ull

OP
to
apl

exl

req

Spa
con
sea

I93(

Chapter 1

INTRODUCTION

1.1 Introduction

Current automated design technology generally requires computation of
gradients of the objective functions and constraints with respect to the design
variables in order to identify a better design. Due to the noisy, chaotic nature of
crash problems, gradients of objectives and constraints with respect to typical
stnictural design variables are difficult, if not impossible, to attain. A gradient-
based optimization approach assumes that the objectives and constraints are
unimodal with respect to the design variables, and will seek extremal portions of
the design space where the gradients are equal to zero. A gradient-based
optimization approach experiences difficulties in multimodal design spaces due
to the fact that many local extrema may exist, making it difficult for such an
approach to escape a local minimum. Typically, many real-life design spaces are
extremely multimodal.

A Genetic Algorithm (GA) is a global search technique that does not
require the computation of gradients for the minimization of objectives or
constraint violations, with the ability to consider stochastic behavior of the design
space. An Injection Island Genetic Algorithm (iiGA) can help increase the
computational efficiency of a simple GA with a divide-and-conquer approach by
searching a discretized design space simultaneously at various levels of

resolution in a heuristic hierarchical fashion (Goodman, 1996).

 

opth
prol
Lair
she:
This

mall

can
slatl
male
there
The
becc

sequ

Fiber reinforced laminated composites have a natural ability to be
optimized. Additional design variables are introduced into the optimization
problem statement when laminated fiber reinforced composites are considered.
Laminated composites are orthotropic in the plane with low ratio of transverse
shear modulus to in-plane modulus and are laminated through-the-thickness.
This allows for a rich design space when compared to standard monolithic
materials. Fiber reinforced laminated composites also have high specific strength
and stiffness when compared to conventional monolithic materials. This increase
in freedom and flexibility is not free; the cost is complex structural behavior at
local and global levels. Fiber reinforced laminated composite lamination schemes
can be specifically constructed to optimally perform under complex loads (large
static loads and large dynamic impact loads). Fiber reinforced composite
material behavior is unique and complex at both the local and global level and
therefore difficult to predict, especially when considering geometric nonlinearities.
The task of manually designing an optimal laminated composite structure
becomes staggeringly hindersome when a large number of discrete lamination

sequences is coupled with complex geometric nonlinear structural behavior.

 

 

3F

slr
sin
brc
line
slrL

don

SIOChg

are sir

1.2 Literature Review

The purpose of this literature review will be to give an overview of the
current level of automated design and composite technology while setting the
stage to define the uniqueness and overall contribution of this research. A huge
amount of research has been performed in the area of structural optimization.
With such a broad area of research, it is impossible to review all aspects of
structural optimization. Since the focus of this dissertation is the optimization of
structures considering nonlinear behavior, this literature review will not include a
broad view of optimization problems where the structure is considered to behave
linearly. Very few researchers have been successful in the area of optimization of
structures while considering nonlinear behavior. The literature review is broken
down into five major parts. The first part will contrast various optimization
approaches. In the second part of the literature review, relevant papers
pertaining to classical optimization will be discussed. The third part of the
literature review will present papers regarding the optimization of robust
structures. The fourth part of the literature review will contain recent applications
of GAs in appropriate areas of structural optimization. Finally, the fifth part of the
literature review will refer to papers dealing with the development of fiber

reinforced laminated composite finite elements.

Optimization approaches include hill climbing, stochastic search, directed
stochastic search and hybrid methods. Hill climbing or gradient-based methods

are single-point search methods that have been applied successfully to many

 

re;
ev;
thu
red
var.

ela

was
The
the r

The

mini:

shape optimization problems (see Soto and Diaz, 1993, and Suzuki and Kikuchi,
1991). Gradient-based methods are less applicable to problems in which the
design space is nonconvex and multimodal (see Sangen, 1990). Random
search methods simply evaluate randomly sampled designs in the search space,
and are therefore generally limited to problems that have small search spaces, if
practical search times are required. A directed random search method, such as
a Genetic Algorithm (GA), is a multiple-point, directed stochastic search method
that can be an effective optimization approach to a broad class of problems
(Goldberg, 1988). The use of GA’s for optimal design requires that a large
number of possible designs be analyzed, even though this number generally still
represents only a miniscule fraction of the total design space. When each
evaluation is computationally intensive, a traditional simple or parallel GA can
thus be difficult to apply. Injection Island Genetic Algorithms (iiGA’s) can help
reduce the computational intensity associated with typical GA’s by searching at
various levels of resolution within the search space (see Eby et al., 1999, Mallot

et al., 1996).

Application of gradient-based topological optimization for crashworthiness
was attempted with success using homogenization (see Mayer, et al., 1996).
The objective of the optimization problem was to maximize the internal energy of
the rear rail of an automobile with a single inequality constraint on total mass.
The objective function was a weighted sum of the internal energy at specified

points in time. The crash analysis was performed with LS—DYNA, an explicit finite

elr

gr:
no.
bOI

Pic

pre
Spe

Spa

element code. A homogenization method was used to compute the elastic and
plastic moduli as a function of the density of the material. The technique begins
searching from a “baseline” design, where gradients of the objective and
constraint are computed with respect to the density of the material for each finite
element. These gradients directed a search algorithm which determined a set of
design variables (material densities) to find a better design in terms of the
objective and constraint. The final design had a 47% increase in internal energy
when compared to the original “baseline” design. The paper makes reference to
a stiffer structure, without specifying the actual peak crush force. No constraints
were placed on: peak crush force, modal response, packaging or
manufacturability. Stochastic variation of gage thickness was not considered.
The paper did show that an optimality condition exists, requiring that each finite
element must have a constant strain energy.

The construction of response surfaces can provide information needed for
gradient-based optimization. Roux et. al (1996) constructed response surfaces of
nonlinear structural analysis considering material, geometric and contact
boundary conditions for a semi-solid tire. The objective of the optimization
problem was to minimize the mass of the design while minimizing strain at a
point. Constraints were placed on the deflection of a point on the design under a
prescribed load. The design variables included the thickness of the tire at
specified regions and the number of “webs” in the design. "Webs” determine the
spacing of material inside the tire which contains air chambers. The authors

began the optimization procedure with a crude finite element model to represent

ot

V3

391

inc

reSp

to ur
been

SUnd

the initial design. Final design iterations in the optimization procedure were
performed using a refined finite element model. This optimization process
represents an attempt to perform a reasonable approximation of the behavior of
the design early in the optimization process while refining the finite element

model later in the optimization run to capture more complex structural behaviors.

A design can be considered robust when small changes in the design
variables and/or loadings produce little or no effect on the behavior of the
objective function and constraints. It is often desirable to desensitize the
objective function and constraints with respect to small changes in the design
variables. Kushel and Rackwitz (1998) formulated a classical gradient-based
approach to reliability-based structural optimization. This was accomplished by
incorporating a smooth differentiable failure probability function into the classical
gradient-based optimization problem statement as a set of constraints. Similariy,
Vietor (1997) approached a stochastic distribution of design variables with a
classical optimization approach as a set of constraints. Gausser and Schuller
(1997) considered optimization of a geometrically nonlinear truss structure with
statistical uncertainties for loading and structural properties through the use of
response surfaces.

The idea of robust optimization through desensitizing engineering systems
to uncertainties, such as manufacturing errors and operational variances has
been the subject of much research (see Belegundu, 1992, Khot, 1991, and

Sundaresan, 1992). Traditionally, the sensitivities of optimal designs are

Opl

Qu

195
Has
Mar

199

GAS
crasl
neuh
the s-
nOnlir
minirr

T080

recorded after the products are manufactured. G. Taguchi (1985) developed a
method to measure sensitivity of a system during the initial design process. The
method was first developed at a time when little computer simulation was
available, so most data was gathered from actual experiments. Taguchi used
orthogonal arrays for approximating the expected design sensitivity to be efficient

with respect to the number of evaluations, due to the lack of compute power.

Many researchers are currently investigating application of GAs for the
optimal design of structures using various numerical methods (Genta, 1995 and
Queipo, 1994). Most papers apply GAs to optimization problems that assume the
structure behaves in a linear fashion (see Chapman, 1996; Eby, 1999; Fabbri,
1997; Flynn, 1995; Foster, 1997; Furura, 1993; Furura, 1995; Genta, 1995;
Haslinger, 1996; Keane, 1995; Kosigo, 1993; Le Riche, 1993; Mallot, 1996;
Mares, 1996; Parrnee, 1997; Punch, 1995; Punch, 1994; Rajah, 1995; Sangren,
1990 and Wolfersdorf, 1997)

Optimization of structures considering nonlinear structural response using
GAs has been applied successfully by Hajela and Lee (1997). Hajela considered
crashworthiness of rotorcraft subfloor structures with a GA in conjunction with a
neural network. A crash model was developed to capture the overall behavior of
the structure through the use of lumped masses, nonlinear beam elements and
nonlinear spring elements. The objective of the optimization problem was to
minimize the accelerations experienced by the occupants for multiple load cases.

The objective was minimized by determining the optimal load-deflection behavior

ar
rel

de

we:
loce
ohh
ohk

SUIT;

SGqL
seho
dent
00m;

of energy absorbing elements while simultaneously finding their optimal location
in the structure. The neural network constructed a response surface of the
objective function (acceleration experienced by the occupant) based on a set of
training data points (in terms of design variables). The surface from the neural
network seeks to approximate the behavior of the structure in a least squares
sense from the data points, and can be viewed as a form of nonlinear regression
analysis. It was shown, for this problem, that the neural network adequately
represented the actual response of the objective function with respect to the
design variables. The load-deflection curves were defined as piecewise linear
functions. An energy balance of the initial kinetic energy of the system and the
internal energy of the structure at rest was performed to produce realistic bounds
on peak force of the load-deflection curves. The second set of design variables
was the topological placement of the energy absorbing components at discrete
locations. The acceleration levels for two load cases were weighted in the
objective function. Multiple independent runs with different weights placed on the
objectives were performed to obtain a Pareto front. A Pareto front is a tradeoff
surface of the objectives (acceleration level for two load cases).

Furuya et al. (1993 and 1995) used a simple GA to optimize the stacking
sequence of a laminated composite plate for buckling load maximization. A
serious problem with GAs occurs when a single analysis is computationally
demanding. In light of this, the objective of the work was to reduce the
computational intensity associated with GAs. One step made in this direction

was the use of a binary tree. A binary tree was used to store all previous

ana
ana
redl
eva
the

pre\

Buc
ope
by a
lhes

lecf

Sim:
coar
find
amo
3W
was
resol
islan.
lSlam

first

analyses performed in the optimization to be retrieved later if an exact duplicate
analysis was reiterated, avoiding costly repeat analyses. A binary tree can
reduce computational costs when a high computational price is paid for each
evaluation (compared to the time used to search the binary tree). A drawback of
the binary tree is the need for a large amount of available memory to store
previous design information. A second step made in the direction of reducing the
computational intensity of GAs was the approximation of certain buckling loads.
Buckling loads were computed exactly for each design string created by genetic
operators. The set of all possible combinations of design strings is approximated
by a linear least squares fit based on bending lamination parameters. The best of
these combinations replaced the nominal design in the population. This
technique searches a small neighborhood of designs for a local optimum.

Punch et al. (1994 and 1995) also used GAs to find optimal composite
structures. They searched for optimal laminated composite beams with a
coarse-grained island injection Genetic Algorithm (iiGA). The objective was to
find stacking sequences for a laminated composites beam that maximized the
amount of mechanical energy absorbed by the beam before fracture. The goal of
approaching the problem with an iiGA was to reduce computational effort. This
was accomplished with an iiGA by searching at various levels of spatial
resolution on separate computational nodes. Structured migration amongst
islands allowed good individuals at low spatial resolution to be injected into
islands of higher spatial resolution. The fitness of the beam was evaluated with

first order zig-zag finite element model developed for composite laminate

 

analysis. The new finite element model accounted for Iayerwise variations of
displacements and stresses by assuming a piecewise linear through-the-
thickness in-plane displacement distribution. A clamped-clamped graphite-epoxy
composite layered beam had a point lead applied at its midspan. The GA was
allowed to place either a 0 or 90 degree ply in each element in the discretized
beam. First, a symmetric beam with a small search space was explored with the
island injection GA. The final results from the island injection GA were confirmed
to be the global optimal solution by enumerating (exhausting) the search space.
Next, larger search spaces that did not assume beam symmetry were explored
with the island injection GA with promising results.

Mallot, et al. (1996) used an island injection GA to optimize an idealized
airfoil. The objective of the problem was to find composite stacking sequences of
sandwich plates that maintain an appropriate opposite twist to compensate for in-
flight aerodynamic loads that cause adverse airfoil twisting, while minimizing
weight under stiffness and ply clustering constraints. Displacements and
stresses were predicted with a finite element model based on a cubic zig-zag
plate theory. The iiGA approach was able to find designs which reduce weight by
66%, increase minimum in-plane stiffness values by 33% and increased the twist
by a magnitude of 77 times the baseline design, but in the opposite direction. In
comparison with a variety of search methods tested, the iiGA approach was the

most efficient in terms of computational time.

10

f2

dis
Re
de'
Mir
the
dell

thro

Out-of-plane loading of a laminated composite can cause appreciable
transverse stresses with respect to other components of stress. Central and edge
delaminations are the primary damage mode in laminated composites. Often,
central delaminations are caused by impact loading while an edge delamination
can often be attributed to the free edge effect (see Liu, 1988). The correct
prediction of transverse stresses is crucial since the interiaminar stresses play a
major role in the energy absorption capability of the laminated composite during
failure.

The first plate theory developed is referred to classical laminate plate
theory and is based on the Kirchoff-Love hypothesis that the in-plane
displacements vary in a linear fashion through-the-thickness while the transverse
displacement remains constant and the transverse shear strain is zero (see
Reddy, 1993). First order shear deformation theory (FSDT) was independently
developed by Reissner and Mindlin, and is typically referred to as Reissner-
Mindlin theory. The theory relaxed the Kirkoff-Love hypothesis which enforced
that the normals to the mid-plane of the plate must remain normal throughout
defamation (see Reddy, 1993). This deformation mode implies a constant
through-thickness transverse shear stress. Shear correction factors are needed
to satisfy equilibrium conditions (see Reddy, 1993). FSDT theory accurately
predicts displacements and stresses of thin and moderately thick isotropic
structures. FSDT theory does not satisfactorily predict displacements and
stresses in some laminated composites (see Cho, 1997). High-order shear

deformation theories (HSDT) have been developed to account for more accurate

11

Ian

80C

was
dhp
COnl
degr
Cho,
Pitt

Vaha

displacements and stresses (especially the transverse stresses) in laminated
composites (Averill, 1992). The through-thickness displacement field is based on
higher order polynomial functions. The higher order polynomials enforce smooth
and continuous variation of displacements. This smooth continuous variation of
displacement causes double valued stresses at inter-laminar interfaces. With
abrupt changes of material properties at ply interfaces, the displacements and
strain field through-the-thickness are not smooth and continuous (see Cho,
1997). Rather, the displacements can “zig-zag” through-the-thickness of the
laminated composite. With this “zig-zag” variation in mind, the goal of layenrvise
theories is to describe each laminated composite as an assembly of individual
layers and then impose one or more continuity condition(s) at the laminate
interfaces. Accurate displacements and stresses are predicted, but the total
number of degrees of freedom is proportional to the total number of layers in the
laminated composite. In view of the strengths and weaknesses of the high order
and Iayerwise theories, an accurate efficient theory is highly sought.

A new class of laminate theories called first order zig-zag theories (FZZT)
was first developed by DiSciuva (1985). Zig-zag theories use a presumed
displacement field for each layer while analytically enforcing interlaminar
continuity of transverse shear stress while maintaining an independent number of
degrees of freedom per layer. Later, many researchers (see DiSciuva, 1993,
Cho, 1993, Averill, 1996 and Yip 1996) improved the accuracy of FZZT, the
primary improvement was achieved by superimposing a piecewise linear

variation of in-plane displacements on a cubic function of the through-thickness

12

coordinate, which Is commonly referred as higher order zig-zag theories (HZZT).
This improved capturing the kinematics associated with warping of unsymmetric
laminates but required C1 continuity of the approximation of transverse
displacement, thus Hennitian-type interpolation functions are required. This
introduced additional rotational degrees of freedom which are gradients of the
through-thickness displacement. These additional degrees of freedom make it
inconvenient, if not impossible, to implement into a commercial finite element
code that allows only six degrees of freedom per node (three displacements and
three rotations).

Averill and Yip (1995) developed a generalized FZZT and a HZZT that
eliminated the requirement of C1 continuity of the through-thickness displacement
with an interdependent element (anisoparametric) scheme and the penalty
method for beams with a two-noded element. These elements still required
additional degrees of freedom associated with the gradient of the transverse
displacement. Yip and Averill (1996) approached this problem by combining the
strengths of discrete-layenlvise theory and zig-zag theory. The new theory
permits the laminate to be subdivided into a set of sublaminates. Each
sublaminate can contain many layers. Each sublaminate is modeled with zig-zag
kinematics. The FZZT in-plane displacement field in each sublaminate is
assumed to be a piecewise linear function through-the-thickness (Cho, 1997).
The finite element was cast as an eight-noded brick with three translations and
two rotations as degrees of freedom. Yip (1996) expanded the HZZT to include

von Karrnan geometric nonlinear Green’s strains with the Total Lagrangian

13

3F

apr

finh

con
doc
aher
any-

990”

wt h
hequr
vahat
devek
whit
finhe

Meir

 

 

approach. Lee (1998) developed C1 continuity zig-zag plate elements with the

full set of nonlinear Almansi strains with an Updated Lagrangian formulation.

1.3 Present Study

The overall objective of the present study was to develop tools and
approaches for optimization of structures using GA’s. Combining a GA with the
finite element method is by now a familiar approach in the optimization of
structures. Using an Injection Island Genetic Algorithm to perform optimization
considering crashworthiness, modal response, mass, manufacturability,
stochastic variability of design variables and load cases has never been
attempted. Accurate, efficient evaluation tools are necessary when performing
any type of optimization, this is especially true when attempting to predict the
geometrically nonlinear response of composite structures.

For reader familiarization, the mechanics of GA’s and Injection Island GA’s
will be reviewed. An optimization problem involving crashworthiness, modal
frequency response, mass, manufacturability, stochastic variability of design
variables and load cases will be approached with an Injection Island GA. The
development of a geometrically nonlinear first order zig-zag finite element model
will be presented. Verification of the geometrically nonlinear first order zig-zag
finite element model will be confirmed with experimental and computational

methods. Lastly, the optimization of composite sandwich panels will be

14

approached with a GA considering geometrically nonlinear response of sandwich

panels using the first order zig-zag finite element model.

15

Chapter 2

THE MECHANICS OF GENETIC ALGORITHMS

2.1 Introduction

The goal of this chapter is to review the mechanics of simple Genetic
Algorithms and island injection Genetic Algorithms. First, the mechanics of
simple Genetic Algorithms are described. Next, the benefits of approaching an
optimization problem with an island injection Genetic Algorithm are discussed.

Finally, the “pros” and “cons” of Genetic Algorithms are presented.

2.2 Simple Genetic Algorithms

Genetic Algorithms (GAs) are a powerful technique for search and
optimization problems, and are particulariy useful when the design space is large
and complex. Many design spaces are discontinuous, nonconvex and
multimodal. These design space characteristics bring strong possibilities that
many local suboptimal solutions exist, which limit the applicability of gradient-type
approaches. The main downfall of a simple GA is the potentially large number of
design evaluations required to obtain a set of satisfactory solutions.

A GA is a search procedure based on the mechanics of natural selection.
Specific knowledge is embedded in a chromosome which represents the overall
design with a set of design variables. The fidelity of each design variable is
determined by the number of bits in the string defining each locus, or position on

the chromosome. These design variables are the building blocks used to

16

0F
in»
org
per
no
larg
help
The

Pope
pope
dlSplg
initial

lando

 

construct a particular design. Typically, each design variable has the potential to
be mapped from the chromosome to a physical attribute of the actual design.

Designs are created and destroyed by the process of decoding
chromosomes. Natural selection is the connection between a chromosome and
the performance of the chromosome and is a process that causes superior
performing chromosomes to have a higher probability of reproduction, while
inducing poorly performing chromosomes to have a lower probability for
reproduction. During reproduction, two genetic operators are commonly modeled
which can produce new chromosomes: crossover and mutation, a third genetic
operator, inversion is sometimes added. The amount of crossover, mutation and
inversion are all controlled by the user. A high crossover rate will produce new
organisms quickly, but will also have a high probability to disregard higher
performing chromosomes. A low crossover rate could stagnate search, producing
no new chromosomes. Many mutation operators exist, and their main purpose is
largely to aid search by introducing new design variables during search. Mutation
helps maintain diversity and reduces the probability of premature convergence.
The inversion operator was not used in this study.

A set of co—existing chromosomes defines a population, while successive
populations are termed generations. A large genetic pool increases diversity in a
population, which in turn improves the ultimate results of the GA. Figure 2.1
displays the structure of a simple GA. The simple GA begins by creating a single
initial population, wherein chromosomes of different design variables are

randomly created. At this point the "goodness" of each design is measured.

17

Biased by the evaluations obtained, the GA uses unary and binary operators on
these designs to create another population. This population probabilistically
maintains the previously ”good" designs while discarding poorly performing
designs. The program evaluates the new population members and continues
with additional rounds of generation and selection. This is repeated until
satisfactory solution(s) are obtained (see Goldberg, 1988). Incorporating these
processes in a computer algorithm will produce an algorithm that solves

problems in a manner reminiscent of natural evolution.

2.2.1 Representation of Design Variables within a Chromosome

A chromosome of “p” potential design variables (x.) can be represented in
vector form as: {x1, x2, x3, xp}. Each design variable is typically referred to as
a field or locus on the chromosome. The size of each field on the chromosome
will define the number of choices for each design variable. If a coarse
representation of a continuous design variable is sought, then the field size of the
continuous design variable should have a small number of choices between the
lower and upper bounds of the design variable. A larger number of choices
between the lower and upper bounds of a continuous design variable can be
defined to refine the representation of the continuous design variable. Each
design variable (field) on the chromosome can have an independent field size. If
a design variable is truly discrete in nature, than the field size of the discrete

design variable should be set at the number of choices in the discrete set.

18

 

2.2.2 Initial Population Generation

First, the GA must perform initialization and evaluation of the initial
population. The GA stochastically creates a population of chromosomes
composed of a randomly chosen allele at each locus. A measure of “goodness”

or fitness is assigned to each individual based on the fitness function.

2.2.3 Selection

The next process the GA must perform is selection. The selection
mechanism for GAs is simply the process that favors the selection of fit
individuals to be placed in a mating pool. The fitness function provides a means
of comparing individuals and the selection process determines how the
individuals are compared. Selection mechanisms determine which chromosomes
are placed in the mating pool. Roulette wheel, stochastic remainder and rank
based are all common selection mechanisms. Tournament selection was used
exclusively in this study due to its low selection pressure. Tournament selection
allows for a group of individuals that are chosen stochastically to gather and
compete in a tournament, the individual with the highest fitness is placed in the
mating pool. A low selection pressure lessens the pressure to converge by
disallowing any “super" individual to dominate the reproduction cycle, promoting
diversity. The mating pool is comprised of the individuals that will create the next
generation through genetic operators. Selection is the GA search catalyst that

works by mimicking natural selection.

19

2.2.4 Crossover
Crossover is a genetic operation performed by a GA that recombines

building blocks from two chromosomes to create new chromosomes. Pairs of
chromosomes that were placed in the mating pool through the selection process
are combined through crossover. Uniform order-based, cycle, partially matched,
one-point, two-point and uniform crossover are commonly used procedures.
One-point crossover was used exclusively in this study and is explained in detail.
Consider the two chromosomes below that have five design variables:

{X1. X2. X3. 54. >55}

{ Y1. Y2. Y3. 14. Its }
One-point crossover only occurs between two loci. One-point crossover occurring
between the third and fourth loci would result in the following new individuals:

{ Y1. Y2. Y3. £4. £5 }

{X1. X2. X3. 14.25}
This new set of designs have simply been created through recombination of the
original two chromosomes by swapping their fourth and fifth design variables

between the two chromosomes.

2.2.5 Mutation
Mutation is another genetic operator that aids GAs in search and
optimization. The mutation operator helps maintain diversity during search by

introducing new alleles. A mutation operator that acts upon a field-by-field basis

20

de

Th

Th

V3

OCT

2.2

W3!

defi

Wher

the f

of the chromosome was used in this study. A field-by-field mutation of this type is
demonstrated in the following example at the first locus on the chromosome:

{51. X2. X3. X4. X5}
The new chromosome would be:

{513 X2. X3. X4. X5}
The field-by-field mutation would randomly change the value of the first design
variable (bounded by the field size). Typically, low rate of field-by-field mutation

occurs for each field on the chromosome.

2.2.6 Fitness Function

For the problems investigated in this study, the “goodness” of each design
was represented with a single scalar value called the fitness function. The formal
definition of the type of optimization problems addressed in this work is:

Find: {x.} i=1 Number of Design Variables

Such that: Minimize F(Xi)

Subject to:

h}.(x,.) s 0 j=1 Number of Inequality Constraints
5,. (xi) = O k=1,...,Number of Equality Constraints
x,’ s x, s x,“ i=1 Number of Design Variables

State Equations

where {x.} is a vector of design variables, F(xr) is the objective function, h 1. (x,) is

the jth inequality constraint, 54):.) is the kth equality constraint while x: and

21

x: represent the lower and upper bounds on the i‘" design variable. The state

equations relate the behavior of the objective function and constraints to the
vector of design variables.

For all problems investigated in this study, the “goodness” of each design
was represented with a single scalar value called the fitness function. Within the
fitness function the objective(s) are either maximized or minimized, for this
research, the fitness, thus the objective(s), are maximized. All constraint
violations are minimized with the penalty method. To evaluate the fitness of a
design, the objectives and constraints are normalized and aggregated. The

formal optimization statement can be rewritten as:

Find: {xm} m = 1,...,Number of Design Variables
Such that: Maximize: Fitness(xm)
Subject to:

x; S xm S x;

State Equations

The ‘fitness” of a candidate design can be defined as:

 

 

. . n. n
,, ob Ob' n ["4 meq. -t_lneq. J n eq —t_eq k
FimeseA+ :Bi 1 — Z Cj[ ’ 1] - ka[—k-—k] 2-1-)

M or»: 1.0 t_inegi M t_eqk

where “n_obj” is the number of objectives to be maximized, “n_eq” is the number
of equality constraints and “n_ineq” is the number of inequality constraints. “A” is
a large coefficient which insures that the fitness is positive. “Br", “C,” and “Dr.” are

weights associated with the i"‘ objective, j‘“ inequality constraint and the k‘’1

22

equality constraint, respectively. C, are penalty coefficients which are set to zero
when the associated inequality constraint is satisfied.

For a candidate design: “Obj.” is a scalar value representing the ith
objective, ineqj is a scalar value of the jth inequality constraint and “eqk” is a
scalar value of the kth equality constraint. “Norm.”, “t_eq,-", and “t_ineqk" are target
values used as scalars that normalize the objectives, equality constraints and
inequality constraints with respect to their target values. “he” and “n,” are penalty
parameters which are chosen such that a violation of a constraint reduces the
overall fitness of the design. Typically take on values of 2 or greater. This
guarantees the fitness of an infeasible design is penalized in at least a quadratic
fashion with respect to the violation of any constraint. This technique ensures the
fitness rapidly decreases when the candidate design violates any constraint.
Often, optimal solutions lie near or at a constraint boundary, so the penalty

method can help maintain solutions that nearly or slightly violate a constraint.

2.3 Pros and Cons of Genetic Algorithms

A good optimization algorithm should not necessitate excessive
computational costs while progressively finding better solutions. Unfortunately,
no single approach is best for all problems. Fortunately, certain optimization
approaches work better for a specific class of problems. Calculus based hill-
climbing methods work well for an explicit class of problems whose objectives
and constraints change in a smooth, well-behaved manner with respect to the

design variables. Most real life design spaces do not have a convex unimodal

23

sir
pr:

ar

obj
req
rep
Objt
prol
mec
that
thou
Spec
mUct
num
Optirr

addit

 

design space. If a design space is known to be convex, a calculus based hill-
climbing approach should be applied. Calculus based hill-climbing will tend to
proceed to local sub-optimal solutions when extended to problems with complex,
noisy design spaces. This handicaps the applicability of calculus based hill-
climbing methods to many nontrivial problems.

To the other extreme, one technique to alleviate these problems is to
simply enumerate the design space. Exhausting the design space is not
practical if the space is large and or if a single computational evaluation requires
a moderate amount of compute time.

GAs do not require the existence or computation of gradients of the
objective and or constraint functions to perform optimization. Rather, a GA
requires that the fitness of each design can be evaluated. The fitness should
represent a design’s “goodness”, or overall performance in terms of the
objectives and constraints. This makes GAs applicable to a large set of
problems. A GA is a multi-point search procedure modeled specifically on the
mechanics of natural selection. A simple GA is a multi-point search procedure
that occurs within isolated surroundings, interacting with nothing else. Even
though a GA requires fitness evaluation for only a fraction of the total design
space, the total number of evaluations to find an optimal solution can still be
much too large for practical applications. This requirement represents the
number one downfall of a GA. This is tme for GA applications in structural
optimization using computational methods such as the finite element method. An

additional downfall to a simple GA is its propensity to converge to suboptimal

24

solutions. Parallel GAs attempt to mitigate both of these issues by allowing
search to happen on a set of multi-point searches (often on separate
computational nodes) with various topological structures allowing structured

migration amongst the search sets.

2.4 Island Injection Genetic Algorithms

An island injection Genetic Algorithm (iiGA) also uses structured migration
of search points amongst a group of search sets, but allows for the exploration of
various levels of resolution of the search space to reduce the overall
computational time required to arrive at a given set of well performing solutions.
iiGAs search at various levels of resolution to quickly find building blocks at low
levels of spatial resolution to inject into higher levels of resolution. iiGAs can
independently search for optimal solutions considering multiple objectives
concurrently while handling stochastic variabilities. An island injection Genetic
Algorithm (iiGA) represents an approach to search at various levels of resolution
within a given space in a structured hierarchical fashion (see Eby, 1996,
Goodman, 1996). This includes first searching at low levels of resolution on
different nodes (islands) and then injecting the high performing individuals into
islands of higher resolution to “fine-tune” the design. If the formal definition of the
fitness is known to change when migrating individuals, the receiving
subpopulation must reevaluate the migrants performance. Figure 2.2 displays a
typical island injection topology that has an island that searches at a low level of

resolution and migrates good designs to islands of higher resolution. The search

25

space is fundamentally divided into hierarchical levels with well defined overlap
(the search space of the parent is contained in the search space of the child).
The rate and structure of migration can be set for each island. Each island could
independently represent the problem in terms of design variables, objectives,
constraints, fitness definition and fitness evaluation. The design variables,
objectives and constraints on each island could be considered stochastic or
deterministic. The rates of each genetic operator could be independently set for
each island. Each island could be evaluated on separate computational nodes.
Many variations on the island topology in Figure 2.2 can easily be performed for
problem dependent custom tailoring. iiGAs embody a divide-and-conquer and
partitioning strategy which has been applied to many problems (see Eby et al.,
1999 and Mallot et al., 1996). The simple GA and iiGA applied in this dissertation

were developed by Goodman (1996).

26

 

Create the Initial Population

1

. Evaluate Fitness

Create the Next Generation

 

 

 

 

 

 

 

 

with Genetic Operators, Selection

 

 

 

    

No

 

Convergence?

 

 

Optimal Solutions

 

 

 

Figure 2.1. Simple Genetic Algorithm Structure.

 

_, ff

. Coarse Design Spacelasland 1)

 

 

 

 

Refined, Design Space (Islands 2 and 3.)

4;. L LMAJ - A...

 

Figure 2.2. A simple iiGA topology that simultaneously searches various levels of

resolution of the design space on three separate computational nodes (islands).

27

Chapter 3

OPTIMIZATION OF CRASHWORTHINESS AND MODAL FREQUENCIES OF
AUTOMOTIVE STRUCTURES

3.1 Introduction

The goal of this chapter is to demonstrate an island injection Genetic
Algorithm’s ability to perform optimization of crash type problems while
concurrently satisfying a set of conflicting constraints considering stochastic
variability of design variables and load cases. The traditional optimization
statement for this problem will be developed in terms of the objective functions,
constraints and design variables with respect to an island injection Genetic
Algorithm. Results display an island injection Genetic Algorithm's ability to
maximize multiple objectives while concurrently satisfying multiple conflicting
constraints. A beam-type structure typically referred to as an automotive lower
front motor compartment rail will be the particular automotive structure
considered during this research. The front motor compartment rail is considered
the main mechanism that governs mechanical transfer of energy during a crash
event. The front motor compartment rail connects the front bumper to additional

beam-type stmctures along the front of dash of an automobile.

3.2 Optimization Problem Statement
In order to define the optimization problem this section describes 1)
optimization design variables, 2) optimization objective functions, and 3)

optimization constraints. The optimization problem can be stated as: find the

28

cross sectional shape and thickness varying along the length of a rail structure to
maximize the energy absorbed by a specific region within the structure subject to
inequality constraints on the first and second modal frequencies, peak crush
force, mass and manufacturability while considering multiple load cases with
stochastic variability of loads and design variables. The objectives were
decomposed by simultaneously (but independently) seeking deterministic optimal
solutions for each load case in hierarchical heuristic levels of refinement.
Concurrently, the objectives were stochastically represented with stochastic
variation of the design variables. An objective function can be considered
stochastic when non-unique objective values are possible for a unique design.
The stochastic nature of the objective functions can be defined by the stochastic

nature of the design variables.

3.3 Optimization Design Variables

Every design variable in structural optimization can be thought of as being
either discrete or continuous in nature while behaving in either a deterministic or
stochastic manner. A design variable behaves stochastically when the mapping
of the intended design variable to the actual design is non-unique. A design
variable acts stochastically when a design is produced multiple times and the
produced designs are not necessarily identical. Each time the design is
produced, a slightly modified design could be created. This inherently occurs in
many processes - each time a design is manufactured, slight differences in the

manufacturing process will produce slightly different design variables. A design

29

 

SE
Sll
8V

lint

fidr
cn‘t
bet

are.

variable behaves deterrninistically when a unique mapping from the intended
design variable to the actual design occurs. Given a set of deterministic design
variables, a unique design will be repeatedly produced, irrelevant of the process
that creates them. Many real life design variables are stochastic in nature.
Stochastic variability of design variables should not be disregarded for problems
that involve noisy chaotic objectives. Nonlinear response of structures is often
very sensitive to small changes in the design variables.

The cross sectional shape of the structure was viewed as a discretized
composition of cross section primitives or structural elements. With this in mind,
the first set of design variables was defined as a deterministic discrete set of
cross sectional primitive shapes. The cross sectional primitives allow variation of
the shape of the cross section along the axis of the stmcture, see Figure 3.1. A
total of five cross sectional primitives biased towards the front of the rail defined
the overall surface of the design, using interpolation from cross section to cross
section with a linear blending function (see Figure 3.2). The surface was
smoothed with a fourth order Bezier surface at the cross sectional primitives to
avoid sharp changes in geometry that would be artificially generated from the
linear blending of the primitive cross sections.

Biasing of the design variables was performed to increase the model
fidelity in critical portions of the stmcture. Refinement of the design variables in
critical regimes of the stmcture was performed to control complex localized
behaviors. Coarsening of the design variables was performed in non-critical

areas of the structure where localized behavior was of a lesser concern. Since a

30

progressive crush is sought, the cmsh of the front of the rail is crucial. Therefore,
the placement of primitive cross sections was biased towards the front of the rail
in the “crush-zone". Figure 3.2 and 3.3 display biasing of the design variables
towards the front of the rail structure.

The second set of design variables was the thickness of each section of
the rail, which was considered a continuous stochastic design variable. The
manufacturing process of the rail stmcture was considered stochastic - each
time a target design was created the thickness of the rail stmcture would change
within ten percent of the target thickness. The stochastic noise in thickness was
considered to play an important role in the objective functions. Thickness
changes occurred between the placement of each cross sectional primitive, see
Figure 3.3. Since the surface of the design was generated with five primitive
cross sections, there were four design variables describing the thickness of the
rail structure. The rail structure was made out of aluminum.

As depicted in Figures 3.1 and 3.2, the cross sectional shape of the rail
was allowed to vary as a function of the rail axis by choosing a set of primitive
cross sections. The choice of primitive cross sections was bounded by a discrete
set of manufacturable primitive cross sections. The iiGA could choose any
combination of cross sectional shapes from a predefined library of
manufacturable primitives. This library of cross section primitives is not arbitrary
— rather, it is knowledge based. Figure 3.4 displays a primitive library that
contains a set of knowledge-based manufacturable cross sections. A design

volume constrained the rail structure to have a smaller width compared to the

31

height of the rail cross section. Each cross section in Figure 3.4 was scaled up
to 1.5 times the dimension of the original perimeter to define the upper and lower
bounds for each cross section design variable.

The representation of the design variables within the chromosome was set
up in an order that staggered the primitive cross sections with the thickness of
each section. The order of the design variables on the chromosome was
constructed to reflect the physical attributes for a given design which has been

known to aid search within a GA (see Holland, 1975, Goldberg, 1988).

3.4 Optimization Objective Functions

The optimization objective function is the maximization of the internal
energy of a specific region of the structure termed the “crush-zone” for multiple
load cases. The internal energy of the crush-zone was computed by LSDYNA,

(see LSTC/KBSZ, 1999). The objective function for a single load case can be

stated as:
Elli
1.11:1; jsi.dei.dc|.-., 31.)
2no 1 1

where Q is the volume of interest called the “crush-zone”, Se is the second
Piola-Kirchoff stress tensor and 8,,- is Green’s strain tensor. The equations of

motion were simulated with LS-DYNA (see LSTC/KBSZ, 1999) to obtain the
objective functions. LS-DYNA is an explicit nonlinear finite element code.
The upper bound of the internal energy of a feasible design was computed

by considering the inequality constraint on the maximum crush force and the

32

lur

St

3.5

we

be

We

maximum distance the rail could displace during crush. This value was used to
normalize the objectives.

Load cases considered were meant to simulate a component of a front
end of an automobile during various crash events. During frontal crash events,
the lower front rail is often designed to be the main mechanism to govern crash
energy. The crash event was modeled by defining an initial velocity for the lower
front rail structure and allowing it to move into a stationary rigid wall. A lumped
mass (with identical initial velocity as the rail) was constrained at the end of the
rail structure to simulate inertial effects that would drive the rail into the rigid
barrier. Different crash events were modeled with various load cases. The first
load case was meant to simulate a rail structure traveling at 35 mph colliding with
a rigid barrier. For the first load case, a lumped mass was placed directly behind
the rail structure as depicted in Figure 3.5. The second load case was meant to
simulate a rail structure traveling at 40 mph colliding with a rigid barrier with
additional moments into the system. The second load case placed an “eccentric”
lumped mass behind the rail at an offset to create additional moments into the

system as shown in Figure 3.6.

3.5 Optimization Constraints

The first and second inequality constraints in this optimization problem
were to ensure that the first two modal frequencies were greater than or equal to
predetermined target values. The cross sectional properties of the structure

were computed and the modal response of the structure was based on a beam

33

analysis performed by ADINA (see ADINA, 1999). It is well documented that an
optimization problem that considers at least two modal frequencies is non-convex
and multimodal with portions of the design space that have undefined gradients
(See Haftka and Gurdal, 1993, and Gurdal, 1999) due to the possibility that order
of modal frequencies may switch as search occurs.

The third inequality constraint was that the maximum force experienced
during crush be less than or equal to the target force. The peak crush force was
computed at the rigid barrier by LSDYNA, (see LSTC/KBSZ, 1999). A definition
of work is the “dot” product of force and distance, while energy can be defined by
integrating work through time. With this in mind, there are two ways to increase
the objective function for a specific load case: the internal energy of a structure
can be maximized by increasing either the force or distance vectors during crush.
The Euclidean length of the force vector was constrained. To this end, we can
conclude that the constraint placed on the maximum crush force will play a major
role in limiting the amount of energy absorbed during crush.

The fourth constraint in the optimization problem was to insure that the
total mass was less than or equal to a target value. Equation 3.1 defines the
objective function of a single load case as the internal energy of a specific
volume of the rail structure. An increase in volume would obviously increase the
total mass of the rail structure. For a given design, the total energy of the system
will remain constant. Kinetic energy will be the only form of energy that is input
into the system. An increase in volume (mass) will slightly increase the initial

(kinetic) energy of the system. The internal energy absorbed by a design is

bounded by the total energy into the system (kinetic energy), which in turn is
bounded by the volume (mass) of the design.

The last constraint measured is the manufacturability of the design.
Manufacturabil'rty constraints were enforced by ensuring good designs had a
maximum “local-taper” less than or equal to a target value. The “local-taper” of a
design was the maximum amount of change in the size of the perimeter of the
rail as a function of the rail axis. This constraint limits the amount of load that the
structure could carry in localized bending. Localized bending is often required to
initiate crush.

All constraints were evaluated prior to a “full” crush analysis. The
maximum force experienced during crush often occurs shortly after the instant of
contact of the rigid barrier with the rail substructure. If all constraints were
satisfied within 30 percent of the target values, then a full crush evaluation was
performed to determine total internal energy of the crush zone and the actual
maximum force experienced during full crush. If a design severely violated any of
the constraints, then the fitness was penalized by assuming the internal energy
of the crush-zone was zero. This technique reduced the total number of crush
evaluations of “poor" designs that harshly violated constraints.

The fitness for each design was computed using equation 2.1. Target
values for all constraints are given in Table 3.1. All coefficients used in the fitness
function are given in Table 3.2. The weights for the fitness function are all equal
except for the manufacturability weight. The optimization process should not be

dominated by manufacturability constraints.

35

3.6 Application of an Island Injection Genetic Algorithm

An iiGA approaches an optimization problem by simultaneously searching
at various levels of resolution within the design space in a heuristic hierarchical
fashion. The following paragraphs present the framework of the iiGA applied
within the context of the current application.

Figure 3.7 shows an iiGA topology that breaks down the total design
space into three levels of design variables and time-based resolutions. The
overall iiGA topology seeks designs which maximize the internal energy for the
two load cases discussed previously. Each island independently searches on
one of nineteen separate computer nodes (termed islands) while simultaneously
sharing information in an ordered fashion. Each computer node had 400MHz
processing speed, on a P586 architecture. Table 3.3 gives genetic operator rates
for all the islands in Figure 3.7.

Each level of resolution was defined by the discretization of the design
variables (rail section thickness and cross sectional primitives) and crush
evaluation time. Islands 0 through 2 have a low level of design variable
representation with a small duration of simulated crush, islands 3 through 8 have
a medium level of design variable representation with a larger crush evaluation
time, while islands 9 through 18 have the highest level of design variable
representation and largest crush evaluation time. The total number of possible
designs grows with an increase in the level of resolution (discretization) of the

design variables. Each cross section depicted in Figure 3.4 was bounded by 35

36

isl

pr

by 127 and 55 by 200 millimeter design areas. Each primitive had an equal
number of choices scaled uniformly between the upper amd lower bounds of the
design areas. Table 3.4 contains the total number of cross sectional primitive and
thickness choices. The lower and upper bounds on the thickness were 2.5 and
6.5 millimeters.

Figure 3.7 also displays the structured mles of migration within the iiGA
topology. The coarse design spaces were constructed to map into the refined
design spaces in a hierarchical fashion. In general, a unique mapping from a high
to low level of design variable resolution cannot be defined, so designs cannot be
shared from a high to low level of design variable resolution. The best design
plus a randomly chosen design was injected according to the iiGA structured
migration rules as depicted in Figure 3.7 at the end of every generation.

The evaluations performed with LS-DYNA are dynamic, thus time
dependent. An evaluation that measures the performance of design over a small
amount of time is computationally cheaper when compared to evaluating the
identical design over a longer period of time. This fact can be taken advantage of
by allowing the iiGA to search with evaluation times that can vary from one set of
islands to another. This can also benefit the development of a design that
progressively crushes. Table 3.5 contains evaluation times for all islands.

Many islands have stochastic design variables and load cases. Islands 0,
2, 3, 8, 9 and 18 evaluate a single deterministic objective function. Islands 0, 3
and 9 evaluate designs based on the first load case. The first load case placed a

lumped mass directly behind the rail with an initial velocity of 35 MPH. Islands 2,

37

8 and 18 evaluate designs based on the second load case. The second load
case placed an eccentric lumped mass at an offset behind the rail with an initial
velocity of 40 MPH. Islands 1, 4 through 7, and 10 through 17 considered
stochastic variability of load case and rail section thickness. For these islands,
the load case was randomly chosen for each evaluation. Each load case had an
equal probability to be chosen for these islands. The thickness of each finite
element was allowed to vary within 10% of the original variable.

The iiGA topology in Figure 3.7 was also dynamic in nature: islands 0
through 2 evaluated 30 generations of designs and then were reinitiated to
monitor the overall progress of the topology, see Figure 3.8. Upon re-initiation,
the islands evaluated the best designs from the highest level of resolution from
the iiGA topology in terms of each load case. This set of islands were meant to
monitor the overall performance of the iiGA topology, not to perform additional
search. No genetic operations were performed on the designs injected into the

monitoring islands.

3.7 Results

The performance of the iiGA topology is represented in Figures 3.7 - 3.8
for each deterministic objective (load case). The overall performance of the iiGA
topology will be measured by viewing the behavior of the objectives and
constraints for a monitoring island in Figure 3.8.

Figure 3.9 displays a progressive increase in the deterministic objective

function of internal energy of the crush zone for the first load case of the best

38

a:

the
as
dis
fun

40,

3.l

ma'

 

fun

 

design found as a function of generation number for island 9 in Figure 3.8.
Figure 3.10 displays the deterministic normalized constraints for the best design
found as a function of generation number for island 9 in Figure 3.8. It is
important to realize that a feasible design has normalized mass, force, and
manufacturability constraints which are less than or equal to unity. A feasible
design also has normalized first and second modal frequency constraints which
are greater than or equal to unity.

Figures 3.9 and 3.10 display a search that allowed for tradeoffs in
constraint violations in order to maximize the deterministic objective of
maximizing the internal energy of the crush zone for the first load case. Many
constraints are violated throughout the evolutionary process depicted in Figure
3.10, but the final design satisfies all the constraints except for the force
constraint, which violates the target force by only 2%.

Figure 3.11 shows the deterministic objective function of internal energy of
the crush zone, for the second load case, increasing for the best design found,
as a function of generation number for island 18 in Figure 3.8. Figure 3.12
displays the deterministic normalized constraints for the best design found, as a
function of generation number for island 18 in Figure 3.8. The run lasted a total of
40 generations and the best design of the run was found at generation 35.

Similar to the results portrayed in Figure 3.9 and 3.10, Figures 3.11 and
3.12 display a search that allowed for tradeoffs in constraint violations in order to
maximize its deterministic objective for island 18. The deterministic objective

function of island 18 was the internal energy of the crush zone for the second

39

load case. Figure 3.11 shows a slight decrease in internal energy of the crush
zone for the second load case from generations 9 to 10 and 19 to 20. This is can
be explained in Figure 3.12 by inspecting where the constraints were violated. A
sudden drop in overall constraint violations occurred over generations 9 to 10
and 19 to 20. This was reflected in the measure of fitness of the designs by
enforcing the constraints with the penalty method.

Figures 3.13 and 3.14 depict the overall performance of the iiGA topology
considering each stochastic objective and constraint. Figures 3.13 and 3.14 show
the objectives increasing along with the behavior of the normalized constraints.
It should be stressed that the monitoring islands in Figure 3.8 were initiated after
islands 0 through 2 completed 30 generations of search, which spanned about
three days. During this initial time, the entire iiGA topology was performing
search seeking good solutions for each load case with stochastic variabilities - in
other words, a huge amount of work was performed before the monitoring islands
were spawned. Figures 3.13 and 3.14 plot the objectives and constraints as a
function of evaluation number due to the fact that the monitoring islands did not
perform any search. The iiGA topology compute time was six and a half days.
Table 3.6 contains the number of evaluations and full crush evaluations. If a
design did not satisfy the constraints, then the objective function was not
evaluated, eliminating the LSDYNA nonlinear finite element evaluation. A full
nonlinear finite element evaluation took about 13 minutes of compute time for a

single load case.

40

Figures 3.15 and 3.16 animate the deformation of the best design found
for the first and second load cases, respectively. The design crushes
progressively from the front to rear of the crush zone for both load cases. Two
additional independent runs were performed to verify that the algorithm
converges to a similar best design. The three best designs found in the three
independent runs are depicted in Figures 3.17 through 3.19. Table 3.7 contains
the thickness of each section for all three designs. All three designs closely
resemble each other in overall shape and thickness - this is an important fact
because it displays the iiGA topology repeatedly converged to similar solutions
even though the algorithm started with different initial search sets. The objective
and constraint function values are listed in Table 3.8 for all three designs. All
three designs have near identical objective and constraint function values. To
verify that the objectives are not highly sensitive to stochastic variability of
thickness, twenty additional independent evaluations were performed on the best
designs depicted in Figures 3.17-3.19. It was found that the standard deviation of
internal energy was small compared to the average internal energy. Table 3.9
contains the average and standard deviation of the internal energy for the design

portrayed in Figures 3.17-3.19.

3.5 Conclusions
The study developed an approach that overcame the shortcomings of
available optimization techniques by allowing concurrent optimization of

structures with an island injection Genetic Algorithm for: 1) crash energy

41

management, 2) modal vibration frequency response, 3) peak crush force, 4)
weight reduction, and 5) manufacturability, while considering stochastic variability
of design variables and loading conditions. Existing automated design technology

has never simultaneously addressed all these competing mechanisms.

42

 

Primitive Cross Sections \§‘

 

 

 

 

 

Front of Rail
R63I 0f Rail Side View of Rail

Primitive Cross Sections N

Top View of Rail

 

, Front of Rail
Rear of Rail

Primitive Cross Sections

Rear of Rail / \g\}

Isometric View of Rail

Figure 3.1 Biasing of cross sectional primitives along the length of a rail structure.

43

Figr
SfCl

Prir

 

 

Front of Rail
Rear of Rail

Side View of Rail

‘ . Front of Rail
Rear Of Rail Top View of Rail

Rear of Rail

Front of Rail

 

Isometric View of Rail

Figure 3.2 Surface of rail structure is developed with a linear blending of the cross
sectional primitives. The surface was smoothed at each of the cross sectional

primitive locations.

Primitive Cross Sections

[f \m

Isometric View of Rail

      
 

First Thickness

Second Thickness
Third Thickness

Fourth Thickness

.55.

V

Figure 3.3. Potential thickness changes occurred at the placement of each cross

sectional primitive.

45

 

 

 

II Zlhlh
- EDIT

Figure 3.4. Cross sectional primitive library that contains a set of knowledge-based

 

 

 

 

 

 

 

manufacturable cross sections.

46

 

Lumped mass Z

Figure 3.5. Finite element model and of rail structure running into a stationary
rigid wall for the first load case. A lumped mass was placed directly behind the rail

structure.

47

Z

Eccentric ltunped mass / z

 

Figure 3.6. Finite element model and of rail structure running into a stationary
rigid wall for the second load case. An eccentric lumped mass was placed behind the
rail structure.

48

L-

 

 

 

 

 

 

   

 

Load Case 1 Stochastic Load Cases Load Case 2
Deterministic Stochastic Design Variables Deterministic
Low
0 F /K‘ Resolution
3.8 ms
_———> g/GAj/l :.Vq4—A Medi“"‘_l
Resolution
8.4 ms
High
Resolution
12.6 ms

 

 

 

 

 

Figure 3.7. An iiGA topology that breaks down the design space into three levels of

design variable and evaluation time based resolutions with multiple load cases.

49

Load Case 1
Deterministic

Stochastic Load Cases
Stochastic Design Variables

Load Case 2
Deterministic

 

—M

@

  

 

 

ere—1r

 

 

 

Weir

 

 

 

 

 

 

a ‘37

Monitoring Islands

 

 

Medium
Resolution
8.4 ms

High
Resolution
12.6 ms

Figure 3.8. An iiGA topology that monitors islands the performance of the best
designs for the two objectives defined by two structural load cases.

50

rm;

 

'—'m=HGH=h—I

anonm

 

 

 

0 5 10 15 Z) 25 ED 35 40

Generation Number

Figure 3.9. Crush zone internal energy as a function of generation number for load

case 1, island 9.

51

1.4

 

 

 

 

 

 

 

 

 

 

N
o ..
r 1.2k . .
m - i"
a [/5\ A
1 "-/ \V/
N: I

l I
Z
c .. ’e
d
C
O A

V V -B-Nhss
n

+Forca
S cal
t +81me
1‘ +Seca'dFm
a
i 02. +mmnbity
n
t
o r 77 r r r r 1 1
O 5 10 15 20 25 30 35 40
Generation Number

Figure 3.10. Normalized constraints as a function of generation number for island 9.
A feasible design has normalized mass, force, and manufacturability constraints
which are less than or equal to unity. A feasible design also has normalized first and

second modal frequency constraints which are greater than or equal to unity.

52

r—‘flflHQc-o-fli—i

wafer-105m

EE
3

 

 

O 5 10 15 20 25 3O 35 40

Generation Number

Figure 3.11. Crush zone internal energy as a function of generation number for load
case 2, island 18.

53

 

moN—owmanoz

 

 

 

 

 

 

 

9
0.7 u U U
C
as
0 t
-B-Nhs
11 est
S +Foroe
t 0.44
-)(—FrrstFreqrarioy
1'
a o.3~ +Seconchremaw/
i 02. +le
n
t 0.1 *
o w w . e e
o 5 10 15 20 as so as 40
Generation Number

Figure 3.12. Normalized constraints as a function of generation number for island
18. A feasible design has normalized mass, force, and manufacturability constraints
which are less than or equal to unity. A feasible design also has normalized first and

second modal frequency constraints which are greater than or equal to unity.

54

 

 

 

 

 

 

 

E 15000; +Load Case 1

n

e -B-Load Case 2

1' 10000 .

8

y 5000 «
kN

o . . . . . . . . . 4
mm 0 55 110 105 220 275 330 385 440 495 550
Evaluation Number

Figure 3.13. Crush zone internal energy as a function of evaluation number for load

cases 1 and 2 from monitoring island.

55

 

 

 

 

 

 

 

 

 

N 1.25
o
r
m
a 1 e
l
i
z
e 0.75 «
d
C CPA) —B—Mass
O 0 5 .
n ' +Force
s
t —>(—First Frequency
r
a 0.25 - +Second Frequency
:1 -e—Manufacturability
t
0 55 110 165 220 275 330 385 440 495 550
Evaluation Number

Figure 3.14. Normalized constraints as a function of evaluation number for monitor
islands. A feasible design has normalized mass, force, and manufacturability
constraints which are less than or equal to unity. A feasible design also has
normalized first and second modal frequency constraints which are greater than or

equal to unity.

56

Side View Top View

 

Rear Front
Rear Front
<—-—>
Crush Zone

 

Crush Zone

time =0.0 ms

— ~

time = 4.0 ms

 

 

time = 8.0 ms

1 ——a

 

time = 13.64 ms

Figure 3.15. Animation of deformation for the best design found (load case 1).

57

Side View Top View

 

Front
Rear
Rear Front
<——>
‘——> Crush Zone

Crush Zone

time = 0.0 ms

time = 4.0 ms

time = 8.0 ms

 

   

time = 11.3 ms

Figure 3.16. Animation of deformation for the best design found (load case 2).

58

 

Figure 3.17. Surface models for best design found, run 1.

59

 

Figure 3.18. Surface models for best design found, run 2.

 

Figure 3.19. Surface models for best design found, run 3.

Table 3.1 . Target values of inequality constralnts for fitness function.

 

 

 

 

 

 

 

 

 

 

 

C1 C2 C3 C4 C5
120kN 5.1 kg 75 Hz 115 Hz 0.1
Table 3.2 Weight values used for fitness function.
A Bl, El'E4 E5
1000000 500.0 100.0

 

 

 

 

 

Table 3.3. Genetic operator rate for each Island In the Island injection
Genetic

 

utatlon
Number Size

1

62

 

Table 3.4. Primitive and thickness design variable lnforrnation.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Number of thickness Number of Primitive
Island Numbers . .
chorces chorces
0 through 2 16 32
3 through 8 64 128
9 through 18 128 256
Table 3.5. Crush evaluation times per island and load case.
Island Numbers Evaluation Time Evaluation Time
(Load Case 1) (Load Case 2)
0 through 2 3.84 3.36
3 through 8 8.4 7.35
9 through 18 12.36 10.82
Table 3.6. Number of full and evaluations.

 

 

Table 3.7. Thickness design variables for each best design.

 

 

 

 

 

First Second Third Fourth
Run Number Thickness Thickness Thickness Thickness
(In!!!) (nun) (MD) (MI!)
1 5.34 5.00 3.98 3.85
5.50 4.75 4.15 3.90
3 5.25 4.85 4.05 4.00

 

 

 

 

 

63

 

 

 

 

Table 3.8. Objective and constraint data for the best design found in
three different independent runs.

 

 

 

 

 

lntemal lntemal
Energy Energy Peak First Second
Run Load Load Mass Modal Modal
Force Manuf.
# Case 1 Case 2 (kg) (kN) Resp. Resp.
(kN- (kN- (HZ) (HZ)
mm) mm)
1 25298.6 24527.1 5.12 123.0 82.1 126.0 0.8
2 25611.2 24201.4 5.05 125.1 84.0 125.8 0.74
3 25010.6 24727.6 5.08 119.5 83.5 124.9 0.70

 

 

 

 

 

 

 

 

Table 3.9. Stochastic variability of lntemal energy for the best
designs found (20

independent evaluations) with respect to

 

 

 

 

 

thickness.
Standard Standard
Average Average
lntemal Dev. lntemal Dev.
Run lntemal lntemal
Energy Ener Energy Ener
Load Case 1 3” Load Case 2 3"
(kN-mm) Load Case 1 (kN-mm) Load Case 1
(kN-mm) (kN-mm)
1 24981.6 405.1 24305.4 601.9
2 25111.2 367.1 24150.3 613.2
3 25310.6 503.5 23891.1 643.1

 

 

 

 

 

 

64

 

Chapter 4

DEVELOPMENT OF A GEOMETRICALLY NONLINEAR ZIG-ZAG FINITE
ELEMENT MODEL

4.1 Introduction

Analysis of laminated composites and sandwich panels often requires a
compromise between accuracy/complexity and efficiency/simplicity. While first-
order shear deformation theory (Yang, et al., 1966) and higher order theories
(see Reddy, 1990) are able to predict the overall structural response of most
relatively thin laminated structures, these theories often pooriy predict the
magnitudes and distributions of local strains and stresses. In sandwich panels,
these theories, as well as more specialized but traditional sandwich theories (see
Vinson, 1999; Zenkert, 1997), fail to properly predict core shear and compressive
stresses, cannot be applied to the full range of core densities and types, and
cannot handle the extremes of skin thicknesses and moduli. At the other end of
the spectrum, fully three-dimensional finite element models typically require
significant through-thickness discretization (e.g., one element per layer),
involving large computing power and an onerous model development effort. The
above issues are further exacerbated whenever a nonlinear analysis is required.

The goal of most new laminate modeling approaches is to increase
accuracy while decreasing modeling effort and/or computational cost. Toward
this goal, a new technical theory and associated finite element model have been
developed for analysis of thick or thin laminated composite and sandwich panels

(Cho and Averill, 2000). The theory employs the sublaminate concept such that

65

each computational layer contains several physical layers. Within each
sublaminate, an accurate approximation of the displacement field accounts for
discrete layer effects without increasing the number of degrees of freedom as the
number of layers is increased. This is accomplished by analytically satisfying the
continuity of transverse shear stresses at layer interfaces. Because the resulting
through-thickness variation of ”in-plane” displacements takes the form of a
piecewise linear function, the theory is often called zig-zag theory. Depending on
the laminate layup and the type of global and/or local response measures that
are sought, the optimal number of sublaminate approximations required for
accurate analysis of thick multilayered composite laminates or sandwich panels
is far less than the number of layers in the laminate. Often, only one sublaminate
is required.

This model has been cast in a form that is very convenient for use in finite
element analysis, namely, eight-noded brick and six-noded wedge three-
dimensional elements have been developed that behave very much like plate
elements, but have nearly all of the advantages of traditional three-dimensional
elements (Cho and Averill, 2000). In particular, these layered solid elements can
be stacked in the thickness direction of a laminate, allowing any desired degree
of through-thickness discretization. Further, these elements may have very large
aspect ratios without exhibiting shear looking or Poisson ratio locking. The three-
dimensional zig-zag sublaminate finite elements contain only the typical five

engineering degrees-of-freedom (i.e., three translations and two rotations) at

66

each node, so they are compatible with the requirements of most commercial
finite element packages.

The first-order zig-zag sublaminate theory and its associated finite
elements have been shown to be very accurate and robust for predicting the
linear stmctural response of a wide range of laminated composite and sandwich
panels. However, the effectiveness of this model in predicting nonlinear response
has not been thoroughly investigated. The zig-zag sublaminate approach seems
to have great promise for predicting progressive failure, sandwich skin buckling,
and other nonlinear structural behaviors. Hence, recent attention has been given
to the development of geometrically nonlinear zig-zag models. Geometrically
nonlinear models of various forms of the zig-zag plate theory have been
developed by several investigators (Yip, 1996; Lee, 1998; DiSciuva, 1996,1997).
The combination of accuracy and simplicity of the first-order zig-zag sublaminate
model makes it more convenient for general purpose use compared to the
available higher-order zig-zag type models. Thus, in the present work the first-
order zig-zag model (Cho and Averill, 2000) is extended to account for geometric
nonlinearity due to moderate rotations of the plate reference surface. Such an
extension is motivated by the need to simulate nonlinear phenomena such as
those mentioned above. Here, the new nonlinear model is developed and

validated.

67

 

4.2 Basic Ingredients Needed For The First-order Zig-Zag Plate Theory

In this section a brief description of the geometrically nonlinear first-order
zig-zag sublaminate plate theory is provided. Many details related to the
mathematical development of the zig-zag coefficients are neglected, and
attention is focused on the extension of the model to account for moderately
large rotations. For a complete description of the geometrically linear theory and
the associated numerical validation studies, the reader is referred to (Cho and
Averill, 2000).

It is assumed that the laminate is composed of N, perfectly bonded

layers. The thickness and material stiffness properties may vary arbitrarily from

layer to layer. The laminate is modeled as M sublaminates. with each
sublaminate containing N," layers, where m is the sublaminate number.

Mathematically, this is represented as:

M

M=ZM an

In order to facilitate the development of the theory, all expressions in the
following derivation pertain to the m-th sublaminate. The sublaminate number
designation m is omitted for brevity. The reference surface of the m—th
sublaminate is assumed to be the bottom surface of the sublaminate, as shown
in Figure 4.1.

The global X3-axis is taken perpendicular to the in-plane X1, X2 coordinate

axes and has its origin at the bottom of the laminate, while x3 is a local

68

sublaminate coordinate in the direction of X3 with its origin at the bottom of the
sublaminate (See Figure 4.1).

The constitutive relations for a three-dimensional stress state in the k-th
layer of the m-th sublaminate can be written as:

(n‘ F (k) (k) (k) at)“ (k) ‘
IS“ C.. C.. Ce 0 0 C.. is.

(k) (k) (k) (k) (k) (1:)
S22 C12 C22 C23 0 0 C26 822

Sac) Cm Cm Cm O 0 Can (k)
< 3(1) }= 13 23 33 m (k) 36 < 83:) q (42.)
See 0 0 0 C. C. 0 7..

SI.“ 0 0 0 Cl? (35‘? 0 73"
Sun Cm Cm Cm O 0 Cm Ly“)
t 12 j _ I6 26 36 66 _ 12 ,

 

 

 

 

 

 

where 85’" are the components of the second PioIa-Kirchhoff stress tensor and

5;“ are the components of the Green strain tensor. For most practical laminated

composites, the 13 coefficients of the material stiffness matrix are obtained from
9 independent material constants and a transformation law (see Jones, 1999).

If strains are small compared to unity but rotations of the plate reference
surface are moderate, then the Green strain tensor can be degenerated to a
simplified form commonly referred to as the von Karrnan strain tensor, defined in

vector format as:

 

 

 

 

 

 

i811 i “1.1 2 3"
(k) (k)
822 “2,2 _1_(u(k)
a”) u“) 2 3'2
3.3
< 3’ i=< l+< 0 i (4.3)
(k) u“) +u(k)
723 2.3 3.2 0
(k) (k) (k)
713 “1.3 +u3,1 O
(k) (k) (k)
i712 a ium +11“ 2 (k) (k)
k “3,1 “3.2 J

69

where uf“ is the displacement of a point in the x; direction and ( )1. implies

differentiation with respect to x).
The sublaminate displacement fields are initially assumed in the following

form:

k—l
uf“ (x,,x2,x3,t) = ub + x3311, + 2(x3 —x3,.)-§,.

i=1

k-l
ug") (x,,x2,x3,t) = vb + x3 -‘I’2 + 2(x3 — x3,)-77,. (4.4)

M
ug") (x,,x2,x3,t) = wb -(1-%)+ w, (Eli-J
where t represents time, h is the thickness of the sublaminate, u, and v, are the
axial displacements in 2:1 and x2 directions, respectively, at x3 =0, LP, and ‘P,
are rotations of the normal at 2:3 =0, and w, and w, are the transverse
deflections of the bottom and top surfaces, respectively, of the m-th sublaminate.
Thus, it is assumed that u,""and ug") vary in a piecewise linear fashion through-
the-thickness of a sublaminate and u§"’ varies lineariy through-the-thickness. The
parameters 5, and 77,. in equation 4.4 are eliminated by enforcing continuity of

transverse shear stress at each ply interface (Cho and Averill, 2000; DiSciuva,
1997). The shear stress continuity conditions at the k-th interface can be

expressed as:

Sr? = Sr”. Sr.“ = $3“) (4.5)
Assuming that infinitesimal strain theory holds for the transverse shear strain
components, the form of the assumed displacement fields in equation 4.5 is such

that the in-plane displacements 14:") and ug") contribute only constant terms to

70

the transverse shear stresses, Sf? and 5;). For consistency, then, the

transverse displacement u§“ should also contribute only constant terms to these

stress components. In order to achieve this consistency, terms in the transverse
shear strain expressions that involve the transverse displacement variables are
evaluated at the midplane of the sublaminate, thereby taking a thickness

averaged value of these terms. Substituting equations 4.3 and 4.4 into equation
4.5 and solving for 5, and 77,, it can be seen that 5, and 77, depend on the ratios

of shear properties between adjacent layers and on the shear deformation in
each sublaminate. If either of these quantities are small, then discrete-layer

effects will also be small.

The rotational variables ‘1’, and ‘1’, in the displacement field are now
eliminated by introducing the variables u, and v, (the in-plane translations at the

top surface of the sublaminate) in order to expedite the development of versatile
finite element models. Thus, rather than describing the in-plane displacement
field by a translation and rotation at one point, it can more conveniently be
described here by the translation at two points.

The displacement fields must be further manipulated in order to develop a

Go finite element. Because the derivatives of transverse deflections

aw” aw' aw” and Efl) appear in the displacement field, their second

(511’5x1’5x2 5x2

 

derivatives will be present in the strain energy functional. C1 continuity of w, and

w, is thus required. It is desirable to alleviate such a requirement by introducing

new rotational degrees of freedom as follows:

71

.6_”)L=_62b,_al=_62“_=61b’ - Ir
6x, 6x, 6x, 6x,

 

(4.6)

These relationships will be constrained in the strain energy via a combination of
the interdependent interpolation and the penalty method. The displacement fields

can now be written using indicial notation as follows:
11,“) = flap-(1)22, a?) = 17¢ -‘Pg;,’, ug") = W” 42,, a =1...4; ,6 = 1,2 (4.7)
where the index ,6 is used to denote the top and bottom of the sublaminate with

1: bottom and 2: top, and (1)23)ng and 0,9 are shape functions in the thickness

direction (see Cho and Averill, 2000). Unless noted otherwise, summation on
repeated indices is implied. The variables 5a,, are:

17,,=u,,, iizfl=v,,, a”,,,=6,,,, 543:6“, (4.8)

The transverse normal strain arising from the assumed displacement field

(4.7) is constant through-the-thickness of a sublaminate and takes the value

(w, — w,)/h. This assumption may generate substantial errors in composites

that have a soft core or layer, such as in sandwich panels, laminates containing
damage, or adhesive layers. It is desirable to improve the distribution of the

transverse normal strain through-the-thickness. The improvement can be

achieved by assuming a constant transverse normal stress, 3,, through-the-

thickness of a sublaminate. 3'3, can be determined using Reissner's mixed

variational principle (Reissner, 1986). From equation 4.2, the transverse normal

strain is obtained as:

72

1 _.
—(k) __ (k) (k) (k) (k) (k) (k)
8 (S33 - C13 811 _ C23 82 — C36 712 ) (4'9)

33 — (k) 2
C33

5,, is then determined analytically from the following relation:
N ’30) k w _w
e: nan—.- b)... (4.10)
r.-. , ’1
art-I)
The constant transverse normal stress 3,, is found to be:
3,, =f,(z7a,,w,,<b;;’,\P;;’,Q;j;) (4.11)
Now, the newly defined transverse normal strain can be written as:
5;,” = f,(17a,,w,,<1>;j;,\11,§,’;’,(2;j;) (4.12)
The functions f, and f2 are defined in (Cho and Averill, 2000) for small

displacement analysis. For the purpose of calculating S”, the von Karrnan

nonlinear strain terms must be included in equation 4.9 to ensure accuracy of the

solution for all types of laminates.

4.3. Virtual Work

The equations of motion, the essential and natural boundary conditions,
as well as the displacement-based finite element model can be developed from
the principle of virtual work. Again, attention will be limited to a single
sublaminate. The work functional is modified here to include the imposition of the
constraints in Equation 4.6 via the penalty method. The total work performed by

internal and external forces on a sublaminate is:

W = W, + W, (4.13)

73

For linear elastic material behavior, the internal work is:

1
W, = 3 LSE dV (4.14)

where E and S are vector representations of the Green's strain and the second
Piola-Kirchhoff stress states. WE is the work of the external forces performed on

the sublaminate:
W, = — jb’udV — It’udl‘ (4.15)

where u is a vector of displacements, b is a vector of body forces and t is a vector

of tractions on the surface I‘.
The work functional defined by equation 4.13 will be subject to the four

constraints in equation 4.6, which can be represented as:
G, (u,u') = O i = l,2,3,4 (4.16)
We now introduce a modified work functional:
WM =W, +W,, +W,, (4.17)

If the constraints (6) are enforced using the penalty method, then:

W, = fill—£— LGfdQ (4.18)

where y is the penalty parameter, Q is the domain upon which the constraints

are to be enforced, 01:921...ga , G2 =a—W'—+6I,, , G, =-%-6,, , and
6x, 8x, 6x,
G4 = 9—“;- 1:
6x,

74

The total virtual work performed on a sublaminate is now defined by the

first infinitesimal variation of the modified work functional:
5W” =6W,+5W,+6'W,,=0 (4.19)
6W. is the virtual work resulting from internal forces:
6W, = (S 5EdV (4.20)
8W5 is the virtual work resulting from external forces:
W, = - [11’de — [tré'udl‘ (4.21)

6Wp is the virtual work resulting from violation of the equality constraints:

W, = fir j[6G,(u,u')G,(u,u')]dV (4.22)

i=1
As the penalty parameter 7 becomes larger, violation of the equality constraints is
minimized in a least squares sense. The governing equations, finite element
model and boundary conditions can be derived directly from equation 4.19. The

details of this procedure are omitted for brevity.

4.4 Finite Element Formulation
The standard form of the finite element approximation is:
u = NH (4.23)
where N is a matrix of element interpolation functions and i is a vector of nodal
degrees of freedom.

Green’s strain can be written as the sum of linear and nonlinear parts:

E = a, + 8”, (4.24)

75

The linear portion of the Green’s strain is:

e, =Lu =LNii=B,fi (4.25)
where L is a matrix of linear differential operators and BL is a matrix containing
gradients of the interpolation functions. All coefficients needed to formulate BL
can be found in (Cho and Averill, 2000). The nonlinear portion of the Green’s
strain is:

5,, =—;-A(l9)l9(u) (4.26)
where 6(a) is a vector of displacement gradients and A(6) is a matrix composed

of the vector 0(u). Green’s strain can then be written as:
E = ms + %A(6)i9(fi) = 1311' (4.27)
and the first infinitesimal variation of Green’s strain is:
an =B,6i+%(6A(0)6(fi)+A(6)5r9(ii)) (4.28)

To simplify the definition of the first variation of Green's strain we note that:

6E .—. B,6i+%(A(66)6(fi)+A(6)66(fi)) (4.29)
and:
619(6): amen (4.30)
so that:
as = 13,511 + A(6)6(Ti)6i (4.31)

The first variation of the Green's strain can now be written as:

76

6E =B,6fi+A(6)6(fi)6fi= (B, +A(6)G)éi =§aa (4.32)
where G is a matrix of nodal shape function gradients.

The virtual work can now be written as:
JWM = air-’(jETSdV + 7 Lp’pfido — [b’NdV — lt’Ndl“)= 0 (4.33)

where p is a vector containing nodal shape functions used to define the penalty
constraints. This relation is valid for an arbitrary variation of the element nodal

degrees of freedom, thus:
r(rr) = [ETSdV + 7 Lp’pfido — [b’NdV — jt’Ndr = 0 (4.34)
In general, an exact solution to equation 4.34 will not be attained, so the residual
r(ii) will be non-zero. The residual is often called the out-of-balance forcing
vector since it represents the amount that the current solution is “out of
equilibrium.”
Equation 4.34 is typically written in the form:
D Ki = f (4.35)
where the direct stiffness matrix (including the penalty term) is written as:
D K = [E’DfidV + y Lp’pdn (4.36)
and the forcing vector is written as:
r = b’NdV + [€th (4.37)

It should be noted that direct iteration of equation 4.35 can be performed
by constructing, inverting and iterating upon the unsymmetric direct stiffness
matrix, but this leads to high computational costs and slow convergence rates.

Instead, we employ the Newton-Raphson solution technique that uses a

77

symmetric tangent stiffness matrix to iterate until a sufficiently accurate solution

is obtained.

4.5 Newton-Raphson Method
To define the tangent stiffness matrix TK, we expand the residual in a
Taylor series expansion about the current equilibrium configuration (denoted by

the subscript o):

 

r(fi)=ro+aE§fi+i§fi’§Léi+-n (4.38)
au 2! er

Higher order terms in the Taylor Series expansion are assumed to be small when
compared to the constant and linear terms. We wish to minimize the residual

r(ii) so equation 4.38 becomes:

 

0 = r, + are: an (4.39)
an

The first variation of nodal degrees of freedom about the current configuration is:

arc ‘1 __r '1
(Si—{afi} r,— (K) r, (4.40)

 

For complex formulations such as the current one, direct application of equation
4.40 can be cumbersome. Alternatively, but equivalently, the tangent stiffness
can be formed by taking the second infinitesimal variation of the modified

functional (see Crisfield, 1991 ). The tangent stiffness can then be defined as:

 

6(5WM)= 55(11): Wiggifi -_- afl’xbn (4.41)

We can express the first infinitesimal variation of the residual roar“) as:

78

(ma) = 611’ { j (éE’S+E’éB)dV+ 7 ((p’pboaa} (4.42)
The first term in (42) can be written as:
an’ [557$ch = 511’ j(5A(e)G)TSdV (4.43)
which can be rearranged and manipulated (due to the symmetry of the second
Piola-Kirchhoff stress tensor) to form:

613’ [G’(5A(e))TSdV = an’ jG’scsA(e)dV (4.44)

A

where S is a matrix of second Piola-Kirchoff stresses which gives rise to the
stress stiffening matrix “K:

an’ jo’éaA(6)dV = as’ (IG’SGdVbn = afi’(°K)aa (4.45)
The stress stiffening matrix represents stmctural stiffening effects due to stresses

that arise during the nonlinear deformation.

The second term in equation 4.42 can be simplified as:
511’ jfi’asarV = an’ [13’0de = 5117([fircfidr/ba = aa’(f<‘)au (4.40)

The stiffness matrix K represents the combined stiffening effects from large and
small displacements.

The last term in equation 4.42 can be represented as:
aa’y(Lp’pdo)aa = éfi’(”K)éi (4.47)

The penalty stiffness matrix P K represents the stiffness required to enforce the
linear equality constraints in a least square sense.

Finally, the tangent stiffness matrix is:

’K="K + K+PK (4.48)

79

The solution at the i-th iteration of the n-th time increment is obtained

using:
[’Kl? {Mir-1. = {F If - [”Kl {i}? (4.49)
and:
his. = Atfilr. + {i}: (4.50)

where {ii},';, and A{ii'}" are the total and incremental displacements at iteration

i+l
i+1 for the n-th time increment. The equations 4.49 and 4.50 are solved iteratively

until A{fi},’.',, is sufficiently small.
4.6 Finite Element Approximations

The first-order zig-zag sublaminate finite element model formulation for
linear strain theory was developed by Cho and Averill (2000). These elements
are based on a formulation that involves sublaminate surface variables, which
facilitates the development of a three-dimensional finite element model. Each
node contains five degrees of freedom — three translations and two rotations.
Both an eight-noded brick and a six-noded wedge finite element have been
developed (see Figure 4.2). For the geometrically nonlinear model formulated
here, a similar approach is used. For brevity, only the eight-noded brick element
is discussed.

The in-plane displacement and rotational degrees of freedom are

approximated by bilinear Lagrange interpolation functions:

Ma
M:-

u, = u,,N. u =

I 1

EN. (4.51)

(I t

.0.

'l-
u.-

80

v, = fie-MN. v = firm, (4.52)

19,, = £67,,N. 19 = i6 .N. (4.53)

a =fijéy,,N, 19,,=_ 5,,N, (4.54)

where N. are bilinear Lagrange interpolation functions and n is the number of
nodes that define either the top or the bottom surface of the element. For the
eight-noded brick element, n=4.

An interdependent interpolation scheme is used for the transverse
deflection degrees of freedom. The interdependent interpolation scheme reduces
the effects of locking and facilitates the satisfaction of the constraints in equation
4.6. The interdependent interpolation concept was first developed by Tessler
and Hughes (1983). This scheme does not require additional degrees of freedom
and it improves element accuracy and efficiency. The interdependent
interpolation scheme couples the transverse deflection degrees of freedom with

the rotational degrees of freedom, as follows:

w, = :(Nfifi + 117,02, + 10,5”) (4.55)

is]
where N. are the bilinear Lagrange interpolation functions, 117, are special

quadratic functions, ,6 ranges from 1 to 2 (1= bottom, 2=top). and n=4 is the

number of nodes on the top or bottom surface.
While the interdependent interpolation scheme alleviates some of the
shear locking in the element, it does not eliminate it completely. A locking-free

element is obtained by enforcing both field and edge consistency conditions, as

81

l0

9E

re

(0

 

shown by Prathap and Somashekar (1988). These concepts have been applied
to the current element (Cho and Averill, 2000) to make it robust for application to
both very thin and very thick laminated plates.

Appendix A defines all matrices needed to form the direct and tangent
stiffness matrices for the geometrically nonlinear first-order zig-zag sublaminate
finite element model accounting for moderately large displacements and
moderate rotations. The current geometrically nonlinear finite element model has
been implemented as a user element subroutine within the commercial finite
element code ABAQUS (Hibbit, et al., 1997). This approach takes advantage of
the efficient iterative solution algorithm native to ABAQUS and allows attention to
be focused on development and testing of the new nonlinear zig-zag sublaminate
theory and its associated finite element model. All numerical results in the
following section were obtained using the present model as a user element

subroutine within ABAQUS.

4.6 Numerical Results

In this section, sample numerical results obtained from the first-order
geometrically nonlinear zig-zag model are compared with available experimental
results (see Zaghloul and Kennedy, 1975; cited in Ochoa and Reddy, 1992) for
bending of square, symmetrically laminated, simply-supported and clamped
composite panels subjected to a uniform transverse pressure load. Since the
geometrically nonlinear response of composite plates is highly dependent upon

the lamination scheme and boundary conditions, it is essential to faithfully

82

r“:

 

 

lar
Tt
tol
me
ple
figr

hrs

19}

Gun

Silua

represent the lamination scheme and boundary conditions in the experiment if
numerical and experimental results are to be compared. The numerical
replication of the experimental setup can often be accomplished to a satisfactory
level of accuracy, but usually cannot be replicated exactly. In order to illustrate
the importance of zig—zag effects in nonlinear analysis of laminated stnrctures,
numerical results are also presented for bending of a square, symmetrically
laminated, simply-supported sandwich composite panel subjected to a uniform
transverse pressure load.

The first laminate studied is a simply supported, single-layer orthotropic
square plate with in-plane dimensions of 12 by 12 inches and a thickness of
0.138 inches. Material properties of this panel are given in Table 4.1. The second
laminate studied is a clamped square cross-ply laminate with a (0/90/90/0) layup.
This panel has the same in-plane dimensions as the first example, but it has a
total thickness of 0.096 inches. Table 4.2 contains the stiffness properties of the
material used in this second panel. Due to symmetry of the panels, only a quarter
plate model was developed, with meshes in the quarter region as given in the
figures (e.g., 4x4x1). Along the simply-supported edges of the plate model in the
first example, simply-supported boundary conditions were imposed at the bottom
of the laminate to best represent the experimental setup (Zaghloul and Kennedy,
1975)

Figures 4.3 and 4.4 contain the experimental and predicted load-deflection
curves for the square simply supported single-layer orthotropic plate and the

square clamped cross-ply laminate (0/90/90/0), respectively. Numerical results

83

 

dr':

arr

are presented for both the first-order zig-zag sublaminate finite element model
and the S4R shell element available in ABAQUS. For the first-order zig-zag
model predictions, results are given for various levels of in-plane and through-
thickness mesh refinement. Both plates stiffen significantly due to the coupling of
bending and membrane stresses. Figures 4.3 and 4.4 show excellent correlation
between results obtained from the geometrically nonlinear first-order zig-zag
sublaminate model, numerical results using ABAQUS, and the experimental
results. The small differences in the solution predicted by the S4R element and
the first-order zig-zag element can most likely be attributed to the slightly different
boundary conditions utilized. In the S4R element, boundary displacement
constraints are imposed at the plate mid-surface, as is common practice for shell-
type finite elements. On the other hand, boundary displacement constraints in the
three-dimensional first-order zig-zag element are imposed at the bottom surface
for the simply-supported case and at both top and bottom surfaces for the
clamped case. Based on the results using the first-order zig-zag sublaminate
model, it is clear that through-thickness refinement does not improve the
predictions of global response for either of these panels. The correlation between
the two numerical solutions is better than the correlation between either of the
numerical solutions and the experimental results. For both cases considered,
discrepancies between the computational results and the experimental results
are likely due to differences in boundary conditions and material properties.

The final numerical example was constructed to compare the predicted

through-thickness variation of displacements and stresses using the

 

96

SU

th
th
st

SL

82
zig
so
thr
tra
lei

BC

geometrically nonlinear zig-zag model. The laminate considered here is a simply
supported sandwich panel with a compliant core, which has been studied
previously (Cho and Averill, 2000). The plate has in-plane dimensions 10 x 10
inches and a total thickness of 1 inch, so the aspect ratio of the plate is ten. The
plate is loaded by a uniform pressure applied on 4 x 4 inch square patch at the
center of the plate. Tables 4.3 and 4.4 contain the material properties and
lamination scheme, respectively, of the sandwich panel. Due to symmetry, only a
quarter of the plate was modeled with a uniform 5x5 finite element mesh in the
plane of the laminate and with a varying number of elements through-the-
thickness. For example, when three elements are used through-the-thickness of
the plate the total thicknesses of the bottom face sheet, the core and the top face
sheet were each represented by one element. As in the first example, simple
supports were imposed along the bottom edge of the panel. The results of the
present model are compared to predictions of the ABAQUS layered solid element
denoted as C3D8I, which uses incompatible modes to overcome locking.

Figure 4.5 displays the predicted load-deflection curves (at center of the
sandwich panel) obtained by the linear and the geometrically nonlinear first-order
zig-zag element and the ABAQUS CSD8I element at the top of the simply
supported sandwich panel. It can be seen that the ABAQUS element requires
three elements through-the-thickness to achieve a converged solution for the
transverse deflection at the top of the panel, while the first-order zig-zag model
requires only a single element through-the-thickness to attain a similar level of

accuracy.

85

 

 

38

ll)

SE

T2
for
of

inc

Figure 4.6 displays the predicted load-deflection curves (at center of the
sandwich panel) obtained by the linear and the geometrically nonlinear first-order
zig-zag element and the ABAQUS C3DBI element at the bottom of the simply
supported sandwich panel. It can be seen that the ABAQUS element requires
three elements through-the-thickness to achieve a converged solution for the
transverse deflection. The deflection at the bottom of the panel decreases as the
number of first-order zig-zag elements is increased through-the—thickness. Figure
4.5 depicts very little change in deflection at the top of the panel as the number of
first-order zig-zag elements is increased through-the-thickness. Together,
Figures 4.5 and 4.6 show an increase in nominal transverse normal strain as the
number of first-order zig-zag elements is increased through-the-thickness. This
occurs due to the compliant core of the sandwich panel. For this problem, three
first-order zig-zag elements through-the-thickness are required to capture the
transverse “squashing" of the core. To determine if the core is “squashing” due to
geometrically nonlinear through-thickness refinement effects, a set of small
displacement analyses were conducted with first-order zig-zag and ABAQUS
C3D8l elements. Table 4.5 contains the deflections at the bottom and top of the
sandwich panel assuming small displacements using the first-order zig-zag and
ABAQUS elements with one, two and three elements through-the-thickness.
Table 4.5 shows the deflection of the top of the panel remaining nearly constant
for the first-order zigzag finite element model, while the deflection at the bottom
of the panel decreases as the number of elements through-the-thickness is

increased. The “squashing” of the core is captured as the number of elements

86

 

thro

88$

thrr
dis

suf

the
cer
thir
anr
pla
in-r

agr

inte
tran
finitr
AB/

bolt

ABA

Withi

through-the-thickness is increased, even though small displacements are
assumed.

For an applied load of 9.5 ksi, Figures 4.7 and 4.8 display the predicted
through-thickness variation of in-plane displacements. The in-plane
displacements from ABAQUS and the first-order zig-zag element agree with
sufficient refinement (three elements through-the-thickness).

Figures 4.9 and 4.10 display the predicted through-thickness variation of
the in-plane Cauchy bending stresses at the nearest integration points to the
center of the plate for an applied load of 9.5 ksi for various levels of through-
thickness refinement. For clarity, Figures 4.11 and 4.12 compare the ABAQUS
and first-order zig-zag model prediction of through-thickness variation of the in-
plane Cauchy bending stresses within the top “skin” of the sandwich panel. The
in-plane bending stresses from ABAQUS and the first-order zig-zag element
agree with sufficient refinement (three elements through-the-thickness).

Figure 4.13 displays the transverse normal Cauchy stress at the
integration points closest to the center of the plate. Due to the assumed constant
transverse normal Cauchy stress, each sublaminate of the first-order zig-zag
finite models will have a constant through-thickness normal Cauchy stress. The
ABAQUS and first-order zig-zag predictions differ significantly in the top and
bottom “skins” of the sandwich panel in Figure 4.13. This can be accounted for
by inspecting the kinematic assumptions of each finite element model. The
ABAQUS finite element model assumes that the displacements vary lineariy

within an element. This enforces the strains to be continuous at each layer

87

interface within an element. This implies that the stresses are discontinuous at
each layer interface due the abrupt changes in material properties within a single
ABAQUS element. Conversely, the first-order zig-zag element assumes that the
stresses are continuous at each layer interface within a sublaminate. For this
reason, the first-order zig-zag element is expected to yield higher accuracy when
compared to the simple linear ABAQUS element.

Figure 4.14 displays the transverse normal strain at the integration points
closest to the center of the plate. Due to the assumed constant transverse
normal Cauchy stress, the transverse normal strains vary linearly within each
layer of the first-order zig-zag finite model. As the number of sublaminates (or
elements) is increased through-the-thickness, a larger amount of nominal
transverse strain is predicted. Figure 4.14 depicts the previously discussed
transverse “squashing” of the core conveyed in Figures 4.5 and 4.6.

Figures 4.15 and 4.16 display the predicted through-thickness variation of
the transverse Cauchy shear stresses for the first-order zig-zag model and
simple linear ABAQUS element for various levels of transverse refinement. The
difference in kinematic assumptions for the first-order zig-zag and linear
ABAQUS finite element explains the discrepancies in the through-thickness
variation of the transverse Cauchy shear stresses within the bottom and top
“skins” of the panel. The linear ABAQUS element assumes a constant state of
strain within an element, thus the stresses have double values at each layer

interface within an element. The first-order zig-zag finite element model assumes

88

PE

de

a piecewise linear variation of in-plane displacements through-the-thickness of
the sublaminate while satisfying continuity of inter-laminar stresses.

Figure 4.17 depicts the predicted through-thickness variation of the in-
plane Cauchy shear stress for the first-order zig-zag model and simple linear
ABAQUS element for various levels of transverse refinement. Both elements
predict similar through-thickness variation of in-plane Cauchy shear stress.

Differences between the first-order zig-zag and linear ABAQUS solutions
can be attributed to the different kinematic assumptions used in the models, with
the zig-zag element being expected to yield higher accuracy than the simple
linear model employed in the CSD8I element. For clarity, Figure 4.18 depicts the

measurement points of the in-plane displacements, stresses and strains.

4.8 Conclusions

A new geometrically nonlinear first-order zig-zag sublaminate plate theory
and finite element model have been formulated, implemented and validated
against existing experimental and numerical results for laminated composite
panels. This model is accurate for a wide range of laminated panel applications
including very thin laminated panels up to very thick laminated and sandwich
panels. Further, this model is cast in a form that facilitates its use within
commercial finite element packages. Namely, the finite element models
developed take the form of a three-dimensional eight-noded brick or six-noded

wedge element with the usual five engineering degrees-of-freedom per node.

89

 

 

 

 

 

 

 

X2

Figure 4.1. Schematic of ply stacking and coordinate system used in first order

sublaminate theory.

 

 

 

   

 

 

Figure 4.2. Element topology for the eight-noded brick and the six-noded wedge finite

elements based on first order zig-zag sublaminate theory.

90

 

:04:ch .-
1...- u \ ‘4‘
Ian!“

 

Figur

plate

refine

 

0.5 ,

 

+ Nonlinear ABAQUS 4x4 S4R "
—x— Nonlinear ngg 4X4X1
+ Nonlinear Zgzag 8x8x1
-e— Nonlinear Zgzag 4x4x2
-3- Linear Zgzag 4X4X1

—1— Linear ABAQUS 34R 4x4

.0
.5

.°
00

 

 

 

\\\
\h

Center Deflection (inches)
hr

0
.5

 

 

Uniform Pressure (psi)

Figure 4.3. Load versus deflection for a square simply-supported single-layer orthotropic
plate under uniform pressure at various levels of in-plane and through-thickness

refinement.

91

 

nu CCU
Ecut‘G-fidc .- t

t!‘

’E

Fig

und

1.2 _

 

+Experiman1ar
+NonlinearABAQUS4x4 84R
,4 +Nonrinearzrgzag4x4x1 /"
+NGWWZIMW ‘r.
+Nonlineafigm4x4xZ /
+1.1an 4X4X1 /
.9.UnearAeAQUS4x4$4R

(cg!)

 

 

 

O
O
n\
,.

Center Deflection

\
I
I

0.2 4 ' /

 

 

0 2000 4000 1011!) 12000 14000

0000 0000
Uniform Preemie (Pa)

Figure 4.4 Load-deflection curves for the clamped square cross ply plate (0/90/90/0)

under uniform pressure at various levels of in-plane and through-thickness refinement.

92

 

1.5

+meaarpeoaseea °
1.35 +mm555¢
+mmssea
+ueaagzagsea

+erirm'2lgm56d c /

+Nrira2igm552

.A
to

lg

\\

p
9

\.

+NrimZIgzag566

’1
'

 

 

 

Deflectkén at Top of Pgnei at Center (inches)
8 a1

8

o
if

 

cummmmmmmummammmmammamam
moi)

 

Figure 4.5. Linear and nonlinear load versus deflection at the top of a simply-supported
square sandwich panel under uniform pressure using first-order zigzag sublaminate

theory.

93

1.5 .

 

 

 

 

1.35. +NoriirmABM11$5xSx1
+Wmm5x5fl

g 12. +NorirnearAeruseec
g +Umer2lsmo5fix1
7:105, a-hbihmllgmfiafiod
: 0‘91 +Noriineabgm516x3 2
E /
E 075. /'j
'6 " f
g 0.0. /2/
§ /
I
c 0.45 ’/
a :-
0
e e. /
8 _/

(115

 

0

 

omo1an1anmmmmmmmsanmeemmmammmm
Preeeuupsl)

Figure 4.6. Linear and nonlinear load versus deflection at the bottom of a simply-
supported square sandwich panel under uniform pressure using first-order zigzag

sublaminate theory.

94

 

 

 

 

 

 

 

+ hbnlinea ABAQB 51516
+ Nonlinea Zgzag 5x5x1
0'75” +l\bnlinea Zgzag 5x56
8
S
0
.E
‘g' 05 .
c
x
O
E
.—
025 .
0 c - / I . . fl .
0.0012100 250502 500502 7.50502 1.00501 125501
Totd X1-dsplacement (inches)

Figure 4.7. Through-thickness variation of total Xl-displacement at the plate midspan
and X1=0.

95

   
  
  

 

 

 

 

 

 

 

—e— Nonlinear ABAQUS 5x5x3
+ Nonlinear Zigzag 5)<5x1
0.75 . .
+ Nonlinear Zgzag 5x516
‘1?
e
.1:
o
.E
g 0.5 .
1:
x
.2
a
.—
025 .
o . . . . .
CIDER!) 2.50602 SIDE-02 7.5(E-02 1.00501 1.25E-01

Total X2-dsplacement (Inches)

Figure 4.8. Through-thickness variation of total X2-displacement at the plate midspan
and X2=0.

96

 

flung—CE: mama—10.2F

Figure

at the 1

 

1111) .
03D 4 4a.

00m 4

 

+ Nonlinear ABAQUS 5x5)<3
+ Nonlinear Zgzag 5x5x1
+ Nonlinw Zigzag 5x56

:3

3

 

 

 

Thickness (inches)
; E
k

43

0.21) .

0.100 . ¢
0 ”m... a
0.000 , .fi . .

-4.0E+05 -ZOE+05 0.06%!) 2% 4.08% 308% &OE+{5
Nonml Cauchy Stress 0,, (psi)

 

 

 

 

 

 

Figure 4.9. Through-thiclmess variation of the normal Cauchy stress (in the Xl-direction)

at the Gauss points nearest to the center of the plate.

97

 

 

 

 

 

 

 

 

 

11“).

0900. %

03!).

0700 +leineaABAQBW
‘é‘ +hbnlineaZgzag5x5x1
.5 0000. .

.E +Nodlnea‘Zgzag5yéxa
£0.51) .

.5

2 0400.

u:

'—

0300.

0.200.

0.100.

1 film-
0000 m

 

.20905 4.0905 00900 1.0905 2.0905 3.0905 4.0905 50905 a0905
WCanohySlreesmw)

Figure 4.10. Through-thickness variation of the normal Cauchy stress (in the X2-
direction) at the Gauss points nearest to the center of the plate.

98

11!!) -

 

 

 

 

 

 

 

 

 

 

 

0975. . _
+inlnea'l-‘BAQB566
'5 +MriineaZigm506x1
.:
2 +Mflinea2igm566
gm _
c
8
2
1:
p.
0925.
0900 . . . a . . . . . .
5905 4.905 a905 .2905 .1905 0900 1905 2905 3.905 4.905 5905
mmmmm

Figure 4.11. Through-thickness variation of the normal Cauchy stress (in the X1-
direction) for the top five layers of the sandwich panel at the Gauss points nearest to the

center of the plate.

99

 

 

 

 

 

 

 

 

 

 

1.000.
+ Nonlinea‘ABAQUS 5x56
+ Nonlinea Zgzag 5x5x1
0.975 . + Nariinea' 29:39 5M
8
.C
c
.5
3 0950
o
c
x
.2
4:
I-
0.925.
0.91) . . r . 9 . . .
4.75905 4.30905 4.00905 025904 -250904 1.2904 500904 075904 1.25905
Mom-I Cauchy are-s 022W)

Figure 4.12. Through-thickness variation of the normal Cauchy stress (in the X2-

direction) for the top five layers of the sandwich panel at the Gauss points nearest to the

center of the plate.

100

1.000 4

 

 

 

1)
0

0.9!) 4

 

m- +Non|inearABAQUS§x§>G
0700. +Nonlinear Zgzag 5x$x1
+ Nonlinear Zigzag 5xS)<3

T5

 

 

 

Thickness (inches)
:5 5
0

3

OLZIH

0100.

0.000 . , . . . k‘ . .
00904 50904 40904 00904 .20904 4.0904 0.0900 1.0904 20904
Nonml CauchyStressmmsl)

 

 

 

 

 

Figure 4.13. Through-thickness variation of the normal Cauchy stress (in the X3-
direction) at the Gauss points nearest to the center of the plate.

101

 

 

 

 

 

 

 

 

 

1.0 . j
09 . T
0.8 .
0.7
3 + NodinearABAQJS (5)156)
0.0 .
3.5 ... Noriinea 39m (5x5x1)
g 05 A 0 0
g 04 ... Nonlinear Zgzag (511510)
E . 4
0.3 .
0.2 .
0.1 . P

 

 

4.00501 4.00901 4.40501 4.20501 4.00901 43.00902 41.001502 4.00502 200602 000900 200602
Tummardnlgu)

Figure 4.14. Through-thickness variation of transverse normal strain (833) at the Gauss

points nearest to the center of the plate.

102

g

 

 

 

 

 

 

 

 

 

 

 

 

0900 : \4
081) 4
0700 .
A
8 + hbriineaABAQJS 5510
.5 0000.
c O
a +MrimeaZgzm5>6x1
3 “5‘” .
2 +l‘bflmeaZgzag 566
x
2 0.41).
S
.—
0300 .
02D .
0100 .
04m 1 I r fi fi
ACBOS -2(E+03 00900 219m 409m GAE-'03 00900

TmCaflySuaressMM)

Figure 4.15. Through-thickness variation of the transverse Cauchy stress (023) at the

Gauss points nearest to the center of the plate.

103

 

 

 

 

 

 

 

 

 

 

 

1.000 . ' h:
v —=———.
0000. J 7%
0&1).
0700 .
‘1?
c +hbdineaABAQJ857676
.5 0W1
: .
c. +hbrl1neaZgzm516x1
0500 0 u ,
+hbri1nea'Zgzag57676
C
.K
2 0400 .
4:
'—
0300.
02134
01(1) . 4
0.000 - :EE‘

 

00900 00900 405100 -2(E+03 00900 20900 40900 00900 00900
Tm Cauchy Shea m 0,, (psi)

Figure 4.16. Through-thickness variation of the transverse Cauchy stress (013) at the

Gauss points nearest to the center of the plate.

104

 

 

 

i,
i

 

 

 

 

 

 

0000.
0700 _
g +NodineaABAQUS§5fi
o 0000. .
é 4.1110111110412929 55x1
0500fi .
g +hb111neaZgzag 5766
.l
.2 am <
5
l-
0000.
0200.
0100
0000 E. , . . .
0.0900 1.(E+04 2.15004 312104 41604 5.1504
ln-leeCauchyShea'm 012 (psi)

Figure 4.17. Through-thickness variation of the in-plane Cauchy stress (cm) at the Gauss

points nearest to the center of the plate.

105

Point of measurement for
displacement (ul) X2

b

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

\ Symmetry Point of measurement for
‘ ’ strains and stresses
I A
K
Simply supported Symmetry
X1
V
4 _ > ’
Srmply supported Point of measurement for
displacement (uz)

Figure 4.18. Measurement points for the displacements, strains, and stresses for the

sandwich panel.

106

Table 4.1. Material properties for square single-layer orthotropic plate.

 

 

 

 

 

E11=3.00E6 ps1 E22=1.30E6 ps1 E33=3.00E6 psr
v12=0.32 vz3=0.32 V13=O'32
G12=0.36E6 psr G23=O.36E6 psi G13=O.36E6 psr

 

 

 

Table 4.2. Material properties for square cross-ply laminated plate (0/90/90/0).

 

 

 

 

 

E, 1=l .82E6 psr E22=1.82E6 psr E33=l .82E6 psr
V12=0'2395 v23=0.2395 V13=O'2395
G12: 0.31E6 pSl G23=0.31E6 pSl G13=O.31E6 pSi

 

 

 

Table 4.3. Material properties for sandwich panel.

 

 

 

 

 

 

 

 

 

 

Material 1 Material 2 Material 3 Core

We) 1 33 25 0.05
E2200" psi) 25 21 1 0.15
E33 (H35 psi) 1 21 1 0.05
0,. 0.01 0.25 0.25 0.01

023 0.25 0.25 0.25 0.15

v ,3 0.25 0.25 0.25 0.15

G 12(106128?) 0.5 8 0.5 0.02
(323,006 psi) 0.5 8 0.5 0.04
G 13(ID6 psi) 0.2 8 0.5 0.02

 

 

 

 

 

 

 

Table 4.4. Lamination scheme for sandwich panel.
No.

 

107

Table 4.5. The effect of transverse refinement considering small displacements.

 

 

 

 

Number of ABAQUS: ABAQUS: FZZ: FZZ:
elements Displacement Displacement Displacement Displaceme
at bottom at top at bottom nt at top
through- (inches) (inches) (inches) (inches)
thickness
1 0.3950 0.3978 1.422 1.484
2 0.9757 1.027 1.425 1.496
3 1.317 1.401 1.375 1.449

 

 

 

 

 

 

 

108

 

Chapter 5

OPTIMIZATION OF SANDWICH PANELS USING GENETIC ALGORITHMS

5.1 Introduction

Optimization problems are encountered on an hourly basis in the everyday
life of an engineer. For example, we will review a typical optimization problem
posed by a naval engineer. The goal of the optimization problem is to find a
sandwich panel design with minimal mass while simultaneously satisfying
stiffness and stress constraints applying state equations that are governed by the
geometrically nonlinear (moderately large displacements and moderate rotations)
equations of motion. The sandwich panel designs found with a simple Genetic
Algorithm will be contrasted to an existing “baseline” sandwich panel design in

terms of its design variables, objectives and constraints.

5.2 Sandwich Panel Optimization Problem Statement

All optimization problems can be defined as finding a set of bounded
design variables to minimize objective(s) while satisfying all constraints. We will
follow in similar fashion by describing the design variables, the objective and
constraints for the sandwich panel optimization problem statement.

Typically, sandwich panel designs are constructed of layered composite
materials on the top and bottom of the panel (referred to as the bottom and top
skin). The center of most sandwich panels is composed of lighter materials

(often called the core). The skins of a classic sandwich panel design are meant

109

to carry the bulk of the bending stresses while the light weight core experiences
most of the shear stresses. With this in mind, the simple GA was allowed to vary
the number of layers in the bottom skin, top skin and core of the sandwich panel.
The bottom and top skins of the sandwich panel were independently required to
each have at least two, but no more than eight plies. The core of the sandwich
panel was required to have at least one, but no more than three layers. The
bounds on the number of layers in the bottom skin, core and top skin of the
sandwich panel were chosen to satisfy manufacturability requirements. The GA
could change each ply orientation ranging from O to 180 degrees with a 15
degree increment. The choice of ply orientation was discretized in such a manner
to meaningfully reflect a realistic manufacturing increment. The GA chose each
ply material from a discrete set of materials given in Table 5.1. For each layer in
the core, the GA could change the core layer thickness ranging from 1/4 to 1%
inches with a ‘/4 inch increment. The thickness of each layer within the core was
discretized to reflect a meaningful manufacturing increment. The GA chose each
core material from a discrete set of materials given in Table 5.2.

The representation of the design variables along the chromosome was set
up in an order that staggered the number of plies in the top, core layers and
number of plies in the bottom with material choice, orientation direction and
thickness to reflect their physical location on the sandwich panel. This choice
was made according to the principal of keeping related physical entities close to
each other on the chromosome, which has been known to aid search within a GA

(see Holland, 1975, Goldberg, 1988). The first discrete variable on the

110

chromosome was the number of plies on the bottom skin of the panel. Next, a
discrete ply material choice and orientation choice was staggered on the
chromosome for each bottom ply. The discrete number of layers in the core was
followed by a core material choice for each possible layer in the core. The next
discrete variable on the chromosome was the number of plies on the top skin of
the panel. A discrete ply material choice and orientation choice was then
staggered on the chromosome for each top ply. The simple GA had a population
size of 50, 70% crossover rate and 1% mutation rate per field.

The objective of the sandwich panel optimization problem was to minimize
mass. Inequality constraints were placed on the maximum transverse
displacement and on the maximum value for a set of failure indices. The first set
of failure indices was defined as the ratio of the inplane stresses to inplane
strengths of each ply and core layer. An additional failure index was defined by
the ratio of the shear strength to shear stress of each layer in the core. A failure
index greater than unity indicates initiation of material failure. It was assumed
that a suitable factor of safety was included in the definition of the loads, so no
additional factor of safety was used in the displacement or stress failure
constraints. The in-plane and shear strengths are listed in Tables 5.1 and 5.2 for
each material choice.

The state equations were the geometrically nonlinear static equations of
motion evaluated with the first order zigzag finite element model developed in
Chapter 4 of this dissertation. The finite element mesh and boundary conditions

are displayed in Figure 5.1. The sandwich panel was considered to be clamped

111

on all edges with no pressure applied on the front half of the panel. A 220psi
transverse pressure load is applied at the middle of the panel at a 45 degree
angle over a six inch band and a 30psi transverse pressure load is applied on the
rear half of the panel. The model was meant to represent a section of boat hull
under pressure loadings that occur during extreme usage. All finite element
meshes of the boat hull referenced in this paper contained a 14 by 6 in-plane
discretization with a single element through the thickness. The finite element
mesh was constructed to align the pressure loads with the finite element
discretization.

The baseline sandwich panel design was a symmetric layup composed of
two plies in the bottom skin, a single core layer and two plies in the top skin. The
orientations of all plies were aligned with the global coordinate system. Starting
from the bottom of the sandwich panel, the first and second plies were composed
of materials CM-1808 and XM-1808, found in Table 5.1. The baseline sandwich
panel core material had properties given by material A550 in Table 5.2, with a
thickness of 1.340 inches.

The baseline sandwich panel weighs 18.47 pounds with a maximum
transverse displacement of 0.73 inches and failure index of 3.14 when evaluated
with the geometrically nonlinear first order zig-zag finite element model
developed in Chapter 4 of this dissertation. The maximum transverse
displacement and failure index occur at the center of the panel within the high
pressure strip. The baseline design is infeasible in terms of the inequality

constraints placed on the maximum failure index.

112

5.3 Results

Two independent simple GA mns were performed. The two best designs
obtained using a simple GA were unsymmetrical layups. Figure 5.2 displays the
mass of the best design as a function of generation number for the first simple
GA run. At first, the mass is much larger than the baseline design, but eventually
converges to a design which is slightly lighter. Figure 5.3 displays the normalized
constraint of the best design found as a function of the generation number for the
first simple GA run. Initially, the stress constraint is severely violated, but it
eventually converges to a feasible final design. Tables 5.3 displays the baseline
and simple GA sandwich panel designs in terms their design variables,
objectives and constraints for two independent runs. The simple GA designs
have neariy identical masses compared to the baseline design while faintly
violating the displacement constraint. The maximum failure indices for the simple
GA sandwich panel designs are below the enforced failure constraint bound. The
two independent runs produced two distinct final designs, indicating the simple
GA prematurely converged to at least one suboptimal design. Each simple GA
mn lasted 5 days on a single 550 MHz cpu; each nonlinear finite element
evaluation took about 3 minutes. As demonstrated on this problem, a simple GA
tends to converge to suboptimal solutions while requiring large overall run times.

Figure 5.4 displays the force deflection response of the simple GA designs
and baseline designs, showing significant stiffening due to geometric

nonlinearities. Whereas the baseline design is stiffer at small loads when

113

compared to the simple GA designs, the stiffnesses for all designs are identical
at large loads. The GA has discovered designs which take advantage of the
geometrically nonlinear behavior of sandwich panels. The ply layup scheme and
the core of the simple GA designs develop smaller stresses by carrying larger
membrane stresses which help reduce stresses associated with bending while
maintaining the appropriate geometrically nonlinear stiffness. Load that is
transferred through bending causes a significant increase in overall stress when
compared to load that is transferred through membrane effects. If the state
equations had ignored geometric nonlinearities. the best designs found by the
simple GA could have been considered to be infeasible designs with respect to

the stiffness inequality constraint.

5.4 Conclusions

The design space of an optimization problem involving geometric
nonlinear response of sandwich panels was searched with a simple GA. The
simple GA found designs that outperformed the original baseline design. The
original baseline design satisfied the displacement inequality constraint, but
harshly violated the inequality constraint placed on the failure index. Both simple
GA designs satisfy all constraints with nearly identical mass to that of the
baseline design. The simple GA designs are more compliant at small loadings
(compared to the baseline design), but are specifically constructed to significantly
stiffen at larger loads, due to geometric nonlinearities. while maintaining feasible

stresses.

114

Location of Location of

 

 

 

maximum maximum 220 psi
failure index for displacement for Pressure
baseline design baseline design Load
A \Clamped _
V
C C
l l
a a
n n
P P
e e
d d
A " L
V V Clamped \ / 7
Back Front
No Pressure 30 psi
Load Pressure
Load

Figure 5.1. Finite element mesh including applied boundary conditions. Location of
maximum transverse displacement and failure index. The sandwich panel is 24 inches

wide and 84 inches long.

115

 

 

504
401
E30.
20.
10
0+- v v v 7 fl
0 5 10 15 20 25 30
mm

Figure 5.2. The mass of the best design found as a function of generation number.

116

4.5

4.0 .

3.5 4

.0
o

 

.N
or

+ Normalized Failure Index

 

 

+ Normalized Displacement

 

Normalized Constraint
N
'o

.A
0|

‘ ”wt-I'—

0.5

 

0.0 , . . , . . .
o 5 10 15 20 25 30
Generation Number

 

Figure 5.3. Normalized constraints of the best design found as a function of generation

number.

117

 

\\

1.4 . +§asline Luear‘
+GA RunJ Linear ‘
+GA Run 2 Linear
1.2 , +Baseline Nonlinear ‘
+GA Run 1 Nonlinear
+GA Run 2 Nonlinear ‘

\\

o

 

 

 

0

0.8 , ,, .

0.6 .

. //
0.4 - 1¢/
0.2 . /

Diephcomont (m)
l
l
l
\
\

 

 

Pereontoftotalload

Figure 5.4. Linear and geometrically nonlinear load deflection response curves for the

baseline and simple GA optimized sandwich panels.

118

Table 5.1. Ply material properties.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. Fiber . Longt. Trans.
N11153:] (iii) (nFiZi) (xiii) V Wight “3:13:88 Strength Strength
(oz/ydz) (ksi) (ksi)
CM-1208 2.17 1.66 1.34 0.3 18.9 37 36.76 14.13
CM-1215 1.86 1.33 1.34 0.3 25.6 37 23 14.13
CM-1608 2.13 2.07 1.22 0.3 22.5 44 24.75 13.17
CM-1615 1.77 1.87 1.22 0.3 29 56 22 13.63
CM-1808 2.59 1.74 1.19 0.3 24.4 43 35.27 13.34
CM-1815 2.2 1.82 1.15 0.3 31.2 57 28.48 15
CM-2308 2.37 2.28 1.36 0.3 29 53 29.9 13.86
CM-2315 1.85 1.79 1.18 0.3 36.7 73 26.33 13.07
CM-3308 2.38 2.52 1.18 0.3 36.7 61 41.19 13.07
CM-3415 1.98 2.01 1.1 0.3 36.7 80 25.46 11.65
CM-3610 2.17 2.27 1.1 0.3 36.7 76 29.86 11.65
XM-1208 1.3 1.33 2.03 0.3 19.2 39 13.59 27.04
XM-1215 1.2 1.02 2.19 0.3 26 47 15.68 27
XM-1708 1.31 1.25 2.05 0.3 24.4 44 14.26 31.82
XM-1715 1.36 1.35 2.25 0.3 31.1 51 16.2 31.49
XM-2408 1.14 1.34 2.05 0.3 31 50 10.58 31.82
XM-2415 1.22 1.62 2.19 0.3 37.8 66 13.53 24.99
XM-1305 1.85 2 2.1 0.3 24 26 8 18
XM-1308 1.85 2 2.1 0.3 24 29 8 18
XM-l708 .5 2.2 2.1 0.3 24 48 13.6 18
XM-1808 .5 2.2 2.1 0.3 24.8 48 13.6 23.4
XM-IBOSB .5 2.2 2.1 0.3 24.8 48 13.6 23.4
XM-2408 .55 2.2 2.2 0.3 30.8 56 14.2 33.2
XM-2415 .5 2.1 3.1 0.3 37.5 71 11.5 39.8
Table 5.2. Core material properties.
. . Shear Compressive
Core Matenal E1 Densrty
Name (psi) V (lb/fi3) 8:325:11 8:2?
A400 4.998e3 0.3 4.0 125.0 80.0
A450 6.5000e3 0.3 4.5 135.0 95.0
A500 8.4440e3 0.3 5.0 142.0 115.0
A550 8.8000e3 0.3 5.7 166.0 142.0
A600 9.2170e3 0.3 6.3 191.0 165.0
A650 9.8850e3 0.3 8.0 241.0 210.0

 

 

 

 

 

 

119

 

 

Table 5.3. Baseline and optimized sandwich panel design variables. All material

names are in reference to Tables 1 and 2.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Baseline Opt. Run 1 Opt. Run 2
# of plies in bottom skin 2 3 3
# of core layer(s) 1 l 1
# of plies in top skin 2 2 4
Orientation of bottom plies (degrees) 0.0, 0.0 30, 150, 60 0, 75, 120
Orientation of top plies (degrees) 0.0, 0.0 45, 75 0’ 15;)6135’
CM-2315, XM-1708,
Bottom ply materials $111883 CM-1208, CM-1208,
CM-3610 CM-1208
Core material(s) A550 A400 A400
Core thickness (inches) 1.340 1.0 0.625
XM- l 8088,
Top ply materials CM-1808, CM-3610, XM-1208,
XM-1808 CM-l808 CM—1608,
CM-23 15

 

Table 5.4. Baseline and GA optimized sandwich panel objectives and constraints.

 

 

 

 

 

 

 

 

 

Ineq. Constraint Baseline Opt. Run 1 Opt. Run 2
Mass (lbs) -- 18.47 17.88 18.94
Transverse
. . 0.73 0.73 0.74 0.74
displacement (Inches)
Max. failure indicy 1.0 3.14 0.92 0.87

 

120

 

 

Chapter 6

SUMMARY AND CONCLUSIONS

This dissertation addressed three major areas: multi-objective optimization
of crashworthiness, geometrically nonlinear analysis of laminated composite
materials and optimization of laminated composite sandwich panels considering
geometrically nonlinear response. Each area of research in this dissertation
contributed unique accomplishments within their field.

The shortcomings of available optimization techniques were overcome by
developing an approach that allows for concurrent multi-objective optimization of
structures with an island injection Genetic Algorithm for". 1) crash energy
management, 2) modal vibration frequency response, 3) peak crush force, 4)
weight reduction, and 5) manufacturability while considering stochastic variability
of design variables and loading conditions. Existing automated design technology
has never simultaneously addressed all these competing mechanisms. Three
independent runs converged to nearly identical designs. This result does not
guarantee that the global optimum was found; rather it indicates that a similar
search of the design space was completed for each run, from different initial
populations, making it at lease possible that a relatively thorough search was
conducted (noting that if radically different designs had resulted, this claim could
not be made).

A new geometrically nonlinear first-order zig-zag sublaminate plate theory

and finite element model was formulated and implemented. The new

121

geometrically nonlinear first-order zig-zag sublaminate plate theory and finite
element model was validated against existing experimental and numerical results
for laminated composite panels. This model is accurate for a wide range of
laminated panel applications including very thin laminated panels up to very thick
laminated and sandwich panels. The model is cast in a form that facilitates its
use within commercial finite element packages. The finite element models
developed take the form of a three-dimensional eight-noded brick or six-noded
wedge element with the usual five engineering degrees-of-freedom per node.
Finally, an optimization problem was considered involving geometric
nonlinear response of sandwich panels with a simple GA. The simple GA
designs were compliant at small loadings, but were specifically constructed to
significantly stiffen at larger loads due to geometric nonlinearities, while
maintaining feasible stresses. This was accomplished by the ply layup scheme
and the core of the simple GA designs, which develop smaller overall stresses by

carrying more load through membrane effects.

122

APPENDIX A

The displacement fields for the first order zig-zag theory can be expressed in

indicial notation as:
(t) _ — (t)
u" — u (15¢ “fl
ugf’ = flaw; a =1...4;,B =1,2 A.(1)

(k) _
13 _

u prfl
where the index ,6 is used to denote the top and bottom of the sublaminate with

1:bottom and 2:top, and (D23, ‘1’? and Opare shape functions in the thickness

direction (see Cho and Averill, 2000). Summation on repeated indices is implied

unless othenrvise noted. The variables 17 are:

afl
17m=ufl, 1725:”, 173fl=t9W 545:0”, A.(2)
Green’s strain can be written as the sum of linear and nonlinear parts:
E=sL+sNL A.(3)
The linear portion of the Green’s strain is:
a, =Lu=LNfi=BLii A-(4)

Since (DZ; ,‘1’3’ and 0’, are functions of the through-thickness direction only,

two dimensional shape functions (equations 4.51 — 4.54) can approximate the

nodal variables 17“”. The linear strains are first approximated as:

123

11 L (15.x; qfl
(k) _— (k)
822 L — u ,Xz‘yafl
(k) _.
833 L T prp x,
(t) _— (k) A. 5
723 1. ”uaqufifl + Wfi,x,Qp 11:3 ( )
2
(k) __ (k)
2’13 L _ uaflwd, + Wfi,lep Li
2
(k) __ - (k) - (k)
712 L - “macho? +uaM‘I’afl

where tensor notations are implied and a comma denotes partial differentiation.

The contribution of the transverse deflection to the transverse shear strains 7%"

(It)

,3 are evaluated at the at the middle of the sublaminate to obtain constant

and )1

shear stresses through-the-thickness of the sublaminate.

A constant transverse normal strain may give rise to substantial errors in
composites which have a soft core or layer (see Cho and Averill, 2000). An
improvement in the through-thickness description of deformation can be found by
assuming a constant transverse normal stress which can be obtained using
Reissner's mixed variational principal (see Reissner, 1986). The newly defined

transverse normal strain for small displacements can be written as:

‘0‘) _

_ (k) — (k)
833 L _ uafi.xl 23.696 + “(1,642 lxzafl + wfll

,, A.(6)

The coefficients ,1“) /l"‘

W. ,2; and 1,. are defined in Cho (1997).

For an eight noded brick element, B L can be defined as:

124

Bg") (1,11) = (1%me

B j“ (1, i2) = <1)??th

B)“ (1, i2) = 0

13;") (1, i4) = cam,

By) (1, i5) = (DEENLX
where:

19;") (1, i6) = (1):) Ni”,

B 2") (1, i7) = (“SM-J

13;") (1, i8) = o

B 2") (1, i9) = (1):]?th

Bg")(1,t10) = $2?)va

B,f*)(2,t1) = WSW,”
Bg‘)(2,i2) = \Pé’i’NL},
3,3“ (2, i3) = o

Bg")(2,t4) = 91$)Ni’y
32") (2, i5) = ‘11:;in
32") (2, i6) = sigma,
39) (2, i7) = Wilme,
32%, i7) = o

Bg")(2, i9) = ‘Pglme,

Effie, :10) = W§:)Ni,y

where:

i=1...4
i1=(i-1)*5+1
i2=(i—1)*5+2
i3=(i—l)*5+3
i4=(i—l)*5+4
i5=(i—1)*5+5
i6=(i—l)*5+21
i7=(i—l)*5+22
i8=(i-l)“'5+23
i9=(i—l)*5+24
i10=(i—1)*5+25

i=1...4
i1=(i-l)*5+1
i2=(i—l)*5+2
i3=(i—l)*5+3
i4=(i-l)*5+4
i5=(i—l)*5+5
i6=(i—1)*5+21
i7=(i—l)*5+22
i8=(i—1)*5+23
i9=(i—l)*5+24
i10=(i—1)*5+25

125

A.(7)

A.(8)

3900,11) = .1531 N,“ + 2;":le

BEWB, i2) = 2(3le + 431%»

32% i3) = -x: N.-

Bl"’(3. 14) = 293,14. + 19?le — 2:.- N1.-
Bi"’(3.i5) = 45:1)le + 43,071,): - It N2:

BEDG: i6) ___ 121% N” + 1(ykl)2Ni’y Where:
3290,17) = 18%,. + 1‘3, N...
32%. i8) = z.- N.-
31.100, 19) = 41,3)2Ni; + 93206,}: 1' Zr N If
31.106, 1’10) = 45.2%,). 1' AfzzNiJ 1’ Zr N2i
Bg")(4,t1) = ‘Pff’z N,
Bg‘)(4,i2) = trig; N,
Bl")<4.i3) = 41$") Iz=h/2 19,,
32") (4,:4) = erg): N,
B,£*)(4,i5) =w§fizN, .
where.

(k) - _ (k) ,
BL (4,:6)—‘¥12,ZN,
(k) - _ (k) ,
BL (4,17)—‘1’22’ZN,

. k ‘
Bik)(4,13)= Q(2 ) Iz=h/2 NW

i=1...4
i1=(i-l)*5+1
i2=(i-l)*5+2
i3=(i-1)*5+3
i4=(i-l)*5+4
i5=(i-l)*5+5
i6=(i—l)*5+21
i7 =(i-l)*5+22
i8=(i—1)*5+23
i9=(i—1)*5+24
i10=(i—1)*5+25

i=1...4
i1=(t'—1)*5+1
i2=(i—l)*5+2
i3=(i-—1)*5+3
i4=(t'—1)*5+4
i5=(i—l)*5+5
i6=(i-1)*5+21
i7=(i-1)*5+22
i8=(i-l)*5+23
i9=(i—l)*5+24
i10=(i-1)*5+25

By) (4, i9) = 4%): N,

Bg")(4,t10) = r32): N,

126

A.(9)

A.(10)

Bj")(5,i1) = ‘Pff: N,

B(k)(5,i2) = ~14") N,
B""(5 i3)- — 42‘ Mr... N.
Btk) - _ (k)
B, (5,14)—\P3LZN,
(k) - _ (k)

B, (5,15)—‘1’4LZN,-

B(")(5,i6) = 51"“ N,
(k)(517)‘1’(k)N-
""(518)n‘.’l.../2N

B(k)(5,i9) = ~14") N,

B(")(5, :10): 4"") N,

B(")(6, i1): ¢(*)N, +Wff’N,,
13g") (6, i2) = <D(2"1)N,y + WS)N,’X
Bg")(6, i3) = o
32") (6,14) = 629%, +‘1’§’;)N,,
By) (6, i5) = ¢§*I)N,y + \Pfi’N,,,
By) (6, i6) = (bilgNW +‘P(k)N,x
3,90 (6,17) = mg‘Z’N, +‘P(’2‘)N,x
(’9 (6 i8) = o
B(")(6,t9)= cth, +w§:)N,,

Ir 2k I:
BE )(6,i10) = (13(42)N,y 4.513%,

127

where:

where:

i=1...4
i1=(i—1)*5+1
i2=(i—1)*5+2
i3=(i-1)"'5+3
i4=(i-l)*5+4
i5=(i—1)*5+5 A.(11)
i6=(i—1)*5+21
i7=(i—1)*5+22
i8=(i—1)*5+23
i9=(i—1)*5+24
i10=(i-1)*5+25

i=1...4
i1=(i—l)*5+1
i2=(t'—1)*5+2
i3=(i—1)*S+3
i4=(i—l)*5+4
i5=(i—l)“'5+5 A.(12)
i6=(i-1)*5+21
i7=(i—1)*5+22
i8=(i—1)*5+23
i9=(i-1)*5+24
i10=(i-1)*5+25

It should be noted that the through-thickness shearing strains represented here
are in their original form. Cho and Averill enforced field consistency, consistent

penalty constraints and edge consistency to define alternate through-thickness

shearing strains which alleviated locking. All coefficients needed to formulate BL
can be found in Cho, 1997.

The nonlinear portion of the Green’s strain is:
1
8”, = EA(6)6(u) A.(13)

where 0(u) is a vector of displacement gradients and A(6) is a matrix

composed of the vector 6(u). For this element formulation, equation A.(13) can

be written in compact vector form as:

 

 

 

 

 

r w
1
I81?) EQIu3pvxiQau3av’r
(t) 1
822 _0 “3,0; 0 “3a,:
.23) 2 " '0’ " ’ a=1,2
4 :4 A. 14
réi’I O I B=1,2 ( ’
m
713 0
k)
trig J NI. 0
k fluJflv‘l «“3042 J z=h/2
where Q , are one-dimensional linear shape functions:
Q,=1--’-‘—3- (1:353— A.(15)
h h

where h is the height of the sublaminate, while b and t denote the bottom and top
of the sublaminate. The one-dimensional shape functions are evaluated at the

mid-plane of the sublaminate to ensure that transverse shear strain expressions

only contribute constant terms to the transverse stresses Sg’and S3”, to

(k)
l

(k)

and u2

maintain consistency with the terms in the in-plane displacements u

(equation 4.4).

128

The van Karrnan strain terms are present in the transverse normal strain due to
the constant transverse normal stress assumption. The nonlinear contribution of

the transverse normal strain for the kth layer is:

_ l —
= (5 we’d“ was“) -C.‘2’7.‘."’N.) A.(1.)

33 NL Cu.) 33M. 11 NL 22 NL
33

.93,” is then determined analytically from Reissner’s mixed variational principle

(equation 4.10) to obtain the following relation:

:3,

._ N . . . . O -
5., =42: llcs’s::’..+c:;’e‘” +C:;’r:;’..ldx. Am)

22 NL
i=1 X3f-I

where N is the number of layers in the sublaminate and the coefficient A1 is:

 

A, = fiI/[OCN ' x,_, ] A.(18)

.1 C,‘;"
Each nonlinear strain term in equation A.(17) is a function of the through the
thickness coordinate (x3). After substituting the one-dimensional through-
thickness shape functions (equations A.(15)) into equation A.(17) and integrating,
the analytical form of the geometrically nonlinear through-thickness constant

normal stress within a sublaminate (equation A.(17)) is:
_ — 1 1 1 A23 A36
S 33,, — -2— wflJ‘WMA 3.1;; + 2 WM w,M1 a, + why“, a, A.(19)

where a comma denotes partial differentiation, while a and ,6 are used in tensor

form, and the coefficients A130,”, A23 0,3. and A36“, are:

N N N
A1364? = ArZ'iitKLB A23ap = [112,233,413 A36ap = 4262K; A.(20)

3'] i=1 is}

129

and:

K;, = locum, A.(21)

x
3i—I

The gradients of the transverse displacement with respect to the in-plane

directions (x1 and x2) are:

4

W“, = 2(w . N,, + (9,111,, + 6MNM) 5 = x, ,x, a =1,2 A.(22)

where a comma denotes partial differentiation, N ,I, are the gradients of the

bilinear Lagrange interpolation functions and 117,; are the gradients of the shape

functions obtained from the interdependent interpolation scheme as found in
Cho, 1997.

The analytical form of the geometrically nonlinear through-thickness
normal strain for the k‘” layer in a sublaminate can be found by substituting

equations A.(15) and A.(17) into equation A.(16) to obtain:

1 1

50:) = __1__ Ewfltxrwav‘rAl 30.5 + 5 width WG-12A23afi

33 ”L Cg) + A36
WB.x.wa.x. aB

A.(23)

1 1
(k) (k)
E rlii Qflwflar, Qa Wag, + E ’23 0,6 W542 Qa me,

 

(k)
+ ’36 Qfiwfi.x.Qawa.x2 z=h/2

where a comma denotes partial differentiation, while a and ,6 are used in tensor

form.

130

Now that the constant transverse normal stress and transverse normal strain
have been defined including the von Kannan strains, we can move on to the first

infinitesimal variation of the nonlinear strain as defined in equation 4.31:
5s”, = A(9)6(‘fi)6‘fi A.(24)

Let 6(u) be defined as a vector containing all the nonlinear coefficients:

7'
0(u) = {WILXI wt,x1 Wb,X2 wt,x2 } A'(25)

A(6) was developed by taking the first variation of the von Kannan strains in

equation A.(14) including the geometrically nonlinear through-thickness normal

strain as defined by equation A.(23). For the k‘h layer in a sublaminate A(6I) is:

A(1,=1)"" -(§2,n w )

a 0.11

A(1,2)‘*’= -(00 .0 . w,_..,)
A(2’3)(k) = (Qb 0 w Nz=hl2

 

z=h/2

 

 

 

 

A(2,4)m =(QQ [Q a Wm) Fri/2
A(3,1)""= l.—.-. (W... A13» + Wu .,.A36 )
—r.(3 (“Q [,Qawa .x1 +r(k)Q anWa x2) =h/2
A(3,2)m=1(rt_) (WM A13” + Wa' ’1A36a)
33
—r|3( (”Q Qawa,x1 +r3ik)Q anWa .32) z=hl2
111(3.3)""=-1.—.,(Wa..AZ3bar+ W M436 )
,, A.(26)
_rn( (“Q (,Qawa .xz +ri")Q anw" ’31] z=hl2

131

A(3,4)(*>_1(_7_)(WMA23M +w,,,A36 )
_((r:(:")Q flaw“; +r,‘,“Q flaw...» ,zm

A(6, 1)‘*’= Q.w .,.,)l, _,,.

A(6,2)"" =(Q,Q.W. .Ll

A(6,3)‘*’ = (0 MGM. .II ,_,,.

A(6,4)"" =(rm w I...

raa,x,

wherea is summed from bottom to top. All other positions of the A(6) matrix are

homogeneous.
The stress stiffening matrix 8 was found from the relationship:
6A’s = 560 A.(27)

For the kth layer within the suiblaminate, S can be defined by:

_

 

 

 

 

 

 

 

 

5(1,1)<*> =,,[5 0, 0,4. [-tA—é—ijb - rff>c,c,,]5,,
33 - zsh/Z
$032)“) = ’7?erer + [—%%L _ rl(:)Qbe]§33]
L 33 rah/2
50,3)“ = 5,0,0, + [Aggy — rjj’Q,Q,]§,,- A.(28)
L 33 - z=hl2
Saar“ =”5.52.0. + (ASE? - r5."0.0.]§.]
- 33 z=h/2
5(2, 2)“) {=5}, (2 (2 +63% — rff’Q,Q,]5,,]
33 z=It/2

 

132

 

 

 

 

 

 

.. A —
S(2,3)(*) {929991—3735— :‘*>Q,Q,)S,]
33 z=lrl2
x A —
8(214)(“ = $12010: +[ :3)“ S’Q‘Qt) 33]
33 z=h/2
. " A23 — -
5053)“) = Szzabe+[ C(kfb _r::)Qbe] 33
- 33 - z=lrl2
5(3,4)<*> = 5,0,0, {Aggy r<*’0,0,]5,,
- 33 = z-h/2
. _ A2 _
S(4,4)(k) = 5,0,0,+( C3)" r:3k)Q,Q,]S,,]
b 33 ssh/2

 

Only the upper portion of the stress stiffening matrix S is displayed due to

symmetry of the second PioIa-Kirchhoff stress tensor. Note that all positions in S

are functions of the newly defined constant geometrically nonlinear transverse

normal stress.

Finally, the matrix of nodal shape function gradients G is:

G(1,i5) = N,
G(1,i5 +1) = N,
G(1, i5 + 2) = N,
G(2, 15 + 20) = N,
G(2, i5 + 21) = N,
G(2, 15 + 22) = N,

G(3, i5) = N,
G(3, 15 +1) = N,
G(3, i5 + 2) = N,
G(4, i5 + 20) = N,
G(4,1'5 + 21)
G(4, i5 + 22)

i=1...4

i5 = 5*(1—1)+3 A129)

1‘2
N
r',x2
A
int)

All other positions of G are zero.

133

We can define B and B as:

13 = B, + éAG A.(30)

"13:13, +AG A.(31)

B, B, G, and 8 form the direct stiffness and tangent stiffness matrices as

defined in equations (4.36) and (4.45).

134

REFERENCES

(1.) ADINA, 1999, “LS-DYNA User's Manual”.

(2.) Averill, R. C. and Reddy, J. N., 1992, “An assessment of Four-Noded
Plate Finite Element methods Based on a Generalized Third-Order
Theory”, Int. Journal for Numerical Methods in Engineering, 33: 1553-
1572.

(3.) Averill, R. C., and Yip, Y., 1996, “Development of Simple, Robust Finite
Elements Based on Refined Theories for Thick Laminated Beams”,
Computers and Structures, 59: 529-546.

(4.) Belegundu A. D. and Zhang S., 1992, “Robustness of Design Through
Minimum Sensitivity" Journal of Mechanical Design, 114: 213-217.

(5.) Chapman, C. and Jakiela, M., 1996, “Genetic Algorithm-Based
Structural Topology Design with Compliance and Topology Simplification
Considerations,” Journal of Mechanical Design, 118:89-98.

(6.) Cho, M., and Parmerter, R., 1993, “Efficent higher Order Composite
Plate Theory for General Lamination Configurations”, AIAAA Joumal,
31 :1299-1306.

(7.) Cho, Y., 1997, New Zig-Zag Composite Plate theories and 3-D Finite
Elements For Static and Explicit Dynamic Analysis of Laminated
Composite and Sandwich Panels, Ph.D. Dissertation, Michigan State
University, E. Lansing, MI.

(8.) Cho, Y3 and Averill, R.C., 2000, “First Order Zig-Zag Sublaminate
Plate Theory and Finite Element Model for Laminated Composite and
Sandwich Panels”, to appear in Composite Structures.

135

(9.) Crisfleld, MA, 1991, Non-Linear Finite Element Analysis of Solids and
Stnrctures, Vol. 1, John Wiley and Sons, Chichester.

(10.) DiSciuva M., 1985, “Development of an Anisotropic Multilayered Shear
Deformable Rectangular Plate Element”, Computers and Structures, 47:
91-105.

(11.) DiSciuva, M., 1993, “A General Quadrilateral Mulit-Layered Plate
Element with Continuous lnteriaminar Stresses”, Computers and
Structures, 47: 91-105.

(12.) DiSciuva, M., 1997, “Geometrically Nonlinear Theory of Multilayered
Plates with lnterlayer Slips,” AIAA Journal, 35:1753-1759.

(13.) Eby, D., Averill, R. C., Punch, W. and Goodman, E., 1999, “The
Optimization of Flywheels using an Injection Island Genetic Algorithm",
Evolutionary Design by Computers, edited by P. J. Bentley, Academic
Press Ltd.

(14.) Fabbri, G., 1997, “A Genetic Algorithm for Fin Profile Optimization,”lnt.
J. of Heat and Mass Transfer, 40:2165-2172.

(15.) Flynn, R., and Sherman, P., 1995, “Multi-Criteria Optimization of
Aircraft Panels: Determining Viable Genetic Algorithm Configurations,”
International J of Intelligent Systems, 10: 987-999.

(16.) Foster, N., and Dulikravich, G., 1997, “T hree-Dimensional
Aerodynamic Shape Optimization Using Genetic and Gradient Search
Algorithms,” Journal of Spacecraft and Rockets, 34:36-42.

(17.) Furuya H., and Haftka R. T., 1993, “Locating Actuators for Vibration
Suppression on Space Trusses by Genetic Algorithms”, Structures and
Controls Optimization, ASME, AD Col. , 38:1-11.

(18.) Furuya H., and Haftka R., 1995, “Placing Actuators in Space
Structures by Genetic Algorithms and Effectiveness lndices", Structural
Optimization, 9:69-75.

136

(19.) Gausser, M., Schueller, G., 1997, “Reliability-Based Optimization of
Structural Systems”, Mathematical Methods of Operations Research, 46:
287-307.

(20.) Genta, G. and Bassani, D., 1995, “Use of Genetic Algorithms for the
Design of Rotors,” Meccanica, 30: 707-717.

(21 .) Goldberg D. 1988, Genetic Algorithms for Search, Optimization and
Machine Learning, New York: Addison-Wesley.

(22.) Goodman, E., 1996, “GALOPPS, The Genetic Algorithm Optimized for
Portability and Parallelism System, User's Guide," Technical Report,
Genetic Algorithms Research and Applications Group (GARAGe),
Michigan State University, East Lansing, July, 1996.

(23.) Haftka, R. T., Gurdal, Z. 1993, Elements of Structural Optimization,
Third Edition, Kluwer Academic Publishers.

(24.) Gurdal, Z., 1999, Optimal Design, Theory and Applications to Materials
and Structures, Technomic Publishing 00., Pennsylvania, Edited by
Vasilev, V. V.

(25.) Hajela P., and Lee, E., 1997, “Topological Optimization of Rotorcraft
Subfloor Structures for Crashworthiness Considerations,” Computers and
Structures, 64: 65-76.

(26.) Haslinger, J. and Jedelsky, D., 1996, “Genetic Algorithms and
Fictitious Domain Based Approaches in Shape Optimization,” Structural
Optimization, 12:257-264.

(27.) Hibbit, Karlsson & Sorensen, 1997 “ABAQUS/Standard User Manual
5.8”.

(28.) Holland J. 1975, Adaptation in Natural and Artificial Systems, Ann
Arbon University of Michigan Press.

(29.) Jones, R. M., 1999, Mechanics of Composite Materials, Second
Edition, Taylor and Francis, Philadelphia.

137

(30.) Keane, A., 1995, “Passive Vibration Control via Unusual Geometries:
The Application of Genetic Algorithm Optimization to Structural Design,"
Joumal of Sound and Vibration, 185: 441-453.

(31.) Khot N., and Veley D., 1991, “Use of Robustness Constraints in the
Optimum Design of Space Structures”, Journal of Intelligent Material
System and Structures, 2:161- 176.

(32.) Kosigo, N., Watson, L., Gurdal, Z. and Haftka, R., 1993, “Genetic
Algorithms with Local Improvement for Composite Laminate Design,”
Struct. & Controls Opt., ASME, Aero. Div., 38:13-28.

(33.) Kuschel, N., Rackwitz, R., 1998 “Two Basic Problems in Reliability-
Bases Structural Optimization”, Mathematical Methods of Operations
Research, 46: 309-333.

(34.) Le Riche, T. and Haftka, R., 1993, “Optimization of Laminate Stacking
Sequence for Buckling Load Maximization by Genetic Algorithm,” AIAA J
31:951-956.

(35.) Lee, C., 1998, “A Continuum-Based Shell Element For Laminated
Composites Under Large Deformation", Ph.D. Dissertation, Michigan State
University, E. Lansing, MI.

(36.) Lee, C., 1998, A Continnum-Based Shell Element for Laminated
Composites under Large Deformation, Ph.D. Dissertation, Michigan State
University, East Lansing, MI.

(37.) Liu, 0., 1988, “Impact-induced Delamination- A View of Bending
Stiffness Mismatching”, Joumal of Composite Materials, 30: 1539-1561.

(38.) LSTC/KBS2, 1999, “LS-DYNA User's Manual".

(39.) Mallot B., Averill R., Goodman E., Ding Y., and Punch lll W., 1996,
“Use of Genetic Algorithms for Optimal Design of Laminated Composite
Cantilever Sandwich Panels with Bending-Twisting Coupling”, presented
at AIAA SDM (Structures, Dynamics and Materials).

138

(40.) Mares, C. and Surace, C., 1996, “An Application of Genetic Algorithms
to Identify Damage in Elastic Structures” Journal of Sound and Vibration,
195:195—215.

(41 .) Mayer, R., Kukuchi, N., and Scott, R., 1996, “Application of Topological
Optimization Techniques to Structural Crashworthiness”, Intemational
Journal for Numerical Methods in Engineering, 39: 1 383-1 403.

(42.) Nakaishi, Y. and Nakagiri, S., 1996, “Optimization of Frame Topology
Using Boundary Cycle and Genetic Algorithm,” JSME lntemational
Joumal, 39: 279-285.

(43.) Ochoa, 0.0. and Reddy, J.N., 1992, Finite Element Analysis of
Composite Laminates, Kluwer Academic Publishers, Netheriands.

(44.) Parrnee, I. and Vekeria, H., 1997, “Co-Operative Evolutionary
Strategies for Single Component Design,” Proc. 7th Int. Conf. on Genetic
Algorithms, T. Baeck, ed., Morgan Kaufmann Publishers, San Francisco,
529-536.

(45.) Prathap, G. and Somashekar, B. R., 1988, “Field and Edge
Consistency Synthesis of a 4-Noded Quadrilateral Plate Bending
Element”, Int. J. Num. Meth. Eng., 26: 1 693-1 708.

(46.) Punch, W., Averill, R., Goodman, B, Lin, S.-C. and Ding, Y., 1995,
“Design Using Genetic Algorithms - Some Results for Laminated
Composite Structures,” IEEE Expert, 10:42-49.

(47.) Punch, W., Averill, R., Goodman, E., Lin, S.-C., Ding, Y. and Yip, Y.,
1994, “Optimal Design of Laminated Composite Structures using Coarse-
Grain Parallel Genetic Algorithms,” Computing Systems in Eng., 5:414-
423.

(48.) Queipo, N., Devarakonda, R. and Humphrey, J., 1994, “Genetic
Algorithms for Therrnosciences Research: Application to the Optimized

Cooling of Electronic Components,” Int. J. of Heat and Mass Transfer,
37:893-908.

139

(49.) Rajan, S., 1995, “Sizing, Shape and Topology Design Optimization of
Trusses using Genetic Algorithms,” Journal of Stnrctural Engineering,
121:1480-1487.

(50.) Reddy, J.N., 1990, “A Review of Refined Theories of Laminated
Composite Plate”, Shock Vibr. Dig., 22:3-17.

(51 .) Reddy, J. N., 1993, An Introduction to the Finite Element Method,
Second Edition, McGraw-Hill, New York.

(52.) Reissner, E., 1986, “On a Mixed Variational Theorem and Shear
Deformation Rectangular Plate Theory“, Int. J. Num. Meth. Eng., 23:93-98.

(53.) Roux, D., Preez, P., Stander, N., 1996, “Design Optimisation of a
Semi-solid Tyre using Response Surface Approximations”, Engineering
Computations, 16: 165-184.

(54.) Sangren, E., Jensen, E., and Walton, J., 1990, “Topological Design of
Structural Components using Genetic Optimization Methods,” Sensitivity
Analysis and Optimization with Numerical Methods, 115:31-43.

(55.) Soto, C. and Diaz, A., 1993, “Optimum Layout and Shape of Plate
Structures using Homogenization,” Topology Design of Structures,
31:407-420.

(56.) Sundaresan S., lshii K., and Houser D., 1992, “Design Optimization
For Robustness Using Performance Simulation Programs”, Engineering
Optimization, 20:163-178.

(57.) Suzuki, K. and Kikuchi, N., 1990, “Shape and Topology Optimization
by a Homogenization Method,” Sensitivity Analysis and Optimization with
Numerical Methods, 1 1 5:15-30.

(58.) Suzuki, K. and Kikuchi, N., 1991, “A Homogenization Method for
Shape and Topology Optimization,” Computational Methods in Applied
Mechanics in Engineering, 93:291-318.

140

(59.) Taguchi, G., 1985, “Quality Engineering in Japan”, Proc. CAD/CAM,
Robotics and Automation International Conference, The University of
Arizona, Tucson, AZ.

(60.) Tessler, A. and Hughes, T.J.R., 1983, “An Improved Treatment of
Transverse Shear in Mindlin-Type Four-Node Quadrilateral Element”,
Computer Method In Applied Mechanics and Eng., 39:311-335.

(61.) Vietor, T., 1997, “Stochastic Optimization for Mechanical Structures”,
Mathematical Methods of Operations Research, 46: 377-408.

(62.) Vinson, JR, 1999, The Behavior of Sandwich Structures of Isotrepic
and Composite Materials, Technomic Publishing Co., Pennsylvania.

(63.) Wolfersdorf, J., Achermann, E. and Weigand, B., May 1997, “Shape
Optimization of Cooling Channels Using Genetic Algorithms,” J. of Heat
Transfer, 119: 380-388.

(64.) Yang, P.C., Norris C.H., and Stavsky, Y., 1966, “Elastic Wave
Propagation in Heterogeneous Plates”, Int. J. Sol. Structures, 2:665—684.

(65.) Yip, Y. , 1996, “ A 3-D Laminated Plate Finite Element with Zig-Zag
Sublaminate Approximations For Composites and Sandwich Panels”,
Ph.D. Dissertation, Michigan State University, E. Lansing, MI.

(66.) Yip, Y. C. and Averill, R., 1995, “On Modeling Local Damage in
Composite Laminates Using Plate Finite Elements with Zig-Zag
Sublaminate Approximations”, ASME lntemational Mechanical
Engineering Congress & Exposition, San Francisco, CA, November 12-17.

(67.) Yip, Y.C., 1996, A 3-D Laminated Plate Finite Element with Zig-Zag
Sublaminate Approximations for Composites and Sandwich Panels, Ph.D.
Dissertation, Michigan State University, East Lansing, MI.

(68.) Zaghloul, SA. and Kennedy, J.B., 1975, “Nonlinear Behavior of
Symmetrically Laminated Plates,” Journal of Applied Mechanics, Vol.
42:234-236.

141

(69.) Zenkert, D., 1997, An Introduction to Sandwich Constnrction, EMAS
Publishing, London.

142

IIllllllllllllllll