Design optimization of composite laminated structures using genetic algorithms and finite element analysis

Design optimization of composite laminated structures using genetic algorithms and finite element analysis

Composite Structures 88 (2009) 443–454 Contents lists available at ScienceDirect Composite Structures journal homepage: www.elsevier.com/locate/comp...

1002KB Sizes 0 Downloads 89 Views

Composite Structures 88 (2009) 443–454

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Design optimization of composite laminated structures using genetic algorithms and finite element analysis F.S. Almeida, A.M. Awruch * Graduate Program in Civil Engineering, Federal University of Rio Grande do Sul, Av. Osvaldo Aranha, 99, 90035-190 Porto Alegre, RS, Brazil

a r t i c l e

i n f o

Article history: Available online 17 May 2008 Keywords: Multiobjective optimization Genetic algorithms Composites laminated structures Finite element analysis

a b s t r a c t A technique for the design optimization of composite laminated structures is presented in this work. The optimization process is performed using a genetic algorithm (GA), associated with the finite element method (FEM) for the structural analysis. The GA is adapted with special operators and variables codification for the specific case of composite laminated structures optimization. Some numerical examples are presented to show the flexibility of this tool to solve different kinds of problems. Two cases of multiobjective optimization of plates under transverse or in-plane load are studied. In these examples the minimization of two objectives, such as weight and deflection or weight and cost, are simultaneously performed and a pareto-optimal set is obtained by shifting the optimization emphasis using a weighting factor. The stiffness maximization of a composite shell under pressure load is presented in the last example, where the geometrically nonlinear behavior of the structure is considered. Some aspects of the optimization performance, such as the apparent reliability and the computational cost, are investigated in each application. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction In the recent decades, structural applications with composite laminated materials have been growing, requiring great effort on the development of analysis and design techniques. An advantage of using fiber-reinforced composites over conventional materials is that the former can be tailored to specific requirements of certain applications [1]. However, the large number of design variables and the complex mechanical behavior associated with such materials turn the structural design much more difficult and laborious than those involving conventional materials. These characteristics have motivated the use of optimization methods in the sense of turn the composite material structural design a more systematic and well defined task, becoming less dependent to the designer sensitivity and achieving the maximum material performance [2]. The earlier works in the field of composite structures optimization employed the same methods already used to optimize conventional material structures. These methods are based on gradients of the objective and constraints functions with respect to the design variables, which are considered to be continuous in the design space. Such works resulted in limited success because composite laminate design falls on a discrete optimization problem, since in practice the variables are restricted to few values imposed by the manufacturing process. Moreover, the composite optimization * Corresponding author. Tel.: +55 51 3308 3587; fax: +55 51 3308 3999. E-mail address: [email protected] (A.M. Awruch). 0263-8223/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2008.05.004

problems typically involve multimodal search spaces, which may lead gradient based methods to converge to locally optimal regions in the design space [3]. Many other optimization techniques have been tested as an alternative to the gradient based methods, having the genetic algorithm (GA) stand out the others because it perfectly adjusts to the characteristics of the composite optimization problem. GAs are probabilistic search methods seeking to mimic the biological reproduction and natural selection process through random but structured operations. The design variables, usually restricted to discrete values, are coded as genes using binary or integer numbers, and grouped together in chromosomes strings that represent an organism (a possible solution in the design space). Instead of working with just one search point in the design space, GA uses a population of designs that, by reproduction and selection operations, evolve through successive generations. Many search points dispersed in the design space prevent the GA to get stuck in locally optimal regions, avoiding a premature convergence of the process. New designs are generated by the reproduction process that consists in the application of the genetic operators in parents selected from the existing population. These genetic operators are counterparts of the natural genetic mechanisms, acting over the chromosomal strings of the organisms [4]. The selection of parents for the reproduction process and the selection of organisms to fill each new generation are both probabilistic. However, the chances of selection of each organisms is proportional to its fitness, as happen in the nature where fittest organisms have more chances to reproduce and to continue in the next generation.

444

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

The organism fitness is obtained directly from an objective function using simple structure information and gradient evaluations are not required. In real designs cases, when the structural geometry is usually complex and the prediction of the structural behavior must be accurate, it is necessary to use numerical tools, such as the FEM, for the structural analysis. These methods are computationally expensive and may turn the optimization processes with GA impracticable when a large amount of analyses is required. Many researchers have proposed modifications to the classical GA structure to take advantage of composite laminates characteristics and minimize the computational cost. Some of these new strategies are applied in this work, consisting essentially in a GA restructuring of the variable codification and the genetic operators. When more than one objective is handled in the design process, a multiobjective optimization problem may be solved considering all the objectives simultaneously, providing a set of optimum designs (pareto-optimal set), depending of the emphasis given to each one of the objectives. The pareto-optimal set may be very useful when the critical objective is not known a priori. Various researchers have studied the problem of multiobjective optimization of laminated structures, but the use of this approach together with GA and the finite element method (FEM) have not been widely explored. In the present work, two examples of multiobjective optimization of composite laminate plates using GA and FEM are studied. A third application deals with the single optimization of a semicylindrical shell considering the geometrically nonlinear behavior of the structure.

2. Composite structure analysis Real problems of composite structure design depend of reliable structural analysis. In the case of composite laminates, the determination of the mechanical behavior is difficult even for simple geometric configurations. It happens because of some complex mechanisms inherent to the material, like coupling between stretching, bending and twisting deformations, depending on the stacking sequence. These coupling effects and the usually complex structural geometric configurations inhibit the use of closed-form solutions in the structural analysis. In consequence, numerical method must be adopted for the accurate prediction of the structural behavior. In this work, a triangular flat plate and shell element with 18 degrees of freedom called DKT (discrete Kirchhoff triangle) is used. This element was developed by Bathe and Ho [5] for the nonlinear analysis of isotropic plates and shells. The original formulation was modified by the introduction of specific constitutive matrices of the membrane and bending parts in order to allow the analysis of symmetric composite laminated structures. Any other element formulation could be adopted, since in GA based optimization no gradient evaluation is needed, and so no additional modification to the analysis tool is necessary. In the solution of geometrically nonlinear problems an incremental iterative scheme referred as generalized displacement control method (GDCM) is used (see [6]). This method is able to solve nonlinear problems with multiple critical points and snap-back points allowing the complete determination of the post buckling behavior of the structure. The only parameter to be set in the GDCM is the basic load increment ki, which corresponds to the ratio between the first load increment and the complete load. This parameter determines the sensibility of the method to the nonlinear characteristic of the problem. Additionally to the structure displacements, the analysis tool must be able to determinate accurately the stress components at the composite layers in order to predict material failure. This is a

very common constraint adopted in most of optimization problems, and it is used in this work too. The Tsai–Wu failure criterion [7] is used for the failure prediction in a ply. A safety factor against failure kf can be obtained in linear problems with the Tsai–Wu failure function using the material strength parameter for traction, compression and shearing at each of the principal material axes. In geometrically nonlinear analysis the structural failure is verified at each load step of the incremental solution method, which is stopped if material failure is detected.

3. Genetic algorithms for composite laminate structure optimization This work uses a GA provided with many modifications with respect to the classical GA structure stated by Goldberg [4]. Although the main concepts and the sequence of operations remain similar to the original formulation, a new scheme for the variable codification and special genetic operators are introduced, increasing the performance of the method in the case of composite structures optimization. An initial population containing P organism is first created in a random process. In order to create successive generations, parents are chosen from the current population based on their fitness. Next, the genetic operators are applied to create children in order to form a children population. An elitist selection scheme is used to obtain the new generation taking organisms from the current population and from the children population just created. This process is repeated until the convergence criterion is met. A description of the variable codification, the genetic operators and the selection scheme for the construction of new generations are given in the following subsections. 3.1. Composite laminate codification A pair of chromosomes is adopted for the representation of each laminate in this work, similar to the scheme implemented by Soremekun [8]. As only symmetric laminates are studied in this work, just half of the layers need to be coded in the gene strings, and so the total number of genes in a chromosome is proportional to half of the maximum admissible number of layers in a specific design problem. Each layer is represented by a pair of genes (with one gene in each chromosome), being the first pair referred to the outermost layer and the succeeding pairs of genes referred to the inner layers. In the first chromosome, named ‘‘orientation chromosome”, is stored information about fiber orientation and about the number of plies contained by each layer of the laminate. The second chromosome, called ‘‘material chromosome”, is used to store the material properties of the layer, containing information like ply thickness, elastic and strength constants. Each design variable must assume one of the admissible discrete values defined in the optimization process. These values are represented by positive integer numbers. They are used to code de variable value in its correspondent gene. Since there are two kinds of variables, the orientation and the material variable, two code alphabets are used in the codification. An empty stack code, represented by the number zero, is used in problems where the variation of the number of layers is considered. This code is introduced in a pair of genes, representing the layer to be deleted, by a specific genetic operator. Another genetic operator acts adding layers to the laminate by changing the value zero of a pair of genes by any other admissible value. The maximum number of layers in the laminate is limited by the number of genes in the chromosome. In Fig. 1 is presented an example of the decodification of a pair of chromosomes in a laminate stacking sequence, based in a given orientation and material gene codification alphabets. These chro-

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

Fig. 1. Decodification of chromosomes in a laminate stacking sequence.

mosomes could represent a laminate with up to 32 plies, since each one of the 8 genes in the strings represent 2 plies, and only one half of the symmetric laminate is coded. However, the presence of the code ‘‘0” in the first pair of genes indicates that this layer does not exist, indicating that the laminate has only 28 plies. 3.2. Genetic operators 3.2.1. Crossover Crossover is an essential GA operator, having the fundamental task of creating new organisms (children) in a reproduction process. It acts combining genetic information taken from a pair of organisms (parents) selected from the current population. The created child will hopefully be better than, or at least equivalent, in fitness to its parents. The crossover operation is applied by first generating a random number to define the crossover point. Then, the gene strings of both material and orientation chromosomes are split at the same point in both parents [8]. The left part of parent 1 and the right part of parent 2 are combined to form a child, as is shown in Fig. 2. The crossover operator is usually applied with some probability, but in this work crossover is always used to create children. 3.2.2. Mutation Mutation is the class of genetic operators responsible to maintain the genetic diversity of the population by introducing new information in the chromosomal strings of each child after it is created by the crossover operation. These operations provide a random search capability to GA, which may be useful to find promising areas in the design space, and prevent crossover to lose its effect due to a standardization of the population. In the classical implementation of this operator, to each gene is given a small probability to switch to any other permissible value, excepting its current value [3]. In spite of the randomness of the mutation process, it is possible to incorporate to this operator some knowledge about the response of composite laminates with respect to the alteration of one or

445

more of its characteristics. This may lead to a less random operation, which is supposed to be more efficient in guiding the evolution towards optimization objectives. New operators named orientation alteration, material alteration, ply addition and ply deletion, are introduced to replace the classical mutation in the present GA [8,9]. A description of these new operators is presented below, followed by an example of their applications as it is shown in Fig. 3. Orientation and material alterations are implemented similarly to the classical mutation, but are independently applied to orientation and material chromosomes, respectively. Different orientation and material operator probabilities (poa and pma) may be adopted, which is useful, since each chromosome may converge at different velocities in most optimization processes. Additionally, these operators are not applied to genes with empty stack code, when it is present in the chromosome. Two distinct operators are used to vary the number of layers in the laminate. The first one, named ‘‘ply addition”, acts in the chromosomal strings, at a given probability ppa, introducing a new layer close to the laminate mid-plane (end of chromosomes), and removing an existing pair of genes with empty stack code. The genes that represent the new layer are randomly created, assuming any value contained in each respective alphabet codification. In an opposite way, the ‘‘ply deletion” operator acts taking out the last pair of genes of the organism chromosomes (innermost layer) and adding a pair of genes with empty stack code at the outermost position. This operator is applied with a given probability ppd. Both operators manipulate at the innermost laminate layer because it has little effect on bend or twist behavior of the plate or shell, causing no abrupt changes in the design.

Fig. 3. (a) Orientation and material alteration; (b) ply addition; (c) ply deletion.

Fig. 2. Crossover operation.

446

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

PN An ¼

Fig. 4. Gene swap operation.

3.2.3. Gene swap The main characteristic of the gene swap operator is the ability to modify laminate stack sequence without changes of the total number of plies with fibers oriented on each permissible direction. This allows GA to change the bending behavior of the laminate without modifying its in-plane mechanical response. The implementation of this operator follows Refs. [8,9], where two pairs of genes are randomly chosen and have their position shifted in the chromosome, resulting on a new stacking sequence. Such operation occurs at a given probability pgs, usually with a larger value than those corresponding to mutation operators probabilities. A description of the gene swap operation is given in Fig. 4. 3.3. Selection scheme There are many ways to obtain the population of successive generations in a GA. In classical algorithms new generations are formed only by children created by genetic operators applied to the current population. This process has many drawbacks since there is no warranty of improvement or maintenance of achieved evolution when all the old organisms are replaced. To solve this problem new selection schemes were created, being one of them the elitism scheme, which consists in transfer some good organisms from the old population to a new generation, preserving desirable genetic information. This papers deals with a multiple elitist scheme proposed by Soremekun [8]. In the implementation, both, parent and child populations of size P are independently ranked from best to worst fitness. These two populations are then combined and ranked together, resulting in a combined population with 2P organisms. Then, the best Ne individuals of the combined population are transferred to the new generation. The best individuals of child population that have not already been used are taken to fill the remainder of the new generation. The number of top elements (Ne) to be transferred to the new generation is a GA parameter to be adjusted at each application. 4. Numerical examples and discussion Three examples of composite laminated structures design using GA optimization and FEM analysis are presented in the following sections. To prove the success of the optimization procedure and to characterize the design space of the problems, all the possible laminate configurations are previously analyzed. Additionally, N optimizations are carried out for each example to obtain the algorithm reliability and computational cost. This is possible to accomplish because results of previous analyses of the design space are stored. This information is used by the GA to evaluate the objective function. The apparent reliability (R) is determined by taking the number of optimizations where the GA finds at least one global optimum (No), divided by the total number of applications of the GA (N). It defines the chances of obtaining the global optimum in a single optimization process. As the structural analysis employing the FEM is usually the most time consuming task in the optimization procedure, the GA cost (An), which is determined by the average number of analyses required in a single optimization process, is given by the following expression:

i i¼1 X g P

N

;

ð1Þ

where X ig is the total number of generations analyzed in the ith optimization procedure. In GA optimization procedures it is very common to occur repeated analyses because one specific design may appear in many generations during the process. When a memory containing information about the performed analyses is used, associated to the GA, these repeated analyses can be avoided, resulting in an important reduction of the computational cost. In these cases, another measurement of the computational cost can be stated tacking the average number of the effectively performed analyses (Ar), obtained dividing the number of analyses that were effectively carried out by the number of GA executions (N). The criterion to stop the optimization process, which was used in all examples presented here, is based in two parameters: the upper limit of the number of generations (NLG) and the maximum number of generations with no improvement of the best design (NSD). Once one of these limits is reached, the optimization process is stopped and the best laminate of the last generation is taken as the optimization result. NLG and NSD are defined in each optimization procedure, depending on the complexity of a specific problem. 4.1. Weight and deflection minimization of a composite laminated plate under transverse load This example deals with the design of a composite laminated square plate, subjected to a uniform pressure load on its surface. Minimization of the structural weight and deflection are the design objectives. The two objectives must be considered in conjunction with the constraints imposed by material failure and maximum values of contiguous plies thickness with the same fiber orientation. These are two opposite objectives, since improvements in one of them leads to depreciation of the other, and they must be considered at the same time in the optimization, requiring a multiobjective approach. In its formulation, the objective function must contain both objectives, which are weighted by a factor which controls the emphasis given to each one of the objectives in the optimization. As a result of the variation of the weighting factor used in the objective function, this problem has a set of optimal solutions (pareto-optimal set) instead of a single solution. The structure geometry, boundary conditions and the mechanical properties of the composite material are presented in Fig. 5. The elastic constants are the Young’s modulus in the fiber direction (E1) and transverse to the fiber direction (E2), the shear modulus (G12) and the Poisson’s ratio (m12), respectively. Strength parameters for traction and compression for longitudinal and transversal directions are given by F1t, F1c, F2t, and F2c, respectively. The remainder parameters are the shear strength (F6) and the specific weight (q). The structure must support a design pressure load of 0.1 MPa with no material failure (the Tsai–Wu failure function must be lower than 1.0 for the whole plate) and thickness of contiguous plies with the same fiber orientation must not be greater than 2 mm or less than 0.75 mm. In the plate design, the laminate is restricted to be symmetric with 8 layers, being represented in the GA by a pair of chromosomes, with 4 genes each one. The fiber orientation angle and the thickness of each layer are the optimization variables. They must assume one of the discrete values given in Table 1, where the codification adopted in the GA is also given. The number of genes combined with the number of possible values of each variable leads the size of the design space (SDS) to be equal to 65536. In this example the fitness evaluation FIT must consider both, weight and deflection reduction, at the same time. It is done by Eq. (2a), where the fitness value is taken as the inverse of the objective function, and the weighting factor a is introduced to allow the

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

447

Fig. 5. Structure geometry, boundary conditions and composite properties.

Table 1 Genes alphabet and possible discrete variable values Orientation gene alphabet

Material gene alphabet

Code

Orientation angle

Code

Ply thickness (mm)

1 2 3 4

1 1 1 1

1 2 3 4

0.75 1.00 1.50 2.00

ply ply ply ply

at at at at

0° 45° +45° 90°

variation of the emphasis given to each objective. The constraints are considered by a penalty formulation of the objective function, where unfeasible designs have a reduction on their fitness proportionally to the magnitude of the constraints violation. The objective function uses the dimensionless variables W* and D*, that represent the total weight and the central deflection of the plate normalized by their maximum and minimum values. This approach is more efficient than a formulation that uses directly the weight and displacement value or even a simple dimensionless value of these variables divided by a reference value (see [10]). The lower and higher weight limits can be easily obtained by taking all the plies thickness equal to 0.75 mm or equal to 2.00 mm, respectively. The maximum and minimum displacements to be used are obtained by the results of optimizations performed with a equal to 0.0 and equal to 1.0, respectively. When a is equal to 0.0 only the displacement is reduced, and the optimization result is a design that have the smallest displacement. When a is taken equal to 1.0 the GA obtains the lightest structure, which may have a large displacement, and it is used as the maximum value in the normalization. The normalization of the variables is given by Eq. (2b), where the coefficient 1.0 is added to avoid nulls values of W* and D*.

(

FIT ¼ ðTvþ1ÞðaW1 þð1aÞD Þ ; FIT ¼

if FF 6 1;

1 ; FFðTvþ1ÞðaW  þð1aÞD Þ

if FF > 1

Tv is referred to the violation of the limit of contiguous plies thickness with the same fiber orientation. Tv is equal to the exceeding value violating the limit fixed to the thickness of contiguous plies with the same fiber orientation. As an example, if the thickness of each one of two contiguous plies with the same fiber orientation is equal to 1.5 mm, the exceeding value violating the limit (i.e. Tv) is 1.0, since this work adopts a limit of 2.0 mm. The constant 1.0 is added to Tv in order to avoid a division by zero in Eq. (2a) when Tv = 0. To allow the analysis of the optimization performance, GA is executed 50 times for each a, which is taken varying from 0.0 to 1.0 with increments equal to 0.05. A population size P = 50 and the elitist scheme parameter Ne = 5 are adopted (see Section 3.3). The genetic operators are used with the probabilities poa = 4%, pma = 2%, and pgs = 80%, while the probability of ply addition (ppa) and ply deletion (ppd) operations are set to zero, since the number of layers must remain fixed during the optimization (see Sections 3.2.2 and 3.2.3). The parameters used as a criterion to stop the process are NLG = 250 and NSD = 125. The distribution of weight and central displacement of all the feasible designs in the problem is shown in Fig. 6. Points A–P in this

ð2aÞ

where

W ¼

W  W min þ 1; W max  W min

D ¼

D  Dmin þ 1: Dmax  Dmin

ð2bÞ

The parameters FF and Tv are introduced in Eq. (3a) to penalize unfeasible designs. The first one represents the maximum value of the failure function evaluated in the structure, while the parameter

Fig. 6. Weight and central displacement of feasible designs.

448

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

figure are the designs that form the pareto-optimal set, which must be obtained by the GA, according to the emphasis given to each of the objectives. Details of the pareto-optimal set are presented in Table 2. Good reliability levels were obtained for the optimization with most of the a values, as can be observed in Table 3, where the column r represents the standard deviation of the apparent reliability R. The loss of reliability in the optimizations using a equal to 0.55 occurs because the GA finds many times the design identified by L instead of the design identified by J, which is the correct solution in this case. It happens because both designs J and L have practically the same fitness value for a equal to 0.55. The same occurs for a equal to 0.70, when the GA finds the design identified by O in many optimization processes instead of the design identified by N, which is the correct solution. In optimizations where a is taken as being equal to 0.80, 0.85 and 0.95, the low level of reliability occurs because the designs obtained in these cases have values of the fitness function very close to the optimal designs identified by O and P. The average number of analyses required (An), the average number of analyses effectively performed (Ar) and the ratio between Ar and the size of the design space (SDS), which is equal to 65,536, are also shown in Table 3 for each value of a. An expressive reduction of the computational cost is observed if the analysis

Table 2 Pareto-optimal designs Optimal design

Laminate

Weight (N)

Deflection (mm)

Weighting factor a

A B C D E F G H I J L M N O P

[902,0, +452,0, 902,0, 451,0]S [902,0, 452,0, 902,0, +450,75]S [902,0, 451,75, 902,0, +451,0]S [902,0, 451,75, 902,0, +450,75]S [902,0, 451,0, 902,0, +451,0]S [902,0, 450,75, 902,0, +451,0]S [902,0, 450,75, 902,0, +450,75]S [902,0, 450,75, 901,75, +451,0]S [902,0, 450,75, 901,75, +450,75]S [902,0, 450,75, 901,0, +451,0]S [902,0, 450,75, 901,0, +450,75]S [902,0, 450,75, 900,75, 450,75]S 1;0 901;0 S ½901;0 2 ; 45 0;75  451;0 S ½901;0 2 ; þ45 0;75 þ 450;75 S ½901;0 2 ; 0

219.7 211.9 204.0 196.2 188.4 180.5 172.7 164.8 157.0 149.1 141.3 133.4 125.6 117.7 109.9

7.9 8.8 9.6 10.8 11.8 13.2 15.0 17.3 19.9 23.3 27.2 32.1 38.5 46.6 59.5

0.0–0.20 – 0.25 – 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75–0.85 0.90–1.0

Table 3 Optimization results with GA for the square plate

a

Optimal designs

R (%)

r (%)

An

Ar

Ar/SDS (%)

0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00

A A A A A C E F G H I J L M N O O O P P P

100 100 100 100 100 98 100 100 100 100 100 84 98 100 80 94 78 68 92 86 92

0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.0 0.0 5.2 2.0 0.0 5.7 3.4 5.9 6.6 3.8 4.9 3.8

7401 7780 7759 7640 7587 8659 8299 7681 7949 7694 7750 8386 8035 8146 7830 8205 7805 7798 8559 8067 8001

3641 3781 3852 3842 3923 4645 4463 4014 4242 4299 4392 4472 4095 4146 3747 3568 3175 3161 3404 3379 3347

5.56 5.77 5.88 5.86 5.99 7.09 6.81 6.12 6.47 6.56 6.70 6.82 6.25 6.33 5.72 5.44 4.84 4.82 5.19 5.16 5.11

of repeated structures is avoided using a memory, since Ar is found to be about 40–57% of An, depending on the value of a. Furthermore, the number of analyses effectively performed Ar is about 4.82–7.09% of the SDS, which means that only a small part of the design space is explored, but the total number of analyses is still huge. A reduction of the average value Ar can be obtained when the parameters of GA are modified. However, in all tests accomplished here this reduction was obtained together with a reduction on the reliability levels, mainly for values of a presenting critical reliability levels. As can be seen in the Tables 2 and 3, the GA is successful in finding most of the pareto-optimal designs, but the designs identified by the points B and D are not obtained. These points are possible solutions of the optimization, but they are located out of a convex curve defined by the other optimal points. Due to this fact the GA does not find the points B and D, since the fitness is evaluated as a convex combination of the objectives. Figs. 7a and b show the difference of the fitness values of the points B and D with respect to their neighbor points in a range of a where the optimal solution changes from point A to C and from point C to E, respectively (see Table 3). The figures show that the fitness of the points B and D are never greater than those of their neighbor points at the same time and so they cannot be obtained by the GA, no matter the value of the weighting parameter. All optimal designs obtained in the pareto optimal set present outer layers with fibers oriented at 90° and the maximum admissible thickness (2 mm). These layers are the most important for the plate bending behavior, and the previous characteristics give the best stiffness properties. The inner layers of the staking sequence of the different designs vary to obtain the different results for the structural weight and deflection, according to the emphasis given for each of the optimization procedures. The lightest design obtained here has less than the half of the weight of the heaviest; however, the last one presented a displacement more than seven times lower than the first one. The chose of the design to be used in a specific application depends on how critical is the weight for such application and on the magnitude of the allowable displacement. 4.2. Cost and weight minimization of an in-plane loaded composite laminate plate This example uses GA to obtain the optimum design of an inplane loaded plate of composite materials. The objective is to find the lightest and cheapest laminate satisfying the design constraints with respect to structural stability and material failure. Two unidirectional composite materials, Kevlar-epoxy and Graphite-epoxy, may be used to construct the stacking sequence, being the first cheaper but heavier and less resistant than the last one. The variables to be manipulated by the GA are the number of layers of the laminate, the material of each one of the layers and the fiber orientation angles of the plies. Plate geometry, load and boundary conditions are presented in Fig. 8. This figure also shows the elastic constants, strength parameters, specific weight, ply thickness and the parameter of cost per unit weight (C) for the Kevlar-epoxy and Graphite-epoxy. The safety factor for material failure kf is given by the Tsai–Wu failure criterion [7], while the safety factor for structural elastic instability kb is obtained solving the eigenvalue problem involving the linear and the geometrically nonlinear stiffness matrices. A design is considered to be feasible when both kf and kb are grater or equal to 1.0. For simplicity, the design cost is considered to be proportional to the amount of each material contained in the different plies of the laminate. This value is obtained by multiplying the portion of the plate weight corresponding to one given material by its cost per unit weight (denoted by C).

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

449

Fig. 7. Differences of the fitness values of points B and D with respect to their neighbor points. (a) When FIT(B) > FIT(A), then FIT(C) > FIT(B); (b) when FIT(D) > FIT(C), then FIT(E) > FIT(D).

Fig. 8. Composite laminate plate under in-plane load.

The laminated is supposed to be symmetric and formed by 6–12 layers. Each layer has 2 plies that may assume the following orientations: 0°2, ±45° and 90°2. They are represented in the orientation gene alphabet by the codes 1, 2 and 3, respectively. The material gene alphabet is formed by the codes 1 and 2, corresponding to Kevlar-epoxy and Graphite-epoxy, respectively. Empty stacks code can be used to reduce the number of layers. Due to the symmetry, only 6 genes are necessary in each of the two chromosomes used to represent the composite material in the GA. Considering the number of variables and the number of possible values of these variables, the size of the designs space (SDS) is 55,944. A distribution of weight and cost of all the feasible designs in the problem is shown in Fig. 9. Points A to F in this figure are the designs that form the pareto-optimal set, which must be obtained by the GA. Details of the points of the pareto-optimal set are presented in Table 4. Due to the simultaneous minimization of cost and weight, the multiobjective approach is used in the formulation of the objective function, given in Eq. (3), which contains both objectives. The weighting factor a is introduced to allow a shifting on the emphasis given to each of the objectives, driving the GA to converge to one of the points in the pareto-optimal set, according to the value attributed to this factor. The constraints of the problem are also considered in the objective function. In this example, a feasible design is determined by the safety factor k*, which is given by the minimum value between kb and kf. A penalization is applied to

the unfeasible designs, driving the GA to feasible areas of the design space. A small ‘‘bonus”, proportional to the safety factor, is incorporated to the objective function in the case of feasible

Fig. 9. Weight and cost of feasible designs.

450

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

Table 4 Designs of the pareto-optimal set

Table 5 Optimization results with GA for the plate with in-plane load

Optimal design

Laminate

Weight (N)

Cost (uc)

kb

kf

A B C D E F

ge ge ge ½45ge ; 90ge 2 ; 02 ; 45 ; 02 S ge ge ke ½90ge 2 ; 02 ; 452 ; 02 S ge ge ke ½90ge 2 ; 02 ; 45 ; 452 S ke ke ; 0 ; 45 ; 90ke ½45ge ; 90ge 2 2 S 2 ke ; 0 ½45ge ; 45ke ; 90ke 4 2 S ke ½45ke 3 ; 904 S

24.49 25.44 26.39 27.34 28.30 29.25

73.46 64.62 55.77 46.93 38.09 29.25

1.04 1.30 1.50 1.64 1.56 1.30

55.12 30.93 31.84 18.00 17.16 16.04

designs. This ‘‘bonus” promotes an additional goal in the optimization process, seeking to maximize the value of the safety factor. As a consequence of this formulation, the result of the optimization process is supposed to be the safest design having the best weighted combination of weight and cost. Parameters W* and C*, used in Eq. (3a), are the dimensionless total weight and cost of the plate, respectively. They are given by Eq. (3b) using the maximum and minimum possible values of weight and cost, which are easily obtained by the extreme combination of materials and number of layers. Instead of working directly with a linear combination of the dimensionless weight and cost, the implemented objective function uses the square of these variables, already weighted by the factor a. This alternative formulation is necessary because the results of optimizations that uses a linear combination of W* and C* as objective function is always one of the extreme points A or F (see Fig. 9). It occurs because the pareto-optimal set is arranged in a straight line in the weight-cost plane. The new formulation provides a curved distribution of optimal points, allowing the GA to find them. The optimal solution is defined as the point which is located at the closest ‘‘distance” from the origin of the weight-cost reference system, being this ‘‘disqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tance” given by ‘‘ ½aðW  Þ2 2 þ ½ð1  aÞðC  Þ2 2 ”.

8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 > > > OBJ ¼ ½aðW  Þ2 2 þ ½ð1  aÞðC  Þ2 2 þ 106 k ; if k P 1; < qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 > > > : OBJ ¼ ðk Þ2 ½aðW  Þ2 2 þ ½ð1  aÞðC  Þ2 2 ;

if k < 1; ð3aÞ

where

W ¼

W  W min þ 1; W max  W min

C ¼

C  C min þ 1: C max  C min

ð3bÞ

To study the GA performance, 25 optimization processes were executed for each a, which is taken varying from 0.0 to 1.0 with increments equal to 0.1. The GA is used with a population size P = 30 and a elitist parameter Ne = 4, together with the following probabilities for the genetic operators: poa = 4%, pma = 2%, ppa = 4%, ppd = 8% and pgs = 80%. The parameters of the criteria used to stop the process are NLG = 300 and NSD = 100. The results of the optimizations performed here show that the GA can obtain one optimal design for every tested a, but not all of the 25 GA executions performed with each value of a are successful. This happens because GA has difficulty to search for the optimal design in regions where the objective function has low gradients, as a consequence of the introduction of the ‘‘bonus” proportional to the safety factor. There are many points (designs) with fitness values very similar to the optimal design fitness value, having the same weight and cost of the optimal design, but with a smaller safety factor. Table 5 shows the apparent reliability (R) and its standard deviation (r), obtained from GA the 25 optimizations for each value of a, considering two situations. In the first column of Table 5 optimizations founding only one optimal design are computed as successful, while in the second column all designs with the same weight and cost of those included in the first column,

a

Pareto-optimal set

Optimal designs R (%)

r (%)

R (%)

r (%)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

F F F E E D C B A A A

100 100 100 100 100 40 20 8 84 92 88

0.0 0.0 0.0 0.0 0.0 9.8 8.0 5.4 7.3 5.4 6.5

100 100 100 100 100 100 100 100 100 100 100

0 0 0 0 0 0 0 0 0 0 0

Near optimal designs

An

Ar

Ar/SDS (%)

3829 3839 3862 3750 3732 4805 4498 5999 4172 4740 4622

1238 1714 1709 1683 1723 2209 2336 2868 2037 2410 2219

2.21 3.06 3.05 3.01 3.08 3.95 4.18 5.13 3.64 4.31 3.97

but having smaller values of safety factor, are considered as valid solutions. These are called ‘‘near optimal designs” and are found in every optimization and for all values of a, showing that GA is efficient in finding near optimal solutions. Table 5 also shows the average number of analyses required (An), the average number of analyses effectively performed (Ar) and the ratio between Ar and the size of the design space (SDS), for each value of a. The algorithm has explored only a little part of the design space, but for optimizations presenting low reliability levels the search was extended to more generations to seek for the best solution, requiring more analyses. Using computer memory to avoid repeated analyses during the optimization gives a significant reduction on the computational cost, as can be seen from the difference between An and Ar in Table 5. Although the variation of the number of layers is allowed in this example, all the designs that form the pareto optimal set are composed by 20 plies. As layers with two plies and the symmetry condition have been used, laminates with total number of plies that are multiple of 4 have been obtained. In this case, the optimization process showed that laminates with 12–16 plies are unfeasible, while laminates with 24 plies are less efficient than those with 20 plies. The transition of the heaviest and less expensive design (F) to the lightest and most expensive design (A) is accomplished modifying gradually the material of the different layers. In each subsequent design contained in the pareto set, the outer layers composed by Kevlar-epoxy were replaced by a graphite–epoxy layer, taking the maximum advantage of the stiffest material. The variation of the weight observed from design A to F is much lower than the variation of the cost. This implies that, unless for weight critical applications, the cheapest design is the best one. However, this conclusion may be strongly affected by the way the cost is evaluated or even if the cost parameter (denoted by C) is modified. 4.3. Stiffness maximization of a composite laminated shell with geometrically nonlinear behavior This example deals with the application of the GA to obtain the design of a composite shallow shell with the maximum stiffness with respect to a pressure load. The nonlinear behavior of the structure is considered in the fitness evaluation by taking into account the critical load level and the maximum displacement of the structure. These values are obtained by the geometrically nonlinear analysis of the problem. The structural optimization is performed by the GA through the manipulation of the laminate stacking sequence, considering a fixed number of layers. Constraints referred to material failure and number of contiguous plies with the same fiber orientation are imposed to the problem. Two cases of optimization are studied in this example, considering the same geometry of the structure, but different magnitudes

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

of the pressure load and thickness of the laminate, what results in different levels of nonlinearity. In the first case the structure is subjected to a pressure load q = 0.25 MPa, and the laminate have a total thickness h = 12.6 mm. The whole design space is previously analyzed, allowing the study of the effect of some GA parameters over the reliability and computational cost of the optimization process by numerous applications of the GA. The second case considers a pressure load q = 0.125 MPa and the laminate thickness h = 6.3 mm. The nonlinear structural behavior observed in this case is stronger than in the first case, arising some difficulties to the analysis requiring a larger number of steps to produce an accurate response for most laminate configurations. A strategy for the automatic refinement of the basic load increment ki, used in the nonlinear analysis by the GDCM [6], is introduced to reduce the computational cost. The effect of this strategy is studied in four optimization processes performed for the second case. In both cases the laminate is symmetric and formed by 14 layers that are represented in the GA by a chromosome containing 7 genes. The material chromosome is not necessary in this example, since the unique layer characteristic that is allowed to be changed in the optimization process is the fibers orientation. The boundary conditions and structural geometry are shown in Fig. 10, together with the elastic and strength parameters of the glass–epoxy composite used in this example. FEM is used to solve the geometrically nonlinear problem. The fitness value of each design obtained by the GA during the optimization process is defined by Eq. (4), where two parameters are used to characterize the structural stiffness and other two parameters are used to consider the constraints violation. The parameters used to evaluate the structural stiffness are the critical load level (NCcrit), determined when the curve pressure  displacement at the central point A reaches the first peak, and the value of the maximum displacement at the same point (Umax), which is taken at the end of the load increment or when material failure is observed. In the case of unfeasible designs, the maximum load level acting on the structure without material failure (NCmax), and the number of violations of the limit of 4 contiguous plies with the same fiber orientation (Vnlc) are used as penalizations.

FIT ¼

! NC crit  NC 2max : U max  ðV nlc þ 1Þ

ð4Þ

The constant value 1.0 is added to Vnlc in order to avoid a division by zero in Eq. (4) when Vnlc = 0. In the first case studied in this example, considering q = 0.25 MPa and h = 12.6 mm, each of the 14 layers that form the laminate contains two plies, oriented at 0°2, ±45° or 90°2. These orientations are represented in the chromosomal string by the codes 1, 2 and 3, respectively. The resulting laminate is formed by 28 plies with an individual thickness of 0.45 mm. Consid-

451

Fig. 11. Critical load level and maximum displacement of feasible designs.

ering the number of genes in the chromosomal string and the number of possible values for each gene, the size of the design space in this example is 2187. The whole design space has been analyzed to allow the evaluation of the GA, and the feasible designs are shown in Fig. 11, where the optimal design is indicated. The optimal design, indicated in Fig. 11, is defined by the stacking sequence [(904, ± 45)2, 902]S and the following value of the parameters were obtained: NCcrit = 0.563, Umax = 0.0272 m and FIT = 20.698. As expected, the optimal design has most of its plies with fibers oriented at 90°, providing the best bending properties. The 45° plies are introduced in order to satisfy the constraint related to the number of contiguous plies with the same fiber orientation. The curve load  displacement at point A for the optimal design is presented in Fig. 12, together with the same curve for a design configuration with a smaller fitness value. In this case the stacking sequence is [±453, 04, ±452]S, NCcrit = 0.228, Umax = 0.0307 m and FIT = 7.427. In the analyses executed for all the points of the designs space the basic load increment ki in the GDCM was taken equal to 0.05. Nine combinations of the GA parameters P and NLG, presented in Table 6, are used to study the influence of these parameters in the optimization performance. Considering that the maximum number of analyses (Amax) to be executed in a optimization process is determined by the product (NLG + 1)  P, it is possible to impose an upper limit to An by adjusting these two parameters. The values of P and NLG are obtained by working with ranges of Amax/SDS equal to 12.5%, 25% and 50%. Using the expression for Amax stated before

Fig. 10. Shallow shell geometry, boundary conditions and material properties.

452

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

Fig. 12. Central displacement  load level for the optimal design and for a low fitness value.

Table 6 Combination of the parameters P and NLG Amax/SDS = 12.5%

Amax/SDS = 25.0%

Amax/SDS = 50.0%

Comb

P (Ne)

NLG (NSD)

Comb

P (Ne)

NLG (NSD)

Comb

P (Ne)

NLG (NSD)

1 2 3

7 (1) 14 (2) 21 (3)

38 (30) 19 (15) 12 (10)

4 5 6

7 (1) 14 (2) 21 (3)

77 (44) 38 (22) 25 (15)

7 8 9

7 (1) 14 (2) 21 (3)

155 (55) 77 (27) 51 (18)

and adopting the population size P equal to 7, 14 and 21 (proportional to 1, 2 and 3 times the number of genes in the chromosome), it is possible to determine the value of NLG for each range of Amax/ SDS. The values of the parameters Ne and NSD are conveniently adjusted in each of the nine combinations, but these parameters are also decisive in the efficiency of the process. The genetic operators probabilities adopted here are poa = 4% and pgs = 80%, being the remainder set of probabilities equal to zero because they are not used in this problem. The GA is executed 50 times for each of the nine combinations of parameters in order to study their influence in the apparent reliability (R) and in the computational cost of the optimization process. Fig. 13 shows the performance of the GA considering each one of the combinations of P and NLG tested in this work. In combina-

Fig. 13. GA performance with different combinations of P and NLG.

tions 1, 2 and 3, referred to Amax/SDS = 12.5%, low computational cost is obtained together with unacceptable levels of reliability. The other combinations have generated good reliability levels (over 95%, excepting combination 4 which has R = 90%) and the increase observed in the computational cost is not proportional to the increase in Amax/SDS to 25% and 50%. Combinations with smaller values of P have shown to be not so expensive in terms of computational cost than the others belonging to the same range of Amax/SDS. A reduction of 30% in the number of analyses to be executed in each optimization is observed as a consequence of the adoption of a strategy to avoid the repetition of analyses that have been performed previously. This aspect is more expressive in optimizations using smaller values of P, as for example in combination 7, where the saving in computational cost is of about 44%. As it is shown in Fig. 13, the optimizations using small populations and many generations have the best performance in terms of computational cost, maintaining good reliability levels. Although the GA performance is highly dependent on the problem to be solved, the tendencies observed in this example may be extended to many other applications, helping to determine suitable values of the parameters P and NLG. Furthermore, if optimizations resulting in designs with fitness values grater or equal to 99% of the optimum design (quasi optimal designs) are considered as being successful, the reliability level is increased to 100% in all the nine combinations, as it is shown in Fig. 13. In fact, GA can find designs very close to the optimum, exploring only few points in the design space. In the second case studied in this example, considering q = 0.125 MPa and h = 6.3 mm, each of the 14 layers that form the laminate contains only one ply, which must be oriented at 45°, +45° and 90°. These orientations are represented in the chromosomal string by the codes 1, 2 and 3, respectively. The individual ply thickness (0.45 mm) and the size of the design space (2187) are the same of the previous case. The GA parameter are adjusted with the values P = 7, Ne = 1, NLG = 155, NSD = 55, poa = 4% and pgs = 80%, according to the best combination found for the previous case. The reduction of the laminate thickness leads to a stronger nonlinear behavior of the structure, even if the load is reduced by the same factor. As a consequence of this fact, many laminate configurations present in the design space of this example result in structures that cannot be accurately analyzed if the basic increment load ki is taken equal to 0.05, as it was adopted in the previous case. There is not a method to determine a priori the best value for the

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

parameter ki. If the adopted value is too big, the nonlinear solution is inaccurate or even may not converge. In the other hand, if the value is smaller than the necessary, the analysis process may became very slow, which is very undesirable when an optimization by GA is used. The challenge in this case is to find a value for ki that provide accurate analyses for all the designs with the minimum computational cost. In this example a strategy for the automatic refinement of the parameter ki is implemented to reduce the cost of the analyses and to guarantee the convergence and the accuracy of the nonlinear solution. The first step in the implementation is carried out monitoring the structural response during the nonlinear solution and identifying anomalous behaviors that can be interpreted as solution instability or error. When this situation occurs, the process is stopped and restarted using a smaller value for ki. The way the structural response is monitored and the factor of reduction of the parameter ki are problem dependent and can be implemented in different forms. In the present example the error is identified when the increment in the displacement of the central point of the shell in a specific load step is greater than 10 times the maximum displacement verified in the previous load step. The refinement of the parameter ki is accomplished by taking one half of its current value. A new structural analysis in the optimization process using GA to look for a new design is always started using the value of ki defined by the user, even if the refinement process was used in the previous analysis. The results of four optimizations executed for this case are presented in Table 7. The optimizations 1 and 2 were executed using the fixed value ki = 0.025, while the optimizations 3 and 4 were executed using the value ki = 0.05, together with the refinement strategy. In the optimizations 1, 3 and 4 the optimal solution for the problem is given by the laminate [904, 45, 902]S, with the following values for the critical load and maximum displacement: NCcrit = 0.451 and Umax = 0.02895 m. This stacking sequence is almost similar to that obtained for the best design of the previous example. It contains many plies with fiber oriented at 90°, in order to maximize the shell stiffness, with a 45° ply introduced to attend the problem constraint. The best design found by the optimization 2 is the laminate [902, 452, 90, 45, 90]S; however, the nonlinear analysis of this structure, using ki = 0.025, results in a wrong solution, invalidating this optimization. The load level  central displacement curves for the designs found in optimizations 1, 3 and 4, together with the curve referred to the right solution (using ki = 0.0125) of the design found in optimization 2, are shown in Fig. 14. In this figure an additional region with snap-back and snap-through points is observed in the curve corresponding to the best design of optimization 2. This behavior occurs in many designs explored during the optimization process, making the nonlinear analysis of such structures much more difficult since these behaviors are correctly reproduced only when small values of the parameter ki are used. The values adopted for the GA parameters resulted in a good performance in terms of computational cost, considering the small size of the design space in this problem. Furthermore, a great reduction in the number of analyses to be executed in each optimization is obtained avoiding repeated analyses using a memory Table 7 Optimization performance with or without the refinement strategy Optimization 1 2 3 4

ki 0.025 0.05 With refin.

a Ttotal and Tan optimization 1.

An

Ar

Ar/SDS (%)

Ttotala

Tana

483 441 504 651

293 297 259 325

13.4 13.6 11.8 14.8

1.000 0.996 0569 0.682

1.000 0.982 0.643 0.616

dimensionless

values

obtained

with

reference

to

the

453

Fig. 14. Load level  central displacement for the best design obtained in optimizations 1, 2, 3 and 4.

which store previous results. In the optimization 4 the number of analyses is reduced in 50%. Some advantages were observed in the use of the refinement strategy. First, it is possible to reduce the average time of the analyses (Tan) from 1.000 in the optimization 1 (when the refinement strategy was not used) to 0.616 in the optimization 4, resulting in 32% of time saving in the total optimization time (Ttotal). The second advantage lies in the guarantee of the correct solution for the nonlinear analysis of the structure, avoiding the use of a wrong structural response by the GA, as occurred in the optimization 2. 5. Final remarks The GA was successfully applied to obtain the optimal design of composite laminate structures such as plates and shells subjected to different load conditions. Two examples of multiobjective optimization were presented, and some aspects related to the formulation of the objective function and its influence on the optimization process were discussed. The performance of the GA in terms of computational cost and reliability was studied, showing that the method is very efficient in finding near optimal solutions, and an important saving in computer time can be obtained by de use of suitable values for the GA parameters and when results of different analyses are stored. In the third example the GA was used together with a nonlinear FEM analysis to maximize the stiffness of a composite shell. Two cases of laminate thickness and load level were considered to explore different levels of nonlinearity of the structural response. The influence of the population size P and the limit of the number of generations NLG over the reliability and the computational cost of the method were investigated in this example. The results demonstrated that when relatively small populations associated with a large limit of the number of generations are used, better performances of the GA are obtained. A strategy for the automatic refinement of the basic load level (ki), used in the nonlinear analysis, was proposed. It showed to be very important for the improvement of the optimization process, since a sensible reduction of the average time of the structural analysis was obtained.

Acknowledgement The authors wish to thank the Brazilian agencies CNPq and CAPES for its financial support.

454

F.S. Almeida, A.M. Awruch / Composite Structures 88 (2009) 443–454

References [1] Walker M, Smith RE. A technique for the multiobjective optimization of laminated composite structures using genetic algorithms and finite element analysis. Compos Struct 2003;62:123–8. [2] Gürdal Z, Haftka RT, Hajela P. Design and optimization of laminated composite materials. Wiley & Sons; 1999. [3] Soremekun GAE, Gürdal Z, Haftka RT, Watson LT. Composite laminate design optimization by genetic algorithm with generalized elitist selection. Comput Struct 2001;79:131–43. [4] Goldberg DE. Genetic algorithms in search, optimization, and machine learning. Addison-Wesley Publishing Company; 1989. [5] Bathe KJ, Ho L. A simple and effective element for analysis of general shell structures. Comput Struct 1981;13:673–81.

[6] Yang Y, Shieh M. Solution method for nonlinear problems with multiple critical points. AIAA J 1990;28(12):2110–6. [7] Daniel IM, Ishai O. Engineering mechanics of composite materials. Oxford University Press; 1994. [8] Soremekun GAE. Genetic algorithms for composite laminate design and optimization. M.Sc. thesis (Master of Science in Engineering Mechanics). Polytechnic Institute and State University, Blacksbourg, Virginia, USA; 1997. [9] Nagendra S, Jestin D, Gürdal Z, Haftka RT, Watson LT. Improved genetic algorithm for the design of stiffened composite panels. Comput Struct 1996;58(3):543–55. [10] Almeida FS. Laminated composite material structures optimization with genetic algorithms. M.Sc. thesis – PPGEC/UFRGS, Porto Alegre, Rio Grande do Sul, Brazil; 2006 [in Portuguese].