See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/41714529

Recent progress in the optimal design of composite structures: industrial

solution procedures on case studies

Article in International Journal for Simulation and Multidisciplinary Design Optimization · December 2008
DOI: 10.1051/ijsmdo/2008038

CITATIONS READS

5 1,156

6 authors, including:

Michaël Bruyneel Caroline Raick

GDTech - Global Design Technology Siemens
106 PUBLICATIONS 1,364 CITATIONS 12 PUBLICATIONS 281 CITATIONS

SEE PROFILE SEE PROFILE

Alain Remouchamps Stéphane Grihon

Siemens Airbus
22 PUBLICATIONS 263 CITATIONS 89 PUBLICATIONS 758 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Michaël Bruyneel on 18 January 2014.

The user has requested enhancement of the downloaded file.

Int. J. Simul. Multidisci. Des. Optim. X, xxx-xxx (2008) Available online at:
© ASMDO 2008 http://www.ijsmdo.org
DOI: xx.xxxx/ijsmdo:xxxxxxx

Recent progress in the optimal design of composite structures: in-

dustrial solution procedures on case studies

Michael Bruyneel1,a, Benoît Colson1, Philippe Jetteur1, Caroline Raick1, Alain Remouchamps1, Stéphane Grihon2
1
SAMTECH s.a., Liège Science Park, 8 Rue des chasseurs-ardennais, 4031, Angleur, Belgium
2
Airbus France, 316 Route de Bayonne, 31060, Toulouse, France

Received xx xx 20xx, Accepted xx xx 20xx

Abstract – In this paper, recent developments carried out in the SAMCEF finite element code and in the BOSS quattro
optimization toolbox are presented. Those developments aim at simulating high non linear effects in laminated composite
structures (post-buckling, collapse, delamination) and at optimising the composites with respect to those structural re-
sponses. The use of Sequential Convex Programming and of Surrogate-Based Optimization methods is discussed on indus-
trial optimization problems.
Key words: optimization methods, composite structures, industrial application, delamination

problem including around 1000 design variables and about

1 Introduction 100000 design restrictions [8,11].

In order to decrease the lead time, virtual testing and numeri-

2 Advanced structural analyses
cal optimization techniques are increasingly used in the dif-
ferent fields of engineering. This is especially the case for This section, briefly describes the analyses involved in the
composite structures, where two main issues have recently optimization problems. The SAMCEF Finite Element code
been identified. Firstly, reliable methods are required to ana- [4] is used to compute the structural responses.
lyse post-buckling [1], collapse and delamination, which is a
specific failure mode of laminated structures [2]. Secondly, 2.1 Stability analysis
designing real composite structures such as aircraft wings
implies handling a large number of design variables and re- Structural stability analysis
strictions, resulting in a large scale optimization problem that Stability analysis consists in solving an eigen-value prob-
is difficult to both define and solve [3]. Moreover, local op- lem of the form (1), where K is the structural stiffness ma-
timization of aircraft components should involve advanced trix, L the geometric stiffness matrix (also termed the initial
non linear analyses simulating delamination, post-buckling stress stiffness matrix), Φ j the jth buckling mode and λj the
and collapse. associated buckling load. The components of vector Φ j are
This paper presents an overview and a synthesis of the so- the structure’s degrees of freedom, usually the displace-
lution procedures recently developed in the SAMCEF finite ments (translations and rotations). The buckling load must
element code [4] and in the BOSS Quattro multidisciplinary be interpreted as the factor by which the external loads must
optimization platform [5] for solving such problems at an be multiplied for the structure to become unstable.
industrial level. The advanced buckling and collapse analysis
capabilities are first presented, together with a discussion on KΦ j − λ j LΦ j = 0 , j = 1,2,... (1)
the related sensitivity analyses, and an approach to evaluate
delamination risks [6]. Then two optimization techniques Typically, Reserve Factors are used in the optimization
used in this paper are recalled: a specific Sequential Convex problem. For example RFj =λj /1.2, where 1.2 is a safety
Programming algorithm [7] and a Surrogate-Based Optimiza- factor.
tion method [8]. Finally, three applications are presented. Sensitivity analysis of the buckling loads
The first one tackles the optimization of laminates with re-
The analytical first order derivative of the buckling load λj
spect to damage tolerance considerations, based on a detailed
is derived in [12] and its implementation in an industrial
(high fidelity) analysis of delamination. The second one aims
software is discussed in [9]. It is given by:
at designing composite structures to withstand buckling and
collapse [9,10]. The last application presents the results ob-
tained for the optimal pre-design of a wing box of the Airbus ∂λ j  ∂K ∂L 
= Φ Tj  −λj Φ j , j = 1,2,... (2)
A350, in which local criteria are used in an optimization ∂xi  ∂xi ∂xi 

a
Corresponding author: michael.bruyneel@samcef.com
 Bruyneel et al.: Recent progress in the optimal design of composite structures: industrial solution procedures on case study.

As shown in this paper and in [9], in order to avoid Sensitivity analysis of the collapse load
excessive oscillations during the optimization process, a The goal here is to compute the value dλ/dxi, where xi is
large number of load factors (and associated modes) must be the ith design variable and λ is the load factor. Starting from
taken into account, to ensure that the structure is sensitive to the equilibrium equations and knowing that the forces de-
local modes everywhere. The sensitivity of RFj is easily ob- pend on the displacement q, the load increment λ and the set
tained from (2). of design variables x, it follows that

2.2 Post-buckling and collapse analyses ∂F ∂F ∂F

dq + dx + dλ = 0 (5)
∂q ∂x ∂λ
Structural non linear analysis and continuation method
Although the stability analysis provides an estimation of The unknown vector is also constrained to be orthogonal
the bifurcation points, it is based on a linear approach, and is to the load-displacement curve rather than being a simple
therefore only an approximation. Moreover, once buckling measure based on the vertical gap ∆λ. This allows a better
has occurred, some thin walled structures can sometimes accuracy of the sensitivity measure. It turns out that dλ/dx
still sustain an increasing loading up to the final collapse [1]. can be computed from the tangent stiffness matrix; the varia-
This is illustrated in Figure 1 in the case of a stiffened com- tion of the internal forces with respect to the design vari-
posite panel. ables is computed by finite differences. Details are given in
[10,14].

2.3 Advanced damage tolerance analysis

When the fracture mechanics approach is applied to the

delamination of laminated composites [2,6], the finite ele-
ment method is used to model the cracked structure and to
compute the energy release rates GI, GII and GIII, related to
each of the three modes (I, II and III) of inter-laminar frac-
ture (crack lips opening, sliding shear and tearing, respec-
tively). Once these values are computed at each node defin-
ing the crack front, they are inserted in a criterion such as (5)
and compared to GIC, GIIC and GIIIC, which characterize the
inter-laminar fracture toughness in the three individual
modes.
GI G G
+ II + III < 1 (5)
Fig. 1. Illustration of the post-buckling analysis of a stiff- G IC G IIC G IIIC
ened composite panel
When the value in (5) reaches unity, it is assumed that the
crack will propagate in the structure. In this paper a specific
Virtual Crack Extension method is used to compute GI, GII
To simulate the large displacements appearing in a struc-
and GIII. This computation is carried out in a linear analysis.
ture beyond its bifurcation point, a geometric non linear
The method is decomposed into two steps. In the first, the
analysis is required, in which, contrary to (1), the stiffness
variation of the total potential energy π at equilibrium is
matrix changes all over the loading. A non linear system of
evaluated:
equilibrium equations (3) is therefore solved, for an imposed
load factor λ or displacements q. When a limit point (e.g. dπ
collapse) must be identified, the classical Newton-Raphson GT = − = G I + G II + G III (6)
dA
method can present problems, and a continuation method
(also called arc length or Riks method) must be used. In this Considering the expression of π in linear elasticity, and
method, q and λ are unknown and linked to an additional taking the derivative with respect to the crack length (sur-
parameter, called the arc length [13], with equation (4). This face) A, it turns out that:
parameter is controlled over the iterations. The complete set 1
of equations can now be now written as: π = q T Kq − g T q
2
F(q, λ ) = 0 (3) dπ 1 T dK dg
=> = q q − qT
c (q , λ ) = 0 (4) dA 2 dA dA
The Reserve Factor associated with the collapse load is where q are the nodal displacements and g is the vector
given by the load increment λ, i.e. RF = λcollapse. At the solu- of applied nodal forces. In a semi-analytical approach,
tion, RF must be larger than or equal to 1, implying that col- dK/dA can be replaced by ∆K/∆A. The total energy release
lapse will not occur at a loading lower than the nominal one. rate GT in (6) is finally obtained from (7), where the index
init corresponds to the initial (actual) crack length, and pert
denotes a perturbed configuration of the crack associated
with a virtual (very small) advance in its local plane. Equa-
tions (6) and (7) clearly show how the sensitivity analysis
 International Journal for Simulation and Multidisciplinary Design Optimization

for structural optimization and the calculation of the energy

release rate in a fracture mechanics problem are related.

dπ (π pert − π init )
GT = − ≅− (7)
dA ∆A
In the second step, the total energy release rate given by
(7) is distributed over the three modes GI, GII and GIII, con-
sidering the relative displacements of the lips and the reac-
tions to the crack opening, sliding and tearing at the crack
front. The VCE method is used in SAMCEF to identify the
most critical cracks in a structure, and to estimate the propa-
gation load, i.e. the amplitude of the load leading to at least Fig. 3. The bottom surface of a wing and its associated su-
one crack propagation, according to criteria such as (5). This per-stiffeners where local computations are carried out
method is illustrated in Figure 2, for mode I.

Sensitivity analysis
The interaction between the local and the global effects is
taken into account by the sensitivity analysis, as illustrated
in Figure 3, through the modification of the internal forces N
with respect to a change in each design variable x. While the
computation of global sensitivities are based on semi-
analytical derivatives, finite differences and massive paral-
lelism are used for computing the derivatives of the local
values, on the set of super-stiffeners. For details, see [8,11].

3 Basic description of the optimization

methods
Several optimization methods are used in this paper, depend-
ing on the application. The selection of an efficient optimi-
zation method depends on the problem to be solved. When
derivatives are available and when the problem includes a
Fig. 2. Principle of the VCE method illustrated on a large number of design variables and design functions, it is
Double Cantilever Beam (DCB) recommended to use a gradient-based method, such as a
Sequential Convex Programming method (SCP). On the
other hand, when derivatives are not available, either be-
A semi-analytical sensitivity analysis of the energy re- cause the problem is non differentiable, or simply because
lease rate with respect to ply thickness and fibres orientation they have not be implemented, the use of a Genetic Algo-
has not yet been developed in SAMCEF, but is somehow rithm (GA) is certainly a good choice, assuming that the
related to the structural stiffness. Moreover, as explained in computational time for a structural analysis is very low and
Section 4, some uncertainties still exist concerning the for- that the problem includes few design variables. When non
mulation of the optimization problem. linear structural analyses are conducted, a Surrogate Based
Optimization method (SBO) is preferred to a Genetic Algo-
2.4 Pre-design criteria for composite aircraft rithm. In this case the number of design variables must how-
structures ever remain relatively small. To conclude therefore, it is
important to note that a gradient based optimization method
Structural analysis will not generally provide the global optimum of the prob-
A wing is divided into a set of super-stiffeners, made of a lem, while a Genetic Algorithm will, for a large population.
portion of the skin and a T-stiffener (Figure 3). The Reserve Moreover, when working with discrete design variables, GA
Factors are computed at the local level, i.e. in each super- is a natural and easy choice. In the following applications,
stiffener. According to [3,11,15], several criteria can appear continuous design variables are considered, and SCP and
in the formulation of the pre-design of a composite aircraft SBO methods are used.
wing. In this case, Reserve Factors reflect buckling, damage
tolerance, reparability and some design rules at the super- 3.1 The general optimization problem
stiffener level. The values of these RFs vary with respect to
design variables, which are the panel thickness, the propor- The optimization problem to be solved takes the following
tions of fibres oriented at 0°, ±45° and 90°, and dimensions form:
of the stiffeners (web and flange height, web thickness, etc).
min FOB( x )
The RFs are also functions of internal forces, N, computed
from a global finite element model of the wing as depicted RF j ( x ) ≥ 1 j = 1,..., m (8)
in Figure 3. The RFs values at the local level (in the super-
stiffeners) are calculated using analytical formulas. xi ≤ xi ≤ xi i = 1,..., n
 Bruyneel et al.: Recent progress in the optimal design of composite structures: industrial solution procedures on case study.

where x is the set of design variables, FOB is the objec-

tive function (weight or total energy release rate in the fol-
lowing applications), and RFj is the jth Reserve Factor. Two
methods are used in what follows. For a detailed presenta-
tion of optimization methods, the reader should refer to
[7,8].

3.2 Sequential Convex Programming - SCP

In this gradient-based approach, the solution of the initial

optimization problem (8) is replaced by the solution of suc-
cessive explicit approximations, as illustrated in Figure 4 for Fig. 5. Illustration of the Surrogate Based Optimization
an unconstrained optimisation problem. A generalization of approach
the MMA approximation, developed in [7] and available in
the BOSS Quattro optimization toolbox, is used. This
method automatically generates monotonous or non mo- 4 case study 1: design with respect to de-
notonous approximations of the structural response with
respect to the design variables, in order to have an accurate lamination
model of the initial optimization problem. An approximation Designing a composite structure with respect to damage
of problem (8) is generated based on the functions value and tolerance is becoming a challenge. Although simplified for-
their first order derivatives (obtained from structural and mulations may be used at a global level for a pre-design (e.g.
sensitivity analyses, respectively). Once the solution of an a complete wing [11]), more sophisticated high fidelity ap-
approximated problem is obtained, it is checked to see proaches should be selected when local detailed structural
whether convergence has been reached. If this is not the components are studied. An ENF (End Notched Flexure)
case, a new approximation is built, and the process is con- specimen is considered here (Figure 6). The function to be
tinued until convergence to a desired accuracy is achieved. minimized is the total energy release rate GT (6) over the
The advantage of this method is that it requires a small crack front, in order to design a laminate less sensitive to
number of iterations to reach the solution (say between 10 delamination. Since the value of GT varies along the crack
and 25) irrespective of the number of design variables. front, minimum, mean and maximum values are therefore
However, derivatives must be available, and generally only minimized, in order to obtain a laminate less prone to crack
local optima can be identified. For constrained optimisation propagation. Only one design variable is selected: it is the
problem, a dual approach is used (see [7]). angle of the following lay-up: [±θ /0/- θ /0/ θ /04/ θ /0/- θ
/0/- θ / θ /d/- θ / θ /0/ θ /0/- θ /04/- θ /0/ θ /0/± θ]. The initial
value of θ is 30°.
The inter-laminar toughnesses GIC, GIIC and GIIIC used in
the criterion (5) depend on the fibre orientations across the
interface [16]. Moreover, a large dispersion in the inter-
laminar toughness values exists. In the simplified approach
adopted here, a constant value of GC is considered. For the
sake of accuracy this variation of GC should be included in
the formulation of the optimization problem, and robust op-
timization methods should be used [17]. The current ap-
Fig. 4. Illustration of the Sequential Convex Programming proach is therefore not rigorous and aims rather at comparing
approach optimization methods. The results must therefore be inter-
preted with care.

3.3 Surrogate Based optimization - SBO

In this method [8], only the function values are used to build
a global approximation of the design domain, as illustrated
in Figure 5. A response surface is generated based on a Neu-
ral Network (NN) or a Radial Basis Function (RBF). The
optimum of this global approximation is obtained with a
genetic algorithm. The new point is then used to build a new
response surface, and the process is repeated until conver- Fig. 6. The ENF test case, and its finite element model
gence. The advantage of this approach lies in the fact that it
can be used when the derivatives are not available. On the
other hand, it results in a prohibitively long computational First the SCP method presented in [7] is used, and the
time when a large number of design variables is used. How- sensitivities are up-to-now computed by finite differences.
ever, parallelism allows a reduction – to some extent – of the The results of the iterative optimization process are
time needed to find a solution. For constrained optimisation presented in Figure 7a. The optimal value of θ (0-deg) is
problems, a penalty is used. obtained after 8 iterations, and 15 structural analyses.
 International Journal for Simulation and Multidisciplinary Design Optimization

fined for the corresponding values in the stiffener. Since

seven panels are used to form the structure, 42 design vari-
ables are defined in the problem.

5.1 Optimization with respect to buckling

The goal is to find the values of the ply thickness that mini-
mize the weight while satisfying RF ≥ 1, according to (1). It
can be seen in Figure 9 that when only the first 12 load fac-
tors are considered in the design problem, large oscillations
appear. Given the local character of the associated buckling
modes, some structural parts may become insensitive to
buckling, and the modes start to move all over the structure
during the iterative optimization process, leading to oscilla-
tions. In contrast, when a large number of buckling loads are
taken into account, the whole structure remains sensitive to
buckling and, insofar as a reliable optimization method is
used, the solution can be obtained in a small number of
iterations (Figure 9). More details can be found in [9].

Fig. 7. The ENF test case, on its finite element model

Surrogate Based Optimization methods [8] – where Neu-

ral Networks or other response surfaces are used to generate
a global approximation of the design problem – may be effi-
cient for treating the damage tolerance problem, as was the
case in [1] for post-buckling optimization of composite
structures. In this case, their main advantage is that analyses
are run on parallel processors, which is not possible with a
sequential gradient-based approach such as SCP. Applying
such a method to the problem depicted in Figure 6 provides
the solution after 63 structural analyses, which is competi-
tive with the SCP method since 4 processors were used for
the parallel solution (Figure 7b). The formulation of this
kind of problems should be further investigated.

5 case study 2: design of a stiffened panel

Fig. 9. Convergence history for 12 and 100 buckling modes
A hat-stiffened curved composite panel subjected to shear
and compression is considered (Figure 8).
5.2 Optimization with respect to buckling and col-
lapse

This problem can be solved when post-buckling behaviour is

considered as well. In this case, post-buckling is allowed to
develop in the structure, and the buckling loads are pre-
scribed to be larger than 0.76. The collapse must appear for
λcollapse ≥ 1. The results are presented in Figure 10, where 100
buckling modes are taken into account. The solution is ob-
Fig. 8. Curved composite panel and its super-stiffeners
tained after 17 iterations. As in the previous section, conver-
gence problems (oscillations) are observed when a small
In each super-stiffener (Figure 8b), three design variables number of buckling modes is used.
are related to the thickness of the plies oriented at 0°, ±45°
and 90° to the skin, and three other design variables are de-
 Bruyneel et al.: Recent progress in the optimal design of composite structures: industrial solution procedures on case study.

2. R. Krueger, Computational fracture mechanics for com-

posites – State of the art and challenges, In NAFEMS
Seminar – Prediction and Modelling of Failure Using
FEA (2006).
3. F. Huttner, M. Grosspietsch, Solving very large scale
structural optimization problems, AIAA J 45(11),
2729-2736 (2007).
4 SAMCEF. Système d’Analyse des Milieux Continus par
Eléments Finis. www.samcef.com
5. Y. Radovcic, A. Remouchamps, BOSS quattro: an open
system for parametric design, Struct. Multidiscipl. Op-
Fig. 10. Convergence history for buckling and collapse opti- tim. 23, 140-152 (2002).
mization 6. M. Bruyneel, J.P. Delsemme, F. Germain, Ph. Jetteur, A
study of complex delaminations in laminated composite
6 case study 3: large scale optimization of structures with SAMCEF, In Fourth International Con-
ference on Advanced Computational Methods in Engi-
a composite aircraft component neering (2008).
Only a part of the composite box structure of a wing is con- 7. M. Bruyneel A general and effective approach for the
sidered in this optimization problem. Based on Section 2.4, optimal design of fibers reinforced composite structures,
the mass is to be minimized, while 71552 restrictions on Compos. Sci. Technol. 66, 1303-1314 (2006).
buckling, damage tolerance and reparability are defined on 8. A. Remouchamps, S. Grihon, B. Colson, C. Raick, M.
the several super-stiffeners composing the wing box. The Bruyneel, Numerical optimization: a design space odys-
problem includes 630 design variables (panel thickness, and sey, In 1st International Conference on Advancements in
size of the stiffeners). As illustrated in Figure 11, the solu- Design Optimization of Materials, Structures and Me-
tion is found after 10 iterations. More details on a full wing chanical Systems (2007).
optimization are available in [9,11]. 9. M. Bruyneel, B. Colson, A. Remouchamps. Discussion
on some convergence problems in buckling optimiza-
tion, Struct. Multidiscipl. Optim. 35, 181-186 (2008).
10. B. Colson, M. Bruyneel, Ph. Jetteur, P. Morelle, A. Re-
mouchamps, Composite panel optimization with non
linear finite element analysis and semi-analytical sensi-
tivities, In NAFEMS Seminar – Simulating Composite
Materials and Structures (2007).
11. L. Krog., M. Bruyneel, A. Remouchamps, C. Fleury,
COMBOX: a distributed computing process for opti-
mum pre-sizing of composite aircraft box structures, In
Fig. 11. Convergence history for the optimization of the 10th SAMTECH Users Conference (2007).
composite wing 12. H. Adelman, R.T. Haftka, A discourse on sensitivity
analysis for discretely-modelled structures, NASA
7 Conclusions Technical Memorandum 104065.
13. E. Riks, C. Rankin, F. Brogan, On the solution of mode
The use of advanced structural analyses and optimization
jumping phenomena in thin walled shell structures,
methods has become challenging in the design of composite
Comp. Meth. Appl. Mech. Engng. 136, 59-92 (1996).
structures. This paper presented three main kinds of recent
14. M. Bruyneel, B. Coslon, J.P. Delsemme., Ph. Jetteur, A.
applications. Depending on the problem, either a gradient
Remouchamps, S. Grihon, Exploiting semi-analytical
based or a surrogate-based optimization method can be used.
sensitivities from linear and non-linear finite element
This comparison should be further investigated, especially
analyses for composite panel optimization, to appear in
when local non linear analyses, such as delamination, post-
Int. J. Struct. Stability & Dynamics (2009).
buckling and collapse, are considered. Indeed, in such cases,
15. D. Weiss, Optimization in aircraft pre-design with siz-
few design variables are used, and both methods could be-
ing criteria, Diploma Thesis N° IMES-ST 06-019,
come competitive.
Swiss Federal Institute of Technology Zurich (2006).
16. O. Allix, D. Lévêque, L. Perret, Identification and fore-
Acknowledgement cast of delamination in composite laminates by an inter-
The support of the European Commission (VIVACE Project) laminar interface model, Compos. Sci. Technol., 58,
and of the Walloon Region of Belgium (OPTISTACK Pro- 671-678 (1998).
ject) is gratefully acknowledged. 17. H.G. Beyer, B. Sendhoff, Robust Optimization – A
Comprehensive Survey, Comput0 Meth. Appl. Mech.
Engng. 196, 3190-3218 (2007).
References
1. L. Lanzi, V. Giavotto, Post-buckling optimization of
composite stiffened panels: computations and experi-
ments, Compos. Struct. 73, 208-230 (2006).

View publication stats