Composites Part B 176 (2019) 107193

Contents lists available at ScienceDirect

Composites Part B
journal homepage: www.elsevier.com/locate/composites

Multi-scale identification of the elastic properties variability for composite

materials through a hybrid optimisation strategy
Lorenzo Cappelli a , Georgios Balokas b , Marco Montemurro a ,∗, Frédéric Dau a ,
Laurent Guillaumat c
a
Arts et Métiers ParisTech, Institut de Mécanique et d’Ingénierie (I2M) de Bordeaux CNRS UMR 5295, F-33400 Talence, France
b Structural Optimisation for Lightweight Design, Hamburg University of Technology, Hamburg, Germany
c Arts et Métiers ParisTech, Laboratoire Angevin de Mécanique, Procédés et innovAtion (LAMPA), F-49100 Angers, France

ARTICLE INFO ABSTRACT

Keywords: The problem of the identification of the variability characterising the elastic properties of the constitutive
Composite materials phases of a composite (at the microscopic scale) is addressed in this work. To this purpose, the information
Homogenisation contained into the probability distribution of the first buckling load of a macroscopic composite specimen is
Buckling
considered, in order to develop a multi-scale identification strategy (MSIS).
Uncertainty quantification
The goal of the proposed MSIS is achieved by solving an inverse problem: the minimisation of the distance
Surrogate model
Inverse problems
between the numerical and the reference buckling response of the plate, at the macroscopic scale, in terms of
Optimisation statistical moments. Furthermore, thermodynamic constraints are considered to ensure the positive definiteness
of the stiffness tensor of each constituent of the composite.
The proposed strategy relies on: (a) a semi-analytical homogenisation method, to perform the microscopic
/ mesoscopic scale transition; (b) the Monte-Carlo technique and an Artificial Neural Network to determine
the material properties variability; (c) a general hybrid optimisation algorithm able to deal with optimisation
problems defined over a domain of variable dimension to perform the solution search. The effectiveness of the
MSIS is proven through two meaningful benchmarks.

1. Introduction resin between adjacent laminae and incomplete cure of resin are only
some examples: environmental factors and uncertain operational aggra-
Composite materials are nowadays widely used into mechanical vate this issue.
components or engineering systems and structures belonging to dif- As outlined in [5], the uncertainties are classified in three main
ferent fields: from automotive to aerospace, from naval to biomedical. categories: aleatory (variability of structural parameters), epistemic
They are mainly employed due to their outstanding strength-to-weight (lack of adequate information about the system) and prejudicial (ab-
and stiffness-to-weight ratios: these features are of paramount impor- sence of stochastic characterisation of the structural system). Composite
tance for lightweight applications, such as aircraft and space vehicles structures are affected by all three forms of uncertainty and the char-
architectures [1]. Composites can be used to build integrated structures acterisation of parameters tuning the variability law becomes of prime
because both stiffness and strength can be tailored point-wise according importance. However, experimental methods commonly used to char-
to the requirements of the problem at hand. This feature allows for
acterise the material properties require a huge number of standard
preserving structural continuity without introducing complex struc-
ASTM tests, if used for uncertainty characterisation, which are de-
tural elements (and the related manufacturing aspects) by opportunely
structive and expensive [6]. Moreover, these tests are only suited to
meeting geometrical and mechanical design requirements.
evaluate mesoscopic uncertainties, in terms of material and geometrical
In the literature, research studies exploiting refined numerical and
properties of the lamina, without providing any information about the
experimental techniques are increasingly used to characterise the me-
variability characterising the properties of the constitutive phases at the
chanical behaviour of composite materials [2–4]. Nevertheless, espe-
cially in large-scale production, a large amount of uncertainty arises microscopic scale.
from unavoidable manufacturing imperfections for both geometrical Standard tests that can be carried out at the mesoscopic scale
and material properties. Intralaminar and/or matrix voids, excess of include the tension test for flat specimens (ASTM D3039 [7]), the

∗ Corresponding author.
E-mail addresses: [email protected], [email protected] (M. Montemurro).

https://doi.org/10.1016/j.compositesb.2019.107193
Received 15 May 2019; Received in revised form 21 June 2019; Accepted 10 July 2019
Available online 12 July 2019
1359-8368/© 2019 Elsevier Ltd. All rights reserved.
L. Cappelli et al. Composites Part B 176 (2019) 107193

three/four points bending test (ASTM D790 [8]), the compression tests (RBO), robust design optimisation (RDO) and model updating. The RBO
(shear loading methods ASTM D3410 [9]) and the shear tests (in-plane technique concerns the solution of an optimisation problem in which
shear tests ASTM D5379 [10]-D7078 [11]-D3518 [12], out-of-plane - the main goal is to design for safety by considering extreme events:
interlaminar shear tests ASTM D2344 [13]-D5379). As far as the mi- common objective functions are defined by the structural weight and
croscopic scale is concerned, only few standard experimental tests can the constraints are both deterministic and probabilistic (e.g. probability
be found in the literature: single fibre tensile test (ASTM D3379 [14]) of failure of the structure) [43–46]. The RDO method is usually imple-
and matrix tensile test (ASTM D638 [15]) to characterise the Young’s mented in order to minimise the influence of stochastic variations on
moduli of the fibre in the longitudinal direction and that of the matrix, the mean design [47]. Finally, the typical goal of the model updating
respectively. In order to characterise the rest of the constitutive phases
technique is to reduce the differences between model prediction and
elastic properties, only non-standard tests are available in the literature:
data from tests [48,49]. In this context, the MSIS can be considered as a
pull-out [16], micro-indentation [17], fragmentation tests [18], etc.
model updating technique that allows identifying the elastic properties
These tests are not able to evaluate the 3D set of the material properties
of the composite (and the related uncertainty) at each scale. This
of the constituents and they are very difficult to be carried out, due to
the fibre diameter size. information can be later used in the framework of both RBO and RDO
In order to get statistically representative results, the aforemen- approaches.
tioned tests must be performed a huge number of times. Of course, The paper is organised as follows. The problem description and
this implies significant costs (and time) and the variability results the MSIS are presented in Section 2. The analytical and the finite
(e.g. mean and standard deviation of material properties) are strongly element (FE) models, developed at each pertinent scale, are shown in
affected by the errors introduced to carry out the experimental cam- Section 3. The uncertainty microscopic quantification with the Monte-
paign, especially for those tests conducted at the microscopic scale. To Carlo technique and the implemented ANN are described in Section 4.
this purpose, Sepahvand et al. developed the inverse stochastic method The sensitivity analyses concerning the meta-model of the considered
based on the general polynomial chaos (gPC) [19–25] to identify benchmarks are presented in Section 5, while the mathematical formu-
uncertain lamina elastic parameters from experimental modal data. lation of the inverse problem is discussed in Section 6. The numerical
Further examples of probabilistic methods are the parametric proba- results provided by the MSIS are given in Section 7. Finally, Section 8
bilistic approach [26] and the Bayesian inference techniques wherein ends the paper with conclusions and perspectives.
all information are included into a prior distribution model [27–30].
However, in the case of composite structures, the uncertainty affecting
the ply elastic behaviour is strictly related to the variability of the 2. Multi-scale identification of the variability of composite elastic
elastic properties of the constitutive phases. To the best of the authors’ properties
knowledge, only few works on the identification of the variability
parameters characterising the material properties of the microscopic
2.1. The multi-scale identification strategy
constituents of the composite are available in the literature [31]. The
majority of researches in this field is devoted to the characterisation of
the material properties uncertainty parameters at the ply-level [32–35]. The goal of the MSIS is the characterisation of the variability related
The research activity here presented focuses on the development to the elastic properties of the microscopic constituents of the compos-
of a multi-scale identification strategy (MSIS) which smartly exploits ite, by using only the information contained into the statistical sample
the data resulting from macroscopic buckling tests to characterise the of the first buckling load of the multilayer plate at the macroscopic
uncertainty of the constitutive phases elastic properties. The proposed scale. In this way, only cheap, standard tests have to be realised,
MSIS has been initially proposed in [36] to identify the elastic prop- with the main advantage of reducing the characterisation time, the
erties of the composite (at each relevant scale), starting from the associated costs and the necessity of specialised skills.
harmonic response of the multilayer composite plate at the macroscopic
The reference macroscopic response can be evaluated either by
scale. Here, the MSIS is extended to the multi-scale characterisation of
means of an extensive experimental campaign of buckling tests or
the variability related to the elastic properties at the microscopic scale
through a wide numerical campaign of tests on a reference configu-
of the composite.
ration of the multilayer plate. To prove the effectiveness of the MSIS,
The MSIS is characterised by some original features. Firstly, it relies
on a particular hybrid optimisation tool used to perform the solution this latter case has been considered in this work.
search, which is an in-house code made by the union of a special To this purpose, the problem of characterising the variability re-
genetic algorithm (GA), i.e. ERASMUS (EvolutionaRy Algorithm for lated to the elastic properties of the fibre and the matrix is stated as
optimiSation of ModUlar Systems) developed by Montemurro [37] a multi-scale constrained inverse problem. Of course, the numerical
(which is able to deal with problems characterised by a number of models involved in the MSIS are characterised by some fundamental
design variables that can change during the optimisation process [38]) hypotheses. As far as the microscopic scale is concerned, the main
and of a gradient-based one, belonging to the MATLAB® fmincon hypotheses are: (a) linear elastic isotropic behaviour for the matrix;
family [39]. Secondly, the MSIS makes use of the Chamis’ micro- (b) linear elastic transversely isotropic behaviour for the fibre; (c) the
mechanical model [40,41] to perform the microscopic/mesoscopic matrix/fibre interface is perfect; (d) the damping capability of both
scale transition. Finally, the MSIS makes use of the Monte Carlo phases is neglected; (e) the uncertainty of the elastic properties is
framework that allows describing the statistical nature of the elastic described by means of a Gaussian probability distribution.
response. To improve the efficiency of the Monte Carlo technique
At the laminate macroscopic scale the following hypotheses hold:
(i.e. to minimise the computational effort related to such a method),
(a) the constitutive lamina has a linear elastic transversely isotropic
an Artificial Neural Network (ANN) [41] is developed as a surrogate
behaviour; (b) the interface between two adjacent plies is perfect; (c)
model: the probability distribution of the first buckling load is predicted
starting from the probability density functions of the elastic properties the damping capability of the lamina is neglected; (d) the kinematics
of the constituent phases. The effectiveness of the MSIS is proven by of the laminate is described by the first-order shear deformation theory
means of two meaningful benchmarks. (FSDT).
Concerning the state-of-the-art of the approaches combining optimi- The general flow chart of the MSIS is illustrated in Fig. 1. The details
sation and uncertainty, three specific research areas can be identified of the optimisation algorithms employed within the MSIS are given
in the literature, as outlined in [42]: reliability-based optimisation in [36,37].

2
L. Cappelli et al. Composites Part B 176 (2019) 107193

Fig. 2. Geometrical parameters of the reference multilayer composite plate (dimensions

are in [mm]).

property 𝑥𝑖 . If 𝑥𝑖𝑗 is the j-th value of 𝑥𝑖 occurring with a probability 𝑝𝑖𝑗 ,

the relative mean value and the variance can be expressed as:
𝑁
( ) ∑𝑖
𝜇 𝑥𝑖 = 𝑥𝑖𝑗 𝑝𝑖𝑗 ,
𝑗=1
( ) ( ( )) (2)
𝜎 2 𝑥𝑖 = 𝜇 𝛽 𝑖 𝑥𝑖 ,
Fig. 1. Flow chart of the MSIS.
( ) ( ( ))2
𝛽 𝑖 𝑥𝑖 = 𝑥𝑖 − 𝜇 𝑥𝑖 .
( )
Usually, the coefficient of variation COV 𝑥𝑖 is introduced as a stan-
2.2. Problem description dard measure of the dispersion of the probability distribution function:

The proposed multi-scale inverse approach for uncertainty char- ( )

( ) 𝜎 𝑥𝑖
acterisation is here applied to a reference multilayer composite plate COV 𝑥𝑖 = ( ) . (3)
made of unidirectional laminae: the relevant geometrical parameters 𝜇 𝑥𝑖
are shown in Fig. 2. Two different benchmarks are investigated to eval- In this work, the reference distribution of the first buckling load of
uate the identification capability of the proposed MSIS. In particular the the structure is determined by means of a multi-scale numerical analy-
geometry of the reference laminate is the same for both cases, the only sis on the reference configuration of the plate for both benchmarks. In
difference being the considered stacking sequence, i.e. particular, the reference material properties of the constitutive phases,
listed in Table 1, are implemented, firstly, to compute the reference
• benchmark 1 (BK1) [0◦ ∕ − 45◦ ∕45◦ ∕90◦ ]S ; distribution of the ply elastic properties. Secondly, the resulting distri-
• benchmark 2 (BK2) [45]8 . bution of the lamina elastic properties is used to compute the reference
For both laminates, the thickness of the elementary lamina is 𝑡ply = distribution of the first buckling load of the composite plate, for each
0.282 mm. The orientation angle of the generic ply is positive according considered benchmark (as described in Section 7.1).
to counter-clockwise rotation around the 𝑧-axis: the 𝑥-axis indicates the
0◦ orientation, as illustrated in Fig. 2. 3. Analytical and numerical models at different scales
The constitutive ply is made of carbon-epoxy fibre-reinforced Hexcel
𝑇 650∕𝐹 584 pre-impregnated tapes: its elastic properties are listed in 3.1. Microscopic/mesoscopic scale transition: the Chamis’ model
Table 1. The mean values are taken from [50,51], while the standard
deviation and the relative shapes of the probability density functions Multi-scale modelling strategies are widely used to assess the be-
haviour of the composite at each relevant scale [52,53]. The transition
are not available experimentally for both the microscopic and the
from the scale of the constitutive phases (microscopic scale) to that
mesoscopic material properties. To this purpose, a Gaussian probability
( ) of the elementary ply (mesoscopic scale) is performed by means of a
density function 𝜒𝑖 = 𝜒𝑖 𝑥𝑖 is selected as a reference distribution for
homogenisation calculation. This phase can be performed either numer-
describing the uncertainty of the generic property 𝑥𝑖 at the scale of the
ically, e.g. by implementing the well-known strain energy homogeni-
constituent phases. The analytical formula of such a distribution is
sation technique of periodic media (SEHTPM) [54], or analytically by
( ( ))2
𝑥𝑖 − 𝜇 𝑥𝑖 using a suitable homogenisation scheme for composites, as the Chamis’
( ) ( ) model [40]. As discussed in [36], the SEHTPM has already been in-
1 2
𝜒𝑖 𝑥𝑖 = ( ) √ 𝑒 2𝜎 𝑥𝑖 , with 𝑥𝑖 ∈ ℜ. (1) tegrated into the MSIS to determine the equivalent elastic behaviour
𝜎 𝑥𝑖 2𝜋
of general periodic materials with complex microstructures. Despite its
In particular, the Gaussian distribution involves two parameters, general nature, the SEHTPM can be quite time consuming (depending
( ) ( )
i.e. the mean value 𝜇 𝑥𝑖 and the variance 𝜎 2 𝑥𝑖 of the i-th material on the problem at hand) since the equivalent elastic properties at the

3
L. Cappelli et al. Composites Part B 176 (2019) 107193

Table 1
Mean value and standard deviation of the elastic properties for the fibre 𝑇 650∕35 − 3𝐾 and the matrix 𝐹 584 (the mean values are taken
from [50,51]).
Fibre Matrix
𝐸1𝑓 [GPa] 𝐸2𝑓 [GPa] 𝑓
𝜈12 𝑓
𝜈23 𝑓
𝐺12 [GPa] 𝐸 𝑚 [GPa] 𝜈𝑚 𝑉𝐹
( )
𝜇 𝑥𝑖 276 17.3 0.25 0.428 11.24 4.14 0.35 0.555
( )
𝜎 𝑥𝑖 27.6 1.73 0.025 0.0428 1.124 0.414 0.035 0.0555

Fig. 3. Loads and boundary conditions (BCs) of the macroscopic FE model.

upper scale are the result of six static FE analyses (i.e. the equivalent properties of the ply are denoted as 𝐸1 , 𝐸2 , 𝐸3 , 𝐺12 , 𝐺13 , 𝐺23 , 𝜈12 , 𝜈13 ,
stiffness tensor of the homogenised material is evaluated column-wise). 𝜈23 .
When dealing with uncertainty quantification, the SEHTPM requires a
strong computational effort to evaluate the propagation of the uncer- 3.2. Mesoscopic/macroscopic scale transition: the finite element model
tainty from the microscopic scale to the mesoscopic one. Therefore, to
reduce the computational cost, an efficient analytical homogenisation The distribution of the first buckling load of the multilayer plate
scheme has been considered in this work, i.e. the aforementioned is the result of an eigenvalue buckling analysis which is carried out
Chamis’ model. Moreover, this choice allows avoiding the integration by considering the distribution of the ply elastic properties evaluated
by means of the Chamis’ model. The FE model is developed into the
of further FE model-related parameters like the mesh size. In particular,
Abaqus® environment [55]: the Abaqus® shell layered element S4R
according to the Chamis’ model, the ply engineering constants can be
having four nodes and six degrees of freedom (DOFs) per node has
determined as follows:
( ) been used to build the FE model of the multilayer plate. The kinematics
𝐸1 = 𝑉𝐹 𝐸1𝑓 + 1 − 𝑉𝐹 𝐸 𝑚 , of the element is described in the framework of the first-order shear
𝐸𝑚 deformation theory (FSDT) [1]. Of course, this type of element is well
𝐸2 = 𝐸3 = ( ),
√ suited to describe the buckling strength of the laminate when its aspect
𝐸𝑚
1 − 𝑉𝐹 1 − 𝑓 ratio is in the range [20, 100]. For the problem at hand the multilayer
𝐸2
plate is characterised by an aspect ratio 𝐴𝑅 = 44.29. Fig. 3 illustrates
𝐸𝑚 the loads and boundary conditions (BCs) for the proposed benchmarks.
2 (1 + 𝜈 𝑚 ) As far as the mesh size is concerned, a sensitivity study of the first
𝐺12 = 𝐺13 = ,
⎛ 𝐸𝑚 ⎞ buckling load of the laminate to the number of elements (not reported
√ ⎜ 2 (1 + 𝜈 ) ⎟
𝑚
1 − 𝑉𝐹 ⎜1 − ⎟ here for the sake of brevity) has been performed: a model with 3654
𝑓
⎜ 𝐺12 ⎟ DOFs is sufficient, to evaluate the first buckling load of the composite
⎝ ⎠ (4)
plate. The mesh of the FE model is illustrated in Fig. 4.
𝐸𝑚
2 (1 + 𝜈 𝑚 )
𝐺23 = ( ) , 4. Probabilistic modelling and uncertainty quantification
⎛ 𝑚 1 + 𝜈𝑓 ⎞
√ ⎜ 𝐸 23 ⎟
1 − 𝑉𝐹 ⎜1 − ⎟ 4.1. Monte Carlo Analysis
⎜ 𝐸2𝑓 (1 + 𝜈 𝑚 ) ⎟
⎝ ⎠
( ) The Monte Carlo (MC) method [41] is the most straightforward and
𝑓
𝜈12 = 𝜈13 = 𝜈 𝑚 + 𝑉𝐹 𝜈12 − 𝜈𝑚 , robust one, among the popular methods used for calculating the re-
𝐸2 sponse variability in stochastic structural mechanics. Based on the law
𝜈23 = − 1.
2𝐺23 of large numbers, MC approximates the statistical moments (e.g. mean,
variance, etc.) of the quantity of interest (QoI), by performing a suffi-
In Eq. (4), 𝐸1𝑓 , 𝐸2𝑓 , 𝐺12
𝑓 𝑓
, 𝜈12 𝑓
, 𝜈23 are the elastic constants of the cient number of model evaluations, while sampling random, indepen-
transversely isotropic fibre, while 𝐸 𝑚 and 𝜈 𝑚 are the Young’s modulus dent variables from the input space. The generated finite sample of
and the Poisson’s ratio of the isotropic matrix. The volume fraction the QoI is then post-processed, to obtain the unbiased statistics of the
of the fibre is indicated as 𝑉𝐹 . Moreover, the homogenised elastic response estimates. In mathematical terms, the first and second moment

4
L. Cappelli et al. Composites Part B 176 (2019) 107193

Fig. 4. Mesh of the macroscopic FE model.

described in Eq. (2) for the discrete case, can be approximated after 𝑁 Fig. 5. Architecture of a single layer feed-forward neural network.
realisations as:

1 ∑
𝑁
𝜇 (𝐫) = 𝑟, space are generated, namely 𝐴 and 𝐵, with 𝑁 being the number of
𝑁 𝑗=1 𝑗
(5) realisations and 𝑘 the stochastic dimension of the problem. After that,
1 ∑[
𝑁
]2 a third matrix 𝐴𝑖𝐵 is formed, identical to 𝐴, except its 𝑖th column which
𝜎 2 (𝐫) = 𝑟 − 𝜇(𝐫) .
𝑁 − 1 𝑗=1 𝑗 is replaced by the 𝑖th column of 𝐵 (𝑖 = 1, … , 𝑘). Finally, the model is
evaluated with respect to the aforementioned input matrices, according
where 𝐫 = {𝑟𝑖 , 𝑖 = 1, … , 𝑁} is the sample of the response QoI (e.g. dis- to the following estimator for the first-order Sobol index, for every
placement, force, bucking load etc.). Although MC can practically input parameter 𝑖:
handle every problem, regardless of the complexity of the response
1 ∑
𝑁
surface topology, the large number of required model evaluations
𝑉𝑖 = Var 𝑥𝑖 (𝐸𝑥∼𝑖 (𝑦|𝑥𝑖 )) ≈ 𝑓 (𝐵)𝑗 (𝑓 (𝐴𝑖𝐵 )𝑗 − 𝑓 (𝐴)𝑗 ). (11)
sets the method prohibitive for high-fidelity models (e.g. FE models), 𝑁 𝑗=1
especially for applications of reliability or uncertainty quantification.
It is noteworthy that there are several other options available,
regarding estimators of this sort [56]. A drawback of GSA, is that
4.2. Variance-based global sensitivity analysis
formulae like Eq. (11) require excessive realisations in order to con-
verge (order of 104 or 105 ). In the context of computationally expensive
In order to understand the cause-and-effect relationship between
simulations, such as FE analyses, a possible remedy is the emulation of
the input variables and the response, a classification of the random
the input–output relationship via a surrogate model, as it is described
parameters in terms of output variability can be achieved through
in the next section.
a global sensitivity analysis (GSA). The total variance of the QoI is
decomposed into parts induced from single input parameters, but also
4.3. Surrogate modelling with Artificial Neural Networks
potential interactions of the latter. Thus, the uncertain parameters can
be qualitatively quantified, and the dominating ones can be later used
Surrogate models (or metamodels) are mathematical functions able
into the models involved into the optimisation process introduced in
to mimic the response of a model, when trained with a relatively
Section 6.
small training set of model evaluations. Afterwards, the demanding
Let 𝑓 (𝑥1 , 𝑥2 , … 𝑥𝑘 ) be a square integrable scalar function over the
model can be substituted from these inexpensive proxy models, for
𝑘-dimensional unit hypercube 𝛺𝑘 model. According to Sobol [56], 𝑓
applications requiring an excessive amount of simulations (e.g. optimi-
can be decomposed into sums of increasing dimensions as follows:
∑ ∑ sation, reliability, GSA etc.) Popular choices, among relevant research
𝑓 = 𝑓0 + 𝑓𝑖 + 𝑓𝑖𝑗 + ⋯ + 𝑓12...𝑘 , (6) studies, are ANNs, Gaussian processes (or Kriging), polynomial chaos
𝑖 𝑗>𝑖 expansions (PCE) and support vector machines (SVM), as outlined
where 𝑓𝑖 = 𝑓𝑖 (𝑥𝑖 ), 𝑓𝑖𝑗 = 𝑓𝑖𝑗 (𝑥𝑖 , 𝑥𝑗 ) etc. After several algebraic manip- in [57].
ulations (the reader is referred to [56] or [41] for details), the final In this work, a surrogate is appropriately trained to emulate the
expression for the variance decomposition is reached: multi-scale modelling strategy described in Section 3. The material
properties of the different phases at the micro-scale, listed in Eq. (15),
∑
𝑘 ∑
𝑘
are used as input, while the output response is the plate first buckling
Var(𝑦) = 𝑉𝑖 + 𝑉𝑖𝑗 + ⋯ + 𝑉12...𝑘 , (7)
𝑖=1 𝑗>𝑖 load. The aim of the surrogate is twofold. Firstly, it is used for the GSA
and the evaluations required by the estimator of Eq. (11). Secondly,
where 𝑉𝑖 = Var 𝑥𝑖 (𝐸𝑥∼𝑖 (𝑦|𝑥𝑖 )), (8)
as explained in Section 1, it is used into the multi-scale identification
𝑉𝑖𝑗 = Var 𝑥𝑖𝑗 (𝐸𝑥∼𝑖𝑗 (𝑦|𝑥𝑖 , 𝑥𝑗 )) − 𝑉𝑖 − 𝑉𝑗 , etc. (9) strategy to boost the optimisation process. Concerning the surrogate
type, ANNs are selected in this study, mostly because, despite their
The 𝑥∼𝑖 notation indicates the set of all variables except 𝑥𝑖 . By dividing
versatility and their good generalisation, they only have few parameters
the term of interest by the unconditional variance Var(𝑦), the first-order
to be tuned within their training procedure, which is beneficial for the
Sobol index is obtained as a fractional contribution:
optimisation algorithm.
𝑉𝑖 In particular, an ANN is a parallel information-processing system,
𝑆𝑖 = . (10)
Var(𝑦) consisting of at least three layers: the input, the output and one (or
In the case of non-analytical models, expressions such as Eq. (8) more) hidden layer. The nodes inside every layer are called neurons
or (9) must be approximated via a sampling (e.g. Monte Carlo) proce- and they are linked by the so-called synapses. When information is
dure. Firstly, two (𝑁, 𝑘) matrices with random samples from the input circulated only in a single direction, the network is called feed-forward.

5
L. Cappelli et al. Composites Part B 176 (2019) 107193

An illustration of a typical single-layer, feed-forward ANN configu- Table 2

Benchmark BK1: bounds of the design variables.
ration is shown in Fig. 5. It is noteworthy that the input neurons
(squares) connect the network to the external environment, without Microscopic parameters Lower bound Upper bound
( )
further processing information, while the hidden layer neurons (circles) 𝜇 𝐸1𝑓 [GPa] 220.8 331.2
( )
process information from a previous layer and feed it to the next one. 𝜎 𝐸1𝑓 [GPa] 22.08 33.12
( )
The learning procedure of an ANN is based on a general uncon- 𝜇 𝑉𝐹 0.444 0.666
( )
strained optimisation problem, where the weight parameters 𝑤𝑖𝑗 as- 𝜎 𝑉𝐹 0.0444 0.0666
signed to every synapse are the design variables, and the objective
function is the sum squared error between the predicted output 𝑡(𝑤𝑖𝑗 )
and the target output 𝑦0 :
characterised by a quasi-isotropic symmetric stack, the sensitivity of the
1∑ first buckling load to the elastic properties 𝐸2𝑓 , 𝐺12
𝑓 𝑓
, 𝜈12 𝑓
, 𝜈23 , 𝐸 𝑚 , 𝜈 𝑚 is
𝐸(𝑤𝑖𝑗 ) = [𝑡(𝑤𝑖𝑗 ) − 𝑦0 ]2 . (12)
2 𝑓
negligible. Accordingly, only 𝐸1 and 𝑉𝐹 affects the laminate behaviour
During the process, the weights are updated through an iterative proce- in terms of first buckling load.
dure, until the desired error level is achieved or the maximum number Conversely, since the multilayer plate of benchmark BK2 is charac-
of cycles is reached: terised by an angle-ply orthotropic symmetric stacking sequence, the
𝑤(𝑡+1) = 𝑤(𝑡) (13) first buckling load is influenced by the following properties: 𝐸1𝑓 , 𝐸2𝑓 ,
𝑖𝑗 𝑖𝑗 + 𝛥𝑤𝑖𝑗 ,
𝐸 𝑚 and 𝑉𝐹 . The sensitivity of the laminate buckling strength to the
where 𝛥𝑤𝑖𝑗 is the correction of the weight at the 𝑡th learning step. other elastic properties remains negligible also for this configuration of
In order to avoid overfitting, a fraction of the sample data is used as the plate.
a validation dataset and the error is monitored over the iterations to According to the aforementioned remarks, the number of param-
stop the training early enough. Regarding the internal process in every eters (characterising the material and geometrical properties uncer-
neuron, each input from the previous neuron is placed into a weighted tainty) to be identified varies with the considered benchmark. These
sum as the following: aspects are discussed in detail in the following Section.
∑
𝑘
𝑧𝑗 = 𝑥𝑖 𝑤𝑖𝑗 + 𝑏, (14) 6. Mathematical formulation of the inverse problem
𝑖=1

which then goes through an activation function where the nonlinearity

6.1. Optimisation variables, objective function and constraints
of the decision boundary is introduced (usually of sigmoid type). The
term 𝑏 in the previous equation is a bias term allowing the neuron to
The multi-scale identification problem is stated as a classical con-
cover a broader range. For more details on ANNs, the interested reader
strained inverse problem: the identification of the elastic properties
is addressed to [41].
variability of the composite constitutive phases can be achieved by
minimising the Euclidean distance between the reference distribution
5. Sensitivity analysis of the meta-model
of the buckling load and that resulting from the numerical simulation.
5.1. Global sensitivity analyses for the two benchmarks As discussed in Section 5, the sensitivity of the buckling load
distribution to the material and geometrical parameters of the micro-
The implementation of the Artificial Neural Network, described in scopic constituents is strongly affected by the stacking sequence of
Section 4.3, allows to apply the variance-based GSA described in Sec- the laminate. Therefore, the number of optimisation variables (i.e. the
tion 4.2, since the computational effort needed to perform the conver- parameters of the distribution law, for each property at the microscopic
gence of the Sobol index is negligible. All the material and geometrical scale, to be identified) depends upon the considered benchmark. As
variables of the constitutive phases with the related uncertainty are a result of the GSA discussed in Section 5, the parameters tuning the
considered here: through the total output variance decomposition, it is distribution of the most relevant elastic and geometrical properties of
possible to identify the dominant microscopic input parameters, from the constitutive phases can be arranged in the vector of design variables
a statistical point of view. 𝝃 𝛼 , (𝛼 = BK1, BK2) as follows:
{ ( ) ( ) ( ) ( )}
According to the hypotheses given in Section 1, a total of eight
𝝃 BK1 = 𝜇 𝐸1𝑓 , 𝜎 𝐸1𝑓 , 𝜇 𝑉𝐹 , 𝜎 𝑉𝐹 , (16)
variables can be identified for the microscopic constituents of the
composite, i.e.
{ } { ( ) ( ) ( ) ( )
𝐱 = 𝐸1𝑓 , 𝐸2𝑓 , 𝐺12
𝑓 𝑓
, 𝜈12 𝑓
, 𝜈23 , 𝐸 𝑚 , 𝜈 𝑚 , 𝑉𝐹 . (15) ( ) ( )}
𝝃 BK2 = 𝜇 𝐸1𝑓 , 𝜎 𝐸1𝑓 , 𝜇 𝐸2𝑓 , 𝜎 𝐸2𝑓 , 𝜇 (𝐸 𝑚 ) , 𝜎 (𝐸 𝑚 ) , 𝜇 𝑉𝐹 , 𝜎 𝑉𝐹 .
The related mean and standard deviation values are summarised (17)
in Table 1, in which, a COV equal to 10% is set, for all the parame-
ters concerning the microscopic scale. The causes at the basis of this Accordingly, benchmark BK1 is characterised by four design variables,
uncertainty are various and often very difficult to be identified. For whilst benchmark BK2 has eight design variables. In both cases, the
example, the uncertainty of the fibre volume fraction is often related elastic properties excluded from the vector of design variables (due to
to the manufacturing process parameters. the negligible sensitivity of the first buckling load to these quantities)
The results of the variance-based GSA for every benchmark are have been set to the reference mean values listed in Table 1.
shown in Figs. 6(a) and 7(a), in terms of the evolution of the Sobol in- Each design variable can vary into a suitable definition domain
dex, defined in Eq. (11), over the number of simulations. It is possible to which depends upon the considered benchmark. Lower and upper
observe that the Sobol index converges after around 15000 simulations, bounds of design variables for benchmarks BK1 and BK2 are given in
for both benchmarks. Tables 2 and 3, respectively.
The pie diagrams shown in Figs. 6(b) and 7(b) highlight a result Moreover, in order to ensure the positive definiteness of the stiffness
of paramount importance: the sensitivity of the first buckling load to tensors of both the lamina (mesoscopic scale) and the constitutive
the material and geometrical properties of the constitutive phases (and phases (microscopic scale) [36], every combination of elastic properties
the related uncertainty as well) is strongly influenced by the nature generated through the Monte-Carlo technique, must satisfy a set of
of the stacking sequence. In particular, for benchmark BK1, which is non-linear constraints 𝐠 (𝝃 𝛼 ) [36]. Of course, these constraints must be

6
L. Cappelli et al. Composites Part B 176 (2019) 107193

Fig. 6. (a) Convergence of the Sobol index and (b) sensitivity analysis results, for the first benchmark (BK1).

Fig. 7. (a) Convergence of the Sobol index and (b) sensitivity analysis results, for the second benchmark (BK2).

Table 3 whilst for the constitutive phases they can be written as

Benchmark BK2: bounds of the design variables.
√
Microscopic parameters Lower bound Upper bound √ 𝑓 𝛼
√ 𝐸 (𝝃 )
( ) √
𝜇 𝐸1𝑓 [GPa] 220.8 331.2 𝑔4 (𝝃 ) = |𝜈12 | − √ 𝑓1
𝛼 𝑓
< 0,
( ) 𝐸2 (𝝃 𝛼 )𝑓
𝜎 𝐸1𝑓 [GPa] 22.08 33.12
( )
𝜇 𝐸2𝑓 [GPa] 13.84 20.76
𝑓
𝑔5 (𝝃 𝛼 ) = |𝜈23 | − 1 < 0, (19)
( )
𝜎 𝐸2𝑓 [GPa] 1.384 2.076 𝑓 𝛼 [ ( )2 ( )2 ]
𝐸 (𝝃 )
𝜇 (𝐸 𝑚 ) [GPa] 3.312 4.968 𝑔6 (𝝃 𝛼 ) = 1𝑓 𝑓
2𝜈23 𝑓
𝜈12 𝑓
+ 2 𝜈12 − 1 < 0.
𝜎 (𝐸 𝑚 ) [GPa] 0.3312 0.4968 𝐸2 (𝝃 𝛼 )
( )
𝜇 𝑉𝐹 0.444 0.666
( )
𝜎 𝑉𝐹 0.0444 0.0666
The objective function 𝛷(𝝃 𝛼 ) is defined as the Euclidean distance
between the reference and the numerical mechanical response, in terms
of the probabilistic parameters 𝜇𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 and 𝜎𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 of the first buck-
imposed at the lamina-level and at the constitutive phases-level. For ling load. In particular, this objective function is a least-square error
the elementary lamina, these constraints read: estimator defined as:
( )2 ( )2
𝜇𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 (𝝃 𝛼 ) − 𝜇𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔
ref 𝜎𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 (𝝃 𝛼 ) − 𝜎𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔
ref

√ 𝛷(𝝃 𝛼 ) = + . (20)
ref
𝜇𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 ref
𝜎𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔
𝐸1 (𝝃 𝛼 )
𝑔1 (𝝃 ) = |𝜈12 (𝝃 )| −
𝛼 𝛼
< 0,
𝐸2 (𝝃 𝛼 ) Finally, the multi-scale inverse problem is stated as a classical CNLPP
√
𝐸2 (𝝃 𝛼 ) as:
𝑔2 (𝝃 ) = |𝜈23 (𝝃 )| −
𝛼 𝛼
< 0,
𝐸3 (𝝃 𝛼 ) (18) min 𝛷 (𝝃 𝛼 ) ,
𝛼
𝝃
𝐸 (𝝃 𝛼 )
𝑔3 (𝝃 ) = 2𝜈12 (𝝃 )𝜈13 (𝝃 )𝜈23 (𝝃 ) 3 𝛼 + ⋯
𝛼 𝛼 𝛼 𝛼
(21)
𝐸1 (𝝃 ) subject to:
𝛼
𝛼
𝛼 2 𝐸2 (𝝃 ) 𝛼 2 𝐸3 (𝝃 ) 𝐸 (𝝃 𝛼 )
+𝜈12 (𝝃 ) + 𝜈23 (𝝃 ) + 𝜈13 (𝝃 𝛼 )2 3 𝛼 − 1 < 0, 𝑔𝑗 (𝝃 𝛼 ) ≤ 0, 𝑗 = 1, … , 6.
𝐸1 (𝝃 𝛼 ) 𝐸2 (𝝃 𝛼 ) 𝐸1 (𝝃 )

7
L. Cappelli et al. Composites Part B 176 (2019) 107193

Fig. 8. Optimisation strategy for the resolution of problem (21).

6.2. The numerical strategy Table 4

Variability parameters of the reference first buckling load.
( ) ( )
𝐼 𝐼
Benchmark 𝜇 𝜎𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 [MPa] COV 𝜎𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔

BK1 83.41 0.12

BK2 59.3 0.098
Problem (21) is a non-convex CNLPP, in terms of constraints and
objective function. The number of parameters, describing the variabil-
ity of material and geometrical properties of the constitutive phases,
depends on the considered benchmark: the first benchmark (BK1) 7. Numerical results
allows to characterise four parameters, while the second benchmark
(BK2) allows to characterise up to eight parameters. Of course, the
non-convexity of problem (21) implies the lack of uniqueness of its 7.1. Buckling response for the reference configuration
solution [36].

Taking into account all these aspects, the CNLPP of Eq. (21) is The multi-scale inverse problem defined in Eq. (21) requires the
solved by means of a hybrid optimisation tool based on the GA ERAS- computation of the objective function 𝛷(𝝃 𝛼 ) of Eq. (20): this function
MUS (EvolutionaRy Algorithm for optimiSation of ModUlar Systems), depends upon the buckling reference response which must be evaluated
which is interfaced with the MATLAB® fmincon algorithm [39], as before starting the optimisation process. Due to the difficulty to get
shown in Fig. 8. The GA ERASMUS has already been used successfully experimental data in terms of variability of the microscopic mate-
to solve different classes of real-world engineering problems [58–66]. rial properties and the related buckling probability distribution at the
macroscopic scale, a numerical test is performed in order to obtain the
The procedure illustrated in Fig. 8 is articulated in two phases. reference data.
The first one represents the global solution search and it is carried
To deal with this task, the reference variability parameters of mate-
out through the GA ERASMUS: the goal is to find potential suboptimal
rial and geometrical properties of the microscopic constituents, listed
solutions which will constitute the starting points for the gradient-based
in Table 1, are considered for each benchmark.
optimisation algorithm. The genotype of the individual is characterised
Firstly, a Monte-Carlo simulation is performed to generate randomly
by one chromosome and four genes for the first benchmark (BK1) and
a statistically representative number of samples. Secondly, for each
eight genes for the second benchmark (BK2).
sample, the homogenisation step is performed by using the Chamis’
The second step is the local optimisation phase and it is performed model, described in Section 3, to get the lamina elastic properties that
by means of the MATLAB® fmincon tool. The selected optimisation are used into the macroscopic FE model, to compute the first buckling
solver is the active-set algorithm, i.e. a Quasi-Newton method, in which load of the plate. After carrying out these operations for the whole
an approximation of the Hessian matrix is used to compute the descent set of samples, it is possible to determine the mean value and the
direction [39]. relative COV of the first buckling load, according to Eqs. (2) and (3),
respectively. The variability parameters of the reference first buckling
Each optimisation algorithm has been interfaced to the ANN, pre-
load distribution are then summarised in Table 4: these quantities have
sented in Section 4, which emulate both the homogenisation phase and
been obtained by performing 1000 realisations.
the eigenvalue buckling analysis. The ANN has been employed in order
Furthermore, a small sub-set of 20 realisations has been used to train
to reduce significantly the computational effort.
the ANN, for each benchmark. In order to check the accuracy of the
In particular, the output of the ANN is the current value of both ANN, 50 samples generated with the Monte-Carlo technique have been
the objective and the constraint functions which are passed to the selected as a validation set and a comparison between them and the
optimisation tool in order to execute the optimisation operations: these results provided by the ANN is carried out, as shown in Fig. 9 (only
operations are repeated until the user-defined convergence criteria are the results related to the first benchmark have been reported for the
met. sake of brevity). As a matter of fact, the results provided by the ANN

8
L. Cappelli et al. Composites Part B 176 (2019) 107193

Table 7
Optimum solution of the multi-scale inverse problem provided by the GA, for
benchmark BK1.
( ) ( ) ( ) ( )
Analysis name 𝜇 𝐸1𝑓 [GPa] 𝜎 𝐸1𝑓 [GPa] 𝜇 𝑉𝐹 𝜎 𝑉𝐹

REF 276 27.6 0.555 0.0555

GE1A 238 27.1 0.628 0.0562
GE1B 237 27.1 0.630 0.0562
GE1C 237 27.1 0.630 0.0564
GE2A 257 25.5 0.589 0.0595
GE2B 257 25.5 0.589 0.0595
GE2C 257 25.4 0.589 0.0597
GE3A 221 25.3 0.665 0.0597
GE3B 223 25.5 0.662 0.0593
GE3C 223 25.5 0.661 0.0593

belonging to each population: the reader is addressed to [37] for more

Fig. 9. Comparison between the validation set of samples and the results provided by
details about these aspects.
the ANN.
Inasmuch as the proposed strategy makes use of a metaheuristic
algorithm, the GA is run three times for each benchmark. The best
Table 5
Optimisation parameters for the genetic algorithm, for benchmarks BK1 and BK2. individual obtained at the end of the genetic calculation is used as a
Parameters BK1 BK2
starting guess for the gradient-based algorithm, in order to execute the
subsequent local optimisation.
N. of individuals 40 80
N. of populations 2 2 In terms of computational effort, the training phase of the ANN
N. of iterations 100 100 needs several seconds to be performed. Then, the hybrid optimisation
Crossover probability. 0.85 0.85 strategy needs 37.8 and 68.2 h for the benchmarks BK1 and BK2,
Mutation probability. 0.025 0.0125 respectively, on an Intel® Xeon® 2.70 GHz CPU with two processors
Isolation time 20 20
and with a RAM of 128 GB.
The results provided by the ERASMUS GA for benchmarks BK1
Table 6 and BK2 are summarised in Tables 7 and 9, respectively, whilst those
Optimisation parameters for the gradient-based algorithm, for benchmarks BK1 and provided by the gradient-based algorithm are reported in Tables 8 and
BK2.
10, respectively. In order to compare the obtained results with the
Parameters BK1 BK2
reference ones, the average of the solutions provided by the gradient-
Solver Active-set Active-set based algorithm is performed for each identified parameter, as it can
Max n. of function evaluation 10000 10000
be seen in Tables 8 and 10, for each benchmark.
Tol. on the objective function 10−15 10−15
Tol. on the gradient norm 10−15 10−15 As it can be easily inferred from the analysis of these results, the
mean value and the standard deviation of the microscopic material
properties are in good agreement with the reference data: the absolute
percentage error ranges from 2.8% to 13.4% for the first benchmark and
are in very good agreement with the samples constituting the validation from 0.2% to 3.2% for the second benchmark. ( ) ( )
set. The discrepancy between the values of 𝜇 𝐸1𝑓 and 𝜇 𝑉𝐹 provided
by the MSIS and the reference ones, for benchmark BK1, is related
7.2. Numerical results of the MSIS for benchmarks BK1 and BK2 to the nature of the laminate stack. Indeed, for this benchmark, the
considered sequence has an isotropic membrane stiffness matrix but
a completely anisotropic bending stiffness matrix. This aspect has a
As discussed in Section 2, two benchmarks are investigated in order strong influence on the solution search for the multi-scale inverse
to show the effectiveness of the proposed MSIS, by varying the stacking problem because the first buckling load is dominated by the bending
sequence of the multilayer plate. stiffness of the laminate. In particular, if the bending stiffness matrix
The parameters tuning the GA and the deterministic algorithms are is not orthotropic, problem (21) becomes strongly non-convex and
summarised in Tables 5 and 6, respectively, according to the main several equivalent optimal solutions exist. Therefore, finding the global
guidelines described in [67]. minimum is anything but trivial in such a case.
The GA calculation is performed with two populations, in which, the These results prove that a particular care should be put in the choice
number of individuals, evolving along the selected maximum number of the stacking sequence, which strongly affects both the number of
of generations, depends on the considered benchmark. Indeed, the best parameters that is possible to identify and the quality of the final result.
practice is to set the number of individuals greater than or equal to ten
times the number of optimisation variables. Accordingly, benchmarks 8. Conclusions
BK1 and BK2 are characterised by two populations composed of 40 and
80 individuals, respectively. The two populations exchange the best In this work the multi-scale identification strategy (MSIS), initially
individual every ten iterations, by using a ring-type operator, whose presented in [36], has been extended to the characterisation of the
probability is automatically computed by the considered GA. Moreover, uncertainty of the geometrical and the elastic properties of the fibre and
as far as the constraint-handling technique is concerned, the Automatic the matrix at the microscopic scale, by using information restrained in
Dynamic Penalisation (ADP) method is used [68]. the macroscopic response of the laminate, i.e. the first buckling load of
It is noteworthy that, the choice of multiple populations, with the multilayer plate.
a small number of individuals, allows finding the global minimum In this context, the multi-scale characterisation problem, is stated
without increasing too much the computational effort. In this way, as a constrained inverse problem. The goal is the minimisation of the
the GA has the possibility to explore the design domain in the most distance between the numerical and the reference variability parame-
effective way, by exchanging information between the best individuals ters that describe the probability distribution of the first buckling load.

9
L. Cappelli et al. Composites Part B 176 (2019) 107193

Table 8
Optimum solution of the multi-scale inverse problem provided by the gradient-based algorithm, for benchmark BK1; the percentage
difference between the solution and the microscopic reference data are given in parentheses.
( ) ( ) ( ) ( )
Analysis name 𝜇 𝐸1𝑓 [GPa] 𝜎 𝐸1𝑓 [GPa] 𝜇 𝑉𝐹 𝜎 𝑉𝐹 𝛷 (𝐱)

REF 276 27.6 0.555 0.0555 0

GR1A 238 27.1 0.625 0.0564 3.01E−05
GR1B 244 27.6 0.616 0.0567 5.58E−07
GR1C 235 27.1 0.632 0.0526 1.62E−05
GR2A 257 25.6 0.588 0.0594 7.33E−06
GR2B 257 25.5 0.589 0.0596 3.00E−07
GR2C 255 25.2 0.594 0.0586 2.88E−07
GR3A 221 25.3 0.665 0.0597 2.48E−07
GR3B 223 27.1 0.663 0.0512 2.11E−06
GR3C 223 25.5 0.661 0.0593 7.37E−08
AVERAGE 239 (−13.4) 26.2(−5) 0.626 (12.8) 0.0571 (2.8)

Table 9
Optimum solution of the multi-scale inverse problem provided by the GA, for benchmark BK2.
( ) ( ) ( ) ( ) ( ) ( )
Analysis name 𝜇 𝐸1𝑓 [GPa] 𝜎 𝐸1𝑓 [GPa] 𝜇 𝐸2𝑓 [GPa] 𝜎 𝐸2𝑓 [GPa] 𝜇 (𝐸 𝑚 ) [GPa] 𝜎 (𝐸 𝑚 ) [GPa] 𝜇 𝑉𝐹 𝜎 𝑉𝐹

REF 276 27.6 17.3 1.73 4.14 0.414 0.555 0.0555

GE1A 255 31.2 15.7 1.50 3.64 0.412 0.620 0.0540
GE1B 311 30.3 20.1 2.02 4.49 0.366 0.482 0.0547
GE1C 313 30.2 20.0 1.59 4.47 0.363 0.482 0.0554
GE2A 304 22.4 18.3 1.51 3.43 0.484 0.568 0.0523
GE2B 304 22.4 18.3 1.51 3.43 0.484 0.568 0.0523
GE2C 258 24.5 16.0 1.93 4.23 0.467 0.581 0.0593
GE3A 233 24.8 16.0 1.59 4.16 0.404 0.615 0.0633
GE3B 233 24.8 16.0 1.59 4.16 0.409 0.615 0.0632
GE3C 233 24.8 16.0 1.60 4.16 0.404 0.615 0.0633

Table 10
Optimum solution of the multi-scale inverse problem provided by the gradient-based algorithm, for benchmark BK2; the percentage difference between the solution and the
microscopic reference data are given in parentheses.
( ) ( ) ( ) ( ) ( ) ( )
Analysis name 𝜇 𝐸1𝑓 [GPa] 𝜎 𝐸1𝑓 [GPa] 𝜇 𝐸2𝑓 [GPa] 𝜎 𝐸2𝑓 [GPa] 𝜇 (𝐸 𝑚 ) [GPa] 𝜎 (𝐸 𝑚 ) [GPa] 𝜇 𝑉𝐹 𝜎 𝑉𝐹 𝛷 (𝐱)

REF 276 27.6 17.3 1.73 4.14 0.414 0.555 0.0555 0

GR1A 288 31.3 17.7 1.74 4.24 0.433 0.520 0.0514 4.50E−04
GR1B 280 28.8 18.8 1.75 4.33 0.401 0.529 0.0539 1.77E−06
GR1C 306 29.4 19.5 1.58 4.38 0.377 0.497 0.0548 4.27E−07
GR2A 288 31.3 17.7 1.74 4.24 0.433 0.520 0.0514 1.71E−09
GR2B 304 22.4 18.3 1.51 3.43 0.484 0.568 0.0523 3.59E−08
GR2C 258 26.1 16.4 1.96 4.14 0.450 0.569 0.0579 2.76E−04
GR3A 236 24.7 16.0 1.61 4.14 0.402 0.610 0.0634 2.05E−05
GR3B 233 24.8 16.0 1.59 4.16 0.409 0.615 0.0632 1.45E−07
GR3C 233 24.7 16.2 1.59 4.14 0.407 0.617 0.0634 1.98E−07
AVERAGE 270 (−2.3) 27 (−2) 17.4 (0.5) 1.67 (−3.2) 4.13 (−0.2) 0.422 (1.9) 0.561 (1) 0.0569 (2.5)

In this case, the solution search is performed by a hybrid optimisation the effectiveness of the proposed MSIS two different stacking sequences
tool, in which, a metaheuristic algorithm and a gradient-based one have have been considered: the first benchmark is characterised by a sym-
been interfaced to solve the related non-convex CNLPP. metric quasi-isotropic stack, while the second one is characterised by a
The MSIS makes use of an analytical homogenisation scheme, symmetric orthotropic one.
i.e. the Chamis’ model, to perform the microscopic/mesoscopic scale
As a consequence, also the obtained results, in terms of the iden-
transition. The elastic properties of the elementary lamina evaluated
tification of the parameters tuning the variability of the elastic and
by means of the Chamis’ model are then used into the FE model of the
geometrical properties of the constitutive phases of the composite, are
multilayer plate to evaluate its first buckling load.
strongly influenced by the nature of the laminate lay-up. In particular,
Moreover, a Monte-Carlo simulation campaign has been performed
for the first benchmark the absolute percentage error ranges from 2.8%
to compute the probability distribution of the first buckling load, start- ( )
to 13.4% for the standard deviation of the fibre volume fraction 𝜎( 𝑉𝐹)
ing from a Gaussian probability distribution of the material properties
of the constituent phases. The obtained samples have been used to and the mean value of the fibre longitudinal elastic modulus 𝜇 𝐸1𝑓 ,
train an ANN which emulates the multi-scale mechanical response of respectively. Conversely, for the second benchmark the absolute per-
the plate: the inputs are the geometrical and elastic properties of the centage error ranges from 0.2% to 3.2% for the mean value of the
microscopic constituents of the composite and the output is the first- matrix elastic modulus 𝜇 (𝐸 𝑚() and
) the standard deviation of the fibre
buckling load of the laminate. Then, the obtained surrogate model has transverse elastic modulus 𝜎 𝐸2𝑓 .
been used into the optimisation process to reduce the computational
Nevertheless, thanks to the proposed MSIS, it is possible to retrieve
effort.
the variability of both longitudinal and transversal effective properties
Before executing the hybrid optimisation process, a sensitivity study of the constitutive phases and this task cannot be easily performed by
has been performed to determine the most relevant microscopic pa- means of standard ASTM tests.
rameters influencing the first buckling load at the macroscopic scale.
In particular, numerical results show that this sensitivity is strongly As far as perspectives of this work are concerned, research is ongo-
affected by the nature of the stacking sequence. Therefore to prove ing in order to include into the MSIS the following aspects:

10
L. Cappelli et al. Composites Part B 176 (2019) 107193

• the validation of the MSIS by means of experimental data. In this [14] ASTM D3379-75(1989)e1, Standard test method for tensile strength and young’s
case, the influence of noise on the results provided by the MSIS modulus for high-modulus single-filament materials (withdrawn 1998). West
Conshohocken, PA: ASTM International; 1975.
should be properly taken into account. To this purpose, suitable
[15] ASTM D638-14, Standard test method for tensile properties of plastics. West
regularisation techniques, as the Tikhonov–Morozov one, which Conshohocken, PA: ASTM International; 2014.
is widely used in different engineering fields [69,70], must be [16] Nairn JA. Analytical fracture mechanics analysis of the pull-out test including
efficiently integrated into the multi-scale identification process to the effects of friction and thermal stresses. Adv Compos Lett 2000;9(6):373–83.
handle noise; [17] Maurin R, Davies P, Baral N, Baley C. Transverse properties of carbon fibres by
nano-indentation and micro-mechanics. Appl Compos Mater 2018;15:61–73.
• the formulation of a dedicated optimisation problem to find a
[18] Feih S, Wonsyld K, Minzari D, Westermann P, Lilholt H. Testing procedure for
suitable stack which maximise the sensitivity of the first buckling the single fiber fragmentation test, vol. 1483(EN). Denmark: Forskningscenter
load to each parameter defined at the microscopic scale of the Risoe. Risoe-R; 2004.
composite; [19] Sepahvand K, Marburg S. On construction of uncertain material parameter
• the extension of the MSIS to the characterisation of the variabil- using generalized polynomial chaos expansion from experimental data. Procedia
IUTAM 2013;6:4–17, IUTAM Symposium on Multiscale Problems in Stochastic
ity of the viscoelastic properties of the microscopic constituent
Mechanics. http://dx.doi.org/10.1016/j.piutam.2013.01.001.
and the evaluation of the variability effects related to further [20] Pajonk O, Rosić BV, Litvinenko A, Matthies HG. A deterministic filter for non-
geometrical parameters, e.g. fibre misalignment, macroscopic ge- Gaussian Bayesian estimation— Applications to dynamical system estimation
ometrical defects, etc. with noisy measurements. Physica D 2012;241(7):775–88. http://dx.doi.org/10.
1016/j.physd.2012.01.001.
Finally, thanks to the versatility of the proposed MSIS, it is possible [21] Desceliers C, Ghanem R, Soize C. Maximum likelihood estimation of stochastic
chaos representations from experimental data. Internat J Numer Methods Engrg
to increase the accuracy in terms of variability parameters by introduc-
2006;66(6):978–1001. http://dx.doi.org/10.1002/nme.1576.
ing more general probability density functions for both the buckling [22] Marzouk YM, Najm HN, Rahn LA. Stochastic spectral methods for efficient
load and the microscopic parameters. In this way, it will be possible to Bayesian solution of inverse problems. J Comput Phys 2007;224(2):560–86.
go beyond the limits of the Gaussian model, in which the shape of the http://dx.doi.org/10.1016/j.jcp.2006.10.010.
distribution is imposed a priori. [23] Soize C. A computational inverse method for identification of non-Gaussian
random fields using the Bayesian approach in very high dimension. Comput
Methods Appl Mech Engrg 2011;200(45):3083–99. http://dx.doi.org/10.1016/j.
Acknowledgements cma.2011.07.005.
[24] Narayanan VAB, Zabaras N. Stochastic inverse heat conduction using a spectral
This research work has been carried out within the project FULL- approach. Internat J Numer Methods Engrg 2004;60(9):1569–93. http://dx.doi.
COMP (FULLy analysis, design, manufacturing, and health monitoring org/10.1002/nme.1015.
[25] Proppe C. Reliability computation with local polynomial chaos approxima-
of COMPosite structures), funded by the European Union Horizon 2020
tions. ZAMM Z Angew Math Mech 2009;89(1):28–37. http://dx.doi.org/10.1002/
Research and Innovation program under the Marie Skłodowska-Curie zamm.200800072.
grant agreement No. 642121. [26] Batou A, Soize C. Stochastic modeling and identification of an uncertain compu-
tational dynamical model with random fields properties and model uncertainties.
References Arch Appl Mech 2012;1–18.
[27] Wang J, Zabaras N. Hierarchical Bayesian models for inverse problems in heat
conduction. Inverse Problems 2004;21(1):183–206. http://dx.doi.org/10.1088/
[1] Jones RM. Mechanics of composite materials. McGraw-Hill; 1975.
0266-5611/21/1/012.
[2] Funari MF, Greco F, Lonetti P, Luciano R, Penna R. An interface approach
[28] Ghanem RG, Doostan A. On the construction and analysis of stochastic models:
based on moving mesh and cohesive modeling in Z-pinned composite laminates.
Characterization and propagation of the errors associated with limited data.
Composites B 2018;135:207–17. http://dx.doi.org/10.1016/j.compositesb.2017.
J Comput Phys 2006;217(1):63–81, Uncertainty Quantification in Simulation
10.018.
Science. http://dx.doi.org/10.1016/j.jcp.2006.01.037.
[3] Funari MF, Greco F, Lonetti P. A moving interface finite element formulation for
[29] Chen C, Duhamel D, Soize C. Probabilistic approach for model and data
layered structures. Composites B 2016;96:325–37. http://dx.doi.org/10.1016/j.
uncertainties and its experimental identification in structural dynamics: Case of
compositesb.2016.04.047.
composite sandwich panels. J Sound Vib 2006;294(1):64–81. http://dx.doi.org/
[4] Funari MF, Greco F, Lonetti P. A numerical model based on ALE formulation
10.1016/j.jsv.2005.10.013.
to predict crack propagation in sandwich structures. Fract. Struct. Integr.
2019;47:277–93, Ten years of ‘Frattura ed Integrità Strutturale’. http://dx.doi. [30] Rosić BV, Litvinenko A, Pajonk O, Matthies HG. Sampling-free linear
org/10.3221/IGF-ESIS.47.21. Bayesian update of polynomial chaos representations. J Comput Phys
[5] Naskar S, Mukhopadhyay T, Sriramula S. Probabilistic micromechanical 2012;231(17):5761–87. http://dx.doi.org/10.1016/j.jcp.2012.04.044.
spatial variability quantification in laminated composites. Composites B [31] Vu-Bac N, Rafiee R, Zhuang X, Lahmer T, Rabczuk T. Uncertainty quantification
2018;151:291–325. http://dx.doi.org/10.1016/j.compositesb.2018.06.002. for multiscale modeling of polymer nanocomposites with correlated parameters.
[6] Sepahvand K, Marburg S. Identification of composite uncertain material pa- Composites B 2015;68:446–64. http://dx.doi.org/10.1016/j.compositesb.2014.
rameters from experimental modal data. Probab Eng Mech 2014;37:148–53. 09.008.
http://dx.doi.org/10.1016/j.probengmech.2014.06.008. [32] Riley EJ, Koudela KL, Narayanan RM. Characterization of the electromag-
[7] ASTM D3039 / D3039M-17, Standard test method for tensile properties of poly- netic parameter uncertainty in single-ply unidirectional carbon-fiber-reinforced-
mer matrix composite materials. West Conshohocken, PA: ASTM International; polymer laminas. Composites B 2019;162:361–8. http://dx.doi.org/10.1016/j.
2017. compositesb.2018.10.089.
[8] ASTM D790-17, Standard test methods for flexural properties of unreinforced [33] Dey S, Mukhopadhyay T, Sahu S, Li G, Rabitz H, Adhikari S. Thermal uncertainty
and reinforced plastics and electrical insulating materials. West Conshohocken, quantification in frequency responses of laminated composite plates. Composites
PA: ASTM International; 2017. B 2015;80:186–97. http://dx.doi.org/10.1016/j.compositesb.2015.06.006.
[9] ASTM D3410 / D3410M-16, Standard test method for compressive properties [34] Dong C. Uncertainties in flexural strength of carbon/glass fibre reinforced hybrid
of polymer matrix composite materials with unsupported gage section by shear epoxy composites. Composites B 2016;98:176–81. http://dx.doi.org/10.1016/j.
loading. West Conshohocken, PA: ASTM International; 2016. compositesb.2016.05.035.
[10] ASTM D5379 / D5379M-12, Standard test method for shear properties of [35] Alazwari MA, Rao SS. Modeling and analysis of composite laminates in the
composite materials by the v-notched beam method. West Conshohocken, PA: presence of uncertainties. Composites B 2019;161:107–20. http://dx.doi.org/10.
ASTM International; 2012. 1016/j.compositesb.2018.10.052.
[11] ASTM D7078 / D7078M-12, Standard test method for shear properties of [36] Cappelli L, Montemurro M, Dau F, Guillaumat L. Characterisation of composite
composite materials by v-notched rail shear method. West Conshohocken, PA: elastic properties by means of a multi-scale two-level inverse approach. Compos
ASTM International; 2012. Struct 2018;204:767–77. http://dx.doi.org/10.1016/j.compstruct.2018.08.007.
[12] ASTM D3518 / D3518M-13, Standard test method for in-plane shear response [37] Montemurro M. A contribution to the development of design strategies for the
of polymer matrix composite materials by tensile test of a 45 Laminate. West optimisation of lightweight structures [Hdr thesis], Université de Bordeaux; 2018,
Conshohocken, PA: ASTM International; 2013. URL http://hdl.handle.net/10985/15155.
[13] ASTM D2344 / D2344M-16, Standard test method for short-beam strength of [38] Montemurro M. Optimal design of advanced engineering modular systems
polymer matrix composite materials and their laminates. West Conshohocken, through a new genetic approach [Ph.D. thesis], France: UPMC, Paris VI; 2012,
PA: ASTM International; 2016. URL http://tel.archives-ouvertes.fr/tel-00955533.

11
L. Cappelli et al. Composites Part B 176 (2019) 107193

[39] Optimization Toolbox User’s Guide. 3 Apple Ill Drive, Natick, MA 01760-2098: [57] Dey S, Mukhopadhyay T, Adhikari S. Metamodel based high-fidelity stochastic
The MathWorks, Inc.; 2017. analysis of composite laminates: A concise review with critical compara-
[40] Chamis CC. Mechanics of composite materials: Past, present, and future. J tive assessment. Compos Struct 2017;171:227–50. http://dx.doi.org/10.1016/j.
Compos Technol Res 1989;11(1):3–14. compstruct.2017.01.061.
[41] Balokas G, Czichon S, Rolfes R. Neural network assisted multiscale analysis for [58] Montemurro M, Catapano A, Doroszewski D. A multi-scale approach for the
the elastic properties prediction of 3D braided composites under uncertainty. simultaneous shape and material optimisation of sandwich panels with cellu-
Compos Struct 2018;183:550–62. http://dx.doi.org/10.1016/j.compstruct.2017. lar core. Composites B 2016;91:458–72. http://dx.doi.org/10.1016/j.compstruct.
06.037. 2014.07.058.
[42] Schuller G, Jensen H. Computational methods in optimization considering un- [59] Montemurro M, Catapano A. On the effective integration of manufacturability
certainties – An overview. Comput Methods Appl Mech Engrg 2008;198(1):2–13.
constraints within the multi-scale methodology for designing variable angle-
http://dx.doi.org/10.1016/j.cma.2008.05.004.
tow laminates. Compos Struct 2017;161:145–59. http://dx.doi.org/10.1016/j.
[43] Enevoldsen I, Sørensen J. Reliability-based optimization in structural engineering.
compstruct.2016.11.018.
Struct Saf 1994;15(3):169–96. http://dx.doi.org/10.1016/0167-4730(94)90039-
[60] Montemurro M, Catapano A. A general b-spline surfaces theoretical frame-
6.
work for optimisation of variable angle-tow laminates. Compos Struct
[44] Gasser M, Schuëller G. Reliability-based optimization of structural systems. Math
Methods Oper Res 1997;46(3):287–307. http://dx.doi.org/10.1007/BF01194858. 2019;209:561–78. http://dx.doi.org/10.1016/j.compstruct.2018.10.094.
[45] Jensen HA. Design and sensitivity analysis of dynamical systems subjected to [61] Costa G, Montemurro M, Pailhès J. A general hybrid optimization strategy for
stochastic loading. Comput Struct 2005;83(14):1062–75. http://dx.doi.org/10. curve fitting in the non-uniform rational basis spline framework. J Optim Theory
1016/j.compstruc.2004.11.016. Appl 2018;176(1):225–51. http://dx.doi.org/10.1007/s10957-017-1192-2.
[46] Papadrakakis M, Lagaros N, Plevris V. Design optimization of steel structures [62] Montemurro M, Izzi MI, El-Yagoubi J, Fanteria D. Least-weight composite
considering uncertainties. Eng Struct 2005;27(9):1408–18. http://dx.doi.org/10. plates with unconventional stacking sequences: Design, analysis and ex-
1016/j.engstruct.2005.04.002. periments. J Compos Mater 2019;51(16):2209–27. http://dx.doi.org/10.1177/
[47] Doltsinis I, Kang Z. Robust design of structures using optimization methods. 0021998318824783.
Comput Methods Appl Mech Engrg 2004;193(23):2221–37. http://dx.doi.org/10. [63] Panettieri E, Montemurro M, Catapano A. Blending constraints for composite
1016/j.cma.2003.12.055. laminates in polar parameters space. Composites B 2019;168:448–57. http://dx.
[48] Farhat C, Hemez F. Updating finite element dynamic models using an element- doi.org/10.1016/j.compositesb.2019.03.040.
by-element sensitivity methodology. AIAA J 1993;31(9):1702–11. http://dx.doi. [64] Bertolino G, Montemurro M, Pasquale GD. Multi-scale shape optimisation of
org/10.2514/3.11833. lattice structures : an evolutionary-based approach. Int J Interact Des Manuf
[49] Hemez FM, Doebling SW. Review and assessment of model updating for non- 2019 [in press]. http://dx.doi.org/10.1007/s12008-019-00580-9.
linear, transient dynamics. Mech Syst Signal Process 2001;15(1):45–74. http: [65] Audoux Y, Montemurro M, Pailhes J. A surrogate model based on non-uniform
//dx.doi.org/10.1006/mssp.2000.1351. rational b-splines hypersurfaces. Procedia CIRP 2018;70:463–8, 28th CIRP Design
[50] Soutis C, Beaumont PWR, editors. Multi-scale modelling of composite material
Conference 2018, 23-25 May 2018, Nantes, France. http://dx.doi.org/10.1016/
systems. the art of predictive damage modelling. In: Woodhead Publishing Series
j.procir.2018.03.234.
in Composites Science and Engineering, New York: Elsevier; 2005.
[66] Garulli T, Catapano A, Montemurro M, Jumel J, Fanteria D. Quasi-trivial stacking
[51] Hexply F584. Hexcell Corporation; 2016.
sequences for the design of thick laminates. Compos Struct 2018;200:614–23.
[52] Feo L, Greco F, Leonetti L, Luciano R. Mixed-mode fracture in lightweight
http://dx.doi.org/10.1016/j.compstruct.2018.05.120.
aggregate concrete by using a moving mesh approach within a multiscale frame-
[67] Montemurro M, Nasser H, Koutsawa Y, Belouettar S, Vincenti A, Vannucci P.
work. Compos Struct 2015;123:88–97. http://dx.doi.org/10.1016/j.compstruct.
2014.12.037. Identification of electromechanical properties of piezoelectric structures through
[53] Bruno D, Greco F, Luciano R, Blasi PN. Nonlinear homogenized properties of evolutionary optimisation techniques. Int J Solids Struct 2012;49(13):1884–92.
defected composite materials. Comput Struct 2014;134:102–11. http://dx.doi. [68] Montemurro M, Vincenti A, Vannucci P. The automatic dynamic penalisation
org/10.1016/j.compstruc.2013.11.018. method (ADP) for handling constraints with genetic algorithms. Comput Methods
[54] Barbero E. Finite element analysis of composite materials. CRC Press, Taylor and Appl Mech Engrg 2013;256:70–87.
Francis Group; 2007. [69] Faghidian S. Inverse determination of the regularized residual stress and eigen-
[55] Système D. Abaqus Analysis user’s guide. Providence, RI, USA: Dassault Systèmes strain fields due to surface peening. J Strain Anal Eng Des 2015;50(2):84–91.
Simulia Corp.; 2013. http://dx.doi.org/10.1177/0309324714558326.
[56] Sobol IM. Global sensitivity indices for nonlinear mathematical models and [70] Faghidian SA. A regularized approach to linear regression of fatigue life mea-
their Monte Carlo estimates. Math Comput Simulation 2001;55(1–3):271–80. surements. Int J Struct Integr 2016;7(1):95–105. http://dx.doi.org/10.1108/IJSI-
http://dx.doi.org/10.1016/S0378-4754(00)00270-6. 12-2014-0071.

12