Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

51, 4, pp. 1039-1052, Warsaw 2013

PLY THICKNESS TOLERANCES IN STACKING SEQUENCE

OPTIMIZATION OF MULTILAYERED LAMINATE PLATES

Jarosław Latalski

Lublin University of Technology, Department of Applied Mechanics, Lublin, Poland

e-mail: j.latalski@pollub.pl

The paper deals with the impact ofmanufacturing tolerances of plies thicknesses on optimal
design of multi-layered laminated plates in compression. It is assumed that the considered
tolerances are represented by the maximum acceptable deviation of every individual ply
thickness from its nominal design value. The robustness of the optimum is achieved dimi-
nishing the buckling load amplitude factor by the product of arbitrary assumed tolerances
and appropriate sensitivities. The discussed optimization problem is solved numerically by
the direct enumerationmethod. The proposed approach is illustrated with examples of the
rectangular multi-layered laminated plate design under uni- and biaxial compression. The
achieved results emphasise the robustness of the proposedmethod compared to the appro-
aches with ignored tolerances.

Key words: laminate composite structures, optimization, manufacturing tolerances, robust
design, structural stability

1. Introduction

Laminated composites are commonly used structural materials, primarily due to their excellent
strength and stiffness to weight ratios. An additional advantage of these materials over conven-
tional isotropic ones is a possibility of tailoring and optimizing their properties to match the
specific requirements of a given application. In a general case, design variables available to the
engineer include multiple material systems, plies orientations, thicknesses and their stacking
sequence.

Numerous research papers deal with optimum design of composite plates, shells and simi-
lar structures. A vast majority of them are prepared on the assumption that design variables
and thus the system performance are not subject to any deviations arising from manufactu-
ring processes, simplifications in operating conditions, etc. Unfortunately, resulting from this
approach optimal solutions might violate the constraints imposed on the system performance
when real designs are considered. The proximate reason is a high sensitivity of the optimized
structure to variations of design variables, and the fact that in optimal light-weight systems at
least one constraint reaches its limit value.These observations seem to beparticularly important
for composite materials which experience larger uncertainties in their properties if compared to
conventional materials. This is due to a number of parameters involved in their fabrication and
manufacturing process.

Considering the above aspects, the expected uncertainties should be recognised and taken
into account in the design code. The original approach to thismatter is represented by thewell-
-known concept of safety (or load) factors.With respect to composites, this has beendiscussed in
numerous papers, e.g. byChryssanthopoulos andPoggi (1995) to estimate limit load of cylinders
made of a lay-upKevlar composite. Although safety factors provide some latitude for variations
between the actual structure and the specific instance considered in the computational analysis,
this approach does not give any further detailed information about the design and on impact



1040 J. Latalski

of uncertainty factors. Therefore, several more advanced approaches are being developed to face
the discussed problem. All of them can be captured into a robust design philosophy, and most
popular ones are stochastic analysis, interval programming and deterministic approach.

The concept of the stochastic based approach to laminate stability estimation is studied in
several papers, e.g. by Kogiso et al. (1997). The authors derive a reliability analysis considering
material properties, orientation angles and applied loads as independent random variables. To
illustrate the method, numerical examples are given for a simply-supported 8-ply fiber lami-
nated plate made of epoxy graphite, where the mean orientation angles are treated as design
variables. Singh et al. (2001) present the application of the stochastic finite element approach
to the eigenvalue problem of buckling arising from the dispersion of laminamaterial properties.
Uncertain material data are modelled by random variables, while the other system parameters
stay nominal. A mean-centred first order perturbation technique is used to find statistics of
buckling loads of cylindrical panels with various boundary conditions.

The stochastic optimization is an accurate and exact approach to the discussed problem
since it guarantees the assumed a priori level of system safety. Nevertheless, taking into account
practicing engineers’ expectations, this approach may have some drawbacks. The main reason
could be the requirement of precise information on probability density functions and incomplete
information about variables that are critical to the structure. Other issue is high complexity of
the problem final formulation and – as a consequence – resulting tough calculations. Because of
the above reasons, the alternative methods of robust design are under parallel development.

Interval programming requires only the bounds of the perturbed parameters to model the
systemuncertainty andno information aboutprobability distributions is necessary.With respect
to compositematerials, thismethod is studiedbye.g. Jiang et al. (2008). Theauthors formulated
a problem to find the laminate plies orientations to maximise the stiffness of a plate if material
properties are perturbed; next, the problem has been solved by a general nonlinear interval
programming method with an arbitrary chosen possibility degree (accuracy). Three numerical
examples are given. First two deal with the uncertain engineering constants of the laminate
layermaterial. In the third example, themethod is extended to incorporate perturbations in the
thicknesses of the laminas as well.

The proposed interval number programming approach makes the uncertainty analysis co-
nvenient and computationally attractive. Comparing to the stochastic programming methods,
the discussed approach can deal with problems of limited information about the uncertainty.
However, examples given in the discussed paper show that the optimal design strongly depends
on the arbitrary predetermined possibility degree level. Therefore, this value must be chosen
appropriately, based on the class of the solved problem, the attitude of the engineer and his
experience.

The third type of approach to the robust design is deterministic one which takes advantage
of complicated and uncertain processes avoidance. The research effort is aimed atmore in-depth
analysis of a considered systembut also at developing approximate simplifiedmodels as well. As
an instance Fong et al. (2010) suggest taking into account second order effects into the design
codeof composite columns.Presented examples prove amoreaccurate reflection of thebehaviour
of the member and a considerable simplification in final calculations for robust design.

An alternative and simplifiedmodel of uncertainties representation in the optimization code
is given by Lee and Park (2001). The authors proposed a new objective function defined as
a linear combination of the mean and standard deviation of the original cost function. The
introduction of a weighting factor into the cost function enables a change in the priority of the
optimum– i.e. to have the lowmean value of the objective function and simultaneously larger its
standard deviation or to shift the solution to the region with a higher cost functionmean value,
but smaller deviations. Among numerous benefits, this method has an advantage in providing



Ply thickness tolerances in stacking sequence optimization ... 1041

full information about the final solution. Nevertheless, the mean value and standard deviations
of design variables still have to be known or assumed.

An interesting paper was published by Kristinsdottir et al. (1996). They propose their own
technique for searching near-optimal safe designs. This is done by changing the right hand side
of the inequality constraints and replacing the initial zero value in the nominal problem by
a small positive number called a safety margin. Next the structure is re-optimized and checked,
if new design values compared to the design variables in the original unperturbed problem have
reached initially specified tolerances. If not, the appropriate safety margin has to be increased
and the optimization procedure is run again. As an example, a hat-stiffened composite panel is
optimized. The discussed method may ensure safe results for complicated problems, but since
the assumed safety margin is not explicitly correlated to the uncertainties of design variables it
leads to excessive near-optimal solutions only. Also, if re-optimization procedures are required,
the choice of the appropriate safety margins to be increased faces some trouble.

Anothermethod for optimum robust design of laminated plates was proposed byWalker and
Hamilton (2005, 2006). The key idea of the method is the triple optimization of the structure.
First the system is optimized without considering deviations in a design variable (angle ply for
presented exemplary designs). Next, the optimization is re-run twice for the perturbed variable
– i.e. for the upper and lower possible design variable value. Following these steps, a comparison
of results of three objective functions is made to get the best solution. As an example, eight-
-layered plates with different side aspect ratios and loadings are optimized. Unfortunately, the
given method is limited to designs, where the only one reinforcement orientation in all plies is
allowed. Therefore, no solutions for multiple amgle-ply oriented laminates are possible.

The deterministic approach to the discussed robust design problem was also proposed by
Gutkowski and Bauer (1999) and next by Gutkowski and Latalski (2003). The authors suggest
to diminish the nominal limit values of state variables in design constraints by positive terms.
These are the products of assumed arbitrary design tolerances and a vector of appropriate state
variables sensitivities. In the paper byGutkowski and Latalski (2003), themethod is illustrated
by composite plates optimum design with tolerances in fibers orientations.

In the current paper, the presented idea is further extended to cases where laminae thickness
tolerances are considered. Updated formulas for buckling load factor in the composite optimiza-
tion problemare derived.Recurrent sensitivity relations are also formulated.Based on these, the
impact of thickness variations on the optimal laminate stacking sequence solution is examined
in detail. The presented method is offered as a deterministic and computationally effective but
an alternative approach to the exact original stochastic problem. Moreover, the proposed me-
thod overcomes certain limitations of other deterministic methods discussed above and directly
correlates manufacturing tolerances to the structural response providing in-depth information
about their actual impact. And finally, the proposedmethodmight be extended to any type of
other possible imperfections.

2. Problem statement

The impact ofmanufacturing imperfectionsdefinedby lamina tolerances on composite structures
design is presented. In this study, themethod is illustrated by optimization of a rectangular and
simply supported sandwich panel subjected to compressive forces. However, this approach could
be also applied to other structures, such as composite beams or shells and other sources of
possible perturbations (material, geometric, etc.).

It is assumed that the value of an individual lamina thickness in the compositemay bevaried
from its nominal dimension tk. This variation corresponds to the accuracy of manufacture and
is represented by a deterministic maximum allowable deviation ∆tk (lower or upper). This



1042 J. Latalski

tolerance is assumed to be set up arbitrary as e.g. a ratio of the nominal lamina thickness tk.
Therefore, if upper and lower deviations stay the same, it is expected that the actual thickness
of every k-th ply stays within a range from (tk−∆tk) to (tk+∆tk).

Since for the considered tolerances only nominal values of design parameters and their ma-
ximum acceptable deviations are known, the exact magnitude of the system performance is
unknown. Therefore, it is estimated by calculating the state variable λ for nominal values of
design parameters andnext diminishing this by apenalty term,which is a product of thicknesses
admissible imperfections ∆t and appropriate sensitivities dλ/dt. Since the actual deviations are
of unknown sign (upper or lower) and also the sensitivities might be of positive or negative sign,
the moduli of successive ∆tkdλ/dtk products are considered in resultant summation

∆λ=
N
∑

k=1

∣

∣

∣

dλ

dtk
∆tk
∣

∣

∣
(2.1)

Finally, the obtained difference of the nominal system performance and the mentioned penalty
term is considered to be the approximation of the imperfect system performance. This approach
ensures the design to stay on the safe side. As it has been shown in previous research papers by
Gutkowski andBauer (1999) andGutkowski and Latalski (2003) the use of the penalty term as
proposedabove results inmore effective solutions if compared to simpleworst case designs,where
the system safety is ensured by considering the structure with all nominal design variables tk
arbitrarily increased (or decreased) by ∆tk.

The imperfect composite material as described above is used for a rectangular and simply
supported on all four edges sandwich panel given in Fig. 1. The plate consists of N plies in
total, each of the equal nominal thickness t= t1 = . . .= tN with fibers orientation denoted by
angle θk (k =1, . . . ,N). The uniform longitudinal stress resultants λNx and λNy are applied
at the edges of the panel, where λ is the amplitude parameter; no shear forces are considered.

Fig. 1. Laminate sandwich plate

In thepresent research, the lamina stacking sequence tomaximise the lowest (over all possible
modes) value of estimated buckling load λ is looked for. The admissible direction of fibers in
every lamina is limited to one of the four angles: 0◦, 90◦ and ±45◦ and they are not subject to
any deviations. It is also assumed that the sandwich plate is symmetric and balanced, i.e. the
number of (+45◦) plies is equal to the number of (−45◦) ones. Moreover, all the assumptions
of classical buckling theory for plates are in force.

3. Optimization problem formulation

Following theabove considerations, thediscussedoptimumdesignproblemwith lamina thickness
tolerances is formulated as follows:



Ply thickness tolerances in stacking sequence optimization ... 1043

1) find the vector θ= [θ1,θ2, . . . ,θN/2]
T,

2) such, that

max
θ

(

λcr(m,n)−

N/2
∑

k=1

∣

∣

∣
∆tk
dλcr
dtk

∣

∣

∣

)

(3.1)

3) with constraints

θk ∈ (0
◦,90◦,45◦,−45◦) where k=1, . . . ,N/2

N/2
∑

k=1

(θk =45
◦)=

N/2
∑

k=1

(θk =−45
◦)

(3.2)

Set of constraints (3.2)1 limits the fiber orientation to admissible values only, while relation
(3.2)2 represents the condition of the plate to be balanced. Plate symmetry requirement is
guaranteed by the fact that the summation index k stays within the range (1, . . . ,N/2) in
relations (3.1) and (3.2). Appearing in (3.1) (m,n) parameters aremode shape indices to ensure
all possible buckling modes are examined.

The plate load factor λ is evaluated according to the classical buckling theory for simply
supported equivalent orthotropic plate subject to in-plane loading. The system stability is as-
sured for the load amplitude λ=λcr(m,n) sustaining the relation as given byAdali and Duffy
(1990), Haftka andWalsh (1992), Nemeth (1986)

1¬λcr(m,n)=
D11m

4+2(D12+2D66)m
2n2(a/b)2+D22n

4(a/b)4

m2Nx+n
2(a/b)2Ny

(3.3)

where m andn are natural numbers corresponding to the number of buckling half-waves (mode
shape) in the x and y directions respectively.

Variables D11, D12, D22 and D66 are flexural stiffnesses and can be expressed in terms of
material invariants Ui (i=1, . . . ,5) and three integrals V0, V1, V3 as follows

D11 =U1V0+U2V1+U3V3 D12 =U4V0−U3V3

D22 =U1V0−U2V1+U3V3 D66 =U5V0−U3V3
(3.4)

Invariants Ui (i=1, . . . ,5) are related to composite material properties

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

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

(3.5)

where

Q11 =
E1

1−ν12ν21
Q12 =

ν12E2
1−ν12ν21

=
ν21E1
1−ν12ν21

Q22 =
E2

1−ν12ν21
Q66 =G12

(3.6)

In the above relations, E1,E2, ν12, ν21 and G12 denoteYoung’smoduli in the fibers direction 1
and lateral direction 2 (see Fig. 1), Poisson’s ratios and shear modulus, respectively.



1044 J. Latalski

The variables V0, V1 and V3 comprise information about fibers orientations, about layers
thicknesses and their stacking sequence

V0 =

h/2
∫

−h/2

z2 dz=
1

3

N
∑

k=1

(z3k −z
3
k−1)=

2

3

N/2
∑

k=1

(z3k −z
3
k−1)

V1 =

h/2
∫

−h/2

z2cos2θ dz=
2

3

N/2
∑

k=1

(z3k −z
3
k−1)cos2θk

V3 =

h/2
∫

−h/2

z2cos4θ dz=
2

3

N/2
∑

k=1

(z3k −z
3
k−1)cos4θk

(3.7)

where h=Nt is the total thickness of the laminate, z is the distance from symmetry plane as
given in Fig. 1.

4. Sensitivity of the buckling load

Following the proposed in this research approach to robust design – as given in ”Problem state-
ment” section – the systemperformance penalty term ∆λ (2.1) is to be calculated. This requires
the information about the sensitivity of the buckling load factor λcr(m,n)— see Eq. (3.3) with
respect to each lamina thickness

dλcr
dt
=
[dλcr
dt1
, . . . ,

dλcr
dtN

]T
=π2

dD11
dt
m4+2

(

dD12
dt
+2dD66

dt

)

m2n2
(

a
b

)2
+ dD22
dt
n4
(

a
b

)4

m2Nx+n2
(

a
b

)2
Ny

(4.1)

According to (3.3) the appropriate derivatives of flexural stiffnesses are to be calculated

dD11
dt
=U1

dV0
dt
+U2
dV1
dt
+U3
dV3
dt

dD12
dt
=U4

dV0
dt
−U3
dV3
dt

dD22
dt
=U1

dV0
dt
−U2
dV1
dt
+U3
dV3
dt

dD66
dt
=U5

dV0
dt
−U3
dV3
dt

(4.2)

To derive the above relations a simple case is analysed first.

4.1. Simple example

Let us consider a 6-layered composite plate fulfilling all the assumptions given in ”Problem
statement” section. The only exception is that this time all the layers are temporarily assumed
to have different thicknesses. Thus the vector of nominal thicknesses is t= [t1, t2, t3]

T since the
symmetry condition is still preserved.

According to (3.7), the integrals Vi for the discussed simple composite are

V0 =

h/2
∫

−h/2

z2 dz=
2

3

3
∑

k=1

(z3k −z
3
k−1)=

2

3
(z31 +z

3
2 −z

3
1 +z

3
3 −z

3
2)=

2

3
z33



Ply thickness tolerances in stacking sequence optimization ... 1045

V1 =

h/2
∫

−h/2

z2cos2θ dz=
2

3

3
∑

k=1

(z3k −z
3
k−1)cos2θk

(4.3)

=
2

3
[z31 cos2θ1+(z

3
2 −z

3
1)cos2θ2+(z

3
3 −z

3
2)cos2θ3]

V3 =

h/2
∫

−h/2

z2cos4θ dz=
2

3

3
∑

k=1

(z3k −z
3
k−1)cos4θk

=
2

3
[z31 cos4θ1+(z

3
2 −z

3
1)cos4θ2+(z

3
3 −z

3
2)cos4θ3]

Therefore, bearing inmind that h/2= z3 = t1+t2+t3, thederivatives
dV0
dt
=
[

dV0
dt1
, dV0
dt2
, dV0
dt3

]T

requested by (4.2) are given by the subsequent relations

dV0
dt1
=
d

dt1

(2

3
z33

)

=
2

3

d

dt1
(t1+ t2+ t3

)3
=2(t1+ t2+ t3)

2 =
h2

2

dV0
dt2
=
d

dt2

(2

3
z33

)

=
2

3

d

dt2
(t1+ t2+ t3)

3 =2(t1+ t2+ t3)
2 =
h2

2

dV0
dt3
=
d

dt3

(2

3
z33

)

=
h2

2

(4.4)

Sensitivity of V1 variable to t1 may be calculated keeping inmind that z1 = t1, z2 = t1+t2
and z3 = t1+ t2+ t3; thus

dV1
dt1
=
2

3

d

dt1

{

t31cos2θ1+[(t1+ t2)
3− t31]cos2θ2+[(t1+ t2+ t3)

3− (t1+ t2)
3]cos2θ3

}

=2
{

t21cos2θ1+[(t1+ t2)
2− t21]cos2θ2+[(t1+ t2+ t3)

2− (t1+ t2)
2]cos2θ3

}

(4.5)

The samemay be done for t2 and t3 parameters

dV1
dt2
=
2

3

d

dt2

{

t31cos2θ1+[(t1+ t2)
3− t31]cos2θ2+[(t1+ t2+ t3)

3− (t1+ t2)
3]cos2θ3

}

=
dV1
dt1
−2t21(cos2θ1− cos2θ2)

dV1
dt3
=
2

3

d

dt3

{

t31cos2θ1+[(t1+ t2)
3− t31]cos2θ2+[(t1+ t2+ t3)

3− (t1+ t2)
3]cos2θ3

}

=2(t1+ t2+ t3)
2cos2θ3 =

dV1
dt2
−2(t1+ t2)

2(cos2θ2− cos2θ3)

(4.6)

Similar recurrent relations may be formulated for V3 derivative – the only difference being
an angle θ factor 2 is to be replaced by 4.

4.2. General relations

Following the above simple case, general relations for dVi/dtk, i=0,1,3, k=1, . . . ,N terms
are formulated.
Set of equations (4.4) remains valid for every ply in the laminate, so

dV0
dtk
=
1

2
h2 for k=1,2, . . . ,N ⇔

dV0
dt
=
[h2

2
, . . . ,

h2

2

]T
(4.7)



1046 J. Latalski

Relations (4.5) and (4.6) for the integral V1 derivativesmay be simplified and generalised taking
into account the initial assumption t1 = t2 = t3 = t

dV1
dt1
=2[t2 cos2θ1+(4−1)t

2cos2θ2+(9−4)t
2cos2θ3] = 2t

2
3
∑

k=1

[k2− (k−1)2]cos2θk

dV1
dt2
=
dV1
dt1
−2t2(cos2θ1− cos2θ2)

dV1
dt3
=
dV1
dt2
−2 ·4t2(cos2θ2− cos2θ3)

(4.8)

Finally, one arrives at the recurrent dependence for individual terms of dV1/dt and dV3/dt
vectors

dV1
dtk
=
dV1
dtk−1

−2(k−1)2t2(cos2θk−1−cos2θk)

dV3
dtk
=
dV3
dtk−1

−2(k−1)2t2(cos4θk−1−cos4θk)



















for k=2, . . . ,N/2

dV1
dt1
=2t2

N/2
∑

k=1

[k2− (k−1)2]cos2θk
dV3
dt1
=2t2

N/2
∑

k=1

[k2− (k−1)2]cos4θk

(4.9)

5. Numerical examples

To illustrate the proposedmethod of incorporating thicknesses uncertainties into a design code,
examples of uniaxial and biaxial plate compression are presented. Both of them come from the
paper by Haftka andWalsh (1992), where they were discussed for nominal design values only.

Computations are performed for the plate consisting of N = 16 plies in total, each of the
equal thickness t = 0.127mm. The axial stress resultant Nx = 175N/m is constant. Material
properties for the graphite-epoxy laminate are considered: E1 = 128.0GPa, E2 = 13.0GPa,
G12 =6.4GPa and ν12 =0.3.

Both the design statement and limitations imposed on admissible fiber orientations as well
result in a combinatorial formulation of the plate optimization task. This class of problemsmay
be solved by any of the integer-programmingmethods, but such an approach does not guarantee
the global optimum to be found.Moreover, since the trigonometric functions appearing in plate
analysis are periodic, multiple equivalent solutions are expected for the given problem. Bearing
this inmind, and the fact that the total number of all possible combinations is not big – e.g. for
a 16 layer laminate it is 48 – all further considered examples are solved by a direct enumeration
approach.

5.1. Axial compression

For uniaxial compression, no lateral stress resultants Ny are present and the buckling load
factor ismaximised for variousplate aspect ratios a/b. It is known(e.g. fromHaftka et al. (1990))
that for low aspect ratios the optimum ply angle is 0◦, and for a/b larger than about 1.4–1.5
the optimum solution has ±45◦ plies only. All these observations deal with nominal designs.
Therefore, a check is performed to see whether there is any transition range of a/b where the
optimumwould include both 0◦ and ±45◦ plies. Alternatively, whether there is any plate where
for different ∆t uncertainties these combined orientations might appear. As demonstrated by
a series of performed numerical tests, all solutions for a/b < 1.445 cases are a combination
of ±45◦ plies only, while for a/b­ 1.446 ratios all the plies are 0◦. More detailed results of the
analysis are gathered in Table 1.



Ply thickness tolerances in stacking sequence optimization ... 1047

Table 1.Laminate plate buckling load in uniaxial compression versus considered allowable plies
thickness tolerances and a/b aspect ratio. Result notation: λn.s.(m,n) –where n.s. is the number
of equivalent solutions, m,n are bucklingmode shape parameters – see Eq. (3.3). Solutions for
all a/b< 1.445 are a combination of ±45◦ plies, while for a/b­ 1.446 all the plies are 0◦

a/b
∆t/t [%]

0.00 1.00 2.00 3.00 5.00

0.5 173.94970 (2,1) 164.51070(2,1) 155.07070 (2,1) 145.63170(2,1) 126.753702,1)
1.0 43.48770(1,1) 41.12770(1,1) 38.76870(1,1) 36.40870(1,1) 31.68870(1,1)
1.4 23.79870(1,1) 22.48370(1,1) 21.16770(1,1) 19.85170(1,1) 17.22070(1,1)
1.44 22.77270(1,1) 21.50970(1,1) 20.24670(1,1) 18.98470(1,1) 16.45870(1,1)
1.445 22.65070(1,1) 21.39470(1,1) 20.13770(1,1) 18.88170(1,1) 16.36870(1,1)
1.446 22.6431(1,1) 21.3881(1,1) 20.1321(1,1) 18.8761(1,1) 16.3651(1,1)
1.45 22.6251(1,1) 21.3721(1,1) 20.1191(1,1) 18.8661(1,1) 16.3591(1,1)
2.0 21.1431(1,1) 20.0961(1,1) 19.0501(1,1) 18.0031(1,1) 15.9101(1,1)

Looking atTable 1, it seems that if a discussed transition range exists, it is extremely narrow,
since even the change in third decimal digit does not alter the optimal solution. Moreover,
no changes in the number of equivalent solutions and buckling mode in relation to thickness
imperfections are observed.

As a next step, the obtained results are compared with each other. Therefore, the relative
decrease ∆λcr in the buckling load is calculated where the nominal solution is considered as the
reference level λcrnom, where ∆λcr =(|λcr−λcrnom|)/λcrnom. The results for different tolerance
levels ∆t/t and selected a/b ratios are given in Fig. 2. It is clear that the observed change is
linear, since there are no combined (0◦,±45◦) plies in any of analysed plates. The rate of the
buckling load decrease is different, although all the lines stay close to each other. The obtained
results indicate a high sensitivity of the optimum solution vs. the assumed lamina thickness
tolerance level.

Fig. 2. Relative change of the buckling load factor (∆λcr [%]) for different a/b aspect ratios and
tolerance levels ∆t/t. The results are calculated with respect to the nominal (ideal) solution

5.2. Biaxial compression

For the purpose of plane stress analysis, the plate dimensions a=50.8cm and b=25.4cm
(b/a = 0.5) are fixed. The lateral loading Ny (see Fig. 1) is varied within a range
(21.875 ¬ Ny ¬ 1225N/m) implying in Ny/Nx ratio changes as 0.125, . . . ,7.0. The results
of the performed analysis are collected in Table 2 and in Fig. 3, where exemplary optimum
configurations for different load ratios and tolerance levels are shown.



1048 J. Latalski

Table 2. Results of the buckling load factor λcr calculations for different levels of lamina
thickness t tolerance and different compression loads ratios. Result notation: λn.s.(m,n) – where
n.s. is the number of equivalent solutions, m,n are the buckling mode shape parameters

Ny/Nx
∆t/t [%]

0.00 1.00 2.00 3.00 5.00

0.125 154.62170 (2,1) 146.23170(2,1) 137.84070 (2,1) 129.00070(1,1) 111.75220 (1,1)
0.15 149.00020 (2,1) 140.43120(1,1) 132.08320 (2,1) 123.93920(1,1) 106.70820 (1,1)
0.20 137.59520 (2,1) 129.95020(2,1) 122.27720 (1,1) 114.47720(1,1) 98.87720(1,1)
0.24 129.72840 (1,1) 122.47240(1,1) 115.21540 (1,1) 107.95940(1,1) 93.44540(1,1)
0.25 127.9466(1,1) 120.8106(1,1) 113.6746(1,1) 106.4806(2,1) 92.0386(2,1)
0.50 94.63212(1,1) 89.20112(1,1) 84.00424(2,1) 78.78212(1,1) 67.95120(1,1)
1.00 62.2182(1,1) 58.5646(2,1) 55.3096(2,1) 51.7016(2,1) 44.8252(2,1)
1.50 46.3552(1,1) 43.6852(1,1) 41.1446(1,1) 38.5654(1,1) 33.3116(1,1)
2.00 36.9374(2,1) 34.8962(1,1) 32.6792(1,1) 30.7342(1,1) 26.5996(1,1)
2.10 35.4942(2,1) 33.4112(1,1) 31.5322(1,1) 29.5442(1,1) 25.4992(2,1)
2.40 31.7602(1,1) 29.9632(1,1) 28.1534(2,1) 26.4082(1,1) 22.8604(2,1)
2.50 30.6952(2,1) 28.8732(1,1) 27.2582(1,1) 25.5492(1,1) 22.0806(2,1)
2.80 27.7281(1,1) 26.3132(2,1) 24.6712(1,1) 23.2042(2,1) 20.0742(1,1)
3.50 22.5521(1,1) 21.4361(1,1) 20.3201(1,1) 19.1052(1,1) 16.5342(1,1)
7.00 11.6651(1,1) 11.0881(1,1) 10.5101(1,1) 9.9331(1,1) 8.7771(2,1)

Fig. 3. Optimal ply stacking sequences for different load aspect ratios Ny/Nx and tolerance ∆t/t levels

Following Table 2, it can be noticed that the change in the thickness tolerance level results
in the buckling load decrease; also changes of the plates buckling mode shape and number of
equivalent solutions are observed. These aremainly related to the number of ±45◦ plies present
(Fig. 3), since the +45◦ and −45◦ plies may be shifted with each other with no impact on the
buckling load – cosine function being even.

As arises from Fig. 3, for a very low ratio of biaxial compression (Ny/Nx = 0.125) the
optimum solution is similar to the uniaxial case presented in Example 1. The difference is that
increasing the tolerance ∆t level makes the 90◦ oriented plies to appear and to stay close to
the symmetry plane. For intermediate Ny/Nx load ratios, mixed solutions are observed. In the



Ply thickness tolerances in stacking sequence optimization ... 1049

case of nominal solutions composed with mostly 90◦ plies (Ny/Nx = 2.0, . . . ,2.4) the increase
in the thickness uncertainty level ∆t makes the existing ±45◦ layers to shift towards the plate
outer surface or even force new ±45◦ layers to appear (i.e. Ny/Nx =2.0 and ∆t/t=0.05). In
such a situation, the already present ±45◦ plies are “not enough” and two additional ones for
each side of the symmetry plane are added (to preserve the balance and symmetry constraints
the optimizer adds not one but two extra ±45◦ plies).
A similar phenomenon is observed for plates having only 90◦ plies in nominal design

(Ny/Nx = 2.8). The increase in thickness uncertainty makes 45
◦ plies to appear and next

to move apart from the symmetry plane. But for high enough load ratios (e.g. Ny/Nx = 7.0
which is not shown in Fig. 3) the robust solution stays the same as the nominal design, and all
the layers are 90◦ oriented for nominal design and all considered tolerance levels.
Figure 3 shows also that for mixed plies solutions the increase in Ny/Nx ratio makes the

existing ±45◦ plies to move towards the symmetry plane. This phenomenon was previously
observed for the nominal designs, and now it is also confirmed for robust ones.
Changes in the buckling load parameter λcr listed in Table 2 for the examined Ny/Nx

ratios are presented in Fig. 4. The depicted relative value ∆λcr/λcrnom is calculated similarly to
Example 1. It is clear that the impact of ∆t parameter stays fully linear for uni-ply plates (e.g.
Ny/Nx = 7.0 with 90

◦ layers only) and quasi-linear for plates composed of mixed plies (e.g.
Ny/Nx =0.125 or Ny/Nx =2.5 both having ±45

◦ and 90◦ fibers).The observed quasi-linearity
is explained by the results of Example 1 and Fig. 2, where the slopes of lines for uni-ply plates
are different, but still close to each other. Therefore, formixed solutions observed in the bi-axial
case and bearing inmind the fact that the individual plies change their position for different ∆t
(see also Fig. 3) a noticeable quasi-linearity is observed.

Fig. 4. Relative change of the buckling load factor λcr for different Ny/Nx load ratios and tolerance
level ∆t/t. Calculated with respect to the nominal (ideal) solution

As observed in Fig. 3 most of the final solutions are of different stacking sequence while
compared to the nominal design (∆t/t=0), so an extra analysis to estimate the benefit of the
robust design seems to be justified. Therefore, the question is: what would be the buckling load
factor for the nominal solution plate subjected to thicknesses tolerances that are not introduced
into the optimization method. This value is denoted in Table 3 by λ′cr and calculated as the
buckling load for the nominal optimum solution stacking sequence diminished by possible tole-
rances and appropriate sensitivities. Next λ′cr is compared to the buckling load amplitude for
optimum solutions with tolerances introduced directly into the design algorithm according to
the approach presented in this paper (i.e. values given in Table 2). In Fig. 5, a relative change
(λcr−λ

′

cr)/λcr in the plate buckling load is presented. This corresponds to actual benefit, which
is possible to achieve if thickness uncertainties are introduced into the design code, compared to
plates designed only for nominal values. It is clear that the best results are achieved for plates



1050 J. Latalski

of mid Ny/Nx load ratios. For low or high ones (e.g. 0.125 or 7.0), there is no actual benefit of
the incorporating thicknesses tolerances in design code. This is because the nominal solutions
and robust solutions (the ones that consider tolerances into the design method) are the same
and contain only plies of one direction.

Table 3. Buckling load factor comparison for plates with (λcr) and without (λ
′

cr) thickness
uncertainties introduced into the design code

Ny/Nx
∆t/t [%]

1.00 2.00 3.00 5.00

0.125 λcr 146.231(2,1) 137.840(2,1) 129.000(1,1) 111.752(1,1)
λ′cr 146.231(2,1) 137.840(2,1) 129.000(1,1) 110.832(1,1)

0.50 λcr 89.201(1,1) 84.004(2,1) 78.782(1,1) 67.951(1,1)
λ′cr 89.174(2,1) 83.515(2,1) 77.857(2,1) 66.539(2,1)

1.00 λcr 58.564(2,1) 55.309(2,1) 51.701(2,1) 44.825(2,1)
λ′cr 58.438(2,1) 54.455(2,1) 50.472(2,1) 42.506(2,1)

2.00 λcr 34.896(1,1) 32.679(1,1) 30.734(1,1) 26.599(1,1)
λ′cr 34.422(2,1) 31.907(2,1) 29.392(2,1) 24.362(2,1)

3.50 λcr 21.436(1,1) 20.320(1,1) 19.105(1,1) 16.534(1,1)
λ′cr 21.436(1,1) 20.320(1,1) 18.911(2,1) 15.603(2,1)

7.00 λcr 11.088(1,1) 10.510(1,1) 9.933(1,1) 8.777(2,1)
λ′cr 11.088(1,1) 10.510(1,1) 9.933(1,1) 8.777(1,1)

Fig. 5. Relative buckling load change for plates without uncertainties introduced into the design
algorithmwith respect to robust ones

6. Conclusions

The problem of optimum design of composite laminates considering ply thicknesses tolerances
has been addressed by the deterministic approach. Updated formulas for the system perfor-
mance and cost function problem are derived based on the classical laminated plate theory as-
sumptions. Advantages of the suggested approach with respect to other deterministic methods
discussed in “Introduction” section are underlined. The presented idea exploits the assumed to-
lerances directly into the design optimization code in a relatively simpleway, thus computational
competitiveness of the method is provided.

The presented theoretical research and further numerical analysis results allow the following
specific conclusions to be formulated.



Ply thickness tolerances in stacking sequence optimization ... 1051

• The algorithm presented in this study offers an efficient and safe approach to incorporate
manufacturing tolerances into optimumdesign. Themethod accounts for realistic parame-
ter variations, so the design is practical rather than amathematical abstraction that is of
quite limited use in real world.

• Multiple local minima of cost function in composite stacking sequence optimization are
confirmed. This applies not only to nominal designs, but also to robust ones. The number
of equivalent optima is strictly related to the plate aspect ratio a/b and load ratio Ny/Nx
as well.

• Considering ply thicknesses perturbations in the laminate plate optimumdesign results in
a different ply stacking sequence when compared to the solutions of the nominal design
problem. This fully justifies the incorporation of thicknesses perturbations in the optimi-
zation design algorithm. It is shown that this approach is necessary for mid-Ny/Nx load
ratios, but has no practical importance for very low or very high load ratios since the
nominal and robust designs in these specific cases are identical.

Acknowledgement

The research leading to these results has received funding from the European Union Seventh Fra-

mework Programme (FP7/2007-2013), FP7-REGPOT-2009-1, under grant agreement No 245479. The

support by Polish Ministry of Science and Higher Education – Grant No. 1471-1/7.PRUE/2010/7 – is

also acknowledged.

References

1. Adali A., Duffy K., 1990, Design of antisymetric hybrid laminates formaximumbuckling load:
I. Optimal fibre orientation,Composite Structures, 14, 1, 49-60

2. Chryssanthopoulos M.K., Poggi C., 1995, Probabilistic imperfection sensitivity analysis of
axially compressed composite cylinders,Engineering Structures, 17, 6, 398-406

3. Fong M., Liu Y.P., Chan S.L., 2010, Second-order analysis and design of imperfect composite
beam columns,Engineering Structures, 32, 6, 1681-1690

4. Gutkowski W., Bauer J., 1999,Manufacturing tolerance incorporated in minimum weight de-
sign of trusses,Engineering Optimization, 31, 3, 393-403

5. GutkowskiW., Latalski J., 2003,Manufacturing tolerances of fiber orientations in optimisation
of laminated plates,Engineering Optimization, 35, 2, 201-213

6. HaftkaR.T.,GürdalZ., LamatM.P, 1990,Elements of StructuralOptimization, 2ndEdition,
Kluwer, Dordrecht

7. HaftkaR.T.,WalshJ.L., 1992,Stacking-sequenceoptimization for buckling of laminatedplates
by integer programming,AIAA Journal, 30, 2, 814-819

8. JiangC., HanX., LiuG.P., 2008,Uncertain optimization of composite laminated plates using a
nonlinear interval number programmingmethod,Computers and Structures, 86, 17/18, 1696-1703

9. Kogiso N., Shao S., Murotsu Y., 1997, Reliability-based optimum design of a symmetric
laminated plate subject to buckling, Structural Optimization, 14, 2/3, 184-192

10. Kristinsdottir B., Zabinsky Z., Tuttle M., Csendes T., Incorporating manufacturing to-
lerances in near-optimal design of composite structures,Engineering Optimization, 26, 1, 1-23

11. Lee K.-H., Park G.-J., 2001, Robust optimization considering tolerances of design variables,
Computers and Structures, 79, 1, 77-86

12. NemethM.P., 1986, Importance of anisotropyonbuckling of compression-loaded symmetric com-
posite plates,AIAA Journal, 24, 11, 1831-1835



1052 J. Latalski

13. SinghB.N.,YadavD., IyengarN.G.R., 2001, Stability analysis of laminated cylindrical panels
with uncertain material properties,Composite Structures, 54, 1, 17-26

14. Walker M., Hamilton R., 2005, A methodology for optimally designing fibre-reinforced la-
minated structures with design variable tolerances for maximum buckling strength, Thin-Walled
Structures, 43, 1, 161-174

15. Walker M., Hamilton R., 2006, A technique for optimally designing fibre-reinforced laminated
plates under in-plane loads for minimum weight with manufacturing uncertainties accounted for,
Engineering with Computers, 21, 4, 282-288

Manuscript received September 24, 2012; accepted for print June 13, 2013