335:...»
h: tn... I.

£55.“ .
yiém .

.. a .. .
{Kidnap

r?

 

 

5;"
“I:
§
3

I)! 3.1

KM.

,2.
.z ‘y if. Sh”...
. L:36... ..
Ida. LBJ“. t :5
1.53.1.2}

sly. ‘
U ..

.rnfit: A... .3...
a .5 a. .
:olv... toe: .

1|..le
flunk War...“ I.
. fr) . “V 41.. 51.3.5.6 1.
In]. And. Kayak)“. r! . l

 

r'HESiS

STIY YLIBRAHI

IHIUWIHUWNIHHHIIIIII‘IH||HH||||| |

mini

LIBRARY
Mlchigan State
Universlty

 

 

 

 

 

 

 

 

 

 

This is to certify that the

thesis entitled

USE OF GENETIC ALGORITHMS FOR
OPTIMAL DESIGN OF LAMINATED
COMPOSITE SANDWICH STRUCTURES WITH
BENDING-TWISTING COUPLING

presented by

Bradley Arthur Malott

has been accepted towards fulfillment
of the requirements for

 

 

Master of Science degree in Mechanics

Mam

Major professor

 

J
Date uly 10, 1996

 

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

 

PLACE IN RETURN BOX to remove thie checkom from your record.
TO AVOID FINES return on or More dete due.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. C]

 

 

MSU le An Afflrmetive Action/Ewe! Opportunity lnetitution
WWW-1

 

 

USE OF GENETIC ALGORITHMS FOR
OPTIMAL DESIGN OF LAMINATED
COMPOSITE SANDWICH STRUCTURES WITH
BENDING-TWISTING COUPLING

By

Bradley Arthur Malott

A THESIS

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

MASTER OF SCIENCE
Department of Material Science and Mechanics

1996

ABSTRACT

USE OF GENETIC ALGORITHMS FOR
OPTIMAL DESIGN OF LAMINATED
COMPOSITE SANDWICH STRUCTURES WITH
BENDING-TWISTING COUPLING

By

Bradley Arthur Malott

The first section of this study compares different theories available for analyzing
laminated composite plates. Symmetric and unsymmetric square sandwich panels sub-
jected to sinusoidally varying transverse surface loadings were analyzed. The LZZHT
model proved to be the most accurate followed by the ZZHT, HSDT, and FSDT theories,
respectively. A genetic algorithm approach was applied to the optimum design of an ideal-
ized laminated sandwich composite aircraft wing using the ZZHT finite element model in
conjunction with the genetic algorithm, GALOPPS v3.0. The optimum stacking sequence
for the top and bottom face sheets was determined which produced the desired twisting
response, minimized weight, and maximized in-plane stiffness while maintaining stress
levels below the failure strength of the material in the fiber and matrix directions.

Four different GA architectures, single node, ring, injection, and hybrid structures
were used and their results compared. The ring structure produced the design with highest
fitness. The highest average fitness was achieved by both the ring. and hybrid structures.

The injection structure converged the fastest.

ACKNOWLEDGEMENTS

I wish to express my gratitude to Dr. Erik Goodman and Dr. William F. Punch III
for the support, insight, guidance and expertise offered throughout my research. I also
would like to express my gratitude to the Genetic Algorithms Research Group for their
patience, advice, and good humor. I also wish to thank Ying Ding for his diligence, hard
work, and assistance during this research.

I also thank all members of the Computational Mechanics Research Group for
their support and assistance. Finally, I extend my sincere gratitude to my advisor Dr.

Ronald C. Averill for his guidance and support through the duration of my research.

iii

Table of Contents

List of Tables ............................................................................................................ vi

List of Figures .......................................................................................................... vii

Chapter 1

Introduction ................................................................................................................ 1
1.1 Introduction ................................................................................................ l
1.2 Problem Definition ....................................................................................... 2
1.3 Genetic Algorithm Design .......... . ................................................................ 3
1.4 Literature Review ...................................................................................... 10

Chapter 2

Description of Laminated Plate Theories .......................................................... 37
2.1 Introduction ............................................................................................. 37
2.2 Equivalent Single-Layer Theories ............................................................ 38
2.3. Zig Zag Theories ...................................................................................... 44

Chapter 3

Assessment of Laminated Plate Theories for Analysis of Composite
Sandwich Panels ......................................................................................... 51
3.1 Introduction .............................................................................................. 51
3.2 Description of Models ............................................................................ 52
3.3 Results .................................................................................................... 57
3.4 Conclusions ................................................................................................ 62

Chapter 4

Optimal Design Laminated Composite Sandwich Structures Using

Genetic Algorithms. ................................................................................................. 77
4.] Introduction .............................................................................................. 77
4.2 Previous Work .......................................................................................... 79

iv

4.3 Model Description ................................................................................... 82

4.4 Optimization Problem ........................................................................... 86
4.5 Genetic Algorithm Topologies ................................................................ 93
4.6 Finite Element Model ............................................................................... 99
4.7 Results and Discussion .......................................................................... 101
4.8 Conclusions ............................................................................................. 109
4.9 Recommendations ................................................................................ 111

Chapter 5

Summary and Conclusions .................................................................................... 130

References ................................................................................................................ 134

List of Tables

Table (3.1) Sandwich Panel Material Properties .......................................................... 54
Table (3.2) Layup l. ...................................................................................................... 55
Table (3.3) Layup 2. ..................................................................................................... 56
Table (3.4) Convergence Data: Transverse center displacement. .............................. 56
Table (4.1) Sandwich Panel Material Properties. ......................................................... 86
Table (4.2) Weighting Coefficient Values. ................................................................... 89
Table (4.3) Mesh Convergence Data. ......................................................................... 100
Table (4.4) Results for Best Design of Each GA Structure. ......................................... 104
Table (4.5) Optimization Result From Various GA Structures. ................................... 113
Table (4.6) Layup Determined by Single Node Structure. ......................................... 125
Table (4.7) Layup Determined by Ring Structure GA. ................................................ 126
Table (4.8) Layup Determined by Injection Structure GA. ........................................ 127
Table (4.9) Layup Determined by Hybrid Structure GA. ............................................ 128

vi

List of Figures

Figure (2.1) - (a) Plate and coordinate system. ...................................................... 39
Figure (2.1) - (b) Undeformed and deformed plate section for CLPT .................... 39
Figure (2.2) - Undeformed and deformed plate section for FSDT. ........................... 42
Figure (3.1) — Sandwich plate geometry .................................................................... 52
Figure (3.2) - Quarter plate geometry and loading. .................................................. 53

Figure (3.3) - Normalized transverse center deflection versus aspect ratio,
ESL theories, layup l ...................................................................... 64

Figure (3.4) - Normalized transverse center deflection versus aspect ratio,
zig zag theories, layup 1. ................................................................. 64

Figure (3.5) - Normal stress vs. normalized thickness for ESL theories,
alh=4 layup 1. .................................................................................. 65

Figure (3.6) - Normal stress vs. normalized thickness for Zig-zag theories,
alh=4 layup 1. .................................................................................. 65

Figure (3.7) - Normal stress vs. normalized thickness for ESL theories,
alh=4, layup 1. .............................................................................. 66

Figure (3.8) - Normal stress vs. normalized thickness for zig-zag theories,
alh=4 layup 1. ................................................................................ 66

Figure (3.9) - Shear stress vs. normalized thickness for ESL theories,
alh=4 layup l. ............................................................................... 67

Figure (3.10) — Shear stress vs. normalized thickness for zig-zag theories,
alh=4 layup l. .................................................................................. 67

Figure (3.11) - Shear stress vs. normalized thickness for zig-zag theories,
alh=4 layup 1. ................................................................................. 68

vii

Figure (3.12) - Shear stress vs. normalized thickness for zig-zag theories,
alh=4 layup 1. ................................................................................ 68

Figure (3.13) - Normal stress vs. normalized thickness for ESL theories,
alh=10, layup 1. 69

Figure (3.14) - Normal stress vs. normalized thickness for zig-zag theories,
alh=10 layup 1. ................................................................................ 69

Figure (3.15) - Normal stress vs. normalized thickness for ESL theories,
alh=10, layup l. ................................................................................ 70

Figure (3.16) - Normal stress vs. normalized thickness for Zig-zag theories,
alh=10, layup 1. ................................................................................ 70

Figure (3.17) - Shear stress vs. normalized thickness for ESL theories,
alh=10, layup l. .............................................................................. 71

Figure (3.18) - Shear stress vs. normalized thickness for Zig-zag theories,
alh=10, layup 1. ................................................................................ 71

Figure (3.19) - Shear stress vs. normalized thickness for zig-zag theories,
alh=10, layup l. ................................................................................ 72

Figure (3.20) - Shear stress vs. normalized thickness for zig-zag theories,
alh=10, layup 1. ................................................................................ 72

Figure (3.21) - Normal stress vs. normalized thickness for ESL theories,
alh=10, layup 2. ............................................................................... 73

Figure (3.22) - Normal stress vs. normalized thickness for zig-zag theories,
alh=10, layup 2. ............................................................................... 73

Figure (3.23) - Normal stress vs. normalized thickness for ESL theories,
alh=10, layup 2. .............................................................................. 74

Figure (3.24) - Normal stress vs. normalized thickness for zig-zag theories,
alh=10, layup 2. ............................................................................. 74

Figure (3.25) - Shear stress vs. normalized thickness for ESL theories,
alh=10, layup 2. ............................................................................ 75

Figure (3.26) — Shear stress vs. normalized thickness for zig-zag theories,
alh=10, layup 2. ............................................................................ 75

Figure (3.27) - Shear stress vs. normalized thickness for zig-zag theories,

viii

alh=10, layup 2. ............................................................................ 76

Figure (3.28) - Shear stress vs. normalized thickness for zig-zag theories,

alh=10, layup 2. ............................................................................ 76
Figure (4.1) - (a) Elliptical spanwise load distribution approximation. ....................... 84
Figure (4.1) - (b) Linear chordwise load distribution approximation. ........................... 84
Figure (4. 1) - (c) Elliptical spanwise load distribution approximation. ........................ 84
Figure (4.2) - Sandwich panel dimensions in meters. ................................................. 85
Figure (4.3) - Schematic of Optimization Process. ..................................................... 87
Figure (4.4) - Chromosome structure. ....................................................................... 92
Figure (4.5) - Ring GA Architecture. ....................................................................... 94
Figure (4.6) - Injection GA Architecture, 7 subpopulations. ....................................... 96
Figure (4.7) - Hybrid GA Architecture. ......... A ............................................................... 99
Figure (4.8) - Fitness versus generation number for single node GA structure. ......... 114
Figure (4.9) - Fitness versus generation number for ring GA structure. .................... 115
Figure (4.10) - Fitness versus generation number for injection GA structure. .............. 116
Figure (4.11) - Fitness versus generation number for hybrid GA structure. ................. 117

Figure (4.12) - Fitness function terms vs. generation number for single node
topology GA. ................................................................................. 118

Figure (4.13) - Fitness function terms vs. generation number for ring structure GA. ...119

Figure (4.14) - Fitness function terms vs. generation number for Injection

structure GA. ................................................................................... 120
Figure (4.15) - Fitness function terms vs. generation number for hybrid

structure GA. ................................................................................... 121
Figure (4.16) - Extensional stiffness matrix versus generation number. ...................... 122
Figure (4.17) - Coupling stiffness matrix versus generation number. ........................ 123

ix

Figure (4.18) - Flexural stiffness matrix versus generation number. ........................... 124

Chapter 1

Introduction

1.1 Introduction

Optimization problems are present all disciplines and thus a wide range of goals
exist. For example, the goals of manufacturing optimization problems could be to maxi—
mize product output, minimize material waste, minimize cycle time, maximize equipment
use, etc., These goals, of course, coincide with the ultimate goal of business-related opti-
mization problems which is to maximize profit. In the engineering field uncountable opti-
mization problems exist. Of these, structural optimization is the concern of the current
study. In structural optimization problems goals are similar to other optimization prob-
lems. Minimum weight, minimum cost, maximum performance, and maximum manufac-
turability are a few common goals. In many aerospace Optimization problems are
concerned with determining the design which minimizes weight while maintaining struc-
tural integrity and function. Various techniques are used to solve design optimization
problems. Traditional techniques use gradient based approaches which have numerous
limitations. Random search techniques are also used, which also have limitations. Guided
random search methods using genetic algorithms are very effective for solving many opti-

mization problems. This technique is used to solve the structural optimization problem in

the current research.

1.2 Problem Definition

During flight an aircraft wing experiences bending and twisting loads which can be
detrimental to the aircraft’s performance. During flight the pressure differential across the
airfoil causes the wing to twist up and back. To compensate for this deflection washout is
built into the wing, providing the best compromise of the lift distribution through a range
of flight speeds. When a wing encounters wind gusts the wing’s lift distribution is
adversely affected. This situation is the motivating factor for the current research. To com-
pensate for such gust loads the wing needs to twist in the opposite direction from which
the load normally causes the wing twist. This type of response can be accomplished using
laminated composite sandwich panels by tailoring the lamination sequence accordingly to
utilize the inherent bending-twisting coupling of such materials.

The goal of the current research is to determine the design of a laminated compos-
ite sandwich panel which achieves the above mentioned “opposite” twisting response
while minimizing weight and maximizing in-plane stiffness (to maximize buckling
strength). To do this the optimum lamination sequence (i.e. number of plies, fiber orienta-
tions, and core thicknesses) must be determined. Sandwich panel dimensions, loading, and
geometry are representative of an idealized aircraft wing. The desired twisting response
was for the wing to twist -1.5 degrees from the horizontal in the direction opposite the
loading. This is four times the amount of twist, but in the opposite direction, exhibited by

a specially isotropic baseline sandwich panel with the same loading.

This optimization problem is solved using the genetic algorithm GALOPPS v.3.0.
All details of the optimization problem are discussed in Chapter 4. To determine the twist-
ing response and stresses in the wing a finite element plate code was used. The theory used

is discussed in Chapters 2 and 3.

1.3 Genetic Algorithm Design

Laminated composite structure design lends itself quite openly to optimization
techniques. Ply thickness and stacking sequence are two common design variables for
laminate design. Numerous optimization techniques are available, one category of which
is calculus-based or gradient-based methods. Gradient-based methods can be classified as
either direct or indirect [16]. Both approaches rely on the gradient of the objective func-
tion providing the necessary information to locate a local optima. One obvious drawback
of this method is the gradient of the function must exist and be obtainable. Determination
of the gradient may be impracticably difficult. Also, these methods seek out local optima,
which may or may not be the global optimum. Gradient based methods are not well suited
to find singular optima common in optimization of laminated composite structures. The
apparent lack of robustness of calculus-based optimization techniques rules out this
method for many applications.

Various algorithms exist which search finite design spaces, evaluating an objective
function at all points in a given domain. This approach is obviously inefficient and imprac-
tical for large domains. Random search techniques perform a random search of the design

space for optima, but are also inefficient. Guided random search techniques are a more

effective means of searching a design space. Genetic algorithms (GA) are guided random
searches which utilize information obtained during the search for direction. This method
mimics nature in the search for an optimum state. This technique is relatively efficientand
quite robust in comparison to other methods. Genetic algorithms are used in this study and

are explained more thoroughly next.

1.3.1 Genetic Operations

Natures means of dealing with existence and reproduction of all species is
described as “survival of the fittest”. Genetic algorithms draw from this phenomena,
employing a similar philosophy when searching a design space for global optima. GAs
perform operations of selection, reproduction, Crossover, mutation, and permutation of the

coded design variables.

1.3.2 Design Coding

In an optimization problem one must first define the problem. Design constraints
and design space must be identified and the objective function determined. Next, a neces-
sary operation is coding of the design variables to produce a string of finite length. These
strings are devised in a way as to represent a particular design. Strings are analogous to
chromosomes in natural systems. Strings are usually coded in either binary or integer form
depending on the situation. As an example, each digit of a string could represent the fiber
orientation of a particular ply in a laminate comprised of five plies:

[01221] (H)

where 0 represents a ply oriented at 0 degrees in some reference coordinate system, 1 rep-
resents 45 degrees, and 2 denotes 90 degrees. This string describes an individual design.
Each digit in a string has a representative counterpart in natural systems which are genes.
Binary coding is useful in creating new design options but does have limitations.
Furura and Haftka [13] provide an explanation of three disadvantages of using binary cod-
ing versus integer coding. First binary coding is complex, which is further intensified by
design solutions which are not a power of two. Another disadvantage is certain designs are
relatively inaccessible since it is not always easy to mutate one design into another similar

design. The third disadvantage is the child string may not resemble the parent.

1.3.3 Initial Population Generation
Creation of “life” is next. With a coding scheme decided upon, an initial popula-
tion must be generated. This is usually done stochastically. The number of individual

strings comprising the initial population is somewhat variable.

1.3.4 Reproduction/Selection

The first operation performed by the genetic algorithm in the reproduction process
is selection. In the selection process individual strings are selected based on their objective
function value, or fitness. Reproduction is essentially a synthetic version of natural selec-
tion. The objective function of each string is evaluated thus assigning a “fitness” to each.
Fitness values aid in determining a string’s probability of being “mated”; the more fit an
individual, the greater the probability it will survive and reproduce.

The most common method of determining which strings will be selected for

mating is called roulette wheel selection (see Goldberg [16]). In roulette wheel selection
each string, or design, is allotted a portion of a roulette wheel, the size of which
corresponds to the relative fitness (or objection function value) of the design. Strings with
relatively high fitness values are assigned a larger portion of the wheel and thus have a
greater probability of being chosen as parents for reproduction once the wheel is spun. As
is the case in nature, the “fittest” survive, and mating the most fit strings should provide fit
“offspring”. Other selection techniques are also used. In addition to roulette wheel
selection, Johnson [19] investigates and compares stochastic remainder selection, binary
tournament selection, and binary tournament with continuously updated sharing selection

techniques.

1.3.5 Crossover

The next genetic operation performed by the GA is crossover. The parent strings
selected for reproduction are mated and characteristics of each are combined in their off-
spring. Different levels of crossover are used. In two-point crossover two digits (genes) of

each string (chromosome) are exchanged. For example, given the two parents:

[01 L11] (1.2)
[11 Q] 1] (1.3)

two offspring are formed by two-point crossover of the underlined two digits to give:
[0 1 0 l l] (1.4)

[11111] (1.5)

The number of points and location(s) where crossover occurs varies. Variations of
averaging, random selection, and uniform crossover methods were compared by Furuya
and Haftka [13]. The probability assigned to the occurrence of crossover is generally quite
high, often near 1.0. Negendra et a1. [35] assign the probability of crossover to be 0.95

which means crossover will occur with a probability of 95% for each generation.

1.3.6 Mutation

Following the crossover operation is mutation. The mutation operator is applied to
introduce new genetic data into a population which helps prevent premature convergence
of the algorithm by maintaining population diversity. In other words, mutation protects

against genetic data being lost. As an example, the following occurs during mutation:

before mutation [0 1 1 1 l] (1.6)

after mutation [0 1 l 0 l] (1.7)

where the digit underlined is replaced by 0 upon mutation.

Mutation is applied with very low probability. LeRiche and Haftka [27] apply
mutation to a string with a probability of 0.01. Variations of the mutation operator were
investigated by Riche and Haftka in efforts to improve the GA. Minor improvements were

achieved by using new variations of the mutation operator.

1.3.7 Permutation

Permutation is another operator. Permutation is actually a type of mutation. This

stochastic operator typically inverts a substring of digits within a string:

before permutation [0 1 I Q l] (1.8)

after permutation [001 1 1] (1.9)

The three underlined digits are inverted. LeRiche and Haftka [27] also investigate varia-
tions of the permutation operator and show a decrease in computational cost from using a

new form of the permutation Operator.

1.3.8 Penalty Method

Genetic algorithms are unable to directly deal with constraints. Consequently a
constrained optimization problem must be transformed into an unconstrained problem.
This can be accomplished via penalty methods. In such methods penalty parameters are
incorporated into the objective function which penalize the fitness for violation of the con-
straints. Penalty methods also serve to drive the optimization process toward optimum
designs. Since many design configurations may satisfy the given constraints, the desired

choices are those closest to or at the global optimum.

1.3.9 Comparison of Genetic Algorithms to Traditional Methods

Genetic algorithms differ intrinsically from traditional search and optimization
techniques. As previously stated, traditional methods utilize derivatives whereas GAs
make use of objective function data. GAs work with a coded set of variables and not the

actual variables as do traditional methods. Another important distinction between GAs

and traditional methods is GAs simultaneously search from numerous points in the design
space and traditional methods search from a single point. To that end, it is apparent tradi-
tional methods may locate a local optima and not a global optima depending on the initia-
tion point of the search. The final fundamental difference is GAs are governed by

stochastic rules, and traditional methods are governed by deterministic rules.

1.3.10 Advantages and Disadvantages of GAs

In the current study GAs are utilized for optimization of a cantilever sandwich
panel. Genetic algorithms are ideas for this type of problem because of the highly non-lin-
ear behavior of the fitness with respect to the design variables. Linearization is not neces-
sary with GAs, while other methods place such. requirements in the optimization problem.
Also derivatives need not be continuous and do not even need to exist when using the GA
approach. GAs are very good to use for large design spaces, however they are relatively
inefficient for small domains.

A unique advantage of GAs is their fostering of innovative and creative designs.
Designs can be categorized as either routine, innovative, or creative. Routine designs are
quite common, being based on existing, well established ideas and constraints. Innovative
and creative designs are much more difficult to achieve. These designs require new con-
straints and goals. GAs can be considered a design search technique which is influenced
only by the design representation and performance which allows existing design para-

digms to be ignored.

10

1.4 Literature Review

Next a summary of literature reviewed and utilized in the current research is pre-
sented. The articles reviewed are categorized according to the content. First a section on
optimization is given followed by a section containing articles pertaining to GA design of
laminated composite plates. Other examples of optimization problems using GAs is listed
next. This section helped to open my mind to new methods for problem representation in
the GA. The next section contains optimization problems solved by other methods, which
served to enforce the use of GAs for solution to the current problem. Various plate theories
for analysis of finite element methods are presented in the next section of this literature
review followed by a section containing articles reviewed pertinent to the aircraft wing

design characteristics.

1.4.] Optimization

Olhoff and Taylor [37] present the basic concepts of structural optimization
emphasizing the fundamentals. They discuss the features and elements involved optimiza-
tion of both discreet and continuum structures, and they discuss the mathematical formula-
tion of such problems.

According to Olhoff and Taylor [37], fundamental structural optimization prob-
lems were of interest to the likes of Galileo (1638), J. Bernoulli (1687, Newton (1687),
Lagrange (1770), and later to Saint-Venant (1864), Maxwell (1869), and Levy (1873).

Olhoff and Taylor [37] group design variables into six categories. Cross-sectional

design variables are the first group. These variables refer to the properties of the structural

11

elements such as plate thickness and second moment of area, etc., which are typically
assumed to be continuous but may also be discreet. Topological and configurational
design variables describe the structural layout by specifying such things as the number
and location of structural members and joints in a design. Shape design variables are vari-
ables which are usually continuous and describe the structures shape. Material design
variables describe material properties and are generally discreet. The final category is the
support or loading design variables which describe the loading and boundary conditions
on the structure. This type of variable can be either discreet or continuous.

An optimal design problem can be termed discreet or continuous depending on the
dependence of the design variables on the spatial coordinates of the problem. Design prob-
lems such as truss structures are discreet where as plates and beams which are often con-
tinuous.

In design problems certain constraints are imposed on all designs to help identify
feasible designs. Olhoff and Taylor [37] divide constraints into two categories; behavioral
and geometric constraints.

Behavioral constraints are, in general, non-linear constraints which can further be
classified as equality or inequality. Equality constraints are, for example, equations of state
and compatibility which govern the structural response. Inequality constraints place
restrictions on the quantities governing the structures response by placing limits on stress
or deflection for example.

Geometric constraints can also be classified as equality or inequality. Geometric
constraints place restrictions on the design variables due to for example manufacturing

limitations or required appearance of the design.

12

To determine the “fitness” or “goodness” of a design an objective function is used.
The objective function is an expression given in terms of the design variables which has a
value for each design or each point in the design space. In an optimization problem the
objective is either maximized or minimized. In single-criteria optimization problems the
objective function represents a single property where as in multi-criteria optimization
problems the objective function contains a weighted sum of two or more properties.

While single-criteria optimization problems have been thoroughly studied, there
are few realistic problems which contain merely a single objective to optimize. Multi-cri-
teria optimization (MCO) problems involve the simultaneous optimization of several crite-
ria and consequently are more difficult to solve. Olhoff and Taylor [37] mention three of
the better known methods for solving MCO problems: weighted sum method, goal attain-
ment or bound method, and the constraint method. A brief discussion of these methods
will be given next. For explanation purposes ci(D) will represent the criteria for design D,

and the intent is to minimize with respect to D the vector {c1,c2,c3,. . . ,cn }.

For the weighted sum method the problem is posed as:

":"(Zwiqj (1.10)

where w, 2 0 are weights specified by the designer based on the experience and judge-

ment of the designer. This method determines D which minimizes the weighted sum of all

criteria.

13

In the constraint method the problem is given as:

min(ck)

D subject to ci-biSO (i=1,2,...n, i ¢k) (1.11)

where the constraint bounds 12,2 0 are again specified by the designer.

The goal attainment or bound method gives the problem as:

min(§)

D subject to wici—fi S. O (i=1,2, . . . n) (1.12)

This method can be used to solve a min—max problem where for example the criteria c,- are

identified with a maximum such as when the problem is stated as:

min ax(Wici
)] (1.13)

D i

In this situation the weighted sum and constraint methods cannot be used.

Stadler and Dauer [47] first present a brief historical review of the field of multicri-
teria optimization, tracing it’s origins to the area of economics where it was discussed by
Adam Smith in his book The Wealth of Nations (1776) in dealing with economics in a
competitive environment. Stadler and Dauer discuss optimization as finding Edgeworth

Pareto optimal points (EP-points), which are named after the two economists Francis

l4

Ysidro Edgeworth and Vrlfredo Pareto. EP-points are optimal or compromise points in a
multicriteria optimization problem.

Stadler and Dauer [47] outlined four steps involved in the solution of a multicrite-
ria optimization problem. First an appropriate mathematical model must be chosen, sec-
ond a design set must be selected, third preferences are determined to allow the best
design be chosen, and finally the optimality concept or condition must be selected.

Stadler and Dauer [47] state the method of constraints approach is a good choice
for solving multicriteria objective problems. This approach is good since nearly any math-
ematical programming code can be used. Many standard approaches for analyzing multi-
objective linear problems are based on the simplex algorithm which generates a set of
acceptable points from which the best design ischosen.

Finally Stadler and Dauer [47] outline examples of multicriteria optimization
problems from a survey they conducted. Two simple examples mentioned were the design
of a simple arch structure, and of an elastic truss structure. More difficult and interesting
optimization problems analyzed by Professor R. Statnikov of the Research Institute for
Mechanical Engineering of the Russian Academy of Sciences in Moscow were also
briefly mentioned. Statnikov applied his optimization methods to air-frame and body
design of the Russian space shuttle Buran. In this problem mass, structural rigidity, rib
spacing, and cross-section shape were the design variables. This optimization effort
gained a 300kg reduction in the tail section of the shuttle alone. Other work by Statnikov
included design of the Russian truck 21]... This problem contained 12 criteria.

Another example mentioned by the authors was conducted by H. Eschenauer, a

professor in the Laboratory for Structural Optimization at the University of Siegen,

15

Germany. In his work Eschenauer solves the general optimal design problem for

laminated conical shells with weight and deformations design criteria.

1.4.2 Application of GAs to laminate design

Design of composite laminates with design variables of ply thickness and fiber ori-
entation is a common continuous optimization problem. With only these two design vari-
ables present, the optimization is an integer problem which can be linearized. Branch
and bound algorithms have been used to solve this type of problem once it has been linear-
ized. Genetic algorithms can optimize nonlinear problems without the need for lineariza-
tion. Inclusion of frequency and strength constraints, such as in buckling load optimization
problems, leads to a nonlinear problem.

Haftka et a1. [24, 25, 27, 33] have used GAs for determining an optimum laminate
stacking sequence which minimizes laminate thickness or weight subject to strength,
buckling, and ply contiguity constraints. The ply contiguity constraint is included to avoid
matrix cracking by limiting the number of contiguous plies having the same fiber orienta-
tion. In these investigations the maximum allowable number of contiguous plies without
penalty to the objective functions four. Another penalty parameter is used which penalizes
the objective function for violation of the critical buckling load constraint. Performance
and accuracy depend on tuning of penalty values which is undesirable.

In these studies the laminates analyzed consisted of a simply supported plates with
compressive loads in longitudinal and lateral directions. The balanced, symmetric,

laminates were comprised of 0, +l-45, and 90 degree plies of equal thickness. To simplify

16
calculations and maintain balanced conditions, stacks of two 0 degree plies (02), two 90
degree plies (902), and a combination of one +45 and one -45 degree ply (44-45) were used

in combination. Flexural stiffnesses were expressed in terms of only two lamination
sequence parameters and material constraints.

An integer coding scheme is used in these investigations. Combinations of stacks
of plies containing 0, +/-45, and 90 degree fiber orientation were used to represent the bal-
anced, symmetric, laminate. Design strings consisted of the above mentioned stacks of

plies with the 02 stack of plies represented by the integer 1, the +/-45 stack represented by
2, and the 902 stack represented by 3. Empty stacks were represented by 0 in the string.

Empty stacks were packed to the outside of the string.

In studies [24] and [25] Haftka et al. utilized genetic algorithms with memory,
storing useful information about past designs in a binary tree. Data stored in the binary
tree consisted of the coded design string, objective function value, laminate thickness, in-
plane lamination parameters, bending lamination parameters, normalized strain failure
load, and normalized buckling load. Since during the evolution process designs are dupli-
cated, the binary tree can be searched to supply necessary data pertaining to a design being
evaluated, thus avoiding expensive duplicate calculations. The use of binary tree storage
was shown to reduce the number of analyses required by 30-40% for the design problems
considered.

A local improvement procedure utilized information stored in the binary tree.
Studies [24] and [25] used local improvement. Local improvement searches for improved
design in the neighborhood of a design generated by the genetic operations. For example,

if two stacks are interchanged, the resulting laminate is not explicitly evaluated but its

17

fitness is estimated based on neighboring designs stored in the binary tree. Neighborhoods
are defined based on relative distance between designs in the binary tree. In [25] one load
case showed an increase in performance by a factor of three from use of local
improvement. However, the average performance of three load cases implementing local
improvement was lower than without local improvement due to load case three which has
a singular optima.

As noted by Haftka et a1. [24, 25, 27], one possible difficulty encountered in lami-
nate design is that of singular optima. In these studies, strain constraints in plies of one
certain orientation can be satisfied by either adding more plies of this particular orientation
or removing these plies completely. As plies of this orientation are gradually removed
strain constraints become more violated until these plies are no longer present and the cor-
responding strain constraint vanishes. The optimum solution contains no plies of the par-
ticular orientation mentioned and thus is termed singular. For certain load cases this
situation typically occurs in plies with a fiber orientation aligned with the loading direc-
tion, usually zero-degree plies.

The genetic algorithms used in these studies difficulty removing these 0-degree
plies from the design because critical strain constraints are violated. To assist the algo-
rithm in doing so the initial population is seeded with designs containing no zero-degree
plies (NZP). This is called NZP seeding and is used in [24, 25, and 27].

Comparison of the basic GA with and without the permutation genetic operator
were made in [27]. In this computational cost was reduced by a factor of two by using per-
mutation. For a nominal increase in price of the search, permutation greatly improved the

practical reliability when compared to mutation. The authors define practical reliability as

18

the probability of reaching a practical optimum. Practical optima were defined as those
designs whose fitness values were within 10% of the “best so far” design. As a means of
evaluating performance, the normalized price is used which is defined as the number of
evaluations of the objective function divided by the practical reliability. Considering only
buckling constraints reveals the need for very high probabilities of permutation (300%)
where as by including strain constraints probabilities of only 100% are necessary [27].

Kogiso et a1. [24] utilize GAs for stacking sequence optimization of a laminated
composite plate for buckling load maximization. Local improvement and a binary tree
were used to improve performance and efficiency of the optimization process. This study
uses an “elitist plan” version of the GA which carries the best design from each generation
into the subsequent generation. Here again, singularities exist in the solution since the
optimum contains no 0 degree plies. This is remedied by NZP seeding.

N agendra et a]. [33] investigated the stacking sequence optimization of a rectangu-
lar, balanced, symmetric, plate with simply supported boundary conditions subjected to
buckling and strain constraints. The buckling constraint was linear in the ply identity
design variables but strain functions were nonlinear. The strain functions were linearized
to allow solution of the problem via sequential linearization using the branch and bound
algorithm. To predict failure due to in-plane compressive loading, a maximum strain crite-
rion (Jones [20]) was used. Results from the branch and bound algorithm were compared
with global optimum designs found using genetic algorithms. The branch and bound pro-
cedure results were within 0.5 percent of the global optima.

LeRiche and Haftka [27] sought to minimize the thickness of a composite laminate

plate subject to in-plane biaxial loading and buckling, strength, and ply contiguity

l9

constraints. Here again, a rectangular, symmetric, balanced laminate is used. In this
optimization, the use of new selection, mutation, and permutation operators are
investigated. A constant penalty is applied, penalizing the objective function for violation
of strength and buckling constraints. The need for large values of the penalty parameter
was alleviated by doing this. The penalty parameter may need to be very large for some
load cases to allow an objective function of a feasible optimum design to be smaller than
the objective function of better feasible optimum designs.

New selection procedures explored the use of various crossover schemes. A thick
crossover technique was used which confines crossover to the “full” portion of thinner
parent laminate being mated, where full means the portion of the string where plies are
present and not empty. Another crossover technique used restricts crossover to the full
portion of the thickest parent and is called thick crossover. Averaging crossover was also
investigated. These three crossover techniques were applied with different orders (number
of points where crossover occurs). Uniform crossover produced the best results and one
point crossover the worst. Each digit or gene of the offspring has equal probability of com-
ing from each parent in uniform crossover.

A new mutation technique was also employed. This new mutation differs from the
traditional mutation technique because of its ability to add plies anywhere in the laminate
and not merely on the outside where empty plies are found. The new permutation tech-
nique employed randomly selects and swaps two stacks from the full part of the string. It
was thought old permutation causes too much change in the string. This new permutation
provided a 56% reduction in search price. Scaling mutation was introduced which pro-

vided a 10% decrease in search price. The scaling mutation operation biases mutation

20

toward adding rather than deleting a stack (of plies) in infeasible designs, and toward
deleting a stack in a feasible design.

The study by Punch et a1. [39] sought to maximize the capacity of a laminated
composite beam to absorb mechanical energy without fracture using a GA with a new
course-grain parallel structure called island injection. This structure allows different levels
of refinement in the representation of the design. The best designs of each level are
“injected” into a node with higher resolution where it is further refined. This method gives
rise to “algorithmic” super linear speedup. A new finite element model developed by the
authors is used for evaluation of each design. This model assumes a piecewise continuous
through-the-thickness in-plane displacement distribution which accurately accounts for
the layerwise variations of stress and displacement in a laminate.

In this study a 24-layer beam made of graphite-epoxy composite layers with
clamped-clamped end conditions and an applied point load at midspan is designed. A thin
layer (about 5.5% of the nominal ply thickness) is placed at the top of each composite
layer, so there are actually 48 layers in the model. Each thin layer may be assigned the
same material properties as the layer immediately below it, or it may be assumed that the
thin layer is compliant, with stiffness properties three orders of magnitude less than the
composite layers. The length-to-thickness ratio of the beam is 50. The length of the beam
is divided into 20 divisions (finite elements), so there are 20x48 = 960 design elements.
The GA must decide whether to place a 0 degree ply (henceforth called material 1) or a 90
degree ply (henceforth called material 2) in each of the 480 structural ply design elements,
and whether or not to place a compliant material in each of the 480 thin layer design ele-

ments.

21

In this application, we used a number of representations but focused on a 480-bit-
string version which fully represents one half of the beam, and which was then mirrored

across the vertical centerline to create a full beam. Thus the search space examined by the

GA was 2480, representing all possible combinations of material and compliant material.
The results of this work were quite promising. The GA generated a number of unique
designs, including designs with “cantilever” beams of compliant layers.

In another application of GAs Nagendra et a1. [34] present a study designing stiff-
ened composite laminate panel with a hole for minimum weight. As in all other studies
mentioned thus far, the rectangular panel is balanced and symmetric and is subjected to in-
plane biaxial loading. Four equally spaced stiffeners with identical laminates were consid-
ered. Both the panel and stiffeners consisted of stacks of 0, +/-45, and 90 degree plies. The
panel and stiffeners possess buckling, strength, and ply contiguity constraints. Also, the
panel must not fail near the hole due to buckling or fail due to stress concentrations near
the hole. A finite element code and a stiffened plate code was used simultaneously with a
sequential approximate optimization procedure, which is a continuous optimization.
Designs based on eigenvalue separation were evaluated. Designs found using the GA
approach were superior to those found by the continuous approach, however, the computa-
tional cost of the GA approach were much higher.

In [35] Nagendra et at. modify the basic GA to reduce computational cost and
improve reliability. The particular subject of this study is the design of stiffened composite
laminates. Designs found by this modified GA were approximately 4 percent lighter than
those found using the unmodified GA. The panel considered was rectangular, balanced,

symmetric, and contained four equal spaced stiffeners. In plane loading was applied in the

22

form of normal and a small shear load which was applied to reduce sensitivity of off-
design conditions. A penalty parameter was employed to enforce ply contiguity con-
straints. Buckling and strength were also constrained. Crossover was modified to improve
performance of the GA. Substring crossover was used which effectively acted as an eight
point crossover compared to the four point traditional crossover. New mutation operators
called stack deletion, stack addition, and orientation mutation improved reliability, per-
forming better than the old mutation.

Arbate [2] examines various aspects of optimal design of laminated plates and

shells subjected to stiffness, strength, buckling load, and natural frequency constraints.

1.4.3 Other GA applications

Johnson [19] implements continuously updated sharing in the simple genetic algo-
rithm to determine the optimal placement of elastic supports on a simply supported plate.
This paper illustrates application of a tournament selection routine with continuously
updated sharing for optimization of elastic support locations on a square plate to maximize
the first system eigenvalues. Sharing serves to aid a GA in possession of multiple optima.
A sharing function is generally a power law function dependent upon the distance between
two individuals. It was found in this study that certain sharing algorithms caused chaotic
behavior when used in tournament selection, but the continuously updated sharing method
was successful. Selection methods investigated included the traditional roulette wheel
selection, stochastic remainder, binary tournament selection, and as discussed, binary
tournament selection with continuously updated sharing. The computational cost of the

GA approach was determined to be much higher than traditional optimization methods for

23

this set of parameters.

Another application of GAs which is unrelated to laminate plates is presented by
Leung and Nevill [28]. In this study GAs are used to generate a stable, light weight, simple
truss structure capable of supporting given forces. In this case, unlike in optimization of
laminated plates, a binary coding of the design variables is used in a two dimensional for-
mat. Consistent results were obtained for the design problem at hand. This was a very
basic study as it was intended to produce preliminary design possibilities.

Sandgren et a1. [46] provide an example of GA use in the automotive industry. A
genetic algorithm was used to both optimize the cross-section and minimize weight of an
automobile bumper subject to simple bending and to minimize the weight of a body panel
with stress and deflection constraints. Both cases involved a two dimensional chromosome
structure, where in most of the previously presented cases of laminate design a one dimen-
sional string of design variables was used. Penalty functions served the purpose of driving
material properties closer to their optimum values. Using the GA approach to design such
structures is feasible however the efficiency is questionable.

In [13] Furura and Haftka determine optimum placement of actuators on large
space structures using GAs. Averaging, uniform, random selection, and multiple-point
crossover were compared. Mutation appeared to be a more effective operator than cross-
over. The comparison of several variations of integer coded and binary coded GAs
revealed integer coding to produce highest performance and lowest cost.

In another application of GAs, Punch et a1. [38] use a genetic algorithm with a K
nearest neighbor (Knn) algorithm to classify soil samples from three environments found

in agricultural experiments. The actual experiment was conducted by researchers at the

24

Center for Microbial Ecology. A large design space of 3x100x96 was necessary for this
experiment which ruled out use of traditional approaches. The Genesis genetic algorithm
program was used with a Knn program which was written by the authors. On test cases
weighted Knn showed improved performance over the standard and normalized Knn. It
was also shown for large data sets the 0/1 weighting did not perform adequately. Parallel
processing was used to reduce run time. A nearly linear speed up was observed by the use
of parallel processing which reduced run time by 40% over the case using a single proces-
sor. The nearly linear speed up was achieved by using multiple processors to calculate fit-

ness which requires more time than all other GA operations.

1.4.4 Design optimization via other techniques
The following articles provide examples of design optimization problems solved by other
techniques.

Graesser et al. [14] present a design methodology for a stiffened laminated com-
posite panel subject to in-plane loads and bending moments. Design variables included ply
orientation angles of both the panel and stiffeners and stiffener geometry. Maximum strain
and minimum strength constraints were imposed in the minimization of structure weight.
The optimization software UWCODA (University of Washington Composite Optimiza-
tion Design Algorithm) was used in this study. This software, originally developed as a
design tool for laminated plates, utilizes classical lamination theory combined with the
global optimization algorithm called “improving hit and run”. This code deals with maxi-
mum strain failure criteria, buckling, and damage tolerance requirements.Effects of fiber

orientation on the objective function (structure weight) is illustrated. One case considered

25

only increments of 15 degrees as possible fiber orientations. The weight of structures with
this constraint imposed were found to be 30% higher than those with no restriction on
fiber orientation angle.

Wang and Costin [50] present an optimum design study of composite structures
including manufacturing constraints. They mathematically define the manufacturing con-
straints which control ply thickness variation, interleaving of plies, and ply orientation
percentage. These constraints were applied simultaneously in the design of a simple wing
structure for minimum weight using ASTROS optimization software. Ply thickness was
treated as a continuous design variable. Ply dropoff rate constraints controlled the rate of
thickness change between zones. Imposition of the ply interleaving constraint had the
smallest affect, increasing structure weight by 2.2%. The largest effect was seen by apply-
ing the ply percentage constraint which increased structure weight by 17.3%. Ply dropoff
rate constraint increased weight by 6.3%. Applying all three constraints produced an
increase of 35.3% in structure weight.

In [32] Miki and Sugiyama present an optimal design method for obtaining sym-
metric, orthotropic, laminated composite plate designs with required in-plane stiffness,
buckling strength, maximum bending stiffness, and natural frequency symmetric, ortho-
tr0pic, laminated composite plates. By minimization of strain energy, maximum stiffness
is determined. Optimum laminate configuration for maximization of buckling stress is an
angle ply laminate. Likewise for maximization of stiffness angle ply laminates are the
optimum configuration.

The multi-start global optimization technique is utilind by Miki and Sugiyama

[32] determination of optimal laminate configuration in multi-layer, symmetric, laminated

26

plates which maximizes stiffness. This study notes the great amount Of work already done
in this area but identifies a need for more effective and efficient techniques in optimization
Of laminated plates. The multi-start technique is shown tO be very efficient for such Opti-
mization problems. TO obtain maximum stiffness, strain energy is minimized. The strain
energy and its gradient are used to search for the optimum. In the multi-start technique tra—
jectories are initialized at numerous random starting points in the design space. Each tra-
jectory converges to a local minimum (since strain energy is being minimized), and is
compared to others. If a local minimum corresponds to the lowest value to date, the proba-
bility Of it being a global minimum is calculated. The constraints are then evaluated for
this case. Optimum ply orientations of centrally loaded plates, uniformly loaded plates,
and plates subject to combined loading were determined. Analysis Of the uniformly loaded
case produced multiple global optima as the number Of plies increased. Reductions Of 9.4-
42.5% in deflection Of both centrally and uniformly loaded cases were predicted for 16-
ply lay-ups when considering optimum configurations.

Fukunaga and Vanderplaats [11] use a mathematical programming method for
strength optimization of laminated composites subject to in-plane loading in a study Opti-
mizing strength Of a laminated composite plate with respect to layer thickness. The modi-
fied direction method Of the ADS (Automated Design Synthesis) program was utilized.
Design variables were ply thickness and orientation angle. Strength constraints for zero
thickness plies were deleted to improve efficiency and reliability Of the approach. Since
optimization with respect to orientation produces nonlinear constraints, the authors trans—
formed the design variables to reduce the nonlinearity of the in-plane stiffness. Results

show a strong dependence on initial points of search in finding the global Optimum instead

27

Of local Optimum.

Fukunaga and Vanderplaats [12] present a stiffness Optimization approach using
lamination parameters by using the example of buckling Optimization Of orthotropic sym-
metric balanced laminated cylindrical shells under combined loadings. The feasible direc-
tion method in the ADS (automated design synthesis) program determined Optimal
laminate configurations. Since stiffness components Of laminated composites are
expressed as linear functions of lamination parameters, their use as design variables is an
efficient method for such Optimization problems. Also, expression of buckling load as a
linear function of lamination parameters is simpler than expressing in terms of ply orienta-
tion angle. When considering the case of compressive axial loading for homogenous lami-
nates, the Optimal configuration was shown to be an isotropic laminate with respect to in-
plane and out of plane stiffnesses. This was also stated as the case by other studies. For
heterogeneous laminates (heterogeneous through the thickness) Optimal configurations for
the same loading case were shown to achieve 27% higher maximum buckling stress.

Kam and Lai [21] present an optimal design study Of optimizatiOn Of a four ply
symmetric laminated plate and of a sandwich plate. The four ply laminate was simply sup-
ported On all edges and subjected to a transverse, central point load. The cantilever sand-
wich plate was loaded by uniform pressure. The purpose of the Optimization was to
determine Optimal fiber orientation and ply thicknesses which minimized structure weight
subject to strength and vibration constraints. A two level approached was used tO do so.
The first level Of the Optimization procedure determined the Optimum fiber orientations
which maximized stiffness and the natural frequency of the first vibration mode. The sec-

ond level Optimized ply thicknesses. A finite element model containing five degrees of

28

freedom per node was utilized to determine strains. Transverse shear effects were included
and results Obtained were shown to differ significantly from those Obtained from classical
laminate theory. A recursive linear programming method was used as the Optimization
algorithm. The main purpose Of this paper was to illustrate performing such an Optimiza-
tion problem and no mention was made of convergence rate or efficiency Of the algorithm.

a method for design Of laminated composite plates for required flexural stiffness is
presented by Miki [31]. The laminates considered were multiple balanced angle-ply sym-
metric laminates. All plies are Of the same material. A design region was determined in
terms of two flexural lamination parameter thus creating a two dimensional parabolic ‘
shaped region. From this region the lamination point is determined for a given stacking
sequence. The lamination point is where all flexural stiffness values are determined. The
converse is also possible, that is, for a given lamination point the stacking sequence can be
found from this method. For the given laminate type it was shown that only two angles are
necessary and sufficient in designing a laminate with required flexural stiffness.

Avalle and Belingardi [3] present a theoretical study Of Optimization Of laminated
plates. The model consisted Of simply supported, rectangular, plates subjected to a
symmetric transverse pressure loading. The laminate was restricted to be symmetric. The
problem was that of determining the Optimal fiber orientation and ply thicknesses which
maximized flexural stiffness. The proposed analytical approach was validated by
comparison tO results from the AN SYS finite element code. Avalle and Belingardi showed
that for a orthotropic laminate (without bending-twisting coupling) the number Of
different fiber orientations necessary to achieve maximum stiffness is restricted to two.

The values Of these orientations are a dependent on ply thickness, stacking sequence, and

29

material properties. The layer thicknesses need to maintain global orthotropy of the
laminate.

Chattopadhyay et al. [8] present a multidisciplinary Optimization Of helicopter
rotor blades which they decomposed into three levels. Baseline performance for the rotors
is compared with a reference rotor. In the first level the total power coefficient is taken as
the Objective function with the rotor thrust coefficient constrained to insure at least the
same lift characteristics as the reference rotor. In level two vibratory and static stresses in
the blade and the dynamic and aeroelastic characteristics Of the rotor are Optimized. Spar
design is determined in level three by using the Optimal structural stiffness Obtained in
level two. In level three the Objective function is the total blade weight.

The nonlinear Optimization algorithm used is based on the method of feasible
directions. A two-point exponential approximation technique was used for Objective func-
tion evaluation in lieu Of exact evaluation Of the objective function which is much more
computationally expensive. Chattopadhyay et al. [8] claim that since gradient information
is retained by this technique it is better than using Taylor series expansions for this pur—
pose. Significant improvements over the reference configuration resulted from this study.
A reduction Of 19.8% in the coefficient Of total power was Obtained, vertical shear was
reduced by 45.8%, and total blade weight was reduced by 3.7%. Also, the Optimum con-
figuration contained a significant reduction in total rotor twist with only approximately 6
degrees of twist from root to tip which is a great improvement over the reference rotor.

The Object of a study by Lake et al. [26] was to provide a basis for understanding
elastic coupling in composite rotor blade structures and to develop a finite element

modeling method capable Of accurately predicting dynamic response of the rotor blades.

30

Experimental results were compared with analytical results for tubular spars of three types
of cross-sectional geometries. These three geometries were square, D-shaped, and
elliptical. Since these geometries are not circular they experience warping. Each spar was
comprised of four plies of T-650/42 graphite fiber with ERLX 1925-2 epoxy resin plane
weave cloth pre-preg. Five global modes were identified for each span type in the range of
0 to 2000Hz. These modes were (in no particular order) first and second lateral bending,
first and second vertical bending, and first torsion mode. Agreement between experimental
and analytical results ranged from 1.6% difference in the best of the elliptical spars to
13.8% in the worst case of the D-shaped spars.

Kam and Lai [23] use a constrained multi-start optimization approach to determine
the optimal design of a laminated composite plate“ for maximum stiffness. The problem is
converted into an unconstrained minimization problem and an unconstrained global opti-
mization method is used. Ply thickness and fiber orientation are treated as design vari-
ables. In the optimization process a shear-deformable laminated composite finite element
is used in the finite element analysis. Only symmetric laminates are considered. Results
from cases with simply supported or with fixed edges subjected to a center point load were
used for comparison with two other optimization routines, one of which was used by the
authors in a previous study and the other was the BCONF minimization routine from the
IMSL mathematical package. The optimization algorithm used was shown to be more effi-
cient and it provided laminates with higher stiffnesses than the other two optimization
methods. Kam and Lai [23] illustrated the affects of aspect ratio, length to thickness ratio,
and number of ply groups on the results.

Adali et al. [1] present two optimization problems. One of which is a lamination

31

scheme optimization for maximum buckling load and the other is for material cost mini-
mization. In both cases, discreet fiber orientations were considered with fixed ply thick-

nesses. The effects bending-twisting stiffnesses D16 and D26 are considered negligible

provided two dimensionless ratios involving these stiffnesses are within a certain range.
To reduce material costs hybrid construction was considered which uses high stiffness,
expensive graphite/epoxy outer plies and low stiffness glass/epoxy inner plies. Results
showed buckling loads of the Optimal hybrid laminates to be only about 10-15% less than
for laminates using all graphite/epoxy plies. The critical buckling load was shown to
quickly decrease as the number of glass/epoxy plies exceeded 1/2 of the total number of
plies. In some cases increasing the total number of plies was shown to increase critical
buckling load. Results also show for the case of luniaxial loading, laminates with aspect
ratios less than 0.6 consist of all 0 degree plies where as for aspect ratios ranging from 0.8

to 2.0 optimal layups consist of all +/- 45 degree ply orientations.

1.4.5 Plate theories and finite element models

In [29] Livine presents an equivalent plate structural model for wing shape optimi-
zation which includes transverse shear effects. Many of the codes available for multidisci-
plinary structural aerodynamic Optimization of wing configurations, such as ELAPS, LS-
CLASS, and ELFINI, are based on the kinematic assumptions of classical plate theory
(CPT). In this study, the author shows use of these assumptions leads to models which are
too stiff and do not accurately capture coupling between chordwise bending, spanwise
bending, and torsion. Natural frequencies of higher mode shapes are overpredicted by

these codes. In this paper a new plate wing modeling theory, which utilizes a first order

32

shear deformation plate theory (FSDT), is compared to CPT, finite element methods
(FEM), and actual measured results.

The model consisted of an all aluminum cantilever wing containing spars and ribs.
The wing was swept back by 30 degrees. A llb vertical force was applied to the tip of the
trailing edge. CPT results were 17% too stiff and natural frequencies were too high, espe-
cially in torsional modes [29]. The FSDT and FEM results correlated nicely with each
other and agree fairly well with measured results. Both plate techniques converged faster
than the FEM technique which is due to the smaller mass and stiffness matrices of the
plate theories. Since the plate theories are based on a Ritz approximation the order of the
polynomial can be varied. When the order was varied from 5 to 8 a change of less than 2%
was achieved for determination of the natural frequency and displacement. The author
feels using FSDT method improves reliability and generality of equivalent plate structural
wing modeling techniques.

Averill and Yip [6] present a high-order zig-zag finite element model for use with
thick laminated and sandwich beams. This element has the nodal topology of a four noded
planar element with beam type kinematics and an interpolation scheme which eliminates
locking. The number of sublaminates is variable for the linear zig-zag theory (LZZT) ele-
ment with increasing accuracy as the number of sublaminates increases. Results are pro-
vided for simply supported thick laminated beams with sinusoidally varying loading. The
LZZT results were compared with first order shear deformation theory (FSDT) and elas-
ticity solutions and were shown to accurately predict in-plane stresses and displacements
with use of only one sublaminate through the thickness. It was also noted that LZZT

results converge to the elasticity solution as the number of sublaminates is increased.

33

Burton and Noor [7] present analytic thennoelasticity solutions for sandwich pan-
els and cylindrical shells. These elasticity solutions are used for comparison with the mod—
els analyzed in the current research. Burton and Noor considered cases of simply
supported, doubly curved, sandwich panels with uniform pressure applied to one side. The
core was treated as a homogenous orthotropic continuum. The affect of variations of mate-
rial parameters of the face sheets and the core was studied by evaluating first-order sensi-
tivity coefficients. The procedure used by Burton and Noor [7] was shown to be very
efficient compared to other methods for obtaining converged normalized transverse dis-
placement. Also, the strain energy for different stress components was calculated and
compared. One set of results showed that for thick face sheets the strain energy in sand-
wich cylinders and panels was very sensitive to 'changes in El and G" of the face sheets
and of G13 Of the core. Geometric parameters such as the ratios of axial length to radius,
total thickness to radius, and face sheet thickness to total thickness were varied to deter-
mine their affect On strain energy components.

In [36] Noor summarizes various aspects of the mechanics of anisotropic plates
and shells. Computational models and deformation types were briefly discussed. A short
discussion on the predictor-corrector approach for analysis of thick anisotropic plates and
shells was also included.

In [10] Cho et al. present a finite element formulation which accounts for both
material and geometric nonlinearities in laminated composite structural components. An
isoparametric degenerated shell element containing independent rotational and displace-
ment degrees of freedom was used. To maintain continuity of shear deformation between

the layers the first shear deformation theory was used and a shear correction factor was

34

included to help approximate transverse shear energy. For treatment of the progressive
failure due to matrix degradation the Newton-Raphson method was used and for geomet-
ric nonlinearities the Lagrange method was used. For filament wound cylinders and lami-
nated plates numerical results were compared with experimental results. Numerical results
appeared to be in good agreement with experimental results however, it was mentioned
that the experimental results accuracy could be improved.

Averill and Reddy [4] present a four noded plate element based on the general third
order theory of Reddy. This theory is employed in the current research. This theory

assumes a cubic expansion of the in-plane displacement field in the thickness coordinate:

u = uo-zi’TwO-i-zlll £6318, (1.14)
v = Vo-Zg—j +z[l -§(%)12]8, (1.15)
W = wo (1.16)

where uo, v0, and wo are displacements in the x, y, and 2 directions, respectively, of a point

on the midsurface of the plate, (x,y,0). Traction free boundary conditions on the top and

bottom surfaces of the element are satisfied.

1.4.6 Aircraft design literature
Raymer [41] presents a conceptual approach to aircraft design, and in doing so he

provides general and historical data for various design parameters. This data proved useful

35

as a guide for selection of geometric and loading parameters used in the current study.
Wing aspect ratios typically are 7.5 for jet transport aircraft, 7.6 for single engine aircraft,
7.8 for twin engine aircraft, and as high as 9.2 for twin turboprop aircraft. Tail aspect
ratios range from 6-10 for sail planes, 3-4 for fighter aircraft, and 3-5 for other aircraft.

Airfoil thickness ratio, t/c (thickness/chord) is another parameter of interest.
Raymer [41] provided a graph containing historical data which he suggests as a good
method for initial selection of thickness ratio when designing an aircraft. For mach num-
bers slightly less than 1.0 this graph shows thickness ratios historically range from 0.06 to
0.14. Also a graph of test data for a NACA 65(216)-415 airfoil section with a chord of 24
inches has a thickness ratio of approximately 0.14. An interesting note is that supercritical
airfoils are approximately 10% thicker than subsonic airfoils. Quite often airfoil thickness
can be 20-60% thicker at the root than at the tip without substantially affecting the drag
due to fuselage effects. Raymer [41] also mentions the thicker root airfoil should not
extend more than approximately 30% Of the total span.

To estimate wing loading some idea of typical aircraft weight must be attained.
Raymer[41] says the take-off weight of a Lockheed S-3A is approximately 52,500 lbs and
that of a 05 is approximately 764,000 lbs. While an aircraft is at cruise speeds, i.e. level
unaccelerated flight, lift equals aircraft weight and thrust equals drag. The highest cruise
weight occurs at the beginning of cruise since at this time the most fuel is present. The air-
craft weight at the beginning of cruise is the take-off weight minus the fuel consumed to
reach cruise. Typical wing loading ratios, WIS = weight/wing lift area, are listed in table
5.5, pp. 84, which shows W/S ranging from 17 for single engine planes, to 26 for twin

engine planes, to as much as 120 for jet transports. As a safety measure a load limit of 5

36

g’s was first considered by the Wright brothers for design purposes. Depending on the type
of aircraft and application different load factors are assumed. Aerodynamic loading in the
chord-wise direction of typical airfoils can be crudely approximated by a constant load for
the first 15% of chord length then by a linearly decreasing load which is zero at the trailing
edge. In the spanwise direction a good approximation is an elliptic load distribution with
the least load at the wing tips.

Twist is often built into a wing to both prevent tip stall and modify the lift distribu-
tion to approximate an ellipse at some desired velocity. Two different types of twist are
aerodynamic and geometric twist. Geometric twist is the change in airfoil incidence angle
with respect to the root airfoil and aerodynamic twist is the angle between the zero-lift
angle of the airfoil and root zero-lift angle. Raymer[4l] recommends not exceeding a twist
of 5 degrees with 3 degrees being typical. Home built aircraft and general aviation type
planes usually have a twist of 2 degrees where transport aircraft use 1 degrees and military

aircraft 0 degrees.

Chapter 2

Description of Laminated Plate Theories.

2.1 Introduction

Composite materials are a good choice for structural materials in many applica-
tions due to their favorable properties of high strength-to-weight and high stiffness-to-
weight ratios. The use of thick-section laminated composites is particularly prevalent in
many military and industrial, aerospace and marine applications. Accompanying the use
of such materials is the structural designer’s need for an accurate means of analyzing
structures composed of these materials. To satisfy this need many different theories have
been developed. Certainly a fully three dimensional analysis is one method capable of
handling the problem. However, this method quickly becomes computationally infeasible
for even simple structures. Equivalent single-layer (ESL) and layerwise theories are two
other approaches to modeling laminated plates. These approaches can furnish adequate
results, depending on the design problem, while also providing a good balance of compu-
tational efficiency and accuracy.

This chapter discusses different plate theories and the associated finite element
models available for analysis of laminated composite plates. Five theories are discussed:

CLPT, FSDT, HSDT, ZZHT, and LZZHT. In the current research FSDT, HSDT, ZZHT,

37

38

and LZZHT are used for analysis of different laminates.

2.2 Equivalent Single-Layer Theories

2.2.1 Classical Laminated Plate Theory (CLPT)

Equivalent single-layer laminate theories reduce a three dimensional continuum
problem to a two dimensional problem. To achieve this reduction non-homogeneous lami-
nated plates are treated as a statically equivalent single layer having certain constitutive
behavior [44]. Classical laminated plate theory (CLPT), based on the kinematic assump-
tions of Kirchhoff plate theory, is the simplest in the hierarchy of ESL theories. This the-
ory is intended for use with thin laminates. CLPT is based on the Kirchhoff assumptions
requiring a straight line normal to the plane of the plate to be inextensible and remain
straight and normal to the plane after deformation [9, 20, 43], see Figure (2.1). Thus trans-

verse shear and normal effects are neglected. This is specified by the following conditions
8:0 8 =0 e,=0 (2.1)

The in-plane displacements are represented as linear functions of the thickness coordinate.

The displacement field for CLPT is given in Equation (2.2) as follows:

39

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

X
P
/ A ----- - —
/
1 /
/// /
Y T Z
Figure (2. 1) - (a) Plate and coordinate system.
no ___..
( ,x
\ i z
, Pt
awo
Z V V37 “’0

 

 

P awe
5x
awo
-z— —>
Bx

Figure (2. 1) - (b) Undeformed and deformed plate section for CLPT.

40

u(x.y.z) = uo(x.y. t)-zaa;wo(x.y)

V(x, y, Z) = V0(X, yr I) - 2%Wo(x, y) (22)

W(X, yr Z) = w0(x: y)

Where u, v, and w are the respective displacements in the x, y, and 2 directions, of a point
(x, y, z), and no, v0, and wo are displacements in the x, y, and z directions of a point on the
midplane of the plate, (x,y,0). as shown in Figure (2.1).

Since laminated composites have low transverse shear stiffness compared to in-
plane stiffnesses, the transverse shear and normal deformations, and thus stresses, are sig-
nificant and cannot be ignored. Classical plate theory neglects these transverse effects and
therefore cannot accurately predict the plate response to complex loads. Consequently, this
theory is Often not a good choice for modeling thick laminated composites unless the load-

ing is predominantly in-plane.

2.2.2 First Order Shear Deformation Theory (FSDT)

The next higher order equivalent single-layer theory is first order shear
deformation theory (FSDT) [48]. This is also called Hencky-Mindlin plate theory [43].
FSDT is intended for use in analyzing thin to moderately thick plates. A finite element
code employing this FSDT theory is evaluated and compared to other theories in the
current research. FSDT extends the kinematics of CLPT by allowing a line initially
straight and normal to the laminate mid-plane to remain straight but no longer normal to
the plane after deformation. This allows a constant state of transverse shear with respect to

the thickness coordinate. Transverse normal strains are assumed to be zero. These

41

conditions are represented as [43]
e, = 0 12yz = eyz(x, y) r»:Jrz = 5x209 y) (2.3)

FSDT assumes the following displacement field [42-44]

u(x.y.z) = uo(x.y)-z$x
V(x.y.z) = vo(x.y)-Z$y (2'4)

W(x, yr Z) = W0(x, y)

where, again, no, v0, and “’0 are displacements in the x, y, and 2 directions of a point on the
mid-plane of the plate 2 = 0, and ox and by are rotations of the transverse normal about

the y and x axes, respectively [42-44]. FSDT requires shear-correction factors to correct
the overpredicted shear strain energy. The finite element code based on FSDT used in this

research calculates the shear correction factors using the approach of Whitney [51]. As

8w
can be seen from Figure (2.2), when ox = 3:0 , FSDT reduces to CLPT.

The finite element model used in the current research which is based on FSDT
contains the “usual” engineering nodal degrees of freedom: displacements u, v, w, and

rotations if), , (13y. This model uses a 9-noded plate element topology.

Finite element models based on FSDT assume the displacement field is C0

continuous which requires the displacement field to be continuous through-the-thickness.

This is an improvement over use of CLPT which requires C1 continuity, in which case the

42

variation of transverse displacements through the thickness must be continuous and
smooth [5, 43]. In other words, CLPT requires transverse displacement and its first

derivatives with respect to the in-plane coordinates to be continuous.

2.2.3 Higher Order Equivalent Shear Deformation Theories (HSDT)

Second, third, and higher order ESL theories have been developed, see for example

[4, 42,, and 44]. These theories present additional degrees of freedom whose physical

‘10 ——u

 

 

 

 

 

 

 

2* wo

 

 

 

 

Figure (2.2) - Undeformed and deformed plate section for FSDT.

43

interpretation is usually not meaningful. With third-order or higher theories the gain in
accuracy is often minimal for laminated composites and is not worth the increased com-
plexity. In the current research a finite element model based on a general third-order the-
ory, developed by Averill and Reddy [4], is evaluated and compared with other theories.

In this HSDT theory the displacement field is of the form:

aw ..
“(xr yr Z) = uo-Z—O+Z 1-5(5)

 

 

3x _ 3 h 4
3w - 2-. (2.5)
v(x, y, Z) = vo-z-é—y-O-i-z 1";(i) ¢y

W(xr yr Z) = W0

where here again, 110, v0, and wo are displacements in the x, y, and 2 directions of a point

on the midplane (x,y,0). and 4),, and o), are rotations of the transverse normal about the y

8w
and x axes, respectively. The nodal degrees of freedom in this theory are u, v, w, —0,

3x

awo
a;

, if), , and by. The element developed employing this theory uses a four noded plate
topology.

Although equivalent single layer theories provide a good compromise between
computational efficiency and accuracy for determination of global response of thin lami-
nates, there are a few drawbacks of ESL theories. First, accuracy of the global response

decreases as laminate thickness increases. Second, such models are unable to accurately

predict strain and stress at layer interfaces and near geometric discontinuities. Finally,

44

transverse stresses cannot be accurately predicted by these theories [44].

2.3 Zig Zag Theories

Layerwise theories were developed to overcome the deficiencies of ESL theories.
The assumed displacement field of these theories provides a more correct kinematic repre-
sentation for thick laminates. Layerwise theories assume a unique displacement field for
each layer. The number of degrees of freedom in the problem may or may not be coupled
to the number of layers in the laminate. If the number of degrees of freedom does depend
on the number of layers, the analysis is similar in computational cost to a full three-dimen-
sional analysis which is the case in full layerwise theories. The number of layers and num-
ber Of degrees of freedom can also be uncoupled which is the case with the zig-zag
theories. To achieve this the transverse stresses at ply interfaces are required to be continu-
ous.

Two finite element models based on the zig-zag theory are used in the current
research. The finite element codes employing these finite element models are called
LZZHT, developed by Yip and Averill [52], and ZZHT. The LZZHT theory uses an eight-
noded brick finite element which allows a laminate to be represented by one or more lay-
ers (or sublaminates) by “stacking” elements through-the-thickness. A four-noded plate
element topology is used in the ZZHT theory. This element cannot be “stacked” in the

thickness direction and thus the sublaminate approach cannot be used.

45

2.3.1 LZZHT Theory

I The element based on the LZZHT theory of Yip and Averill [52] can be stacked in
the thickness direction. This is called a sublaminate approach. In the sublaminate
approach each sublaminate (or layer of elements) may contain several physical layers.
Within each sublaminate an independent displacement field is assumed in which the
through-the-thickness variation of in-plane displacements is represented by a piecewise
linear (or Zig-zag) function superimposed upon a cubic polynomial in the local thickness
coordinate [6, 52]. The transverse deflection is assumed to vary quadratically with the

thickness coordinate. The displacement field for the nth layer of the mth sublaminate is

[52]:

“(m n)

(x y.z)= 20km ‘ u).('")+ 2(2- 2.15m)

i=1]

(x m- 2 z W 2 (2 we (26>

i=1

.0" ">(., y. z)- wb("')(l - g).w,<m>(§). wgm>(§)(1-§)

“(m r!)

where subscripts t and b refer to the top and bottom surfaces of the sublaminate, n is the

layer number, h the sublaminate thickness, E.- and f], are the zig-zag functions, and ux, uy,

and 112 are displacements in the x,y, and 2 directions, respectively. The displacement field
variable uz contains the non-conforming bubble function to eliminate the Poisson’s ratio

stiffening effect which occurs in continuum elements when only one layer of elements is

46

used through the thickness of a region [52]. This stiffening effect arises from bending

energy terms due to the approximation function’s inability to linearly interpolate the trans-

verse normal strain through the thickness of the plate [52]. In the computation of trans-

verse shear strain W3 is ignored since the purpose of this term is only to eliminate

Poisson’s stiffening effects. This term is regarded as a nodeless variable degree of freedom

which is condensed out at the element level.

The linear elastic constitutive relations for the jth layer of a laminate are:

 

 

 

 

 

"Ox;(j) -C11 C12 C13 0 0 C130) 15x3”)

o” C22 C23 0 0 C26 er)

a“ = C33 0 0 C36 8“ (2.7)
Tyz C44 C45 0 ‘ 'sz

1;“ symmetric C55 0 7x2

EX); - C66. yo;

 

Using the small strain assumption, the infinitesimal strain-displacement relations are [30]:

e _Bu 8 _3v 8 _BW
”-33’ ”’33, 22-52’
(2.8)
_3w 3v _8w Bu 8v Bu
”2‘37"? “waits; ”Nata—y

The E.- and 1‘],- degrees of freedom in Equation (2.6) can be eliminated by

47

enforcing transverse shear stress continuity at each layer interface. The conditions for

shear stress continuity between layers 0) and (j+1) are:

(1') (j + 1)
xz XZ
(j) = 0+1)
yz yz

(2.9)

Using the strain displacement relations of Equations (2.8) along with the linear elastic

constitutive equations in Equation (2.7), Equation (2.9) can be solved giving expressions

for E.- and iii in terms of the remaining degrees of freedom (see [52]).

Further simplification is possible by satisfying the applied transverse shear traction

conditions on the top and bottom surface of the sublaminate as follows:

(1)

 

sz z=0 - Txb Tyz 2:0 _ Tyb
(N) (N) (2.10)

 

where z = 0 is the bottom surface and z = h is the top surface of the sublaminate, N is the
number of layers in the sublaminate, and Txb, Tyb, 1,“, and Ty, are applied shear tractions
(or interlaminar shear stresses when more than one sublaminate is used) at the bottom and
top surfaces of the sublaminate, respectively.

The final form of the displacement field is obtained by introducing the surface vari-

ables:

48

 

 

b = “inzgo i“0 “r “imbgh
2.11
v — um| = i7 v u(N)| ( )
b y 2:0 0 t y z=h
e _aw, e _aw,,
Xb-Fx- yb" (212)

The rotational degrees of freedom are expressed in the above form to allow the

assumed displacement field to remain C0 continuous. The constraints in Equation (2.12)
were enforced by penalty methods in the element formulation [52]. The displacement field
for each sublaminate is expressed in final form in terms of the operative degrees of free-

dom at the top and bottom surfaces of the sublaminate:

ub, ut, vb, vt, wb, wt, Oxb, 9,“, Oyb, 9)... 1,“), 1,“, Tyb, ‘12).t (2.13)

As mentioned previously, the number of degrees of freedom is not coupled to the
number of layers in each sublaminate. In the finite element model using this theory each
node contains seven degrees of freedom u, v, w, 9 and T. The 1 degrees of freedom are
easily dealt with when one sublaminate is used by setting them equal to the applied shear
traction on the top and bottom surfaces of the laminate. When more than one sublaminate
is used 1 is not specified at the internal nodes but is still specified on the top and bottom

surfaces of the laminate in the same manner as when one sublaminate is used.

49

The finite element model based on this theory contains the above mentioned seven
degrees of freedom per node. The element topology is an eight-noded brick, and thus the

total number of degrees of freedom per element is 56.

2.3.2 ZZHT Theory

Another third-order zig-zag theory used in the current research is called ZZHT.
The element developed from this theory contains seven degrees of freedom and a four
noded plate element topology for a total of 28 degrees Of freedom per element. This ele-
ment is computationally more efficient than the code employing the LZZHT element due
to the number of degrees of freedom per element for ZZHT being one-half that of the
LZZHT element. The ZZHT theory is more accurate than any Of the ESL theories. Good
accuracy coupled with acceptable efficiency makes the ZZHT theory a likely choice for
use in the Optimization portion of this research.

The displacement field for this theory uses a cubic polynomial superimposed onto
piecewise linear zig zag functions for the in-plane displacements, as done in the LZZHT

theory, however the transverse displacement is assumed constant:

u = ”’0 (2.14)

Interlaminar transverse shear stress continuity at ply interfaces is enforced to elim-
inate the associated zig zag degrees of freedom as was done in the LZBIT theory (see
Equation (2.9)). The displacement field is cast in final form in terms of the following

degrees of freedom:

50

U0, V0, W0, 9X’ By, ¢X’ (by (2.15)

where 110, v0, and wo are displacements as defined previously, Oxand 9,. are the same as in
the LZZHT theory, and (1),, and (by are the generalized rotations as used by Averill and

Reddy [4].

Chapter 3

Assessment of Laminated Plate Theories for Analysis of

Composite Sandwich Panels.

3.1 Introduction

The first portion of this research compares different theories for the analysis of
laminated composite plates. Various thick sandwich panels are analyzed using finite ele-
ment models based on‘four different laminated plate theories. Results are compared with
the elasticity solution of Burton and Noor [7]. The four theories compared are those previ-
ously discussed in Chapter 2. Two equivalent single layer (ESL) theories are used. One is
based on first order shear deformation theory (FSDT) [48] using shear correction factors
of Whitney [51] and a the second based on a modified version Of Reddy’s general third
order theory (HSDT). The remaining two theories used for this comparison are third -
order zig zag theories (ZZHT and LZZHT). This comparison of theories also aided in the
selection of a theory for use in the optimization portion of this research. An appropriate
theory is selected based on solution accuracy and computational efficiency which is
extremely important in optimization processes where large numbers of calculations are

performed.

51

52

3.2 Description of Models

The above mentioned four theories were compared for square sandwich panels
with simply supported edges and subjected to a sinusoidal surface loading as shown in

Figures (3.1) and (3.2). The simply supported boundary conditions were in the form:

u=w=0;0’ =0 on y=0,b 0<x<a O<§<h (3.1)

v=w=0;o =0 on x=0,a O<y<b 0<§<h (3.2)

where u, v, and w are displacements in the x, y and z direction, respectively, and C is a

    
    

.¢Q
9‘ / /
mmet lines ‘ \
)3“ ry \9 \

  

 

 

Figure (3.1) - Sandwich plate geometry

53

 

Figure (3.2) - Quarter plate geometry and loading.

thickness coordinate which is zero at the bottom surface of the laminate and 1 at the top.
Due to the symmetry of the problem it was only necessary to model one quarter of the
plate (see Figure 3.2).

Two lamination schemes were used in the present investigation. Material prOper-
ties used for both cases are given in Table 3.1. The first case analyzed consisted of a sym-
metric sandwich panel having 5 layers top and bottom surrounding a thicker core layer.
The lamination scheme for this case, called layup l, is given in Table 3.2. For this lamina-
tion sequence two different aspect ratios were analyzed, alh=4 and alh=10, where a is the
length dimension of one side of the full plate and h is the thickness. The lamination
sequence for the second case analyzed, layup 2, is given in Table 3.3. In this case the sand-
wich panel has a non-symmetric layup with 5 layers top and bottom surrounding a thick

COI'C.

54

Table (3.1) Sandwich Panel Material Properties.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Property Material 1 Material 2 Material 3 Material 4
E22 2.5e7 2.1e7 1.0e6 1.5e5
E33 1.0e6 2.1e7 1.0e6 5.0e4
v12 0.01 0.25 0.25 0.01
v23 0.25 0.25 . 0.25 0.15
v13 0.25 0.159 0.10 0.15
G12 0.5e6 8.0e6 0.5e6 2.17e4
G23 0.5e6 4.0e6 0.2e6 4.2e4
G13 0.2e6 8.0e6 0.5e6 2.17e4

 

 

55

Table (3.2) Layup l.

 

 

 

 

 

 

 

 

 

 

 

 

Layer No. Material Thickness Orientation
1 1 0.01 0
2 2 0.025 0
3 3 0.015 0
4 1 0.020 0
5 3 0.030 0
6 4 (core) 1 0.800 0
7 3 0.030 0
8 1 0.020 0
9 3 0.015 0
10 2 0.025 0
1 1 1 0.010 0

 

 

 

 

 

56

Table (33) Layup 2.

 

 

 

 

 

 

 

 

 

 

 

 

Layer No. Material Thickness Orientation
=1=1=T=W=—fl6‘=l
2 3 0.01500 90
3 1 0.02000 90
4 3 0.03000 0
5 1 0.01000 0
6 4 0.40000 0
7 4 0.40000 90
8 3 0.03000 0
9 1 0.02000 0
10 2 0.02500 90
l l 1 0.01000 90
12 3 0.01500 90

 

 

 

 

 

 

Table (3.4) Convergence Data: 'h'ansverse center displacement.

 

 

 

 

 

D333; FSDT HSDT ZZHT LZZHT
W505? .m
16 -O.3529E-3 -0.1659E-3 -0.2144E-3 -0.2065E-3
64 -0.3518E-3 -0.l700E-3 —0.2118E-3 ozone-3

 

 

 

 

 

 

 

57

3.3 Results

In all cases the square domain was discretized into a mesh containing 8 elements in
the x-direction and 8 elements in the y—direction. This mesh was chosen based on a con-
vergence study of four meshes with mesh densities of 4, 16, and 64 elements (see Table
(3.4)). Mesh density is defined as the number of elements per layer in the mesh. The mesh
for the FSDT theory contains 9-noded plate elements, while the HSDT and ZZHT theories
uses a mesh of 4-noded plate elements, and the LZZHT theories uses a mesh of 8-noded
brick elements.

Results produced by the finite element model based on first order shear deforma-
tion theory are denoted FSDT [53] in the following discussion. Results from the third
order, equivalent single theory of Reddy [46] are labeled HSDT, and those from the third
order zig zag theories are labeled ZZHT and LZZHT. Two different meshes are used with
the LZZHT finite element code, denoted as LZZHT-l, which contains one layer of ele-
ments (or one sublaminate) through the thickness of the plate, and LZZHT-3 which con-
tains 3 layers of elements (or sublaminates) through the thickness. LZZHT-3 uses one

layer of elements for each of the top and bottom face panels and one for the core.

3.3.1 Layup 1:

A laminate having the lamination scheme defined in Table (3.2) was first analyzed.
Figures (3.3) and (3.4) display results for center deflection normalized with respect to the
elasticity solution of Burton and Noor [7] versus plate aspect ratio. Large errors are

present when the plate is thick. For the equivalent single layer theories, the predicted

58

center deflections are accurate for aspect ratios greater than 100. FSDT results are within
4% of the elasticity solution for an aspect ratio of 100 and within 1% at alh=500 and a/
h=1000. HSDT results are within 2.1% of the elasticity solution for an aspect ratio of 50
and within 1% for aspect ratios of 100 and greater. At an aspect ratio of 4, ZZHT and
LZZHT-1 results are within 3.8% and LZZHT-3 results are within 2.9% of the elasticity
solution. For aspect ratios greater than 4 the center deflections predicted by all three zig-
zag theories are within 1% of the elasticity solution.

Figures (3.5) through (3.10) contain results from the four different theories for in-
plane stresses for an aspect ratio of alh=4. Since the core experiences very little in-plane
stress, results are presented for the top face sheet only. As can be deduced from these
graphs, FSDT results are the least accurate and only slight improvement is gained using
the third-order HSDT theory. Both zig-zag theories predict in-plane stresses very well.
LZZHT-1 and ZZHT results are comparable in accuracy to each other and are very close
to the results obtained from the elasticity solution. LZZHT-3 results are the most accurate,
virtually identical to the elasticity results. These results illustrate the improvement in accu-
racy as the mesh is refined from 1 element through the thickness in the LZZHT-l model to
3 elements through the thickness using the LZZHT-3 model. The elements from ZZHT

theory cannot be refined through the thickness.

The maximum error in 0',“ of the FSDT theory is 118% and the maximum error

from the HSDT theory is 71%, (see Figure (3.5)). Both zig-zag theories predict on very

well as shown in Figure (3.6). The ZZHT theory error is the largest of the zig-zag theories
with a maximum of 18.5%. Maximum error of the LZZHT-1 model is 13.9%, and an error

of 4.3% was found for the LZZHT-3 model.

59

Results for (1'yy in Figures (3.7) and (3.8) are similar to the results of 0’“. The

FSDT theory produced the least accurate results with a maximum error of l651%. The
HSDT theory maximum error was 27%. In the top ply of the face sheet the error from the
LZZHT-1 model was the largest Of the zig-zag theories at 9.7% followed by 7.4% error by
the ZZHT theory and finally 2.3% from the LZZHT-3 model.

Figures (3.9) and (3.10) give in-plane shear stress results. The maximum error of
the FSDT theory was 43%, 19% for the HSDT theory, 10% for the LZZHT-1 model, 2.5%
for the ZZHT theory, and 1.4% for the LZZHT-3 model.

Transverse shear stresses are compared in Figures (3.11) and (3.12) for the case a/
h=4. Only results from the zig-zag theories are compared since the ESL theories do not
sufficiently model transverse stresses. As before, results from the ZZHT and LZZHT-1
models are nearly identical in accuracy. It is important for a laminate theory to provide
accurate results in the core region of sandwich panels. Core layers are typically much
softer than the face sheets and failure due to transverse shear effects can be predicted if
accurate transverse shear stresses can be determined. Both the ZZHT and LZZT theories
provide good results in the core region of the plate with errors of 4% for the ZZHT and

LZZHT-1 models and 1% for the LZZHT-3 model for Iyz. Results for sz in the core

region are better yet with errors of 1% for the ZZHT and LZZHT-1 models and 0.3% for

the LZZHT-3 model.

The maximum transverse shear occurs in the face sheets and not in the core. The
maximum error in Tyz in the top face sheet was 7.1% for the LZZHT-1 model, 6.9% for
the ZZHT model, and 0.2% for the LZZHT-3 model. The maximum error in In for the

ZZHT and LZZHT-1 models was 42% and for the LZZHT-3 model it was 12%.

60

Figure (3.13) through (3.20) contain in-plane stress results for the plate with an
aspect ratio of alh=10. Accuracy of the results from all models show improvement over

the case with alh=4 as expected, however, results are still quite poor for the two ESL theo-

ries. The maximum error in the prediction of on for the FSDT theory was 279%, and for

the HSDT theory it was 313% (see Figure (3.13)). All zigzag theories predicted on very
well with a 3% error from the LZZHT-1 model, a 0.8% error from the LZZHT-3 model,
and a 0.7% error from the ZZHT theory, (see Figure (3.18)

Results for 6”,, given in Figures (3.15) and (3.16), are similar to the results of On.

The maximum error in the predictions of O’yy by the FSDT theory was 23% and for the
HSDT theory it was 13%. Here again, the zigzag theories did an excellent job predicting
O” with a maximum error of 3% from the ZZHT and LZZHT-1 models and a maximum

error of 0.7% for the LZZHT-3 model.
In-plane shear stress predictions by the FSDT model contained a maximum error

of 15%, and the HSDT model’s maximum error was 2.9% (see Figure (3.17)). The zigzag
theories again did an excellent job predicting an. Predictions from the ZZHT model,

given in Figures (3.18), contained a maximum error of 2.8%. The maximum error of the

LZZHT-1 model is 0.6%, and for the LZZHT-3 model the maximum error is 0.1%.
Figures (3.19) and (3.20) contain results for transverse shear stresses In and In as
predicted by the zigzag theories. For both stresses the maximum error occurred in the face
sheets. The maximum error in the predictions of In by the ZZHT model was 18%, maxi-
mum error from the LZZHT-1 model was 17%, and maximum error from the LZZHT-3

model was 4.8%. Maximum errors in In are as follows: 21% for the ZZHT model, 20%

for

’13

the
M

1.2

fa

tic

tra

 

61

for the LZZHT-1 model, and 11% for the LZZHT-3 model.

3.3.2 Layup 2:

A laminate having the lamination scheme defined in Table (3.3) and an aspect ratio
of alh=10 was analyzed next. Layup 2 is a non-symmetric, random layup compared to
layup 1 which has fiber orientation angles of zero degrees and a symmetric lamination
scheme. Figures (3.21) - (3.28) contain the results comparing stresses from the different
theories for this case. The same conclusions were arrived at as were for the case of layup 1
with aspect ratio alh=10. Results from FSDT and HSDT codes are the least accurate and
LZZHT-3 the most accurate. Transverse shear stresses are very accurate through the entire
plate thickness for the LZZHT-3 model. Results from the LZZHT-1 and ZZHT models
were again comparable, with a maximum error of 7% for LZZHT-1 model and 6% for
ZZHT compared to the elasticity solution, in the core region at the core/bottom face sheet

interface for Tn. At the top face sheet/core interface results for “tn from the LZZHT-1

model are a maximum of 6.5% less than the elasticity solution with results from the ZZHT
model being under predicted by 5.6%. This is the maximum error in the core region.
Results for Tn are very good in both face sheets. In the bottom face sheet In at the inter-
face between the bottom two plies is predicted to be 36% greater than the elasticity solu-
tion by the LZZHT-l model and 33% greater by the ZZHT model. The maximum
transverse shear stress in the core region of all cases analyzed is predicted very well by

both LZZHT-1 and ZZHT theories.

62

3.3.3 Computational Efficiency

Computational efficiency of the finite element codes employing the ZZHT and
LZZHT-1 models were compared since both use the same mesh and produce similar
results. All theories are programmed in FORTRAN77. Both the ZZHT and LZZHT-1
models used the same mesh size (64 elements) but the ZZHT model contains only 4 nodes
per element giving 28 degrees of freedom per element where as the LZZHT models con-
tains 8 nodes per element and 56 degrees of freedom. The mesh used with the LZZHT-3
model contained three layers of 64 elements or a total of 192 elements.

The average system time required for one run of the code employing ZZHT finite
element model was 0.21 seconds while one run required 0.68 seconds using the LZZHT
model. This is 224% longer than the ZZHT model. User time for the ZZHT model was
8.58 seconds but for the LZZHT-1 model user time was 383 seconds, 4,364% longer than
ZZHT. The LZZHT-3 model takes slightly longer than LZZHT-1. Based on these results
the most efficient of the zig zag code is the ZZHT. This code is also very accurate for in-
plane stresses and transverse shear stresses in the core region. This code was used for the

optimization portion of this research.

3.4 Conclusions

The LZZHT-3 model is the most accurate in predicting stresses in sandwich panels
and FSDT is the least accurate. HSDT shows a slight improvement in accuracy over FSDT
for predicting in-plane stresses, however, neither theory can predict transverse stresses. An

improvement in accuracy results from these models as the aspect ratio increases. For the

63

cases examined the ZZHT and LZZHT-1 models produced comparable results which were
also in close agreement with the elasticity solution for in-plane stresses. For transverse
stresses both the ZZHT and LZZHT-l models predicted the maximum shear stress in the
core very well but were not as accurate in the face sheets. The results illustrated an
increase in accuracy is obtained by increasing the number of sublaminates used through
the thickness. This is possible with the LZZHT theory. Also the results obtained from
analysis of layup 2 indicate the zigzag theories perform very well for unsymmetric and
symmetric layups.

i The ZZHT model is superior to the LZZHT-1 model in computational efficiency.
The best compromise of solution accuracy and computation efficiency is given by the
ZZHT finite element model, thus it will be used in the optimization portion of this

research.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

FSDT .
HSDT 1
3. €
IE 1
3 :
E I
3 €
m. ......... ‘1
.A ............. - . - . - 44 RAIL--. A . AAAA J
10. 20. 50. 100. 200. 500. 1000.
‘1
h
Figure (3.3) - Normalized transverse center deflection versus aspect ratio,
ESL theories, layup 1.
r a I I I - 3
1.035 f 1
: LZZHT-1 ‘
1.03 ;-_ ....... LZZHT-3 1,
’ . ----- ZZHT 3
.§‘ 1.025 ,- 1
.2 . 1
a ’ i
2% 1.02 :- -‘
3 1.015 E
1.01 f
l d
1.005L ----------- -—-—--~—-—---—-—---—-i
. . I . I 1" - ;:'#3
10. 20. 50. 100. 200. 500. 1000.
2
h

Figure (3.4) - Normalized transverse center deflection versus aspect ratio, zig zag

theories, layup 1.

65

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

-50 -202 20

 

 

 

 

Elasticity
- ----- - ZZHT
---------- LZZHT-1 —
— — — LZZHT-3

 

 

 

 

 

 

 

 

 

 

o4 - L A A A - A A 4 A A a A A
.

 

Figure (3.6) - Normal stress vs. normalized thickness for zig-zag theories, alh=4 layup l.

66

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Elasticity
————— FSDT
- - — - HSDT _
'3 T.
l
0"330 A -2202 A 10 I 20

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Elasticity
- ----- — ZZHT
---------- ZZHT-1 ,
— — — LZZHT-3
0430* 2120 A-ioé 0 U 21.02 A A ‘20 A I
Cy),

Figure (3.8) - Normal stress vs. normalized thickness for zig-zag theories, alh=4,
layup 1.

67

 

0.5 V'fiI—i i fifi'V'fi "ifi*fifv vvvvvvvvv
I

A

 

 

 

 

 

 

 

 

 

 

 

 

__ Elasticity
————— FSDT
- — — - HSDT '
0'4 2 4 6 8A A 10 A A A1A2A A 14
1,,

Figure (3.9) - Shear stress vs. normalized thickness for ESL theories, alh=4 layup 1.

 

0.5 ] ,--,--L,----T---,3-II-

 

 

 

 

 

S
__ Elasticrty

 

 

 

 

 

 

 

- ----- - ZZHT —
---------- LZZHT-1
- - — - LZZHT-3
0'4 0 2A *4 A6ATAAsA A 10A++A1A2AI ArirT
rxy

Figure (3. 10) - Shear stress vs. normalized thickness for zig-zag theories, alh=4
layup 1.

68

 

0,5 1—-—-—====e_ - e+ -

_ . - I," - -
0.4 _ 7

 

 

 

 

 

 

 

 

 

 

 

 

0.3 . __ Elasticity ‘
_ - ----- - ZZHT J
0'2 ---------- LZZHT-l
0.1? — — - LZZHT-3 . .
0 1
fi l
-0.1 t ;
-0.2 t
-0.3 r
I
-o.4 J I; a
_ . fl - . - A - - - - L I L - - - L
0 5 0.2 0.4 0.6 0.8 l
Tyz

Figure (3.11) - Shear stress vs. normalized thickness for zig-zag theories, alh=4 layup 1.

 

0.5

0.4
0.3 . __ Elasticity

 

 

 

 

V

0.2 .......... LZZHT-1 i
0.1 . _ — — LZZHT-3 J

 

 

 

 

--0.1 ~ .
-0.2 ~
-0 .3 _
-0.4
-0.5

 

 

 

 

 

 

Figure (3.12) - Shear stress vs. normalized thickness for zig-zag theories, a/hfl layup 1.

69

 

0.5- 4
T

‘ __ Elasticity
‘ ————— FSDT

\ I -—-- HSDT
\\|
L -

 

 

 

 

 

 

 

 

 

 

 

 

--M-.A-\P --r--2---
03100 -50 0 50 100

GXX

 

Figure (3.13) - Normal stress vs. normalized thickness for ESL theories, alh=10,

 

 

 

layup 1.
0.5 - 1 ~ . - - . - I .....
Elasticity
— ----- - ZZHT
---------- LZZHT-1 J

 

 

 

 

 

 

 

 

A A 4 AAAAA A A A A A A A A A A k AAAAAA A

 

Figure (3.14) - Normal stress vs. normalized thickness for zig-zag theories,
alh=10 layup 1.

7O

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(L5 . . .. - . I - 2. j
\ \
\\
t\
\ lflmflkky
( ----- FSDT ._
--—- HSDT
‘ v —
\\\1
‘ \
0'4 -107 -50 A A 50 A A 100

Figure (3.15) - Normal stress vs. normalized thickness for ESL theories, alh=10,

0.5

0.4 ‘

layup 1.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

__ Elasticity

— ----- — ZZHT

---------- LZZHT-l '—

— — — LZZHT-3
-100A -50 A A A AsAoA A A A100

Figure (3.16) - Normal stress vs. normalized thickness for zig-zag theories, alh=10

layup l.

 

71

 

0.5- -._--LI--I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

l
I
1
fi In ______, Elasticity
ll ----- FSDT
; — — — - HSDT '—
l.
l
.5
LA--. - . -2-.A---.--
0'4 10 20 30 40 50
Tn,

Figure (3.17) - Shear stress vs. normalized thickness for ESL theories, alh=10 layup 1.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.5 m.
____Elasticity
- ----- -ZZHT "
---------- LZZHT-1
1 ——-LZZHT-3
0'4 A AAlAOTAAAZAoAAAHBAoI Aé4AoAT+A50
In

Figure (3. 18) - Shear stress vs. normalized thickness for zig-zag theories, a/h=10
layup 1.

0.5
0.4

0.3-

0.2

0.1 ~

-0.1 -
-0.2 ~
-0.3 -

-0.4
-0.5

72

 

 

Elasticity
- ----- - ZZHT
---------- LZZHT-1
— — — LZZHT-3

T

 

 

 

 

 

 

 

 

 

Figure (3.19) - Shear stress vs. normalized thickness for zig-zag theories, alh=10

 

 

 

 

 

 

 

 

 

 

 

 

 

layup 1.
0,5 ———.===; - a ~ 4 ~
0.4
0,3 » Elasticity
- ----- - ZZHT
0'2 ---------- LZZHT-l
0-1- — — — LZZHT-3
0
-0.1 -
-0.2 -
-0.3 -
-0.4
_ _.._-===—A-AA444AAA-A4A4
0'5 0 0.25 0.5 0.75
AtXZ

Figure (3.20) - Shear stress vs. normalized thickness for zig-zag theories, alh=10

layup 1.

73

 

 

 

 

 

 

 

 

 

 

 

 

 

0.5vvvv+vvvvvv..-e--..---.---
fil
U
j ‘ ‘
\Q. 1
\\
\ \
\\
‘ ,
___Elasticity 1
----- FSDT
-—--HSDT
0'4 -100 A A-sAo A .20 0

 

Figure (3.21) - Normal stress vs. normalized thickness for ESL theories, alh=10
layup 2.

 

().5..........,....,...,A..I

 

 

 

 

 

 

 

 

 

0.4 ‘ -20

 

Figure (3.22) - Normal stress vs. normalized thickness for zig-zag theories, alh=10
layup 2.

74

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.5 v 1
__ Elasticity —
————— FSDT
- — — - HSDT
0"1.100 A 750 A 0 A 50 A A 100
ny

Figure (3.23) - Normal stress vs. normalized thickness for ESL theories, alh=10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

layup 2.
0.5 Y +\ ' 2 2 m
\ __ Elasticity _
— ----- — ZZHT
---------- LZZHT-1
— — — LZZHT-3
0'4-1AooA A A .50 0 A 50 A T00
610'

Figure (3.24) - Normal stress vs. normalized thickness for zig-zag theories, alh=10

layup 2.

75

0.5 TAITIA --L--L, - A-

 

 

 

 

 

 

 

 

----- FSDT

 

 

 

 

 

 

 

0’4 A10 A20 30 A 40 AAso

Figure (3.25) - Shear stress vs. normalized thickness for ESL theories, alh=10 layup 2.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.5 I -2. A - -
i
___Elasticity
_ ..... —zzrrr
---------- LZZHT-l
———LZZHT-3
0'4 AA 10 AionAA I30 AAA4AoA A50

Figure (3.26) - Shear stress vs. normalized thickness for zig-zag theories, alh=10
layup 2.

76

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

asfifip . vvvvvvvvvv L-----,I--A
0.4 f
0 ,3. __ Elasticity l. ,
————— ZZHT ..

”-2‘ ---------- LZZHT-1 ii
0.1- — — — LZZHT-3 5.

0 -1
-0.1 »
-03.
-03-
-0.4
-0.5 0.25 Egg—A735 1A

1,,

Figure (3.27) - Shear stress vs. normalized thickness for zig-zag theories, alh=10

 

 

 

 

 

 

 

 
 

 

 

 

 

 

 

 

 

layup 2.
0.5 ' _' ‘ ‘__* i_ - _. _~ v ~ ~ 1 ~ . a 1 a I
0.4 \—
0.3 . Elasticity I
- ----- - ZZHT x.
0'2 l ---------- LZZHT-l t
0.1 b — — — LZZHT-3
0 'r
-0.1 ’ El 4
-0.2 ' 5'
1 el
.03 _ 5,
-0.4 __ i”
o s Méfid . 1 -
. 0.25 0.5 0.75 l 1.25 1.5 1.75 2
1X2

Figure (3.28) - Shear stress vs. normalized thickness for zig-zag theories, alh=10
layup 2.

Chapter 4

Optimal Design Laminated Composite Sandwich

Structures Using Genetic Algorithms.

4.1 Introduction

Laminated composite structure design lends itself quite openly to optimization
techniques. Ply thickness and stacking sequence are two Often optimized design variables
when considering a laminate design. Numerous optimization techniques are available, one
category of which is calculus-based or gradient-based methods. Gradient—based methods
can be classified as either direct or indirect [16]. Both approaches rely on the gradient of
the objective function providing necessary information to locate a local optimum. One
obvious drawback of this method is the gradient of the function must exist and be obtain-
able. Also, these methods seek out local optima, which may or may not be the global opti-
mum. Gradient based methods are clearly not well suited for finding singular optima
common in optimization of laminated composite structures.

Various algorithms exist which search finite design spaces. Evaluating an objective
function at all points in a given domain is obviously inefficient for large domains. Random
search techniques perform a random search Of the design space for optima, but are also

inefficient. Directed random search techniques are a more effective means of searching a

77

78

design space. Such methods may, for example, mimic nature in their search for an Opti-
mum state. Genetic algorithms (GA) are guided random searches which utilize informa-
tion obtained during the search for direction. This technique is relatively effective and
robust by comparison to other methods.

Nature’s means of dealing with existence and reproduction of all species is
described as “survival of the fittest”. Genetic algorithms draw from this phenomenon,
employing a similar philosophy when searching a design space for global optima. G-As
maintain one or more populations containing a number of possible designs, and perform
operations of selection, reproduction, crossover, mutation, and permutation of the coded
design variables.

Genetic algorithms differ intrinsicallyfrom traditional search and Optimization
techniques. As previously stated, traditional methods utilize derivatives whereas GAs
make use Of objective function data. Also, GAs work with a coded set of variables and not
the actual variables as do traditional methods. Another important distinction is that GAs
simultaneously search from numerous points in the design space while traditional methods
search from a single point. Thus it is apparent traditional methods may locate a local max-
imum and not the global maximum, depending on the initiation point of the search. GAS
can provide many “near optimal” designs which is advantageous to the designer. Finally,
GAs are governed by stochastic rules, and traditional methods are governed by determinis-

tic rules.

79

4.2 Previous Work

The design of composite laminates with design variables of ply thickness and fiber
orientation is a common continuous Optimization problem. With only these two design
variables present the optimization is an integer problem which can be linearized. Branch
and bound algorithms have been used to solve this type of problem once it has been linear-
ized as done by Haftka et al. [33]. In their study results for optimization of stacking
sequence Obtained via branch and bound methods were found to be within 0.5% of those
found using GAS. Genetic algorithms can optimize nonlinear problems without the need
for linearization. Inclusion of frequency and strength constraints, such as in buckling load
optimization problems, leads to a nonlinear problem.

Haftka et al. [24, 25, 27, 33, 39] have used GAs for determining an optimum lami-
nate stacking sequence which minimizes laminate thickness or weight subject to strength,
buckling, and ply contiguity constraints. In these studies, balanced, symmetric, simply
supported laminates comprised of 0, +/-45, and 90 degree plies of equal thickness were
subject to in-plane biaxial loading. To improve efficiency of the algorithm Haftka et al.
[24, 25, 27] utilized a genetic algorithm with memory, storing useful information about
past designs in a binary tree. Use of binary tree storage was shown to reduce the number
of analyses required by 30-40%. Another attempt at improved efficiency was a procedure
called local improvement which utilized information stored in the binary tree [24, 25].

GAs are gaining popularity and becoming more frequently used in a variety of
optimization problems. Applications include minimizing the weight of a stiffened com-

posite laminate with a hole [34], optimizing elastic support locations on a square plate to

80

maximize the first system eigenvalues [19], Optimizing a simple truss structure capable of
supporting given forces [28], optimizing the cross-section and minimizing the weight of
an automobile bumper and minimizing the weight of an automobile body panel with stress
and deflection constraints [46], and optimizing the placement of accouters on large space
structures [13].

Previously work has been done by Punch et al. [39] on laminated composite
beams. The objective was to design the beam for maximum energy absorption, such as is
required for tank armor or automobile bumpers. This study reported work designing a 24-
layer beam made of graphite-epoxy composite layers with clamped-clamped end condi-
tions and an applied point load at midspan. A thin layer (about 5.5% of the nominal ply
thickness) is placed at the top of each composite layer, so there are actually 48 layers in
the model. Each thin layer may be assigned the same material properties as the layer
immediately below it, or it may be assumed that the thin layer is compliant, with stiffness
properties three orders of magnitude less than the composite layers. The length-to-thick-
ness ratio of the beam is 50. The length of the beam is divided into 20 subsections (finite
elements), so there are 20x48 = 960 design elements. The GA must decide whether to
place a 0 degree ply (henceforth called material 1) or a 90 degree ply (henceforth called
material 2) in each of the 480 structural ply design elements, and whether or not to place a
compliant material in each of the 480 thin layer design elements.

In this application, a number of representations were used but the focus was on a
480-bit-string version which fully represents one half of the beam, and which was then

mirrored across the vertical centerline to create a full beam. Thus the search space exam-

ined by the GA was 2480, representing all possible combinations of material and compliant

 

81

material. The results of this work were quite promising. The GA generated a number of

unique designs, including designs with “cantilever” beams of compliant layers.

In aerospace structures, the use of laminated composite materials offers many
advantages over conventional materials such as aluminum. One particular advantage is the
ability to tailor the lamination scheme (or ply stacking sequence) to achieve desired struc-
tural response to a given loading situation. For example, during certain maneuvers, an air-
craft wing may be subjected to aerodynamic loadings which cause the wing to twist in an
aerodynamically undesirable manner. By taking advantage of the stretching-bending-
twisting coupling of laminated composite materials, more desirable responses to such
loads are possible.

During flight, aircraft wings experience aerodynamic loadings which cause bend-
ing and twisting to occur. The twisting load is due to a pressure differential across the air-
foil in which pressure is greater at the leading edge than at the trailing edge. Many aircraft
are designed and built with up to five degrees of initial wing twist or washout. Washout
means the plane of the wing twists forward as it extends from the fuselage. This helps
compensate for gust loads the wing experiences causing the wing to twist up and back,
adversely affecting lift. Washout is built into an aircraft to reshape the spanwise lift distri-
bution to approximate an ellipse and to prevent tip stall [41]. At a given lift coefficient the
lift distribution can be Optimized by correct choice of initial wing twist. In this study, for a
given load, a desired wing twist which is Opposite in direction from that caused by aerody-
namic loads is sought. This can be achieved using composite sandwich structures in which

inherent bending-twisting coupling allows tailoring of the twisting response.

82

4.3 Model Description

The model used is a rectangular, cantilever sandwich panel, intended as an ideal-
ized aircraft wing, see Figure (4.1). So as not to detract attention from the main objective
of the current research, no attempt was made to exactly model actual aircraft wing load-

ings or geometry.

4.3.1 Loading

The intent was to choose a loading which produced a bending and twisting
response in the rectangular sandwich panel roughly similar to that of an aircraft wing
while in service. This type of loading allows illustration of the ability to achieve a desired
twist in a rectangular sandwich plate by optimizing the plate configuration via use of
genetic algorithms. In the spanwise direction an elliptic lift distribution was used with the
maximum load being at the wing root as illustrated in Figure (4.1)c. This is considered a
good approximation of actual loading [18, 40, 41, 45, 49] in the spanwise direction. In the
chordwise direction the lift distribution depends on numerous factors such as airfoil shape
and angle of attack. A load decreasing linearly from leading to trailing edge was used (see
Figure (4.1) b) since this type of loading will produce a twisting response in the wing sim-
ilar to the response produced by actual wing loading, as desired. The total lift distribution

is given in Equation (4.1).

2 1/2
P(x.y) = 4.00420) = p,(1 {X1 -:—,) (4.1)

83
P0 was determined so as to give the desired total loading on the wing.

To determine the actual wing loading the aircraft weight was assumed to be 13,349
N (3,0001bf) and a design factor of 3.5g was used, where g is the acceleration due to grav-

ity. Thus the load each wing must carry is Ws=23,362 N, which is 1/2 the total design

weight. P0 can now be determined by integrating Equation (4.1) as follows:

2 1/2
_ -2 _x_
w, _ j; [2 190(1 CXI 02) dydx (4.2)

and since Ws is known, P0 is found to be:

 

8W
P0 = ‘ = 13, 021N (4.3)
can

4.3.2 Geometry

The intent was not to exactly model the geometry of an aircraft wing but to use as
a gage the gross dimensions of typical aircraft wings. Dimensions of the sandwich panel
were chosen based on historical data provided by Raymer [41]. Wing aspect ratio, as used
in industry, is defined as the ratio of wing span (8) to chord length (c), with wing span
being the total tip to tip dimension of the wing. Typical wing aspect ratios are 7.5 for jet
transport aircraft, 7.6 for single engine aircraft, 7.8 for twin engine aircraft, and as high as
9.2 for twin turboprop aircraft. Tail aspect ratios range from 6-10 for sail planes, 3-4 for
fighter aircraft, and 3-5 for other aircraft.

The wing dimensions used in this study were chosen to be similar to common,

 

84

non-commercial aircraft. The wingspan was considered to be 9.75 meters (32 ft) with an
aspect ratio s/c=8, or for one wing an aspect ratio s/2c=4. Fuselage width, assumed to be
0.122 meters (4 ft), was subtracted from the wing span dimension thus making the dimen-
sions of one wing a=4.27 meters (14 ft) by c=1.07 meters (3.51 ft), where a is the length of
the single wing (see Figure (4.2)).

The National Advisory Committee for Aeronautics (NACA), which is the

 

 

 

 

 

 

 

 

 

 

‘ \ \ C120’)
T l ‘
leading , 7 I F _ trailing
edge 4 edge
chordwise direction
——->
Figure (4.1) - (a) Elliptical spanwise Figure (fl-1? - (13) Linear ehordWise
load distribution approximation. load 915919990" approxrmatron.
(110‘)
root T ' I tip
+

 

Spanwise direction

Figure (4.1) - (c) Elliptical spanwise
load distribution approximation.

 

 

 

 

85

precursor to NASA, devised various numbering schemes for airfoil classification. The first
was a four digit system which was developed around 1932 [18]. An example of an airfoil
using this classification is the NACA2412. Five digit and other classification systems
arose in the late 19305 for representation of more complicated airfoils. These numbering
schemes contain information about an airfoil’s shape such as thickness to chord ratio. The
NACA652-415 airfoil is of the same series first used for the P-51 Mustang aircraft [45]. As
is indicated by the last two digits in the code, this airfoil has a thickness to chord ratio of
0.15. This is the ratio of thickness to chord used for the sandwich panel in this study, and

for this thickness to chord ratio the corresponding actual thickness of the wing is

t=0.1605m, see Figure (4.2).

 

Figure (4.2) - Sandwich panel dimensions in meters.

 

86

4.3.3 Material Properties

The wing consists of a sandwich panel comprised of two materials. Material prop-
erties are contained in Table (4.1). The outer reinforcing layers have the properties of an
orthotropic graphite/epoxy material and the core region exhibits the properties of an alu-

minum honeycomb material consistent with the theory of Gibson and Ashby [15].

Table (4.1) Sandwich Panel Material Properties.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Property Material 1 Material 2
(face sheets) (core)
E11(Pa) 0.1390E+12 0.1587E+06
E22(Pa) 0.9860E+10 0. 1587E+06
E33(Pa) 0.9860E+10 0.6900E+09
v12 0.300E+00 0.9900E+00
v23 0.373E+00 0.000E+00 1
v13 0.2128E-01 0.300e+00
012(Pa) 0.5240E+10 0.3933E+05
623(Pa) 0.3250E+10 0.1500E+09
013(Pa) 0.5240E-1-10 0.1500E-1-09

 

4.4 Optimization Problem

Structural optimization problems have been around for centuries. In aerospace

applications the goal of structural optimization is to minimize structure weight while

 

87

maintaining sufficient structural integrity. Various works on optimization of laminated
composite plates have been published (see for example [2, 3, 23, 25, 27, 34, 35, 50]).
Nearly all these cases require the laminates to be symmetric. In the current research, no
such restrictions are placed on the layup.

In the current research, a rectangular laminated composite sandwich panel (see
Figure (4.1) serves as an idealization of an aircraft wing and the lamination scheme is
sought which minimizes weight, maximizes in-plane stiffness, and achieves a desired
twisting response of the wing. To obtain an Optimal design the genetic algorithm
GALOPPS, developed by Eric Goodman at Michigan State University, is used in conjunc-
tion with the finite element code employing the ZZHT theory discussed in Chapters 2 and
3 (see Figure (4.3)). This finite element code-provides the necessary displacements and
stresses for each design.

This optimization problem is presented as a maximization. This is necessary for

 

Input:
(convergence condition, crossover rate, mutation rate, etc.)

L,

GALOPPS v3.0 information

1 analysis of the

(optimization engine) ‘ Cantilever Plate

----------'

output:
(Best Design)

 

Figure (4.3) - Schematic of Optimization Process.

 

 

 

88

problems Optimized by genetic algorithms which seek the design having a maximum fit-
ness. However, all optimization problems can be posed in both ways (as a maximization or
minimization) [17]. The optimization problem in the current research is classified as a

constrained multicriterion optimization problem which is posed as follows:

N I
maximize[f(i) = 4000 - chfnCt) - EEO-0]

subject to
OS 1516, 05x2516,

x
E (4.4)
2 9

NlF-t

5x35 0.148me4SO.169m,

and

<6,“ = 1.24x109Pa

01,0101

92,max< 02'" = 4.1x107Pa

where f is the fitness, cn are weighting coefficients, fn are design objectives, 1: are the
design variables, P,- are penalty functions, O’Lmax is the maximum stress in the fiber direc-
tion, and 02mm, is the maximum stress in the transverse direction. Values for cn are given
in Table (4.2). The design variables in Equation (4.1) are as follows: x1 and x2 are the
number of layers in the top and bottom face sheets, respectively, x3 is the fiber orientation

angle, and x4 is the core thickness.

 

89

Table (4.2) Weighting Coefficient Values.

 

c1 c2 c3 d1 d2 d3 k2 k3
50 50 -1000 10 20 100 500 1000

 

 

 

 

 

 

 

 

 

 

The design Objectivesfn in Equation (4.4) are comprised of wing twist, weight, and

in-plane stiffness as follows:

 

 

~ _ Wbaseline 2
f1(x) - (TOO—j (4.5)
not) 40:91:” (4.6)
~ t(~)-t esire 2
f3(X) = ( xtd .dd d) (4.7)

where W is the panel weight, a is the minimum of the two in-plane stiffness values, A“
and A22 from the extensional stiffness matrix, and t is the wing twist. These three quanti-

ties are normalized by values Obtained from a baseline case. This normalization provides
the same order of magnitude for each design objective. The baseline case consisted of a
symmetric sandwich panel having face sheets oriented such to produce a specially ortho-
tropic laminate. The lamination scheme for both the top and bottom face sheets was [0/90/

+45l-45]s and the core thickness was 0.158 meters.

The wing twist is defined as:

t(3't) = (WA -wB) (4.8)

where WA and “’3 are transverse deflections of the tip at the leading edge and trailing edge
(see Figure (4.1)). The desired twist, Idem-m1 , was chosen based on the twist of the baseline

panel which was 0.007 meters. The intent was to achieve a twist of —4 times the baseline

twist or tdesired = -0.028 m, which gives an angle of -1.5’ with respect to the horizontal.

Negative twist indicates the leading edge of the wing will deflect less vertically upward
than the trailing edge indicating the wing is twisting in the direction opposite that caused
by the aerodynamic loading.

The normalizing values used in Equations (4.5) and (4.6) are the weight and in-
plane stiffness, respectively, of the baseline sandwich panel. For the baseline laminate

Wbaseline = 60.03 kg, and abaseun = An = A22 = 1.223E+08 Pa.

Second-order polynomials were found to work best for representation of design
objectives in the fitness function. To achieve the desired amount of twist a quadratic func-
tion was used to represent the normalized difference in twist. The vertex of the function is
located at the desired twist value thus driving the twist toward the desired value by reward-
ing the fitness of the designs whose twist is near the desired value. As a matter of consis-
tency, quadratic functions were used for the normalized weight and in-plane stiffness
objectives of Equations (4.5) and (4.6), respectively.

Fitness is penalized when stress levels in a design exceed the material strengths in

91

the fiber and transverse directions. A penalty is also applied to the fitness for exceeding the

given maximum deflection. The penalty functions, P,- , in Equation (4.4) are as follows:

0.0 (we 5 W0)
P (i) = C(32) - w 2 (4.9)
‘ dill—”l

otherwise
“’0

0.0 (01mm s ol'")

 

 

 

132(3) =1 0 3: —o' 2 (4.10)
k2+d2[ l’max( 3‘ 1‘] otherwise
91
0.0 (62. m, s 62’“)
P3(5'c) =1 0 g _o . 2 (4.11)
k3+d3[ 2,max( 1):. 2 otherwise
92

 

where we is the transverse deflection of the wing tip at mid-chord, wo is the transverse tip
deflection of the baseline case, and d 1’ d2, d3, k2, and k3 are weighting coefficients whose

values are listed in Table (4.2).

The weighting coefficient values in Table (4.2) were chosen in part based on the
relative importance of each objective and penalty. Another factor influencing the choice of
weighting factor values is the ease with which the constraint is violated. In this study the

constraint on 62 was most easily violated. This is an important constraint since designs in

violation of this constraint are not valid. Consequently a large weighting factor was used

for this function thus assessing a large penalty to discourage the GA from choosing

92

designs in violation of this constraint. After many iterations of fitness function revisions
the values in Table (4.2) were chosen. These values enabled the GA to determine designs
which produced the desired amount of twist, had the maximum in-plane stiffness and min-

imum weight, and not did not exceed maximum allowable stress values.

4.4.1 Design Coding
Each design is represented by a coded string of binary numbers called a chromo-
some. The chromosome structure used in the current research is given in Figure (4.4). The

chromosome consists of a string of binary numbers in which different portions represent

 

 

 

 

 

 

 

Sandwich Cross-Section:
No. of layers
on top
(1 - 16)
No. of layers
K on bottom
(1 - 16)
Chromosome:
orientation of each layer
(6x16=96 bits each)
6 bits per layer 6 bits per layer
/ \ I \
Wllllllllllllflllllllllllll
, 4 bits , \ 6 bits / . 4 bits ,
# of layers Core thickness # of layers
(1 - 16) (0.148 -O.169) (1 - 16)
(bottom) (tap)
Figure (4.4) - Chromosome structure.

 

 

 

 

93

different physical characteristics of the design. The string contains 206 total bits giving a

2206

design space of . The first and last four bits of the chromosome represent the number

of layers on the bottom and top of the laminate, respectively, providing a range of l to 24 =
16 layers.

The next 96 bits inward represents the fiber orientation angle of each ply in the top
and bottom face sheets. Fiber orientation angles are given a range of 180° (-90" to 90°).

Each angle is represented by a 6-bit portion of the total 96 bits allotted for the top or bot-

‘ tom face sheets. With six bits representing each angle the 180° increment is divided into 26
= 64 divisions, thus producing an angular resolution of 2.8125°.

The middle six bits in the chromosome structure determine the core thickness
which is confined to a range of 0.148 to 0.169 meters. This range was chosen to maintain

the thickness to chord ratio of the panel near 0.15 as discussed previously.

4.5 Genetic Algorithm Topologies

A comparison is made between four different GA topologies or structures. The
four different t0pologies are the single node, ring, injection, and hybrid architectures. In
all four cases results were obtained for five runs of each topology. Each of the five runs
contained a different random seed to generate a different initial population. The total pop-
ulation size for each case was 140 and each run terminated after 150 generations which
provided sufficient time for results to converge. A generation is defined as the completion
of a cycle in which two “parent” strings are selected from the population to be mated.

They then reproduce, thus producing two offspring who are now a part of the population.

94

The GA uses a tournament selection method in which individuals compete (as the name
tournament suggests) for the chance to mate. All four GA topologies use a two point
crossover with a probability of 0.75. Mutation rate was 0.001. A description of the differ-

ent topologies follows.

4.5.1 Single Node Architecture GA
The least complicated GA structure is the single node. While this structure is sim-
ple in appearance it is no less robust than the other structures. The single node GA con-

tains the entire population at a single node where all genetic operations are performed.

4.5.2 Ring Architecture GA

The ring structure GA contains seven nodes arranged in a ring configuration as

 

node 4
node 3
node 5
7 sub-populations
od 2
high resolution 0 n e
node 6
node 1
node 0
Figure (4.5) - Ring GA Architecture.

 

 

 

95

shown in Figure (4.5). Each node is an individual, self sufficient, sub-population similar to
the single node case. A population size of 20 individuals at each of the seven nodes gives
the same total population size as the other cases of 140. With the same population size for
all topologies a direct comparison between structures can be made. The same total number
of generations (150) is used for the ring structure with 5 generations per cycle and 30
cycles totaling 150 generations.

On the node level, the ring structure works similar to the single node case since
each node simultaneously searches for the best fit individual, performing all the genetic
operations of selection, reproduction, crossover, and mutation. The difference lies in the
fact that each suprpulation performs 5 generations after which the best fit individual
(chromosome) is copied to the adjacent subpopulation.This chain of events constitutes one
cycle. This continues for 30 cycles for a total of 150 generations. The chromosome
exchange occurs in one direction only, in this case from node 0 to 1, node 1 to 2, node 2 to
3, etc., as indicated by the arrows in Figure (4.5). This structure was devised to improve
GA performance by accelerating convergence. One inherent characteristic of this structure
is the tendency of a single design with high fitness to dominate the entire population in

later generations.

4.5.3 Injection Architecture GA

In previous work with parallel genetic algorithms for design, Punch, et al. [39] exam-
ined a number of exchange topologies in attempts to discover the most appropriate topolo-
gies for the laminated composite beam design example. In so doing, a new topology called

injection island parallelism, or iiGA, was created. A similar injection structure is used in

96

 

Top level node 0 node 1 node 2 node 3
4 bit (low)
resolution

Middle level
5 bit (moderate)
resolution

node 5

 

Bottom level

6 bit (full) resolution
node 6

Figure (4.6) - Injection GA Architecture, 7 subpopulations.

 

 

 

the current research.

TWO interesting aspects of an iiGA are its migration rules and the heterogeneous nature
of its populations. In most GA parallel processing approaches, multiple subpopulations
are used to generate solutions, and these solutions are exchanged between subpopulations
in hopes of acquiring better solutions. Two common properties of most such approaches
are first all suprpulations use the same GA representation of the problem, and second the
exchanges are typically two-way meaning each subpopulation gives and receives solutions
from other subpopulations.

The iiGA used in [39] as well as in the current research differs from both of these situ-

ations. First, the subpopulations are organized in a hierarchical exchange which is one way

 

97

and tree-like; that is, solutions are passed from children to parent subpopulations, but no
exchanges occur between subpopulations on the same level (see Figure (4.6)). Second,
various subpopulations are allowed to work with different representations of the same
problem. The form of the representation typically follows the hierarchical exchange topol-
ogy, with “child” subpopulations using a less refined or more abstract representation of the
problem, which is then passed to a more detailed representation in the parent for refine-
ment. This allows subpopulations to search smaller, but more abstract, search spaces for
areas of possible interest, which can then be refined by parent subpopulations using a
more detailed and more focused search. The top level or leaf nodes are autonomous, so
their search experiences no influence from local optima discovered in other nodes on the
way to the global optimum.

In previous work, the iiGA approach was applied to a laminated composite beam
design problem [39]. The result of the iiGA approach to this problem showed dramatic
improvement both in the energy absorption characteristics of the beams and the computa-
tional time required to find these solutions. Lessons learned from using iiGAs for the
beam design in [39] are extended to the more complicated domain of composite sandwich
panel design in the current study.

The 7-subpopulation tree shown in Figure (4.6) illustrates the topology of the iiGA
used in the current research. The injection structure is capable of representing a design
with different levels of refinement to decrease computational costs and improve efficiency.
Three levels of refinement are used for representation of the fiber orientation angle. As
mentioned previously, 6 bits are allotted for representation of the fiber orientation angle of

each ply. In the top level of the injection structure (low resolution) orientation angles are

98

represented by only 4 bits, where in the middle level (moderate resolution) a 5 bit repre-
sentation is used, and finally in the bottom level the full resolution of 6 bits is used.

The iiGA has the same population size of 140 as do the other structures which is
obtained from a population size of 20 per each of the seven nodes. Each subpopulation
simultaneously performs all genetic operations, increasing the fitness of the best individ-
ual. As with the ring structure, the iiGA performs 30 cycles consisting of 5 generations
each for a total of 150 generations. After each cycle the four subpopulations in the top
level inject the best fit individual into the two middle subpopulations and the middle sub-
populations inject the best fit individuals into the bottom level subpopulation where the
designs are further refined. The top level designs receive no data from outside the popula-
tion and the bottom level does not inject data'into any other population. Sharing of data is

one way as indicated by the arrows in Figure (4.6).

4.5.4 Hybrid Architecture GA

The final GA structure used is a hybrid structure which combines the ring and
injection topologies (see Figure (4.7)). This structure was developed in an attempt to
achieve improved performance over the ring and injection GAs. In the injection structure
the moderate and low resolution subpopulations very quickly converge to a low fitness and
thereafter do not contribute to the performance of the GA. The intent is to take advantage
of the initial rapid increase in fitness of the injection GA and then finish with the ring
structure to achieve higher fitness.

This hybrid approach uses the same 7 subpopulation injection structure shown in

99

Figure (4.6) for the first one half of the total generations with the orientation angle resolu-
tions exactly the same as the iiGA discussed previously. For the last half of the generations
the GA structure is converted into the same ring structure discussed previously with each

node having full resolution.

4.6 Finite Element Model

The wing was modeled as a rectangular sandwich panel which was fixed at one end
and free on all other edges. The rectangular domain was discretized into a mesh of ele-
ments, the size of which was chosen based on a convergence study. The convergence study

determined the coarsest mesh possible which gives sufficiently converged results (see

 

 

"0°18 1 2 3 node 4

resolution

4
resolution 4 5

5
resolution

6 6

First Stage Second Stage
First 1/2 generations. Final 1/2 generations

Figure (4.7) - Hybrid GA Architecture.

 

 

 

100
Table (4.3)). It is imperative to use a mesh with the fewest number of elements since com-
putational expense is directly affected by the number of elements. In an optimization prob-
lem tens of thousands of calculations are performed so even a seemingly small savings in
computational time can amount to a large savings over the course of an optimization run.

The quantities we and W]; in Table (4.3) are the transverse deflections of the free end (or

tip) of the plate at the leading edge and mid-chord respectively.

Table (4.3) Mesh Convergence Data.

 

 

 

 

Mesh O,“ 6” Oxy WC W3
(x109) (x10'8) (1110”)
2x4 0.5291 0.0980 ‘ -O.1906 0.2280 0.2314
2x8 0.5774 0.1 105 .0.1370 0.2335 0.2300
2x12 0.5737 0.1 165 0.1343 0.2341 0.2306

 

 

 

 

 

 

 

A mesh containing two 4-noded plate elements in the chordwise direction (the
short dimension of the wing) and eight elements in the spanwise direction was found to be
sufficient. This mesh provided in-plane stresses which were converged to within 5% and
transverse tip displacement which were converged to within 3%.

The displacements and stresses used by the optimization algorithm were deter-
mined by the FORTRAN finite elemcnt code employing the ZZHT model discussed in

Chapters 2 and 3.

101

4.7 Results and Discussion

4.7.1 GA Results

Five runs of each GA topology were performed and results compared. Table (4.5)
contains results for each run of each structure. To obtain different results for each run the
different random seed was changed in the input file. Changing the random seed assures a
different initial population is generated for each run. Table (4.5) contains the fitness of the
best fit individual (design) from each run of each structure, the subpopulation number
from which this individual came, and the average of the best fitness values.

Each structure achieved nearly the same level of fitness with the maximum differ-
ence in average fitness being 0.3%. This indicates the ability of all structures to solve the
optimization problem sufficiently. The individual with the highest absolute fitness of
4,140.25 came from node 6 of run 4 of the ring structure. The ring and hybrid structures
produced the same average fitness of 4,130.80 indicating the performance of these two
structures is superior to the injection and single node structures. The injection structure
gave the lowest performance which is due to the application specific nature of the injection
structure. For problems with much larger design spaces the injection structure should out-
perform at least the single node structure and likely would out-perform the ring structure
as well. While one GA topology can be termed the best for this study, all structures per-
formed sufficiently well.

Figures (4.8) through (4.11) contain results of fitness versus generation number for
each of the topologies. In all cases results are converged long before the maximum number

of generations (150) is reached. In the case of the single node structure the fitness is

.. aup

 

102

converged to within 2% of the maximum by generation 10 and to within 1% by generation
number 20. Figure (4.8) contains the single node results.

The ring structure results are given in Figure (4.9). Seven sets of data points are
plotted, one for each of the seven subpopulations or nodes. Node 5 initially increases in fit-
ness the most rapidly while node 1 increases the slowest. Node 4 produced the individual
(design) with the highest absolute fitness. At node 4 fitness increases to within 4% of the
maximum by generation 10 and to within 2% of the maximum by generation 20 which is a
slightly slower convergence rate than the single node case but on the same order.

Figure (4.10) shows results from each of the seven nodes of the injection structure.
As discussed previously nodes 0 through 3 use the lowest resolution representation of the
design, nodes 4 and 5 use moderate resolution, and node 6 uses the full resolution (see
Figure (4.6)). The design fitness from node 6 increases to within 4% of the maximum by
generation 10 and to 0.7% by generation 20 making the injection structure the quickest of
the four structures to converge.

The hybrid structure’s convergence performance (see Figure (4.11)) is similar to
the injection structure since for the first 75 generations the hybrid structure uses the injec—
tion structure (see Figure (4.7)). From generation 75 to 150 the hybrid structure uses the
ring architecture. Once the structure changes from injection to ring at generation 75 the fit-
ness of all the nodes converge quickly as can be seen in Figure (4.11). As expected the
lower resolution nodes begin to converge to a low fitness value when the injection struc-
ture is used. This is evident for the hybrid structure (see Figure (4.11)) more so than for
the injection structure results of Figure (4.10) which is due to a different initial population.

In both the injection and hybrid structures the final fitness converges to similar maximum

103

values and the fitness convergence rate of node 6 is very similar. For the hybrid structure
the fitness is within 7% of the maximum by generation 10 and within 1% by generation
20.

The various quantities in the fitness function (Equation (4.4)) are presented in Fig-
ures (4.12) through (4.15). These figures illustrate how the GA works to achieve maxi-
mum fitness while satisfying all constraints. Seven quantities are plotted for each GA
structure. These quantities show how the design is evolving physically as the GA conducts
its search of the design space.

Fitness is the first quantity plotted. The fitness is plotted for the subpopulation
from which the best fit individual came (see Table (4.5)). This plot exhibits the increase in
fitness as the generation number increases, which is shown in Figures (4.8) through (4.11).

The next quantity plotted, labeled wt*, is the normalized weight objective multi-

plied by the weighting factor and is calculated as follows:

W - 2
w = c,(————;;E;;") (4.12)

where c], Wbaselm, and W are the previously defined quantities. These graphs show how
the structure weight evolves during the optimization process with final weights for each
structure given in Table (4.4). The single node structure determined the lightest weight

design at 66.19 kg and the hybrid structure determined the heaviest design at 71.76 kg.

 

104

Table (4.4) Results for Best Design of Each GA Structure.

 

 

 

 

 

 

GA S Wt* Sufi»? twistak
mm” (actual weight kg) (an,a22,(*10'8)) (actual twist rn)

{fight 45.35 86.49 0.0044
(66.19) (2.1 10,2.133) (0.0279)

ring 44.00 96.84 0.0249
(68.22) (2363,2390) (0.0273)

injection 45.04 82.58 0.0380
(66.65) (2201,2015) (0.0209)

hybrid 41.83 97.56 0.0120
(71.76) (2427,2380) (0.0277)

 

 

 

 

 

 

The next quantity stifP" is the stiffness objective term in the fitness function, calcu-

lated from the following equation:

stiff’l‘ = 62(fl (4.13)

abaseline

where all quantities are as previously defined. The in-plane stiffness of the wing structure
is maximized to maximize the buckling strength of the structure. As these figures illus-
trate, stiffness is maximized as the generation number (and fitness) increases. In all cases
in-plane stiffness converges to similar values. The hybrid structure determined the design
having the highest in-plane stiffness and the injection structure determined the design with

the lowest in-plane stiffness (see Table (4.4)). The difference in stiffness values for each

 

105

GA structure is due to the fact the total fitness is a weighted balance of all design objec-
tives which allows a design with the similar fitness to have different values for each design
objective.

The next quantity plotted, twist”, is the relative difference in twist of the wing,
where twist is defined in Equation (4.8). The twist* quantity is calculated from the follow-

ing equation:

 

t (x) - ’desired) (4.14)

twist* = (

tdesired
As the graph of this quantity indicates twist? converges to zero as the generation number
(and fitness) increases indicating the wing structure has achieved the desired amount of
twist. This was the most difficult objective to achieve so the largest weighting was given to
it. The single node structure produced a design having the “best” twist, deviating only
0.4% from the desired twist. The largest deviation, 3.8%, from the desired twist was from
the design produced by the injection structure.

The final three quantities plotted are the constraint or penalty terms: delta“,
sigma1*, and sigma2*. The delta* quantity is the transverse tip deflection penalty calcu-
lated from Equation (4.9). The plot of this quantity shows all best designs of each genera-
tion do not exceed the maximum allowable tip deflection. The stress constraints are
plotted as the quantities sigma] * and sigma2* where sigma1* is the stress constraint in the
fiber direction calculated from Equation (4.10), and sigma2* is the stress constraint in the

transverse direction calculated from Equation (4.11). It should be noted that these three

 

106

constraints are satisfied for the best design of each generation throughout the 150 genera-
tions as indicated by the final three graphs of figures (4.12)-(4.15). Thus no penalty is
assessed to the fitness for these designs. In the remainder of each population there may
exist many designs which violate these constraints however their fitness is not high enough
to make them the “best” individual of a generation. The design space does not exclude
such designs and thus the population diversity is maintained.

To ensure the designs found were near global optima an optimization was per-
formed with an initial population of 1,000. A larger initial population increases the num-
ber of designs from which the GA can choose which increases the chances of the best
design being found by the GA. From this optimization run fitness changed less drastically
in early generations and the fitness of the best individual was 4,123.94 which is on the
same order as the results obtained from using a population size of 140 individuals. This
validates the use of the smaller population which is preferable to use since it conserves

computational resources.

4.7.2 Stiffness Matrix Results

Fiber orientation angles, number of plies, and core thickness were design variables
in this optimization project. Over the course of the optimization these variables are
changed by the genetic algorithm in the search for the optimum design. Altering these
variables affects the structure’s response since the overall laminate stiffness matrix, often
referred to as the [ABD] matrix [9,20]. In preceding sections the evolution of the design
was traced by examining the weight, stiffness, and twist for each generation. Another

means of illustrating the evolution of the design is to examine the various terms of the

 

107

[ABD] stiffness matrix. Figures (4.16) through (4.18) display the behavior of the stiffness
matrix as the genetic algorithm seeks out the optimum design. Results are given for run 4,
node 6, of the ring ga structure which produced the individual (design) with the highest fit-

11688.

The [ABD] matrix as well as each of the [A], [B]=[B]T, and [D] matrices are sym-
metric so only the upper diagonal half of each matrix needs to be considered. In Figure
(4.16) [A] matrix values are given for the best individual of each generation. This portion
of the [ABD] matrix is called the extensional stiffness matrix and it relates in-plane stress
resultants to in-plane mid-surface strains (shear and normal) [48]. One design objective
was to maximize the in-plane stiffness or the A11 and A22 terms of the [A] matrix. These
terms relate the in-plane normal stress resultants to in-plane normal strains. Maximizing
these two terms increases the buckling strength of the wing (plate). As Figure (4. 16) indi-
cates these terms are maximized during the optimization, as desired, making these two
terms the largest in the [A] matrix, A11=2.363e+08 and A22=2.390e+08. The A12 and A66
terms exhibited similar behavior to each other during the optimization. Tension-shear cou-

pling terms A16 and A26 are non—zero indicating coupling occurs between in-plane normal
stress resultants and in—plane shear strain. Their final values are A16 =-7.007e+07 and
A26=-2.276e+07.

The [B] matrix portion of the [ABD] stiffness matrix is called the bending-stretch-
ing coupling stiffness. Values for the [B] matrix are given in Figure (4.17). The [B] matrix,
in part, determines the twisting response produced by the wing, which was the main objec-
tive of this optimization. Bending-shearing and extension-twist coupling is present since

the 816 311d 826 terms 31°C non-zero [48], B16=-8.873C+05 and B26='2°159e+05' The Bzz

108
term had the largest absolute value at B22=-7.317e+06.

The [D] matrix, called the flexural stiffness matrix, relates stress couples to curva-
tures [48]. The behavior of this matrix is very similar to the behavior of the [A] matrix

with the D11 and D22 terms being the largest at l.321e+06 and 1.341e+06, respectively.
Values for the coupling terms are D16=-3.935e+05 and D26=-1.275e+05. These terms
determine the bending—twisting coupling in the wing making them instrumental in achiev-
ing the desired amount of wing twist since a bending load was applied to the wing and a

predetermined twist response was sought. As was the case with the [A] matrix, the D12
and D66 terms behaved similarly throughout the optimization, achieving final values of

2.294e+05 and 2.725e+05 respectively.

4.7.3 Layups

The lamination schemes corresponding to the best design (the design with the
highest fitness) as determined by each of the GA structures are contained in Tables (4.6)
through (4.9). As indicated by Table (4.5) the ring structure produced the design with the
highest fitness at 4,140.25 and the injection structure the lowest at 4,096.00. The design
produced by the single node was the lightest of the four cases weighing 66.19 kg while the
heaviest design came from the hybrid GA weighing 71.76 kg. The hybrid GA also pro-
duced the design having the highest in-plane stiffness value with All being 2.4273e+08
Pa.

The layups produced by each GA structure are quite different. Genetic algorithms
typically locate several optimal solutions to a design problem which is a desirable situa-

tion in many instances. This allows the designer many options to choose from since one

 

109

design may be more desirable for reasons not accounted for in the fitness function such as
aesthetics. All layups produced have no more than two contiguous plies with the same
fiber orientation. This helps to alleviate matrix cracking which is a concern when a lami-

nate contains several adjacent plies with the same fiber orientation.

4.8 Conclusions

The goal of this optimization problem was to illustrate the ability of a laminated
sandwich structure (an aircraft wing) to achieve a specific twisting response which was
opposite the direction the loading was applied. This objective was obtained while structure
weight was minimized and in-plane stiffness maximized, subject to constraints on tip
deflection and maximum stress. Optimum designs were determined satisfying all desired
design criteria using genetic algorithms.

Four different GA structures were used to determine the best design: single node,
ring, injection, and hybrid structures. Five runs were performed for each structure for
determination of the maximum fitness and for evaluation of the average performance of
each structure so the four structures could be compared. The maximum fitness of 4,140.25
was produced by run 4, node 6, of the ring structure and the lowest fitness of 4,096.00
came from the injection structure. Both the ring and hybrid architectures produced the
highest average fitness of 4,130.8 and the lowest average fitness came from the injection
structure. The maximum difference in average fitness was only 0.3%, thus all architectures
performed comparably in regards to average fitness. Another measure of performance was

the convergence rate of each structure. The injection GA converged the most rapidly of the

 

110

four structures and the ring GA the slowest. Again, the four structures exhibited a very
small difference in convergence rates with the injection converging to 0.7% of the maxi-
mum fitness by generation 20, the hybrid and single node structures converging to within
1% and the ring structure within 2% by generation 20.

The single node structure produced the design with the lowest weight of 66.19 kg
as well as the design having the “best” twist, -0.0279 m, which is only 0.4% less than the
desired twist of -0.028 m. The design having the maximum in-plane stiffness,

A11=2.427e+08 Pa, was produced by the hybrid structure.

The population size of 140 used in all cases was found to be sufficient by compar-
ing results with a case using a population size of 1,000. Maximum fitness from the case
using the larger population size was 4,123.94 which is 0.4% less than the results obtained
by the ring structure using a population size of 140. This indicates the designs found are
near the global optima.

The lamination sequence of the best design from each GA structure was included
in the results section. The four lamination schemes are quite different from each other,
however, they are all valid designs which satisfy all conditions and constraints of the prob-
lem. Each design provides a compromise between stiffness and weight. Equal weighting
was given to both objectives. If maximum stiffness was more important than minimum
weight, a larger weighting could be given to the stiffness objective to drive the design
toward higher stiffness. The converse is also possible. Having a variety of optimal designs
to choose from allows the designer greater flexibility. Genetic algorithms, in many cases,
can provide several near optimal solutions to an optimal design problem, provided the

problem lends itself accordingly.

 

Q
E
'a
5
L

 

111

The various terms in the [ABD] stiffness matrix were examined over the course of

the Optimization. The A and D matrices behaved similarly with A“, D“, A22, and D22

terms having maximum values. The presence of extension-shear coupling was illustrated

by non-zero A16 and A26 terms and bending-twisting coupling was present since the D16
and D26 terms were non-zero. 'l‘wisting-stretching and bending-shear coupling is indicated
by non-zero B16 and B26 terms. The B22 term had the largest absolute value in the B

matrix. The 16 and 26 terms of the D matrix were instrumental in achieving the desired
“opposite” twisting response in the wing.

One concern which was not directly addressed in this problem is design manufac-
turability. Due to the highly unsymmetric lamination scheme necessary to achieve the
desired twisting response in the wing, warping will occur in the manufacturing process as
the laminate cures. Manufacturing of this structure will need to be done in such a way to
minimize (or eliminate) the warping due to curing. The shape of the wing after curing can
be determined by applying the appropriate thermal loading which will produce the same
response in the wing. With the post-cured shape known it should be possible to manufac-
ture the laminate in a shape which will correct for the warping by using an appropriately

shaped mandrel.

4.9 Recommendations

The motivation for this design problem came from the desire for aircraft wings to
adjust their angle of attack when experiencing gust loads. Gust loading causes the wing to

twist up and back which is detrimental to the wings lift. The desired response to such

 

112

loadings is for the wing to twist in the opposite direction. This research shows such a
response is possible with laminated sandwich structures. Specific design data was
idealized, i.e. wing shape, sweep back, aircraft size and weight, actual wing loadings,
desired wing twist, etc., to devise a general case which illustrates the ability to solve such
a problem using an experimental finite element plate model with genetic algorithms. This
problem can be extended and refined as desired to account for more realistic designs when
specific design information is known. Also a more refined failure criteria could be
incorporated such as the Tsai-Hill, Tsai-Wu, Chamis, and Hoffman interactive failure
criteria [48]. Design manufacturability and cost constraints could also be included in the

fitness function.

 

113

Table (4.5) Optimization Result From Various GA Structures.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Subpopulation Average
Run Number Best Fitness number Best
Fitness
= J =: 4: =
1 4.121654e+03 1
5mm 3 4.131813e+03 1 4.1215e+3
(popsrze 140,
‘50 gene’a‘i‘ms) 4 4.114829e+03 1
5 4.118105e+03 1
l 4.132970e+03 6
Ring Structure 2 4.119689e-1-03 3
(7 “F’des' 3 4.138185e+03 3 4.1308e-1-3
popsrze 20,
150 generations) 4 4.14025164'03 6
5 4.123152e+03 0
1 4.126204e+03 6
Injection 2 4.096004e+03 6
Structure
(7 nodes, 3 4.126848e+03 6 4.1165e-i-3
popsize 20
150 generations) 4 41074648403 5
5 4.125922e+03 6
l 4.118339e+03 5
Hybrid Structure 2 4.121376e+03 6
(7 “9d“: 3 4.1379l4e+03 4 4.1308e+3
popsrze 20
150 generations) 4 4.136995e+03 0
1 5 4.139252e+03 0

 

 

 

 

 

l_‘-l"-Il'l

 

 

thess

114

 

4200.0 . t

 

 

I

4000.0

1

3800.0

 

3600.0

'x‘

 

l J l

 

 

3400.0 ‘
0 50 100
Generation Number

Figure (4.8) - Fitness versus generation number for single node GA structure.

150

 

Fitness

115

 

4200.0 f t w 1

3700.0

 

 

 

 

 

150

3200.0 cl -— node 3 _
- -— node 4
~— node 5
--- node 6
2700.0 - 4 ‘ ' -
0 50 100
Generation Number

Figure (4.9) - Fitness versus generation number for ring GA structure.

 

thess

116

 

 

 

 

 

 

150.0

4200.0 , I .
g F-——" -- -- -- ---.---..--_ 2
, -.- J
"i i [L] 1
1 .l,’
1 i:
3700.0 ,5; J, -
.1
l
5 ' --- node 0
‘ — node 1
--- node 2
3200.0 3 - - 4 node 3 _
9- - —0 node 4
< G—o node 5
—-- node 6
2700.0 ‘ J
0.0 50.0 100.0
Generation Number

Figure (4.10) - Fitness versus generation number for injection GA structure.

Ffiness

117

 

 

 

 

 

 

Generation Number

Figure (4.11) - Fitness versus generation number for hybrid GA structure.

150

4200.0 - i . t
,3. = -A _‘———':‘—- '7 - —
. ---A-’ ‘ I
,- ' .5.._.._&”.:. £1" £1
r ............ Lu-'
#3 ’3’-
l
1' 5"”?
J . ,3"
3700.0 " [’51 d
1 ix,
4.2!
r‘
’1' 5': --—- node 0
:3], 9 -------- a node 1
‘ . node 2
'.' --— node3
200.0 . ~
3 f 4r - - a node 4
.; 9—9 node 5
,' —-- node 6
2700.0 1 - 1
0 50 100

fitness

118

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4500 -
4000-
3500 .
0 50 100 150
60 . 100 . .
yam—M a 50W
40 1 0 - 4
0 50 100 150 0 50 100 150
1 1 -
L .m
4 0 = o
E 8
-1 . . -1 - .
o 50 100 150 0 50 100 150
1 . . 1 . fi
3 0 g o
g: .9
U) (D
_1 a a _1 . 1
0 50 100 150 0 so 100 150

Figure (4.12) - Fitness function terms vs. generation number for single node
topology GA.

 

119

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3000 r
0 50 100 150
50 ~ 100 ~
E45 .6
40 a a 50 a a
0 50 100 150 0 50 100 150
1 . 2 1 . -
3— in
.9 0’ — = o
2 [SJ/“V— g
—1 a 4 -1 - .
0 50 100 150 0 50 100 150
1 . . 1 . .
g o (g o
2’ .9
U) U)
-1 . - -1 - .
0 50 100 150 0 50 100 150

generation no. generation 110.

Figure (4.13) - Fitness function terms vs. generation number for ring structure
GA.

 

 

120

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4500
4000 {fr
3500 ‘
0 50 100 150
50 . 100 . .
. 5 l ‘ (—\ E H
E 45 9‘5 50
40 ‘ ‘ 0 ‘ ‘
0 50 100 150 0 50 100 150
1 ~ r ' 1 7 '
:— it
.9 0 . r_ z o
2 f 8
-1 4 a -1 a a
0 50 100 150 0 50 100 150
1 . f 1 . .
=- a
3 0 s o
9 .9
m (D
—1 t ‘ -1 ‘ ‘
0 50 100 150 0 50 100 150
generation r10. generation 110.

Figure (4.14) - Fitness function terms vs. generation number for Injection
structure GA.

 

 

fitness

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

generation no.

Figure (4.15) - Fitness function terms vs. generation number for hybrid structure

5000 . -
3000 ‘
0 50 100 150
50 .
E 45 J 1
1
40 a
0 50 100 150
o .
-0.5'
-1 . .
0 50 100 150
1 . .
3 0
.9
in
-1 a a ‘
0 50 100 150

121

GA.

 

 

L

 

 

150

 

 

 

 

 

150

 

 

 

 

 

50 100
generation no.

150

 

 

Stiffness Value (Pa)

3.0e+08

122

 

2.0e+08

1 .0e+08

I

 

 

 

 

 

0.064.001 ‘
‘~-‘ I-d-1\/"A---VA'”-A---t\/-\_A_"‘b'---A-’ \k/l-A’-
~ __1 r ‘

113-1 /~A-’\--d’

\ ..

“1 _ —.\' ........... — ----------------
II.“ I \rl

-1.0e+08 L4 1 1
0 50 100 150

Generation Number

Figure (4.16) - Extensional stiffness matrix versus generation number.

 

 

Stiffness Value (Pa)

123

 

 

 

 

 

5.0e+06 r T
.’ 1
1.7 " 1-1 ‘ 9--- _-/°-_'l -‘
J '3 Izv-‘h 1| 1‘ l"! V b-I 1‘6’__*--+-- ‘
l | [(9 ‘ ‘ \9- '9\-
0.00400 F. 1 1. ‘ , F . "Gm ...-
I 1-415 “"1-..‘
19' ' c«norm
{#12, ll J e—ob12
1 ' r. _-_
l 1‘: 1 11 “-32
I l 1 ll.
.5.Oe+06 ‘- kJ 1 I" ..1 A ......... Am
..1L__
\~____r\.__ ___
\
LVA—J-‘
4.06407 L 1 . L
0 50 100 150
GenerationNumber

Figure (4.17) - Coupling stiffness matrix versus generation number.

 

 

Stiffness Value (Pa)

124

2.0e+06 4 t . 1 .

 

....................

 

1.09106

 

 

 

-1.0e+06 . J . 1 .
0.0 50.0 100.0 150.0

Generation Number

 

Figure (4.18) - Flexural stiffness matrix versus generation number.

 

 

 

125

Table (4.6) Layup Determined by Single Node Structure.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Ply No. Mtl. Orientation Thickness
l l -l.406250e+01 1.27e-04
2 l -1.968750e+01 1.27e-04
3 l -2.250000e+01 1.27e-04
4 1 -2.531250e+01 1.27e-04
5 1 -2.812500e+01 1.27e-04
6 l -2.53l250e+01 1.27e-04
7 l -1.687500e+01 1.27e-04
8 1 -l.406250e+01 1.27e-04
9 l 7.593750e+01 l .27e-04

10(core) 2 0.000000e+00 1.5076938e-01
l l l 8.156250e+01 l .27e-04
12 l -3.656250e+01 l .27e-04
13 l 9.000000e+01 1.27e-04
14 l -8.437500e+00 1.27e-04
15 1 8.718750e+01 1.27e-04
16 1 -2.812500e+01 1.27e-04
17 1 8.437500e+01 l .27e-04
18 l -8.437500e+01 1.27e-04
19 1 8.718750e+01 1.27e-04

20 l 7.593750e+01 1.27e-04
21 l -1.125000e+01 1.27e-04
22 l 7.875000e+01 l .27e-04
23 l 0.000000e+00 1.27e-04
24 l -7.593750e+01 l .27e-04
25 1 -8.718750e+01 1.27e-04
26 1 -2.250000e+01 1.27e-04

 

 

 

 

 

 

Table (4.7) Layup Determined by Ring Structure GA.

126

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Ply No. Mtl. Orientation Thickness
1 1 -l.40625e+01 1.27e-04
2 1 -8.71875e+01 1.27e-04
3 l 8.15625e+01 1.27e-04
4 l -2.53l25e+01 1.27e-04
5 1 8.71875e+01 1.27e-04
6 l -2.25000e+01 1 .27e-04
7 l -8.71875e+01 1.27e-04
8 1 -l.40625e+01 1.27e-04
9 1 -8.71875e+01 1.27e-04
10 1 -8.43750e+00 1.27e-04
l 1 1 -8.15625e+01 l .27e-04
12 1 8.71875e+01 1.27e-04
13 1 -7.87500e+01 1.27e-04
l4 1 -l.96875e+01 1.27e-04
15 l -8.43750e+01 1.27e-04
16 l -2.25000e+01 1.27e-04

l7(core) 2 0.00000e+00 l .4809438e-01
18 l -5.62500e+00 1.27e-04
19 l 9.00000e-1-01 1.27e-04
20 1 -l.68750e+01 1.27e-04
21 l -2.25000e+01 1.27e-04
22 l -2.25000e+01 1.27e-04
23 l -2.25000e+01 1.27e-04
24 l -8.7l875e+01 1.27e-04
25 1 -2.53125e+01 1.27e-04
26 l -3.09375e+01 1.27e-04
27 1 -3.37500e+01 1.27e-04

 

 

 

127

 

 

 

 

 

 

 

Table (4.7 (cont’d).
Ply No. Mtl Orientation Thickness
=_-=l———
28 1 -2.53125e+01 1.27e-04
29 1 -8.43750e+01 1.27e-O4

 

Table (4.8) Layup Determined by Injection Structure GA.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Ply No. Mtl. Orientation Thickness
1 1 -2.812500e+01 1.27e-04
2 l 7.312500e+01 1.27e-O4
3 l -7.875000e+01 1.27e-04
4 l 8.437500e+01 1.27e-04
5 l -3.375000e+01 1.27e-04
6 1 -2.531250e+01 1.27e-04
7 1 -8.156250e+01 1.27e-04
8 1 -1.406250e+01 1.27e-04
9 l -4.SOOOOOe+01 1.27e-04
10 1 8.437500e+01 1.27e-04
l l l -2.250000e+01 1.27e-04
12 1 -l.968750e+01 1.27e-04
13 1 -3.656250e+01 1.27e-04
l4 1 -l.125000e+01 1.27e-04

15(core) 2 0.000000e+00 1.490975e-01
16 1 -l.968750e+01 1.27e-04
17 1 -l.968750e+01 1.27e-04
18 l -7.031250e+01 1.27e-04
19 l -2.250000e+01 1.271304

20 1 -2.250000e+01 1.27e-04

 

 

 

 

 

128

Table (4.8) Layup Determined by Injection Structure GA.

 

 

 

 

 

 

 

 

 

 

 

 

Ply No. Mt]. Orientation Thickness
2] 1 -].406250e+01 1.27e-04
22 l -3.093750e+01 1.27e-04
23 l -7.875000e+0] 1.27e-04
24 1 -8.437500e+01 1.27ev04
25 l -].]25000e+01 1.27e-04
26 1 78750001340] 1.27e.04
27 1 -8.437500e+01 1.27e-04

 

 

' ' I . “Ir-.31!

 

Table (4.9) Layup Determined by Hybrid Structure GA.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Ply No. Mtl.. Orientation Thickness
1 l -1.968750e+01 1.27e-04
2 l -3.375000e+01 1.27e-04
3 ] 8.718750e-1-01 1.27e-04
4 l 8.156250e+01 1.27e-04
5 1 -2.250000e+01 1.27e-04
6 l -8.7]8750e+01 1.27e-04
7 l ~6.187500e+01 1.27e-04
8 ] -3.656250e+01 1.27e-04
9 1 -2.812500e+01 1.27e-04
10 l -8.437500e+01 1.27e-04
ll 1 -].687500e+01 1.27e-04
12 l 7.875000e+01 1.27e-04
l3 1 -3.375000e+01 1.27e-04
14 1 -2.250000e+01 1.27e-04

 

 

129

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Table (4.9) (cont’d)
Ply No. Mt]. Orientation Thickness
= 4
15 ] -2.531250e+01 1.27e-04
16 1 -8.7l8750e+01 1.27e-04
l7(core) 2 0.000000e-1-00 1.5076938e-01
l8 1 8.437500e-1-01 1.27e-04
19 1 -3.375000e+01 1.27e-04
20 l 4.781250e+01 1.27e-04
21 1 -3.656250e+01 1.27e-04
22 l -2.531250e+01 1.27e-04
23 1 -2.812500e+01 1.27e-04
24 1 -].125000e+01 1.27e-04
25 1 -3.375000e+01 1.27e-04
26 l -8.437500e+01 1.27e-04
27 1 -8.156250e+01 1.27e-04
28 1 -3.656250e+01 1.27e-04
29 l -8.437500e+00 1.27e-04
30 1 -2.531250e+01 1.27e-04
31 1 -8.7]8750e+01 1.27e-04
32 l -3.375000e+01 1.27e-04

 

 

 

 

 

 

Chapter 5

Summary and Conclusions.

This document consists of three major parts dealing with the research performed.
The first section, discussed in Chapter 2, dealt with describing different theories available
for analyzing laminated composite plates. Classical laminated plate theory, first order
shear deformation theory, a third order shear deformation theory, and two zig-zag theories
were described and their details discussed along with the associated finite element models.
Finite element models associated with the first order shear deformation theory (FSDT),
third order shear deformation theory (HSDT), and two third-order zig-zag theories, ZZHT
and LZZHT, were used to analyze symmetric and unsymmetric square sandwich panels
subject to sinusoidally varying transverse surface loadings. Results from this comparison
are contained in Chapter 3.

The FSDT theory proved the least accurate and only slight improvement was
gained from using the third order shear deformation theory. As plate thickness became
small in comparison to the other dimensions, accuracy of these theories improved, how-
ever the zig-zag theories were far superior in all cases analyzed. The third order zig-zag
theories used two different displacement fields and different finite element models. The
ZZHT theory used a four-noded plate element and the LZZHT theory an eight-noded
brick element, both had 7 degrees of freedom per node. A sublaminate approach was pos-

sible with the LZZHT model because the element was formulated in such a way as to

130

 

 

131

allow elements to be stacked in the thickness direction. This was achieved by including
transverse shear stress degrees of freedom in the element. The elements in the ZZHT the-
ory cannot be stacked. The ability of the LZZHT model to be refined through the thickness
of the plate (by stacking elements) makes this by far the most accurate plate theory ana-
lyzed. Results from this theory were shown to be nearly identical to the elasticity solu-
tions. Even transverse shear stresses were in very close agreement with the elasticity
solution when using only three sublaminates through the thickness. Results from the
ZZHT theory were comparable in accuracy to the LZZHT model when only one sublami-
nate was used. Computational efficiency of the ZZHT model was far greater than the
LZZHT model for the same number of elements since the LZZHT elements have eight
nodes per element and the ZZHT elements have 4. The best combination of solution accu-
racy and computational efficiency belongs to the ZZHT model and thus it was chosen to
be used in the design optimization portion of this research which is discussed in Chapter 4.
Computational efficiency is very important in an optimization problem because many cal-
culations must be performed. Any savings in computation time can amount to a large sav-
ings over the course of an optimization.

Chapter 4 discusses the results of a design optimization. The motivation for this
design problem came from the desire for aircraft wings to adjust their angle of attack when
experiencing gust loads. Gust loading causes the wing to twist up and back which is detri-
mental to the wing’s lift. The wing needs to respond to such loadings by twisting in the
opposite direction. This research shows such a response is possible with laminated sand-
wich structures. The ZZHT finite element model was used with the genetic algorithm,

GALOPPS v3.0.

132

The desired amount of twist was achieved while maximizing in-plane stiffness and
minimizing structure weight. To do this the number of plies in the top and bottom face
sheets, the fiber orientation of each ply, and the core thickness were treated as design vari-
ables. Maximum tip deflection was constrained by penalizing the design fitness for
exceeding a specified maximum value. The maximum stress in the fiber and matrix direc-
tions were constrained to be less than the strength of the face sheet material. This was
achieved by applying a strong penalty to the fitness of designs exceeding the maximum
allowable stresses.

Also in Chapter 4 four different GA architectures were used to solve the optimiza-
tion problem and the performance of each compared. The ring structure generated the
design with the highest fitness (4,140.25) and the injection structure generated the design
with the lowest fitness (4,096.00). The average fitness of four runs for each of the four
structures was compared with the ring and hybrid structures having the highest average fit-
ness (4,130.8) and the injection structure the lowest at 4,116.5. This is a 0.3% difference
which indicates all GA structures performed sufficiently well for this design problem. The
injection structure seemed to converge the fastest followed by the hybrid, single, and ring
structures.

The single node structure produced the design with the lowest weight of 66.19 kg
as well as the design having the “best” twist, -0.0279 m, which is only 0.44% less than the

desired twist. The design having the maximum in-plane stiffness, A11=2.427e8 Pa, was
produced by the hybrid structure.

The designs generated by the four structures are valid designs, however they are

slightly conservative because the ZZHT finite element model slightly overpredicts

133

maximum in-plane stresses. The LZZHT model is far too computationally expensive for
large Optimization problems so the ZZHT model was the best choice of the theories
available.

As mentioned in Chapter 4, the wing design data, i.e. wing shape, sweep back, air-
craft size and weight, actual wing loadings, desired wing twist, etc., was idealized to
devise a general case which illustrates the ability to solve such a problem using an experi-
mental finite element plate model with genetic algorithms. This problem can be extended
and refined as desired to account for more realistic designs once specific design informa-
tion is known. Also interactive failure criteria could be incorporated into the problem in
lieu of the maximum stress failure criteria used here.

One concern which was not directly addressed by this research is design manufac-
turability. Due to the unsymmetric lamination sequence necessary to obtain the desired
response, the wing will warp during curing. This could be corrected by manufacturing

techniques. Alternatively, the Bi j, A16, and A26 terms could be forced to be zero using

constraint or penalty methods. This would reduce coupling effects which cause warping

during curing. Cost constraints could be incorporated into the design process as well.

 

LIST OF REFERENCES

 

134

List of References

Adali, S., Richter, S., Verijenko, V. E., and Summers, E. B., 1995, “Optimal
Design of Laminates with Discreet Ply Angles for Maximum Buckling Load and
Minimum Cost”, Composite Structures, Vol. 32, pp. 409-415.

Arbate, S., 1994, “Optimal Design of Laminated Plates and Shells”, Composite
Structures, Vol. 29, pp. 269-286.

Avalle, M., and Belingardi, G., 1995, “A Theoretical Approach to the Optimization
of Flexural Stiffness of Symmetric Laminates”, Composite Structures, Vol. 3], pp.
75-86.

Averill, RC, and Reddy, .1. N., 1992, “An Assessment of Four-Noded Plate Finite
Elements Based on a Generalized Third-Order Theory”, Int. Journal for Numerical
Methods in Engineering, Vol. 33, pp. 1553-1572.

Averill, R. C., 1994, “Static and Dynamic Response of Moderately Thick Lami-
nated Beams With Damage”, Composites Engineering, Vol. 4, no. 4. pp. 381-395.

Averill, R. C., and Yip, Y. C., 1996, “An Efficient Thick Beam Theory and Finite
Element Model with Zig-Zag Sublaminate Approximations”, to be published in
AIAA journal.

Burton, W. 8., and Noor, A. K., 1994, “Three-Dimensional Solutions for Thermo-
mechanical Stresses in Sandwich Panels and Shells”, Journal of Engineering
Mechanics, Vol. 120, pp. 2044-2071.

Chattopadhyay, A., McCarthy, T. R., and Pagaldipti, N., 1995, “Multilevel Decom-
position Procedure for Efficient Design Optimization of Helicopter Rotor Blades”,
AIAA Journal, Vol. 33, No. 2, pp. 223-230.

Chawla, K.K., 1987, “Composite Materials, Science and Engineering”, Springer-
Verlag, New York.

 

10.

11.

12.

13.

14.

15.

16.

17.

18.

19.

135

Cho, W-M., Youn, S-K, and Lee, Y-S., 1995, “A Numerical And Experimental
Study on The Nonlinear Behaviour of Laminated Composite Structural Compo-
nents”, Computers and Structures, Vol. 57, NO. 6, pp. 1071-1077.

Fukunaga, H., and Vanderplaats, G.N., 1991, “Strength Optimization of Laminated
Composites With Respect to Layer Thickness and/or Layer Orientation Angle”,
Computers and Structures, Vol. 40, No. 6, pp. 1429-1439.

Fukunaga, H., and Vanderplaats, G.N., 1991, “Stiffness Optimization of Orthotro-
pic Laminated Composites Using Lamination Parameters”, AIAA Journal, Vol. 29,
NO. 4, pp. 641-646.

Furura, 11., Haftka, R.T., 1993, “Locating Actuators for Vibration Suppression on
Space Trusses by Genetic Algorithms”, Structures and Controls Optimization,
ASME, AD Col. 38, PP. 1-11.

Graesser, D.L., Zabinsky, Z. B., Tuttle, M.E., and Kim, G.I., 1993, “Optimal
Design of a Composite Structure”, Composite Structures, Vol. 24, pp273-28l.

Gibson, L. J ., and Ashby, M. F., 1988 Cellular Solids, Structure and Properties,
Pergamon Press, New York.

Goldberg, DE, 1989, “Genetic Algorithms in Search, Optimization, and Machine
Learning”, Addison-Wesley Publishing Company Inc., New York.

Haftka, R. T, and Gurdal, Z., 1992, “Elements of Structural Optimization”, third
edition, Kluwer Academic Publishers, Boston.

Houghton, E. l... and Carpenter, P. W., 1993, “Aerodynamics for Engineering stu-
dents”, 4th ed., John Wiley & Sons, inc., New York.

Johnson, E.A., 1994, “The Implementation of Continuously Updated Sharing in
The Simple Genetic Algorithm Code, and it’s Application to the Optimal Place-
ment of Elastic Supports on a Simply Supported Plate”, Proceedings of the AIAA/
ASME/ASCE/AHS Structures, Structural Dynamics, and Materials Conference,

 

20.

21.

22.

23.

24.

25.

26.

27.

28.

136

AIAA, New York, Vol. 4, pp. 2276-2286.

Jones, R.M., 1975, “Mechanics of Composite Materials”, Hemisphere Publishing
Corporation, New York.

Kam, T.Y., and Lai, M.D., 1989, “Multilevel Optimal Design of Laminated Com-
posite Plate Structures”, Computers and Structures, Vol. 31, NO. 2, pp. 197-202.

Kam, T.Y., and Snyman, J .A., 1991, “Optimal Design of Laminated Composite
Plates Using a Global Optimization Technique”, Composite Structures, Vol. 19,
pp. 351-370.

Kam, T. Y., and Lai, F. M., 1995, “Maximum Stiffness Design of Laminated Com-
posite Plates via a Constrained global Optimization Approach”, Composite Struc-
tures, Vol. 32, pp. 391-398.

Kogiso, N., Watson, L. T., Gurdal, Z., and Haftka, R. T., 1993, “Genetic Algo-
rithms with Local Improvement for Composite Laminate Design”, Structures &
Controls Optimization, ASME, Aerospace Division, New York, Vol. 38, pp. 13-28.

Kogiso, N., Watson, L. T., Gurdal, Z., Haftka, R. T., and Nagendra, S., 1994, “Min-
imum Thickness Design of Composite Laminates Subject to Buckling and
Strength Constraints By Genetic ALgorithms”, Proceedings of the AIAA/ASME/
AS CE/AHS/AS C Structures, Structural Dynamics, and Materials Conference,
AIAA, New York, Vol.4, pp. 2257-2275.

Lake, R. C., Izadpanah, A. P., and Baucom, R. M., Jan 1993, “Experimental and
Analytical Investigation of Dynamic Characteristics of Extension-Nist-Coupled
Composite Tubular Spars”, NASA TP-3225 ARL-TR-30.

Le Riche, T., and Haftka, R. T., 1993, “Optimization Of Laminate Stacking
Sequence for Buckling Load Maximization by Genetic Algorithm”, AIAA Journal,
Vol. 3]. No. 5. pp. 951-956.

Leung, M., Nevill, G. E., 1994, “Genetic Algorithms for Preliminary 2-D Struc-
tural Design”, Proceedings of the AlAA/ASME/ASCE/AHS/AS C Structures, Struc-
tural Dynamics, and Materials Conference, AIAA, New York, Vol. 4, pp. 2287-

29.

30.

3].

32.

33.

34.

35.

36.

37.

38.

137

2291.

Livine, E., 1994, “Equivalent Plate Structural Modeling for “ring Shape Optimiza-
tion Including Transverse Shear”, AIAA Journal, Vol. 32, NO. 6, pp. 1278-1288.

Mase, G. E., and Mase, G. T., 1992, “Continuum Mechanics for Engineers”,CRC
Press, Ann Arbor.

Miki, M., 1985, “Design of Laminated Fibrous Composite Plates with Required
Flexural Stiffness”, Recent Advances in Composites in the United States and
Japan, ASTM stp 864, pp. 387-400.

Miki M., and Sugiyama, Y., 1993, “Optimum Design of Laminated Composite
Plates Using Lamination Parameters”, AIAA Journal, Vol. 30, No. 5, pp. 921-922.

Nagendra, S., Haftka, R. T., and Gurdal, Z., 1992, “Stacking Sequence Optimiza-
tion Of Simply Supported Laminates With Stability and Strain Constraints”, AIAA
Journal, Vol. 30, No. 8, pp. 2132-2137.

Nagendra, S., Haftka, R.T., Gurdal, Z., 1993, “Design of 3 Blade Stiffened Com-
posite Panel by a Genetic Algorithm”, Proceedings of The 34th Structures, Struc-
tural Dynamics and Materials Conference, La Jolla, CA. pp. ??.

Nagendra, S., Jestin, D., Gurdal, Z., Haftka, R.T., and Watson, LT, 1994,
“Improved Genetic Algorithm for the Design of Stiffened Composite Panels”, Pro-
ceedings of the International Symposium on Mathematical Programming.

Noor, A. K., 1992, “Mechanics of Anisotropic Plates and Shells-A New Look at an
Old Subject” Computers and Structures, Vol. 44, No. 3, pp. 499-514.

Olhoff, N., Taylor, J .E., 1983, “On Structural Optimization”, Journal of Applied
Mechanics, Vol. 50, pp.l 139-1 15 l .

Punch, W. F., Goodman, E. D., Pei, M., Chia-Shun, L., Hovland, P., and Enbody,
R., 1993, “Further Research on Feature Selection and Classification Using Genetic
Algorithms”, Procedings of the Fifth International Conference on Genetic

 

39.

40.

4].

42.

43.

45.

46.

47.

48.

138

Algorithims, Urbana Champaign, pp. 557-564.

Punch, W.F.III, Goodman, E.D., Lin, S.-C.,,Ding, Y., Averill, R.C., and Yip, Y.C.,
1994, “Optimal Design of Laminated Composite Structures Using Course-Grain

Parallel Genetic Algorithms”, Computing Systems in Engineering, Vol. 5. pp. 415-
423.

Rauscher, M., 1953, “Introduction to Aeronautical Dynamics”, John Wiley &
Sons, inc., New York.

Raymer, D. P., 1992, “Aircraft Design: A Conceptual Approach”, AIAA Educa-
tion Series, Published by AIAA Inc., Washington DC.

Reddy, J. N., 1990, “A Review of Refined Theories of Laminated Composite
Plates”, The Shock and Vibration Digest, Vol. 22, no. 7, pp. 3-17.

Reddy, J. N ., 1993, “An Introduction to the Finite Element Method, Second Ed.”,
McGraw-Hill, New York.

Reddy, J. N., and Robbins, D., H., 1994, “Theories and Computational Models for
Composite Laminates”, Applied Mechanics Reviews, Vol. 47, no 6, part 1, pp. 147-
169.

Smith, H. C., 1992, “The illustrated guide to Aerodynamics” 2nd ed., TAB Books,
inc., Blue Ridge Summit.

Sandgren, E., Jensen, E., Welton, J .W., 1990, “Topological Design of Structural
Components Using Genetic Optimization Methods”, Sensitivity analysis and Opti-
mization with Numerical methods, ASME Advanced Materials Division, Vol. 115,
pp. 31—43.

Stadler, W., and Dauer, J., 1992, “Multicriteria Optimization in Engineering: A
Tutorial and Survey”, Structural Optimization: Status and Promise, Progress in
Astronautics and Aeronautics, Vol. 150, AIAA Inc., pp. 209-249.

Vinson, J. R., and Sierakowski, R. L., 1986, “The Behavior of Structures Com-

 

49.

50.

51.

52.

139

posed of Composite Materials”, Martinus Nijhoff Publishers, Boston.

VonMises, R., 1945, “Theory of Flight”, McGraw-Hill Book Co. inc., New York.

Wang, BR, and Costin, D.P., 199], ”Optimum Design of a Composite Structure
With Three Types of Manufacturing Constraints”, AIAA Journal, Vol. 30, No. 6,
pp. 1667-1669.

Whitney, J .M., 1973, “Shear Correction Factors for Orthotropic Laminates Under
Static Load”, Journal of Applied Mechanics, Vol.40, pp. 302-304.

Yip, Y. C., 1996 “A 3-D Laminated Plate Finite Element Mth Zig-Zag Sublami-
nate Approximations For Composite and Sandwich Panels”, doctoral disseration,
Department of Material Science and Mechanichs, Michigan State University.