Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

54, 2, pp. 423-436, Warsaw 2016
DOI: 10.15632/jtam-pl.54.2.423

ANALYTICAL DISCRETE STACKING SEQUENCE OPTIMIZATION OF

RECTANGULAR COMPOSITE PLATES SUBJECTED TO BUCKLING

AND FPF CONSTRAINTS

Aleksander Muc, Małgorzata Chwał

Cracow University of Technology, Kraków, Poland

e-mail: olekmuc@mech.pk.edu.pl; mchwal@pk.edu.pl

Aspecial type of newdiscrete designvariables is introduced in order to find optimal stacking
sequences for laminated structures.Using theproposednewdesignvariables,wedemonstrate
how to find analytical (i.e. without anynumerical optimization algorithm) optimal solutions
for laminates made of plies having three different fibre orientations: 0◦, ±45◦, 90◦. It is
proved that the definition of design variables enables us to distinguish two types of optimal
solutions, i.e. unimodal and bimodal ones. The form of optimal stacking sequences affects
the multiplicity (bimodal problems) or uniqueness (unimodal problems) of the solutions.
The decoding procedure between membrane and flexural design variables is also proposed.
The results demonstrate the effectiveness, simplicity and advantages of the use of design
variables, especially in the sense of the accuracy, repeatability of results and convergence of
the method.

Keywords: stacking sequence optimization, composite plates, buckling, first-ply-failure

1. Introduction

Thesuperiormechanical properties of compositematerials suchashigh stiffness,weight ratio and
anisotropic properties that can be tailored through variation of fiber orientations and stacking
sequence give the designer an added degree of flexibility. However, this additional tool should be
used by engineers in a proper manner, i.e. it requires an application of optimization methods.
Formanydesign problems using the 2D approach (beams, plates, shells),most notably for those
where stiffness requires domination, there are multiple designs with similar performance. These
designs may have very different stacking sequences but very similar or almost identical values
of stiffnesses A and D. In such cases, it is important to produce all or most of the design
alternatives.

The effectiveness of optimal design, especially for composite structures, is stronglydependent
on the proper choice of two elements: i) definition of design variables, ii) application of an
appropriate optimization algorithm.

The simplest definition of design variables depends on direct application of real continuous
variables (i.e. fibre orientations θl and thicknesses tl in the l-th ply, l=1,2, . . . ,N). Now, such
an approach is commonly used in finite element codes, such as e.g. ANSYS, ABAQUS etc. This
method is not very convenient formany engineering or analytical applications and is replaced by
the introduction of so-called lamination parameters (Miki, 1986; Fukunaga and Vanderplaats,
1991). The lamination parameters are usually determined for thin walled composite structures
(i.e. beams, plates or shells) with the use of the Love-Kirchhoff kinematical hypothesis. For an
arbitrary laminate, the structural stiffness is characterized by12 independentparameters instead
of 2×N −1 variables for the previous real continuous variables (θi and ti). In fact, analytical



424 A.Muc, M. Chwał

studies deal mainly with the use of four of them, i.e. ξ
{A}
[1]
, ξ
{A}
[2]
, ξ
{D}
[1]
, ξ
{D}
[2]
corresponding

to laminates in which the stiffnesses Bij, A16, A26, D16 and D26 are assumed to be equal
to zero.

The past few decades have seen an increased interest in general-purpose “black-box” optimi-
zation algorithms that exploit limited knowledge concerning the optimization problem onwhich
they are run. In particular, the two most popular black-box optimization strategies, evolutio-
nary algorithms and simulated annealing, mimic processes in natural selection and statistical
mechanics, respectively. Commonly, they are based on the standard use of the lamination para-
meters. In this case it is impossible to finda unique laminate configuration, and, in addition, the
correctness and accuracy of solutions can be verified by the comparison with other numerical
result only.

According to the definitions and classifications of IEEENNC (1996) the four types of algo-
rithms constitute evolutionary computation methods (Genetic Algorithms – GA, Evolutionary
Programming – EP, Evolution Strategies – ES andGenetic Programming – GP) but in general
they are based on the Darwinian concept of evolution. In fact, now hybrids of the four me-
thodologies are becoming most popular. The distinguishing feature of traditional Darwinistic
evolution is selection, the survival of the fittest members of each generation. For composite
materials it is much better to look beyond that view in the sense of the above algorithms
and to explore a new view of evolution that includes natural selection plus self-organization –
see Kauffmann (1993). It is worth to emphasize that Grosset et al. (2006) formulated almost
the same conclusions and stated that it was necessary to abandon partially the Darwinistic
theory of evolution and finally introduced to the analysis of composites the concept of esti-
mation of distribution algorithms. They used a statistical framework to formalize the search
mechanisms.

However, it has become important to understand the relationship between how well an
algorithm performs and the optimization problem on which it is run. In this paper, we present
an analysis that contributes toward such an understanding by addressing questions like the
following: how we can best match design variables and algorithms to the problems, i.e., how
best we can relax the black-box nature of lamination parameters and the algorithms and have
them exploit some knowledge concerning the optimization problem?

In the present paper, we intend to solve the problem of optimal design of bi-axially com-
pressed rectangular multilayered composite plates having discrete fibre orientations in each in-
dividual ply and subjected to buckling andFPF constraints. Contrary to the identical problems
discussed in the literature, we look deeper into the physical problem considered herein.We de-
monstrate that the appropriate definition of design variables allow us to obtain unique, exact
and accurate solutions. Using the proposed new design variables we show how to find analytical
(i.e. without any numerical optimization algorithm) optimal solutions for laminates made of
plies having three different fibre orientations: 0◦, ±45◦, 90◦. For a higher number of different
discrete fibre orientations, we propose the application of an effective numerical algorithm based
on the evolution strategy, see Muc and Muc-Wierzgoń (2012). We explain also the discussed
in the literature problem of multiplicity of optimal solutions. We prove that the multiplicity
is an artificial result since it is caused by wrong coding of laminate configurations and wrong
interpretation of optimal solutions.

The numerical results, presented in the paper, are obtained for a graphite/epoxy resin
material having the following mechanical properties: E1 = 127.59GPa, E2 = 13.03GPa,
G12 = 6.41GPa, ν12 = 0.3, and the thickness of an individual ply in the laminate is equal
to 1.27mm. The ultimate allowable strains are following: εlocal1allowable =0.008, ε

local
2allowable =0.029,

γlocal12allowable =0.015. A safety factor equal to 1.5 is used to calculate the strain allowables.



Analytical discrete stacking sequence optimization... 425

2. Fundamental relations for 2D multilayered composite structures

Fig. 1. Global and local coordinate systems (x1 ≡x′,y1 ≡ y,z≡ z′)

Usually thin plies in the laminate (Fig. 1) can be considered to be under a plane stress with
all the stress components in the out-of-plane direction z being approximately zero. In the 3D
case the generalized Hooke law (the local coordinate system x1x2z associated with fibres) is
reduced to

σ′i(x1,x2,z)=C
′
ijε
′
j(x1,x2,z) i,j =1,2,6 (2.1)

where

C′11 =Q11 =
E1

1−ν12ν21
C′12 =Q12 =

ν12E2
1−ν12ν21

C′22 =Q22 =
E2

1−ν12ν21
C′66 =Q66 =G12

and σ′ denotes the tensor of in-plane stress components, and ε′ the tensor of in-plane strain
components. Let us consider that the ply material axes are rotated by an angle θ with respect
to the global reference system (xyz) – Fig. 1. In the global system, writing the Hooke law

σi(x,y,z) =Cijεj(x,y,z) i,j =1,2,6 (2.2)

and using the Tsai-Pagano invariant formulation, all components of the stiffness matrix C can
be written in the invariant form

C11 =U1+U2cos2θ+U3cos4θ C12 =U4−U3cos4θ

C22 =U1−U2cos2θ+U3cos4θ C16 =
1

2
U2 sin2θ+U3 sin4θ

C26 =
1

2
U2 sin2θ−U3 sin4θ C66 =U5−U3cos4θ

(2.3)

where

U1 =
1

8
(3Q11+3Q22+2Q12+4Q66) U2 =

1

2
(Q11−Q22)

U3 =
1

8
(Q11+Q22−2Q12−4Q66) U4 =

1

8
(Q11+Q22+6Q12−4Q66)

U5 =
1

2
(U1−U4)



426 A.Muc, M. Chwał

Theabove relations aredeveloped for a singleply (lamina).The laminate canbebuilt ofN layers,
see Fig. 2, so that the stresses in the l-th ply are related to the strains in the following way

σ
(l)
i (x,y,z) =C

(l)
ij εj(x,y,z) i,j =1,2,6 l=1,2, . . . ,N (2.4)

assuming that all ply strains are equal to the laminate strains. The stiffness matrix coeffi-

cients C
(l)
ij are defined with the use of Eq. (2.3) where the fibre orientations θ are replaced by

the symbol θl referring to the fibre orientations of the l-th layer. From the assumption that
the strains vary linearly through the laminate thickness, i.e. εi(x,y,z) = ε

0
i(x,y) + zκi(x,y)

(i=1,2,6) one can find that relation (2.4) can be rewritten in the following form

Ni(x,y)=Aijε
0
i(x,y)+Bijκi(x,y)

Mi(x,y)=Bijε
0
i(x,y)+Dijκi(x,y)

i,j =1,2,6 (2.5)

where the in-plane stress resultants Ni(x,y) and the stress couples Mi(x,y) are expressed as

Ni(x,y)=

t/2
∫

−t/2

σ
(l)
i dz=

N
∑

l=1

σ
(l)
i (zl−zl−1)

Mi(x,y)=

t/2
∫

−t/2

σ
(l)
i z dz=

1

2

N
∑

l=1

σ
(l)
i (z

2
i −z2i−1)

i=1,2,6 l=1, . . . ,N (2.6)

where A, B, D are the extensional, coupling and bending stiffnesses, respectively, defined as
follows

Aij =
N
∑

l=1

C
(l)
ij (zl−zl−1) Bij =

1

2

N
∑

l=1

C
(l)
ij (z

2
l −z2l−1)

Dij =
1

3

N
∑

l=1

C
(l)
ij (z

3
l −z3l−1)

i,j=1,2,6 (2.7)

where zl and zl−1 are the location coordinates of the top and the bottom surface of the lamina l.
ε0i(x,y) are the components of the in-plane (membrane) strains, andκi(x,y) are the components
of the vector of curvature (i=1,2,6).

3. Definition of design variables

In the 2D approach, topological variables defining the connectivity of particular structural ele-
ments in the structure (in the paper it denotes the stacking sequence of the individual layers in
the laminate) are understood in the sense of the sequence of layers having prescribed discrete
fibre orientations θl in each individual ply. Commonly, it is assumed that the thicknesses of
individual plies are identical, i.e. tl = t/N – see Fig. 2. In order to assure great flexibility and
generality in the formulation of various optimisation problems, different types of the above-
mentioned discrete design variables must be represented in a similar unified manner, i.e. each
design variable must be coded as a finite string of digits. Let us note that the angle-ply anti-
symmetric laminates are considered only, however, it can be easily extended for an arbitrary
laminate configuration. Using the classical method of coding, 1 represents 0◦2, 2 –±45◦, 3 – 90◦2.
Each design variable s representing the fibre orientation (i.e. 1, 2 and 3) is coded as a binary
number and is called as a gene. The sequence {1,2,1,3} is called a chromosome. Such a repre-
sentation is not very convenient for optimisation problems since there is a lot of design variables



Analytical discrete stacking sequence optimization... 427

Fig. 2. Cross-section of the laminate (N =6)

(increasing with the total number of plies N) and, in addition, various stacking sequences are
described by the identical values of the A, B, D matrices. Therefore, we propose to adopt he-

rein a special type of integer variables x
{A,D}
r (r = 1,2,3) introduced by Muc (1997) that are

completely different than those introduced byMiki (1986), Fukunaga and Vanderplaats (1991).
The new design variables represent triangles in the design space, see Fig. 3. However, there is

Fig. 3. Graphical representation of design variables

no uniquemapping between the spaces xAr and x
D
r . For the assumed laminate configuration, the

B matrix is identically equal to zero, whereas the stiffnesses A and D can be written in the
following way

A11 = t(U1−U3)+
4t

N
U2(x

A
1 −xA3 )+

8t

N
U3(x

A
1 +x

A
3 )

A12 = t(U4+U3)−
8t

N
U3(x

A
1 +x

A
3 )

A22 = t(U1−U3)−
4t

N
U2(x

A
1 −xA3 )+

8t

N
U3(x

A
1 +x

A
3 )

A66 = t(U5+U3)−
8t

N
U3(x

A
1 +x

A
3 )

(3.1)

and

D11 =
t3

12
(U1−U3)+

1

12

(4t

N

)3
U2(x

D
1 −xD3 )+

1

6

(4t

N

)3
U3(x

D
1 +x

D
3 )



428 A.Muc, M. Chwał

D12 =
t3

12
(U4+U3)−

1

6

(4t

N

)3
U3(x

D
1 +x

D
3 )

(3.2)

D22 =
t3

12
(U1−U3)−

1

12

(4t

N

)3
U2(x

D
1 −xD3 )+

1

6

(4t

N

)3
U3(x

D
1 +x

D
3 )

D66 =
t3

12
(U5+U3)−

1

6

(4t

N

)3
U3(x

D
1 +x

D
3 )

where

x{A,D}r =

N/4
∑

k=1

{1, [3k(k−1)+1]}cos(2θ)Ξ(αr)

Ξ(αr)=

{

1 where αr = θ

0 where αr 6= θ

αr =90
◦r−1
2

r=1,2,3

(3.3)

Let us note that using the above notation all terms in the stiffness matrix are uniquely repre-
sented by the set of four integer variables {xA1 ,xA3 ,xD1 ,xD3 }. The terms having the index r = 2
are identically equal to 0 since they correspond to the plies having fibres oriented at 45◦. the
integer numbersxA1 ,x

A
3 represent directly the number of plies with fibers oriented at 0

◦ and 90◦,
respectively, since the following relation is always fulfilled

xA1 +x
A
3 +x

A
2 =
N

4
(3.4)

where xA2 = N45/4, and N45 denotes the total number of plies oriented at 45
◦. With the aid

of Eqs. (3.1)-(3.4) it is possible to define feasible regions for our new design variables – the
triangles presented in Fig. 3. Two sets of variables {xA1 ,xA3} and {xD1 ,xD3 } are not independent,
however, using their definition (Eq. (3.3)) it is possible to evaluate the ranges of their variations
demonstrated in Fig. 3 in the form of quadrilaterals

(xA1 )
3 =

xA
1
∑

k=1

[3k(k−1)+1]¬xD1 ¬
N/4
∑

k=N/4+1−xA
1

[3k(k−1)+1]=
(N

4

)3
−
(N

4
−xA1

)3

(xA3 )
3 =

xA
3
∑

k=1

[3k(k−1)+1]¬xD3 ¬
N/4
∑

k=N/4+1−xA
3

[3k(k−1)+1]=
(N

4

)3
−
(N

4
−xA3

)3

xD1 +x
D
3 ¬
(N

4

)3

(3.5)

Knowing the values of the {xD1 ,xD3 } variables, one can derive from relations (3.1), (3.2) the
upper and lower bounds of the {xA1 ,xA3} variables that have to be integer numbers belonging to
the triangular domain shown in Fig. 3. However, in the optimisation procedure they are treated
as continuous variables since they are always normalised (by the division of them by N/4 and
(N/4)3, respectively) and theybelong to the interval [0,1]. Thenormalised variables are denoted
by the bar over the symbols, i.e. as {xA1 ,xA3}, {xD1 ,xD3 }. The optimal solutions are completely
independent of the total number of plies N and, in addition, the definition of the appropriate
terms in the stiffness matrices A and D has an identical form, although they are functions of
different variables. Such a definition may be also very useful for the pseudorandom generation
of the design variables.
If the optimal solutions are found (in the sense of four variables {xA1 ,xA3}, {xD1 ,xD3 }) the

decoding procedure is required to represent the above-mentioned variables by the appropriate



Analytical discrete stacking sequence optimization... 429

stacking sequences.Thedecodingprocedure canbe easily conductedwith theuse of the symbolic
package Mathematica. The fundamental two operations are given below

a= Table[3∗l∗ (l−1)+1,l,N/4];
f= Subsets[a,L];

First of all, a list of values of the expression “3l(l − 1)+ 1” is generated when the natural
number l runs from 1 to N/4 – see Eq. (2.3). Assuming N = 16, table “a” takes the form:
{1,7,19,37}. Then a finite number of subsets “f” having exactly “L” elements is constructed
from the list “a”. For instance, for L=2 the subsets “f” are following: {1,7}, {1,19}, {1,37},
{7,19}, {7,37}, {19,37}. If L=2 defines the number of plies oriented at 0◦ (i.e. xA1 =2) then
each of the subsets represents the laminates (in fact the upper part of the laminate where the
symbol 0◦ corresponds to the pair of plies having the 0◦ orientation): [0◦,0◦,÷,÷], [0◦,÷,0◦,÷],
[0◦,÷,÷,0◦], [÷,0◦,0◦,÷], [÷,0◦,÷,0◦], [÷,÷,0◦,0◦].According to definition (3.3), all laminates
have the same value xA1 = 2 but different values of x

D
1 equal to: 8, 20, 38, 26, 44, 56, respec-

tively. In addition, the first laminate can be characterized by two different values of xD3 , i.e.:
[0◦,0◦,45◦,90◦]−xD3 =37 or [0◦,0◦,90◦,45◦]−xD3 =19. In both cases, xA2 =xA3 =1.The exam-
ple demonstrates evidently the non-uniqueness of the mapping presented in Fig. 3. However, as
four design variables are known (e.g. xA1 = 2, x

A
3 = 1 and x

D
1 = 8, x

D
3 = 19) the lamination

sequence can be derived uniquely – it corresponds to [0◦,0◦,90◦,45◦]. Therefore, to decode the
lamination sequence from the set of the normalized design variables {xA1 ,xA3}, {xD1 ,xD3 }, it is ne-
cessary to conduct the following operations: Step1) to round the values {xA1 ∗(N/4),xA3 ∗(N/4)}
to the nearest integers, Step2) to find the subsets “f” corresponding to xA1 ≈xA1 ∗ (N/4), com-
pute the sum of each of the elements in the subsets “f” and find the nearest integer to the real
numberxD1 ∗(N/4)3, Step3) to repeat step 2 for the valuexA3 ≈xA3 ∗(N/4), creating new subsets
for the layers oriented at 90◦, but selecting empty spaces only (noted as÷ in the laminate in the
above example), Step4) to fill the rest of empty spaces in the laminate by the layers oriented
at 45◦. To verify the convergence, it is possible to increase the total number of layers N.

The decoding procedure is simple using symbolic packages. In many cases (one of themwill
be discussed further), the optimal design is not represented by all design variables, and the
decodingmethod has to be slightly modified.

4. Buckling and the First-Ply-Failure of plates

Manyexperimental results on thebucklingof compositematerial plates have beenpresentedover
the last years.Theyare summarized inMuc (1988),MucandGurba (2001). In general, they tend
to indicate that the theory for composites is rather in good agreement with experiments. The
experiments demonstrate evidently that for thin-walled flat composite plates, the loss of stability
is not equivalent to the catastrophic failure of structures. The catastrophic failure in form of the
limit carrying capacitymayoccur for thicker plates and it is associatedwith theFirst-Ply-Failure
(FPF). For plates with a cutout, delaminations or stiffened, the failuremodemay be in form of
bifurcation buckling, but the final damage is usually associated with other modes of failure. For
flat plates, the comparison between experimental and theoretical results is conducted with the
use of linear prebuckling theory according to experimental observations. For more complicated
plated structures, the nonlinear prebuckling and postbuckling analysis is required since the final
failure mode is associated, e.g., with local buckling modes or failure of a core for sandwiches.
Therefore, for flatbi-axially compressed rectangular plates, it is assumed that a criticalmultiplier
of loading corresponding to the global loss of stability λb can be expressed in the following form



430 A.Muc, M. Chwał

λb(s)=
(mπ/a)2

Px(1+kβ2m)
[D11+(D12+2D66)β

2
m+D22β

4
m] βm =

na

mb
k=
Py
Px
(4.1)

where a, b are geometrical plate dimensions, andm,n are numbers of half-waves in two perpen-
dicular directions corresponding to the plate co-ordinate system, andPx is the axial compressive
force in the x direction. s denotes the vector of design variables having 2 independent normali-
zed real variables {xAr ,xDr } defined in the interval [0,1] representing 3 different fibre orientations
and various stacking sequences in the laminate. Using the notation introduced in Eq. (4.1), the
critical multiplier of loading can be written as follows

λb =Ωm(Z1+Z2x
D
1 +Z3x

D
3 ) (4.2)

where

xDr =
( 4

N

)3
xDr Ωm =

(nπ/b)2

Px(1+kβ2m)β
2
m

t3

12

Z1 =U1(1+β
2
m)
2+U3(6β

2
m−1−β4m) Z2 =U2(1−β4m)+2U3(1−6β2m+β4m)

Z3 =2U3(1−6β2m+β4m)−U2(1−β4m) r=1,2,3

It is also worth to point out also that the validity of relation (4.1) is strictly limited by the
values of t/[Min(a,b)] ratio. For the ratio higher than 0.05, it is necessary to include transverse
shear effects employing, for instance, theMindlin hypothesis.

It iswell-known thatbuckling loads (4.2) are straight lines in the convexdesign space (Fig. 3).
Since for the knowncompositematerialsU1 is always greater thanU3, the coefficientZ1 is always
positive, whereasZ2 andZ3 maybepositive, equal to zero or negative. According to the classical
theory ofmathematics for the prescribedmode of buckling (i.e.m andn) in the feasible domain
of design variables xD1 , x

D
3 , the maximal value of the parameter λb defined by Eq. (4.2) may

exist at the vertices of the triangle (points A, B, C – unimodal solutions) or along the lines
AB and BC (e.g. points D and E, bimodal solutions – degenerated solutions since one of the
variables xD1 or x

D
3 is equal to zero) – Fig. 4.

Fig. 4. Positions of the optimal buckling loads for three discrete fibre orientations

4.1. Unimodal solutions

The variations of the coefficients Z1,Z2 andZ3 with the geometrical ratios a/b are demon-
strated in Fig. 5. Their values correspond directly to the location of optimal fibre orientations
with respect to the a/b ratio, i.e. xD1 = 1, x

D
3 = 0 – orientation 0

◦; xD1 = 0, x
D
3 = 0 – orienta-

tion 45◦, xD1 =0, x
D
3 =1 – orientation 90

◦.



Analytical discrete stacking sequence optimization... 431

Fig. 5. Variations of the coefficients in Eq. (4.2)

With the use of the introduced design variables, the analytical unimodal optima can also be
easily derived. They can be written in the following way:
— equivalent to the pointA,Z2 > 0 andZ3 < 0

xD1 =
(N

4

)3
xD3 =0 if βm ¬

√

−3ϕ+
√

8ϕ2+1

1−ϕ
(4.3)

— equivalent to the pointB,Z2 < 0 andZ3 < 0

xD1 =0 x
D
3 =0 if

√

−3ϕ+
√

8ϕ2+1

1−ϕ
¬βm ¬

√

3ϕ+
√

8ϕ2+1

1−ϕ
(4.4)

— equivalent to the pointC,Z2 < 0 andZ3 > 0

xD1 =0 x
D
3 =
(N

4

)3
if βm ­

√

3ϕ+
√

8ϕ2+1

1+ϕ
(4.5)

whereϕ=2U3/U2.
In the above relations, the estimations are computed from the equalities:

— equivalent toZ1 =0

λb
[

xD1 =
(N

4

)3
,xD3 =0,m,n

]

=λb[x
D
1 =0,x

D
3 =0,m,n] (4.6)

— equivalent toZ3 =0

λb[x
D
1 =0,x

D
3 =0,m,n] =λb

[

xD1 =0,x
D
3 =
(N

4

)3
,m,n

]

(4.7)

where the first relation describes the equality of buckling loads for the laminates oriented at 0◦

and 45◦, and the second, the equality of buckling loads for the laminates oriented at 45◦ and 90◦.

4.2. Bimodal solutions

The bimodal solutions correspond to situations as the buckling load is identical for two
neighbouring buckling modes (e.g. m and m+1) – see also Muc (1988) (continuous angle-ply
orientations). Let us note that the identical procedure is conducted in the buckling analysis of
isotropic structures although the solutions are not called the bimodal ones. For isotropic plates
constructing the classical “chain curve” for buckling loads versus the a/b ratio for different wave
numbersm, it is possible to find such values of thea/b ratio that correspond to twoneighbouring



432 A.Muc, M. Chwał

Fig. 6. Optimal unimodal and bimodal solutions

buckling modesm andm+1 (e.g. for uniaxially compressed plates: a/b=
√
2 as m=1 and 2,

a/b=
√
6 asm=2 and 3 etc.). Figure 6 shows a part of the chain curve for laminated plates as

well as the optimal unimodal solutions – the lower bound envelope of buckling loads for different
values ofm. The dot represents the classical bimodal solutions.

For laminated plates, the critical multiplier λb (Eqs. (4.1) or (4.2)) is not only a function
of the geometrical a/b ratio but also of the design variables s. Therefore, for the identical a/b
ratio, there may exists a finite number (discrete variables) of bimodal solutions. Some of them
may correspond to the higher values of buckling loads than for unimodal solutions. They will
be called the optimal bimodal solutions (Fig. 6).

The bimodal solution presented as the dot in Fig. 6 demonstrates that the unimodal optimal
solutions exist for two different modes of buckling m = 2 and m = 1, however, each of them
corresponds to different fibre orientations (±45◦ and90◦, respectively). For uniaxial compression
(k=0), themaximal buckling load occurs always for plies having fibres oriented at α2 =±45◦.
For each bucklingmode (considering separately themodesm=1 andm=2), the buckling load
at the vertex B is higher than at the vertices A and C, and it takes a lower value (for m = 1
and a/b < 1.445) than the corresponding buckling load at the vertex B for m = 2 – see Eq.
(4.4). If a/b=1.445, the optimal unimodal fibre orientations switch from±45◦ to 90◦ form=1.
However, the latter case does not satisfy the buckling criterion – it is not the lower boundwith
respect tobucklingmodes since unimodalbuckling loads forplies having the orientation 90◦ have
lower values form=2 than those for the lower bucklingmode (m=1). Of course, the unimodal
solutions for ±45◦ and m=2 cannot be treated as optimal ones because this orientation gives
lower buckling loads form=1. Thus, the bimodal constraint becomes active.

For discrete design variables, it is possible to find two potential candidates for the optimum
considering each wave numbers of buckling, i.e. (m,n) and (m+1,n). Such type of optimal
solutions is represented by the points D and E in Fig. 4. The identical analysis can be carried
out in the similar manner for the (m,n) and (m,n+1) buckling modes. Of course, the points
D and E are not single candidates for the bimodal solutions. The equalisation of the buckling
load coefficients (Eq. (4.1)) for the neighbouring bucklingmodes (i.e.m, n andm+1, n) leads
to a finite number of solutions represented by a straight line in the design space (xD1 ,x

D
3 ). Using

Eqs. (4.6) and (4.7), the bimodal solutions can be found at the boundaries of the triangle only
since those values offer the highest buckling load among them. The proof of this conclusion is a
trivial one since for each bucklingmode the straight lines described by Eq. (4.1) create a family
of parallel lines (parametrized by the value of buckling load) whose maximum occurs at the
opposite vertices of the triangle – one vertex where themaximum occurs at the pointB for the
mode (m+1,n), and the secondmaximum located at the pointC for the mode (m,n).



Analytical discrete stacking sequence optimization... 433

The bimodal solutions can be found analytically from the following relations

xD1 =0

λb(m,n)=Ωm[Z1(m)−xD3 U2(1−β4m)+2xD3 U3(1−6β2m+β4m)]
=λb(m+1,n)=Ωm+1Z1(m+1)

(4.8)

The relations are valid for a/b > 0.7. In the opposite case, the bimodal solutions are located
along the line xD3 =0.

4.3. FPF constraints

For the assumed laminate configuration and considering the membrane state only in the
global coordinate system, the strain tensor is reduced to two nonzero components written in the
following way (see Eq. (2.5)):

ε01 =
λFPFPx(A22−kA12)
A11A22−A212

ε02 =
λFPFPx(kA11−A12)
A11A22−A212

ε06 =0 (4.9)

The strains in the local coordinate systemof the ply having the orientationαr take the following
form

ε′1 = ε
0
1cos

2αr+ε
0
2 sin

2αr ε
′
2 = ε

0
1 sin

2αr +ε
0
2cos

2αr ε
′
6 =(ε

0
2−ε01)sin2αr

(4.10)

For each individual ply, the above mentioned local strains are compared with the allowable
strains along the fibres εlocal1allowable, in the direction perpendicular to the fibres ε

local
2allowable andwith

the shear strains γlocal12allowable, respectively. Thus, for each discrete fibre orientation αr we have
three inequality FPF constraints. However, one can easily find that each FPF relation may be
presented in the identical form

2t(U1+U4)(tQ66+F
A
2 )− (FA1 )2 =λFPFPx(p1+p2FA1 +p3FA2 )

1

εlocal
allowable

FA1 =
4t

N
U2

rmax+1
∑

r=1

xAr cos(2αr) F
A
2 =
8t

N
U3

rmax+1
∑

r=1

xAr cos
2(2αr)

(4.11)

where εlocalallowable denotes the appropriate allowable strains for the plies, and pi (i = 1,2,3) are
appropriate constants derived from equations (4.9) and (4.10). The curves described by Eq.
(4.11) represent ellipses in the design space {xA1 ,xA3}. It is obvious that with the aid of any of
numerical packages (Mathematica, Maple, Matlab, Mathcad) it is possible to compute values
of λFPF (Eq.14) for all values of x

A
1 , x

A
3 as well as the prescribed mechanical and geometrical

properties of the plate and the loading parameter k. Such a procedure can be easily conducted
for all components of strains, assuming different allowable values for tension and compression.
Finally, the results are collected in a table (called as FPFTable) that is parametrized by the
values of xA1 , x

A
3 and the failure mode (Eq. (4.11)).

5. Optimization problem

The optimisation problem is formulated as follows

max
s
(min
m,n
λb) (5.1)



434 A.Muc, M. Chwał

where λb denotes a critical multiplier of loading corresponding to the global loss of stability and
s denotes the set of design variables {xD1 ,xD3 }.
The analysed problem may be subjected to various subsidiary constraints written in the

following form:
—Bimodal constraints

λb[m,n]¬λb[m+1,n] and λb[m,n]¬λb[m,n+1] (5.2)

—FPF constraints

λFPF ¬λb (5.3)

Constraint (5.2) presents two conditions for each wave number in buckling independently, and
the values of buckling coefficients λb are derived from Eqs. (4.1) or (4.2).
The optimization analysis is carried out for the prescribed mechanical constants, the total

number of layers N, the loading parameter k and the geometrical ratio a/b. At the beginning,
the buckling loads are computed for all vertices of the triangle (Fig. 3); they correspond to the
fibres oriented in all plies at 0◦, 45◦ and 90◦, respectively. Those values are collected in a table
(called as UNITable) and parametrized by the valuesm and n (selected in the prescribed range,
let say from 1 to 5). Then, the optimization is divided into four procedures described below.

I. Unimodal optimum. There are two integer numbersm and n such that for all vertices of the
triangle, the values in the UNITable have the global minimum with respect to them. In
addition, condition (3.5)1 is satisfied (comparison with the values in the FPFTable). The
optimal values can be computed from relations (4.3)-(4.5). The laminate stacking sequence
is easily determined.

II. Bimodal optimum. There are two integer numbersm and n such that for all vertices of the
triangle the values in theUNITable have not the globalminimumwith respect to them, i.e.
the bimodal constraints become active and if condition (5.2) is not satisfied (comparison
with the values in the FPFTable) then the optimal normalised design variables {xD1 ,xD3 }
are computed from relations (4.8). It is necessary to decode the results to obtain the
optimum in form of the laminate stacking sequence. To demonstrate it, let assume that
the plate ismade ofN =64 layers, a/b=0.5 and is biaxially compressed, i.e. k=2,n=1.
The optimum occurs for xD1 =1963.307, x

D
2 =3899.693, x

D
3 =0. The domain of possible

variations of thenumber“L” is createdwith thehelpof inequalities (3.5), and let it beequal
to 7. For each generated subset “f”, the sum of the elements in the subset is computed in
the loop and comparedwith the rounded to two natural numbers real values of the optimal
solution (1963 or 1964). Having, for instance, the subset {91,127,169,217,331,397,631}
(the sum is equal to 1963) it is possible to recognize immediately that from the whole
laminate representedbythe table “a” (N/4=16) theelementshaving thenonzerovalues in
the whole table, i.e. {÷,÷,÷,÷,÷,91,127,169,217,÷,331,397,÷,÷,631,÷} are replaced
by the pairs oriented at 0◦, whereas the symbol “÷” by the pairs oriented at ±45◦. It
is worth to add that it is possible to find other 35 laminate configurations (subsets “f”)
that give the same sum 1963. There is multiplicity of the solutions since the optimum is
represented by two variables only.

III. Bimodal and FPF optimum. Similarly as in the previous case II, the bimodal optimum
becomes active, however, the FPF load is lower than the optimum bimodal solution (con-
dition (5.3) is satisfied) and the optimum exists inside the triangles plotted in Fig. 3. The
optimum can be found from the relation

λb[m,n] =λb[m+1,n] =λFPF (5.4)



Analytical discrete stacking sequence optimization... 435

Using the appropriate relation for λb (4.2) and λFPF (4.9), it is possible to express the
above equality constraint conditions in the following form:xD1 = p(x

A
1 ,x
A
3 ),x

D
3 = q(x

A
1 ,x
A
3 )

where p and q are analytical algebraic functions. Inserting those results to the definition
of buckling loads (4.2), it is possible to find themaximal buckling loadwith respect to the
values xA1 ,x

A
3 searching for the maximum by building the table for all possible variations

of the values xA1 , x
A
3 ((N/4)[(N/4)−3]+2 possible values inside the half of the triangles)

or using theMathematica procedure “Maximize”. For the optimal values of xA1 , x
A
3 design

variables, it is possible to derive the optimal values ofxD1 ,x
D
3 andfind the optimal stacking

sequences with the aid of the decoding procedure presented in Section 3. The discussed
problem of both FPF and bimodal active constraints occurs, e.g. for plates: a/b = 4,
k=0.25.

IV. FPF optimum. If the FPF is the dominant failuremode (thick plates), the optimum can be
found directly from relation (4.9) searching for themaximum similarity as in the previous
case. Since the optimum is a function of themembrane parameters, one can observe again
the multiplicity of optimal laminate configurations.

6. Concluding remarks and further works

In the paper, the proposal of a new discrete design variables that are used in the buckling and
FPF optimisation problems is shown. It is demonstrated that in the buckling and FPF analysis
of plates, four types of solutions may exist, i.e.:

• Unimodal buckling solutions – unique determination of optimal stacking sequences.
• Bimodal buckling solutions – multiplicity of optimal stacking sequences (bending state
only).

• FPFandbimodal buckling solutions – unique determination of optimal stacking sequences
(flexural andmembrane design variables).

• FPF –multiplicity of optimal stacking sequences (membrane state only).

Inouropinion, thenewset of designvariables demonstrates a lot of theadvantages in comparison
with the existing ones, i.e.:

• It allows us to derive unique analytical solutions.
• It explains the uniqueness and multiplicity of possible laminate configurations being a
representation of laminate discrete configurations.

• It shows a simple symbolic method for the derivation of multiple optimal configurations
from the used set of discrete design variables (the so-called decoding method).

The analysed example (buckling and FPF of bi-axially compressed plates) is relatively simple
since the optimisation problem is characterized by flexural andmembrane properties only, given
byanalytical formulae.Thenext step is connectedwith simultaneous application of theproposed
design variables to optimisation problems characterised by analytical solutions that include both
bending,membrane and coupling effects. Finally, we intend to include theproposedmethodology
into numerical FE codes in the similar manner as described byMuc and Gurba (2001).

Acknowledgement

The National Science Center (Poland) UMO-2013/09/B/ST8/00178 is gratefully acknowledged for

the financial support.



436 A.Muc, M. Chwał

References

1. Fukunaga H., Vanderplaats G.N., 1991, Stiffness optimization of orthotropic laminated com-
posites using lamination parameters,AIAA Journal, 29, 641-648

2. Grosset L., LeRiche R., Haftka R., 2006, A double-distribution statistical algorithm for
composite laminate optimization, Structural and Multidisciplinary Optimization, 31, 49-59

3. IEEENeuralNetworkCouncil,Glossary ofEvolutionaryComputation terms, StandingCommittee
on Standards, IEEE, Piscatway, NY, 1996

4. KauffmanS.A., 1993,TheOrigins ofOrder: Self-Organization and Selection inEvolution,Oxford
University Press, NewYork

5. Miki M., 1986, Optimum design of fibrous laminated plates subjected to axial compression,Pro-
ceedings of 3rd Japan-US Composite Materials Conference, Tokyo, 673-81

6. Muc A., 1988, Optimal fibre orientations for simply -supported plates under compression,Com-
posite Structures, 9, 161-172

7. MucA., 1997,A new set of design variables in stacking sequence optimization under buckling and
strength constraints,Proceedings of WCSMO-2, 2, 693-698

8. Muc A., Gurba W., 2001, Genetic algorithms and finite element analysis in optimization of
composite structures,Composite Structures, 54, 275-281

9. Muc A., Muc-Wierzgoń M., 2012. An evolution strategy in structural optimization problems
for plates and shells,Composite Structures, 94, 1461-1470

Manuscript received May 17, 2015; accepted for print August 13, 2015