Lorenzo Cappelli

Marco Montemurro
email: marco.montemurro@ensam.eu

Frédéric Dau

Laurent Guillaumat

Multi-scale identification of the viscoelastic behaviour of composite materials through a non-destructive test

Keywords: Composite materials, Homogenisation, Harmonic analysis, Inverse problems, Optimisation, Viscoelasticity List of acronyms ADP: Automatic Dynamic Penalisation ANM: Asymptotic Numerical Method BCs: Boundary Conditions CNLPP: Constrained non-linear programming problem DOF: Degree Of Freedom Fast Fourier Transform FRF: Frequency Response Function Inverse Iteration algorithm ISIM: Iterative Shift-Inverter Method LB: Lower Bound MSIS: Multi-Scale Identification Strategy NLAM: Non-Linear Arnoldi's method

diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Composite materials have a long story: they were employed by Israelites for the manufacturing of strengthen mud bricks and by the Egyptians to produce plywood [START_REF] Jones | Mechanics of composite materials[END_REF]. Nowadays, because of the introduction of constitutive phases with adequate mechanical properties in terms of strength, stiffness and heat transfer behaviour, high-performance engineering structures are more and more made of composite materials.

Composite structures are subjected to a wide variety of loading conditions including both static and dynamic loads: in order to decrease the design costs related to experimental tests, accurate numerical simulations coupled to cheaper tests are required to predict the complex behaviour of the structure.

In order to carry out numerical simulations, reliable material properties have to be defined, under both static and dynamic conditions. In particular, for fibre-reinforced composites, the macroscopic behaviour depends upon that of the constitutive phases at the lower scales. In particular, for multilayer plates three characteristic scales can be identified. At the macroscopic scale the laminate is usually modelled as a homogeneous anisotropic plate whose structural response depends upon the constitutive stiffness, mass and damping matrices (membrane, bending and membrane/bending coupling behaviours). The mesoscopic scale focuses on the lamina-level: each constitutive ply is characterised by some geometric and material parameters, i.e. the orientation angle, the thickness, the position and the material properties of the ply. Of course, each one of the previous parameters is involved in the definition of the laminate stiffness and mass matrices and has a strong influence on the macroscopic response of the multilayer plate. Finally, the microscopic scale is that of the constitutive phases (e.g. fibres and matrix for fibre-reinforced composites): at this scale each phase is characterised by a given material behaviour and by a set of geometrical parameters, e.g. fibres volume fraction, fibre shape, fibres arrangement, etc.

The characterisation of the composite material properties at each relevant scale is a rather complex problem. Concerning the identification of the elastic properties, ASTM standard tests can be found in literature: mesoscopic destructive tests (as tension [START_REF]D3039M-17[END_REF], three/four points bending [START_REF]Standard Test Methods for Flexural Properties of Unreinforced and Reinforced Plastics and Electrical Insulating Materials[END_REF] and shear tests [START_REF]Standard Test Method for Shear Properties of Composite Materials by the V-Notched Beam Method[END_REF]) and few microscopic ASTM standard tests are available (to characterise the matrix [START_REF]Standard Test Method for Tensile[END_REF] and the fibre [START_REF]Standard Test Method for Tensile Strength and Young's Modulus for High-Modulus Single-Filament Materials[END_REF] longitudinal Young's modulus, respectively). Unfortunately, the above mentioned tests are not able to provide the required 3D set of material elastic properties.

Moreover, composite materials show also a dynamical behaviour that is strongly affected by damping properties [START_REF] Mahmoudi | Experimental and numerical investigation of the damping of flaxepoxy composite plates[END_REF]. Accordingly, a proper characterisation of the damping capability of the material, at each relevant scale, is a challenging task [START_REF] Suarez | The influence of fiber length and fiber orientation on damping and stiffness of polymer composite materials[END_REF]. This problem is more difficult than that of the elastic properties characterisation essentially because of the non-linear nature of the viscoelastic matrix behaviour, in terms of time response, which influences the damping capability of the composite at all characteristic scales.

From an experimental point of view, two methods are commonly used: the direct method and the indirect one. On the one hand, the direct techniques are based on the measurement of the dissipated energy per load cycle, which can be evaluated from the area of an hysteresis loop [START_REF] Krasnobrizha | Hysteresis behaviour modelling of woven composite using a collaborative elastoplastic damage model with fractional derivatives[END_REF]. On the other hand, the indirect methods allow estimating the dissipated energy from the analysis of the spectrum response: free vibration-decay, resonant-dwell, bandwidth and impedance methods are some of the experimental techniques used for damping characterisation [START_REF] Kostopoulos | A new method for the determination of viscoelastic properties of composite laminates: a mixed analytical-experimental approach[END_REF].

However, it is possible to describe the damping behaviour of a composite structure by directly looking at the viscoelastic properties [START_REF] Jayendiran | Micromechanical modeling and experimental characterization on viscoelastic behavior of 13 active composites[END_REF][START_REF] Swain | Viscoelastic modelling and dynamic characteristics of cnts-cfrp-2dwf composite shell structures[END_REF]: nowadays, a common and useful method is the so-called Dynamical Mechanical Analysis (DMA) [START_REF] Melo | Time and temperature dependence of the viscoelastic properties of CFRP by dynamic mechanical analysis[END_REF][START_REF] Finegan | Modeling and characterization of damping in carbon nanofiber/polypropylene composites[END_REF][START_REF] Chandra | Experimental evaluation of damping of fiber-reinforced composites[END_REF]. DMA is an indirect method to characterise the material properties of reinforced polymers in terms of thermal, elastic and viscoelastic behaviours [START_REF] Abedi | Viscoelastic characterization of out-of-autoclave composite laminates: Experimental and finite element studies[END_REF].

The DMA test is performed by applying harmonic loads to the specimen. By measuring the sample response, it is possible to compute an apparent modulus that can be used to estimate the viscoelastic material parameters of the specimen. In the case of a composite multilayer plate wherein the lamina has an isotropic transverse behaviour [START_REF] Suarez | The influence of fiber length and fiber orientation on damping and stiffness of polymer composite materials[END_REF], the identification process has to be carried out three different times, e.g. by considering a symmetric angle-ply stack, to determine the longitudinal E L , transversal E T and shear G LT moduli.

Unfortunately, when high modulus composite materials are investigated, the DMA technique provides less accurate results [START_REF] Abedi | Viscoelastic characterization of out-of-autoclave composite laminates: Experimental and finite element studies[END_REF] compared to the ASTM three-points bending test [START_REF]Standard Test Methods for Flexural Properties of Unreinforced and Reinforced Plastics and Electrical Insulating Materials[END_REF]. Indeed, the DMA test provides an apparent modulus giving only an average approximation of the plate flexural stiffness response which group both structural and material aspects. However, the DMA test does not allow to extract information about microscopic properties and only fibre-reinforced polymers can be tested.

From the engineer's viewpoint it is more interesting to look for those tests which allow to identify material properties at all relevant scales and which are not limited by the size of the composite sample or by the geometrical and material properties of the constitutive phases composing it. The formulation of a suitable inverse problem for material proper-ties characterisation is a widely studied topic in the literature [START_REF] Meng | An objective meta-modeling approach for indentation-based material characterization[END_REF][START_REF] Yap | A non-destructive experimental-cum-numerical methodology for the characterization of 3d-printed materialspolycarbonate-acrylonitrile butadiene styrene (pc-abs)[END_REF][START_REF] Ghorbal | Mechanical characterization of brittle materials using instrumented indentation with knoop indenter[END_REF]. In this background, Barkanov et al. [20] proposed an inverse technique based on modal analysis and on the response surface method to characterise the nonlinear behaviour of the viscoelastic core layer in sandwich panels. Elkhaldi et al. [START_REF] Elkhaldi | A gradient method for viscoelastic behaviour identification of damped sandwich structures[END_REF] worked on the viscoelastic parameters identification for a sandwich panel where a generalised Maxwell model is considered and a gradient algorithm is used to solve the associated inverse problem. Cortés et al. [START_REF] Cortés | An approximate numerical method for the complex eigenproblem in systems characterised by a structural damping matrix[END_REF] developed an identification strategy to characterise the parameters of the fractional derivative model representing the viscoelastic behaviour of a sandwich beam. The goal is the minimisation of the error between the predicted Frequency Response Function (FRF) and the measured one. Ledi et al. [START_REF] Ledi | An identification method for frequency dependent material properties of viscoelastic sandwich structures[END_REF] proposed an identification method for frequency dependent material properties of viscoelastic sandwich beams able to take into account for the property of the interface between layers.

As it can be inferred from the aforementioned works, the damping capability related to the viscoelastic behaviour of the matrix can be characterised by exploiting the information restrained in the dynamic response of the structure. In these works, sandwich beams/plates manufactured by interposing a viscoelastic layer between metallic ones were considered because this configuration is well suited to reduce noise and vibration. However, in multilayer plates the macroscopic damping capability is mainly related to the viscoelastic behaviour of the matrix at the microscopic scale.

This work focuses on the damping capability of multilayer plates made of unidirectional fibre-reinforced laminae. In particular, this study aims to generalise the multi-scale identification strategy (MSIS) developed in a previous study [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF] to the case of the viscoelastic behaviour of composites. The MSIS relies on the information restrained in a non-destructive harmonic test conducted at the macroscopic scale. The idea is to exploit this information to characterise the viscoelastic behaviour of the constitutive phases at the microscopic scale.

In the context of the MSIS, the multi-scale identification problem is stated as a constrained non-linear programming problem (CNLPP). The goal is to minimise the distance between a reference harmonic response (that can be obtained either experimentally or numerically) and the numerical one. This function is subject to some requirements involved at different scales: (a) on the positive definiteness of the stiffness tensor of the constitutive phases (microscopic scale); (b) on the damped natural frequencies of the composite (macroscopic scale); (c) on a non-negative internal work and a non-negative dissipation rate as far as the viscoelastic model is concerned.

Nevertheless, the identification of the viscoelastic behaviour of the constitutive phases (mostly due to the matrix) at the microscopic scale is characterised by two difficulties: (a) the equivalent viscoelastic properties of the constitutive lamina at the mesoscopic scale depend upon the frequency ; (b) since the ply material properties depend upon the frequency, the problem of determining the structure damped natural frequencies becomes a non-linear eigenvalue problem, thus a suitable iterative method must be foreseen to perform the related modal analysis.

Therefore, the harmonic response of the laminate, at the macroscopic scale, is strongly affected by the matrix viscoelastic behaviour.

In this context, the MSIS presented in [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF] is generalised here to the characterisation of the parameters of the law tuning the viscoelastic behaviour of the constitutive phases, for a given frequency range. The proposed MSIS relies upon the following features: (a) an hybrid optimisation tool called HERO (Hybrid EvolutionaRy-based Optimisation) algorithm, see [START_REF] Montemurro | A contribution to the development of design strategies for the optimisation of lightweight structures[END_REF]; (b) an extension of the numerical homogenisation method based on the strain energy of periodic media and on volume-averaged stresses and strains [START_REF] Barbero | Finite element analysis of composite materials[END_REF] to the case of viscoelastic materials; (c) the Arnoldi's method [START_REF] Voss | An Arnoldi Method for Nonlinear Eigenvalue Problems[END_REF][START_REF] Hamdaoui | Comparison of non-linear eigensolvers for modal analysis of frequency dependent laminated visco-elastic sandwich plates[END_REF] to solve the non-linear modal analysis for materials with frequency-dependent viscoelastic properties. The effectiveness of the MSIS for viscoelastic materials is proven on a meaningful benchmark taken from the literature.

The paper is structured as follows. The problem and the MSIS are introduced in Section 2. The mathematical formulation of the inverse problem and the related numerical aspects are discussed in Section 3. The finite element (FE) models of the composite at both microscopic and macroscopic scales are presented in Section 4. The numerical results of the MSIS are illustrated and discussed in Section 5. Finally, Section 6 ends the paper with some conclusions and perspectives.

Multi-scale identification of composite viscoelastic properties

Problem description

In this work, the MSIS is applied to a rectangular composite plate made of unidirectional viscoelastic plies, whose geometrical parameters are shown in Figure 1. The constitutive ply is made of carbon-epoxy fibre Hexcel T 650/F 584 pre-impregnated tapes, whose fibre volume fraction is V f = 0.555. Concerning the available material properties of the constitutive phases, only the elastic properties can be found in litterature (they are taken from [START_REF]Multi-scale modelling of composite material systems. The art of predictive damage modelling[END_REF]). The parameters tuning the viscoelastic response of the F 584 matrix are set a priori because they are not available in [START_REF]Multi-scale modelling of composite material systems. The art of predictive damage modelling[END_REF]. They are set to reasonable values to give all the necessary microscopic material parameters defining the matrix behaviour. The reference material properties for both the fibre and the matrix are reported in Table 2. It is noteworthy that the viscoelastic behaviour of the F 584 matrix is described by means of the Bagley-Torvik model (briefly discussed in Section 2.2 and taken from [START_REF] Bagley | On the fractional calculus model of viscoelastic behavior[END_REF]).

The reference laminate is constituted of sixteen identical plies with the stacking sequence [0 • / -45 • /45 • /90 • /45 • /90 • / -45 • /0 • ] S . The average thickness of the ply is t ply = 0.28225 mm and the orientation angle of the lamina is defined positive according to counter-clockwise rotation around the z-axis.

The goal of this study is to provide a numerical validation of the MSIS for viscoelastic materials. To this purpose, the reference response of the structure is determined by means of modal and harmonic analyses performed on the reference configuration of the laminate described above. As described in Section 5.1, the reference material properties of the constitutive phases are implemented at the microscopic scale in order to determine the reference effective viscoelastic properties of the single lamina. Due to the viscoelastic behaviour of the matrix, the lamina elastic properties evolve according to the frequency. This variation is determined by generalising, to the viscoelastic case, the well-known homogenisation technique for periodic media based on the strain energy [START_REF] Barbero | Finite element analysis of composite materials[END_REF], as detailed in Section 4.1. Finally, the reference harmonic response and the reference natural frequencies of the laminate are determined, at the macroscopic scale, on a FE model of the multilayer plate making use of the reference properties provided by the homogenisation method.

In order to easily follow the flow of information throughout the manuscript, Table 1 lists all the adopted models, techniques and methods which are needed to achieve the ambitious goal of the multi-scale identification of the viscoelastic behaviour of the composite.

Model

The Bagley-Torvik viscoelastic material model

Composite structures show a dynamical behaviour that is significantly influenced by the damping capability of the matrix. The time-dependent response of materials can be classified into elastic (crystalline materials), viscous and viscoelastic. A viscoelastic material can be characterised by either a linear or a non-linear time-strain relationship and it can be in the form of a liquid (unrecoverable viscous flow) or a solid (fully recoverable viscous deformation).

From a numerical point of view, different linear viscoelastic material models are available in literature. These models are usually implemented to fit experimental data (usually creep and relaxation tests) [START_REF] Barbero | Finite element analysis of composite materials[END_REF][START_REF] Krasnobrizha | Modelisation des mecanismes d'hysteresis des composites tisses a l'aide d'un modele collaboratif elasto-plastique endommageable a derivees fractionnaires[END_REF]. The most common models are: (a) Maxwell model; (b) Kelvin-Voigt model; (c) Zener model; (d) power law-based models; (e) Prony seriesbased models; (f) generalised Kelvin model. These models differ essentially in terms of the number of parameters required to get the best possible fitting of experimental data. However, these laws are usually applied in the α-region of polymer creep (the characteristic time of the load application varies from seconds to years) [START_REF] Barbero | Finite element analysis of composite materials[END_REF]. Nevertheless, when looking at the frequency range characterising the application presented in this study (see Section 5.1), one can state that harmonic excitation falls in the β-region of polymer creep [START_REF] Barbero | Finite element analysis of composite materials[END_REF]; accordingly, a different viscoelastic model must be considered.

The effectiveness of a mathematical model in describing the viscoelastic behaviour of a given material can be seen as the ability of fitting a set of data points by using the least number of parameters tuning the model. Moreover, from a computational viewpoint and for optimisation purposes, the interest is always to have a limited number of parameters to be identified without reducing the accuracy of the model. Among the most effective mathematical representations of the viscoelasticity, the models based on fractional derivatives have been widely studied in the last three decades [START_REF] Krasnobrizha | Modelisation des mecanismes d'hysteresis des composites tisses a l'aide d'un modele collaboratif elasto-plastique endommageable a derivees fractionnaires[END_REF]. For such models, the general constitutive law reads:

σ (t) + b m D β m σ (t) = E m 0 ε (t) + E m 1 D α m σ (t) , where α m , β m , E m 0 , E m 1 , b m ∈ R. (1)
In Eq. ( 1), D α is the fractional derivative operator which represents a generalisation of the concept of derivative of a function. Consider a function

f ∈ L 1 ([a, b]) | a, b ∈ R. If α ∈ R * +
, the fractional derivative of order α is defined as the fractional integral of order n -α derived n times (Riemann-Liouville definition -RL):

(D α a f ) RL (x) = 1 Γ (1 -α) d dx x a f (x) (x -t) α dt, where 0 < α < 1.
(

) 2
The function Γ (z) , with z ∈ R * + , is the extension of the factorial function to real numbers. The reader is addressed to [START_REF] Krasnobrizha | Hysteresis behaviour modelling of woven composite using a collaborative elastoplastic damage model with fractional derivatives[END_REF][START_REF] Krasnobrizha | Modelisation des mecanismes d'hysteresis des composites tisses a l'aide d'un modele collaboratif elasto-plastique endommageable a derivees fractionnaires[END_REF] for further information concerning the mathematical formulation.

The fractional derivative model has been generalised to materials presenting a reticular complex molecular structure, because the choice of 0 < α < 1 well reproduces the relaxation function of different kind of polymers, as widely discussed in [START_REF] Schiessel | Applications to problems in polymer physics and rheology[END_REF]. The Fast Fourier Transform (FFT, whose operator is indicated as F ) of the fractional derivative operator can be computed as:

F [(D α a f ) RL (x)] = (Ωi) α F [f (x)] , ∀α ∈ R * + | 0 < α < 1 ∧ ∀x ∈ R + , (3)
where Ω is the frequency and i the imaginary unit. It is noteworthy that working in frequency domain allows dealing with a very simple mathematical expression of the fractional derivative operator. Applying the FFT to the Eq. ( 1), one obtains

σ (Ω) = F [σ (t)] = F [E m (t)] F [ε (t)] = Êm (Ω) ε (Ω) , (4)
where the relaxation modulus Êm (Ωi) reads [START_REF] Bagley | On the fractional calculus model of viscoelastic behavior[END_REF] Êm

(Ω) = E m 0 + E m 1 (Ωi) α m 1 + b m (Ωi) β m , (5)
that represents the Bagley-Torvik viscoelastic model. It can be observed that, when Ω → +∞, the relaxation modulus of the matrix Êm (Ω) → E m 1 b m (iΩ) α m -β m (vitreous domain). Otherwise, when Ω → 0, the relaxation modulus tends to the elastic constant E m 0 . However, a reliable viscoelastic model must be characterised by a non-negative internal work and a non-negative rate of dissipation, as highlighted in [START_REF] Bagley | On the fractional calculus model of viscoelastic behavior[END_REF]. In order to ensure these thermodynamical properties, the material parameters must satisfy the following relationships:

0 < α m = β m < 1, E m 0 > 0 and E m 1 > E m 0 b m . (6)
Therefore, only four parameters are needed to describe the viscoelastic behaviour of the considered polymeric matrix.

Fibre properties

E f 1 [MPa] E f 2 [MPa] ν f 12 ν f 23 G f 12 [MPa] 276000 17300 0.25 0.428 11240 Matrix properties E m 0 [MPa] E m 1 [MPa] b m α m = β m ν m 4140 30 0.0053 0.5 0.35
Table 2: Reference material properties for the fibre T 650/35 -3K and the epoxy matrix F 584 (taken from [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF][START_REF]Multi-scale modelling of composite material systems. The art of predictive damage modelling[END_REF][START_REF] Corporation | [END_REF]).

The multi-scale identification strategy

The goal of the MSIS is to find the optimum value of the parameters tuning the viscoelastic behaviour of the composite, at each scale, by smartly exploiting the information restrained into the harmonic response, measured in some precise locations, of the multilayer plate.

The reference response can be obtained either by a non-destructive experimental harmonic test (e.g. performed with shaker, hammer or solenoidal excitation system) or by carrying out a numerical harmonic test on the reference structure. This work deals with the latter case : the reference configuration of the multilayer plate as well as the reference dynamical results are presented in Section 5.

The MSIS aims to identify the parameters defining both the elastic behaviour of the fibre and the viscoelastic behaviour of the matrix by using the information available into the dynamical response of the composite at the macroscopic scale. The proposed approach relies on some hypotheses. As far as the microscopic scale is concerned, the following hypotheses are considered:

• the matrix has a viscoelastic isotropic behaviour, described by the Bagley-Torvik model, with a constant Poisson's ratio according to [START_REF] Luciano | Analytical expressions for the relaxation moduli of linear viscoelastic composites with periodic microstructure[END_REF];

• the fibre has an elastic transversely isotropic behaviour;

• the fibre-matrix interface is perfect (i.e. perfect bonding condition between the two phases).

Regarding the mesoscopic and macroscopic scales, the following hypotheses apply:

• the constitutive lamina has a viscoelastic orthotropic behaviour with only six parameters, due to the plane of symmetries characterising the considered Representative Volume Element (RVE), as shown in Figure 3a;

• perfect bonding condition at the interface between two consecutive plies;

• the first-order shear deformation theory (FSDT) is considered to describe the kinematics of the multilayer plate.

The general flow-chart of the one-shot MSIS for viscoelastic materials is shown in Figure 2.

Mathematical formulation of the multi-scale inverse problem

Optimisation variables, objective function and constraints

The multi-scale identification problem considered in this work is formulated as a classical inverse problem: the goal is the determination of the material properties of the constitutive phases of the composite by minimising the euclidean distance between the reference harmonic response at macroscopic scale and that provided by the numerical simulation.

According to the hypotheses introduced in Section 2, the effective elastic properties of the ply depend upon the frequency (due to the viscoelastic behaviour of the matrix) and must be computed through a suitable numerical homogenisation procedure. To achieve this task, a dedicated FE model of the RVE at the microscopic scale is built to perform the numerical homogenisation: the results are the frequency-dependent elastic constants of the lamina. In this analysis, the fibre volume fraction V f is set up a priori, because usually it is a value that the manufacturer can provide reliably, as highlighted also in [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF].

Considering the behaviour of the constituent phases, as discussed in Section 2.3, the design variables of the inverse problem, i.e. the material parameters to be identified, are five elastic constants for the fibre, supposed to hava an isotropic transverse behaviour and five parameters tuning the law of the isotropic viscoelastic matrix. In particular, as discussed in [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF], the fibre Poisson's coefficient ν f 23 is not considered among the design variables due to negligible sensitivity of the harmonic and modal responses of the plate to this parameter. Indeed, the reference plate (illustrated in Figure 1) is very thin and ν f 23 does not significantly influences the macroscopic dynamical response of the structure. Accordingly, ν f 23 has been set to the reference value given in Table 2. Therefore, the nine material properties involved at the composite microscopic scale can be collected in the vector of design variable x, as:

x = E f 1 , E f 2 , G f 12 , ν f 12 , E m 0 , E m 1 , b m , α m , ν m . (7)
In order to guarantee the positive definiteness of the fibre and matrix stiffness tensors [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF] and the thermodynamic requirements related to the viscoelastic behaviour of the matrix (see Section 2.2), the following set of non-linear constraints must be considered:

g 1 (x) = |ν f 12 | - E f 1 E f 2 < 0, g 2 (x) = E f 1 E f 2 2ν f 23 ν f 12 2 + 2ν f 12 2 -1 < 0, g 3 (x) = -E m 1 + E m 0 b m < 0, g 4 (x) = -α m < 0, g 5 (x) = -E m 0 < 0, g 6 (x) = ν m - 1 2 < 0, g 7 (x) = -ν m -1 < 0. (8)
In order to get a numerical harmonic spectrum really close to the reference one (and also to match the reference damped natural frequencies), a set of constraints on the laminate (damped) eigenfrequencies must be integrated into the problem formulation:

g 7+i (x) = f in -f ref in f ref in -ǫ i ≤ 0, i = 1, ..., n f . (9)
In Eq. ( 9), n f is the number of natural damped frequencies falling in the selected frequency spectrum range (as discussed in Section 5.1), whilst f in and f ref in are the i-th computed and reference damped eigenfrequency, respectively. ǫ i is a user-defined tolerance that establish the relative error for each eigenfrequency: here, a maximum relative error equal to 0.005 is chosen. The microscopic material parameters vary within the design space defined in Table 3. The objective function Φ(x) is defined by introducing an Euclidean distance between the reference and the numerical harmonic responses, both for real and imaginary parts. In particular, an error estimator of the least-squares type has been chosen:

Material properties

Φ(x) = Np q=1 Ns r=1 2 f r -f ref r f ref r 2 + ℜ H r,q (x) -H ref r,q ℜ H ref r,q 2 + ℑ H r,q (x) -H ref r,q ℑ H ref r,q 2 . (10)
In Eq. ( 10), f r is the r-th sampled frequency, while H r,q is FFT of the FRF determined at the q-th sample point of the multilayer plate and evaluated at the r-th sampled frequency. Of course, f ref r , H ref r,q are the same quantities evaluated on the reference configuration of the laminate. ℜ (• • • ) and ℑ (• • • ) represent the real and imaginary part, whilst N s and N p are the number of sampled frequencies and of sample points over the laminate plate (where the FRF is computed/measured), respectively. These quantities are detailed in Sections 4 and 5.

Finally, the multi-scale inverse problem can be stated as a classical CNLPP:

min x Φ (x) ,
subject to:

g j (x) ≤ 0, j = 1, ..., 7 + n f . (11)

The numerical strategy

Problem [START_REF] Jayendiran | Micromechanical modeling and experimental characterization on viscoelastic behavior of 13 active composites[END_REF] is highly non linear and non-convex in terms of both constraint and objective functions, see Eqs. ( 8), ( 9) and [START_REF] Kostopoulos | A new method for the determination of viscoelastic properties of composite laminates: a mixed analytical-experimental approach[END_REF].

For inverse problems, the uniqueness of solution is not a priori guaranteed: the set of parameters matching a given observed state may not be unique. Nevertheless, no proved theoretical rules exist in literature [START_REF] Sun | Inverse Problems in Groundwater Modelling[END_REF][START_REF] Tarantola | Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation[END_REF], to define the number of data points N p for a given number of unknowns (n) that have to be identified. Often, the inverse problem is stated as a CNLPP and it can be viewed as an over-determined system of equations [START_REF] Sun | Inverse Problems in Groundwater Modelling[END_REF][START_REF] Tarantola | Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation[END_REF]. Since more observation points than parameters exist (N p is usually much greater than n) there are more equations than unknowns. If an optimal point exists, of course it may be not unique, thus implying the existence of many combinations of parameters that result to be equivalent optimal solutions for the CNLPP at hand. Considering all these aspects and according to the practice always employed in literature, in this work a number of observed states (i.e. sample points N p ) greater than two times the number of design variables n has been considered. As explained in the next Section, the number of sample points has been inferred by means of a numerical sensitivity analysis of the FRF of the plate with respect to parameter N p : as a results N p = 329 has been chosen to properly perform the optimisation calculations.

Taking into account all of the aforementioned points, a hybrid optimization tool composed of the genetic algorithm (GA) ERASMUS (EvolutionaRy Algorithm for optimiSation of ModUlar Systems) developed by Montemurro [START_REF] Montemurro | A contribution to the development of design strategies for the optimisation of lightweight structures[END_REF], interfaced with the MATLAB fmincon algorithm [START_REF]Apple Ill Drive[END_REF], has been used. The GA ERASMUS has already been successfully applied to solve different kinds of real-world engineering problems, see for instance [START_REF] Montemurro | A multi-scale approach for the simultaneous shape and material optimisation of sandwich panels with cellular core[END_REF][START_REF] Costa | A General Hybrid Optimization Strategy for Curve Fitting in the Non-Uniform Rational Basis Spline Framework[END_REF][START_REF] Montemurro | Variational analysis and aerospace engineering: mathematical challenges for the aerospace of the future[END_REF][START_REF] Montemurro | On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates[END_REF][START_REF] Montemurro | Least-weight composite plates with unconventional stacking sequences: Design, analysis and experiments[END_REF][START_REF] Panettieri | Blending constraints for composite laminates in polar parameters space[END_REF][START_REF] Montemurro | A general multi--scale two-level optimisation strategy for designing composite stiffened panels[END_REF].

As shown in Figure 2, the optimisation procedure for problem [START_REF] Jayendiran | Micromechanical modeling and experimental characterization on viscoelastic behavior of 13 active composites[END_REF] is split in two phases. During the first phase, solely the GA ERASMUS is used to perform the solution search. Due to the strong non-linearity of problem [START_REF] Jayendiran | Micromechanical modeling and experimental characterization on viscoelastic behavior of 13 active composites[END_REF], the aim of the genetic calculation is to provide a potential sub-optimal point in the design space, which constitutes the initial guess for the subsequent phase, i.e. the local optimisation, where the MATLAB fmincon tool is employed to finalise the solution search. The optimisation algorithm is the active-set which is a Quasi-Newton method making use of an approximation of the Hessian matrix to estimate the descent direction. For more details on the active-set algorithm see [START_REF]Apple Ill Drive[END_REF].

For the resolution of the multi-scale inverse problem, both optimisation algorithms have been interfaced with the FE models of the multilayer plate at two different scales: microscopic (constitutive phases-level) and macroscopic (laminate-level). As shown in Figure 2, for each individual at each generation, the optimisation tool performs three different type of FE analyses:

1. an homogenisation analysis to determine the frequency-dependent equivalent elastic properties of the lamina (microscopic / macroscopic scale transition); 2. a non-linear modal analysis (by means of a suitable in-house coded solver) to extract the n f natural frequencies; 3. a non-linear harmonic analysis for the evaluation of the FRF of the laminate.

Then, the optimization algorithm elaborates the results provided by the two FE analyses in order to execute the optimization operations on the basis of the current value of both objective and constraint functions. These operations are repeated until the algorithm satisfies the user-defined convergence criterion. The details of the FE analyses are given in Section 4.

The number of design variables and that of constraint functions is nine and n f + 7, respectively. The generic individual of the GA ERASMUS represents a potential solution for the problem at hand. The genotype of the individual for problem [START_REF] Jayendiran | Micromechanical modeling and experimental characterization on viscoelastic behavior of 13 active composites[END_REF] is characterised by only one chromosome composed of nine genes, each one coding a component of the vector of design variables, see Eq. ( 7).

Finite element models at different scales

The finite element model at the microscopic scale and the homogenisation strategy

The microscopic / macroscopic scale transition is carried out through a homogenisation step performed on the RVE shown in Figure 3a.

The frequency-dependent elastic properties of the ply are obtained by means of the strain energy homogenisation technique (SEHT) of periodic media [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF][START_REF] Barbero | Finite element analysis of composite materials[END_REF]. This technique, originally introduced for elastic heterogeneous materials, can be generalised to different kinds of composites showing a general non-linear behaviour, e.g. fabrics, lattice structures, etc. This technique has already been successfully used in other works, see [START_REF] Montemurro | A multi-scale approach for the simultaneous shape and material optimisation of sandwich panels with cellular core[END_REF][START_REF] Montemurro | Identification of electromechanical properties of piezoelectric structures through evolutionary optimisation techniques[END_REF][START_REF] Catapano | A multi-scale approach for the optimum design of sandwich plates with honeycomb core. Part I: homogenisation of core properties[END_REF][START_REF] Catapano | A multi-scale approach for the optimum design of sandwich plates with honeycomb core. Part II: the optimisation strategy[END_REF]: in this work, the SEHT is generalised to the case of viscoelastic materials subjected to harmonic loads. The SEHT is based on the main hypothesis that the RVE of the periodic heterogeneous material and the corresponding homogenised volume undergo the same deformation having, therefore, the same strain energy. Consequently, at the ply scale, an equivalent homogeneous anisotropic material replaces the heterogeneous medium, by using the frequency-dependent stiffness tensor resulting from the homogenisation phase.

In this study, the real random micro-structure of the lamina (which is usually characterised by misalignments of the fibres, porosity, damaged zones, etc.) is not taken into account and the topology of the RVE is described by a perfect hexagonal array, as illustrated in Figure 3a. The FE model of the RVE has been realised into the commercial FE code ANSYS R . A 20-nodes solid element (SOLID186) with three DOFs per node has been used. The model together with its structured mesh is shown in Figure 3b. A sensitivity study (not reported here for the sake of brevity) on the proposed FE model with respect to the mesh size has been conducted: it was observed that a mesh having 19551 degrees of freedom (DOFs) is sufficient to properly evaluate the set of frequency-dependent homogenised elastic properties of the lamina.

The RVE is submitted to an average strain field εij (tensor notation), to evaluate the stiffness matrix components Cij . The six components of the average strain tensor are applied, by means of the classical periodic boundary conditions (PBCs) as follows [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF][START_REF] Barbero | Finite element analysis of composite materials[END_REF]:

u i (a 1 , x 2 , x 3 ) -u i (-a 1 , x 2 , x 3 ) = 2a 1 εi1 , -a 2 ≤ x 2 ≤ a 2 , -a 3 ≤ x 3 ≤ a 3 , u i (x 1 , a 2 , x 3 ) -u i (x 1 , -a 2 , x 3 ) = 2a 2 εi2 , -a 1 ≤ x 1 ≤ a 1 , -a 3 ≤ x 3 ≤ a 3 , u i (x 1 , x 2 , a 3 ) -u i (x 1 , x 2 , -a 3 ) = 2a 3 εi3 , -a 1 ≤ x 1 ≤ a 1 , -a 2 ≤ x 2 ≤ a 2 ,
where i = 1, 2, 3.

(12)

As stated above, the RVE is subjected to harmonic excitations, in order to compute the frequency-dependent elastic properties at the upper scale. Consequently, the internal stresses σ α and strains ε α , α = 1, ..., 6 (Voigt's notation) are varying harmonically with different amplitudes and phases, for each RVE internal point x p :

σ α (f 0 , x p , t) = |σ α (f 0 , x p ) | exp [(2πf 0 t + ϕ σα (f 0 , x p )) i] , (13)
ε α (f 0 , x p , t) = |ε α (f 0 , x p ) | exp [(2πf 0 t + ϕ εα (f 0 , x p )) i] . (14)
The equivalent stresses and strains at the lamina-level can be evaluated from the corresponding fields by considering an average over the RVE volume, i.e.

σα (f 0 , t) = σ α (f 0 , x p , t) = exp [(2πf 0 t) i] V RVE V RVE |σ α (f 0 , x p ) | exp [ϕ σα (f 0 , x p ) i] dV, (15)
εα (f 0 , t) = ε α (f 0 , x p , t) = exp [(2πf 0 t) i] V RVE V RVE |ε α (f 0 , x p ) | exp [ϕ εα (f 0 , x p ) i] dV. ( 16
)
where V RVE = 8a 1 a 2 a 3 according to Figure 3. The internal RVE stresses and strains of Eq. ( 13) and Eq. ( 14)can also be written in the Laplace-Carson (L is the related operator) space:

σ * α (f 0 , x p , f ) = L [σ α (f 0 , x p , t)] = |σ α (f 0 , x p ) | exp [ϕ σα (f 0 , x p ) i] 2π (f -f 0 ) i , (17)
ε * α (f 0 , x p , f ) = L [ε α (f 0 , x p , t)] = |ε α (f 0 , x p ) | exp [ϕ εα (f 0 , x p ) i] 2π (f -f 0 ) i . ( 18
)
By following the same logical steps, also the average stresses and strains components of Eqs. ( 15) and ( 16) can be written in the Laplace-Carson L space:

σ * α (f 0 , f ) = L [σ α (f 0 , t)] = 1 V RVE 2π (f -f 0 ) i V RVE |σ α (f 0 , x p ) | exp [ϕ σα (f 0 , x p ) i] dV, (19) ε
* α (f 0 , f ) = L [ε α (f 0 , t)] = 1 V RVE 2π (f -f 0 ) i V RVE |ε α (f 0 , x p ) | exp [ϕ εα (f 0 , x p ) i] dV. (20)
In order to perform the numerical homogenisation process, harmonic analyses are required and the imposed strain field, applied through the PBCs of Eq. ( 12), assumes the following expression:

ε * α (f 0 , f ) = L [ε α (f 0 , t)] = |ε α (f 0 ) | V 2π (f -f 0 ) i , ∀ α = 1, ..., 6, with ε * β (f 0 , f ) = 0∀ β = 1, ..., 6 and β = α. (21)
Finally, the frequency-dependent components of the homogenised stiffness complex tensor Cαβ of the lamina can be evaluated as:

Cαβ (f 0 ) = σ * α (f 0 , f ) ε * β (f 0 , f ) = 1 V RVE |ε β (f 0 ) | V RVE |σ α (f 0 , x p ) | exp [ϕ σα (f 0 , x p ) i] dV, ( 22
)
ℜ Cαβ (f 0 ) = 1 V RVE |ε β (f 0 ) | V RVE |σ α (f 0 , x p ) | cos [ϕ σα (f 0 , x p )] dV, (23)
ℑ Cαβ (f 0 ) = 1 V RVE |ε β (f 0 ) | V RVE |σ α (f 0 , x p ) | sin [ϕ σα (f 0 , x p )] dV. (24)
Even in the case of complex stiffness tensor, the compliance matrix at a given frequency f 0 can be determined as: S (f 0 ) = C-1 (f 0 ). Finally, the frequency-dependent lamina (complex) elastic properties can be computed from the components of the compliance matrix [START_REF] Barbero | Finite element analysis of composite materials[END_REF].

To give an idea of the homogenisation of the frequency-dependent elastic properties of the lamina at the mesoscopic scale, an analysis is performed by considering the material properties of the constitutive phases listed in Table 4, in the frequency range f ∈ [100, 6000] Hz. Figure 4 illustrates the viscoelastic behaviour of the matrix, represented through the Bagley-Torvik model: the storage and the loss moduli are respectively the real and the imaginary part of the matrix Young's modulus. The trend of the engineering moduli of the lamina vs. the frequency is give in Figures 5 and6.

Fibre properties

E f 1 [MPa] E f 2 [MPa] ν f 12 ν f 23 G f 12 [MPa] 275622 20435 0.32 0.451 10693 Matrix properties E m 0 [MPa] E m 1 [MPa] b m α m = β m ν m 3000 30 0.0053 0.5 0.33
As it can be easily inferred from these figures, the lamina Poisson's ratios, namely ν 12 , ν 13 and ν 23 , can be considered constant with the frequency. This result is of paramount importance to reduce the required computational effort for the multi-scale identification process.

The finite element model at the macroscopic scale

The FE model of the multilayer plate is built within ANSYS R environment [START_REF] Ansys | ANSYS Mechanical APDL Basic Analysis Guide[END_REF] by using SHELL281 layered shell elements with eight nodes and six DOFs per node: the plate kinematics is described by the first-order shear deformation theory (FSDT) [START_REF] Jones | Mechanics of composite materials[END_REF].

The laminate together with the applied excitation load and boundary conditions (BCs) is illustrated in Figure 7: the model is characterised by 1974 DOFs and the mesh size has been chosen after a convergence study (not reported here for the sake of brevity).

The choice of shell elements is due to the aspect ratio of the considered laminate (AR = 22.14), which is in the range [START_REF] Barkanov | Characterisation of viscoelastic layers in sandwich panels via an inverse technique[END_REF]100] where the FSDT gives satisfactory results.

As far as the computation of both the constraint functions of Eqs. ( 8), ( 9) and the objective function of Eq. ( 10) two macroscopic FE analyses are run for each point of the design space. Firstly, a non-linear modal analysis is performed, to extract the first n f natural frequencies and, secondly, a non-linear harmonic analysis is carried out to calculate the harmonic response of the plate for each sampling harmonic frequency of the chosen spectrum (Section 5.1). The spectrum harmonic response is obtained by measuring the displacement u z in the generic node of the macroscopic FE model mesh, at every sampled frequency f r , as shown in Figure 8. The harmonic response can be obtained by solving the following problem [START_REF] Hamdaoui | Variability analysis of frequency dependent visco-elastic three-layered beams[END_REF][START_REF] Daya | A numerical method for nonlinear eigenvalue problems application to vibrations of viscoelastic structures[END_REF][START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. part ii: Homotopy and complex nonlinear eigenvalue method[END_REF]:

[K sys (Ω) -λM sys ] {U sys } = f ext , (25)
where K sys (Ω), M sys , {U sys } and f ext are the stiffness matrix (which depends upon the pulsation Ω = ℜ λ

1 2
due to the viscoelastic matrix behaviour), the mass matrix, the nodal displacements vector and the external nodal forces, respectively.

The harmonic response, for each sample point, is obtained by evaluating the ratio of the FFT of the displacement along the z-axis u z q (f r ) to the nodal force along the same direction F z (f r ):

H r,q = u z q (f r ) F z (f r ) . (26)
It is noteworthy that the problem of determining the structure natural frequencies becomes non-linear due to the viscoelastic behaviour of the ply. The following non-linear eigenvalue problem must be faced:

[K sys (Ω) -λM sys ] {U sys } = {0} , det [K sys (Ω) -λM sys ] = 0, (λ, {U sys }) ∈ C × C n . (27)
Unfortunately, the non-linear eigenvalue problem of Eq. ( 27) cannot be solved by means of FE codes because it requires a dedicated algorithm / solver. Some research works are explicitely devoted to the implementation of a suitable algorithm for solving non-linear eigenvalue problems. As discussed in [START_REF] Hamdaoui | Comparison of non-linear eigensolvers for modal analysis of frequency dependent laminated visco-elastic sandwich plates[END_REF], several algorithms are available in the literature: the asymptotic numerical method (ANM) [START_REF] Daya | A numerical method for nonlinear eigenvalue problems application to vibrations of viscoelastic structures[END_REF], the inverse iteration algorithm (IIA) [START_REF] Schreiber | Nonlinear Eigenvalue Problems: Newton-type Methods and Nonlinear Rayleigh Functionals[END_REF], the iterative shift-inverter method (ISIM) [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF], the non-linear Jacobi-Davidson method (NLJDM) [START_REF] Voss | A Jacobi-Davidson method for nonlinear and nonsymmetric eigenproblems[END_REF] and the non-linear Arnoldi's method (NLAM) [START_REF] Voss | An Arnoldi Method for Nonlinear Eigenvalue Problems[END_REF]. Each method is characterised by its own advantages and drawbacks: in this work, the NLAM has been implemented into the MATLAB environment and interfaced with the FE model of the multilayer plate implemented into the ANSYS software.

In particular, the stiffness matrix of the FE model is recovered from ANSYS R and exported (and elaborated) into the MATLAB R software. The frequency-dependent stiffness matrix of the FE model can be expressed as:

K sys (Ω) = ℜ [K sys (Ω)] + iℑ [K sys (Ω)] . (28)
The stiffness matrix of the structure depends upon the frequency because of the viscoelastic behaviour of the ply (due to the matrix), see Eq. ( 22). However, this means that the numerical homogenisation method, discussed in Section 4.1, must be performed for each sampled frequency in the considered range. This would require a strong computational effort (and simulation time) which is not accetable for optimisation purposes. This issue can be easily overcome by looking at the volume-averged stress tensor of Eq. ( 19) which can be expressed in the following form:

σ (Ω) = 1 V RVE V f V RVE σ f (x p ) dV + 1 V RVE (1-V f )V RVE σ m (x p ) dV. (29)
Since the Poisson's ratio of the matrix ν m does not depend upon the frequency, the previous expression can be rewritten as:

σ (Ω, ε) = 1 V RVE V f V RVE σ f (x p , ε) dV + ... + Êm (Ω) V RVE (1-V f )V RVE 1 1 + ν m ε m (x p , ε) + ν m 1 -2ν m tr [ε m (x p , ε)] dV, (30)
therefore, the previous equation can be rearranged in a more compact form as follows

σ (Ω, ε) = M ( ε) + Êm (Ω) R ( ε) , M ( ε) = 1 V RVE V f V RVE σ f (x p , ε) dV, R ( ε) = 1 V RVE (1-V f )V RVE 1 1 + ν m ε m (x p , ε) + ν m 1 -2ν m tr [ε m (x p , ε)] dV. (31)
As a consequence of Eq.( 32), the equivalent frequency-dependent stiffness matrix of the homogeneous anisotropic material of the ply can be evaluated as (Voigt's notation):

Cjk = σj εk = M jk + Êm (Ω) R jk , where j, k = 1, ..., 6. (32)
It is straightforward to verify that the global stiffness matrix of the FE model at the macroscopic scale can be decomposed as follows:

K sys (Ω) = K sys 0 + Êm (Ω) Ksys , (33)
where K sys 0 is the term related to the constant part of the ply stiffness tensor, whilst Ksys is the contribution related to the frequency-dependent part.

Eqs. ( 32) and ( 33) allow obtaining the frequency-dependent stiffness matrices of the material and of the structure, respectively, directly within the ANSYS R software by means of only two homogenisation analyses carried out at two arbitrary frequencies (in this case the lower and the upper bounds of the considered frequency spectrum), instead of performing an homogenisation calculation (recall that each homogenisation corresponds to six FE analyses) for each sampled frequency in the considered range. Accordingly, the computational costs of the whole optimisation process is significantly reduced (the FE analyses at both microscopic and macroscopic scales must be carried out for each point in the design space).

Numerical results

Harmonic response for the reference configuration

Before running the optimisation procedure, the reference response must be calculated. To this purpose, firstly the numerical harmonic homogenisation process is performed on a RVE characterised by the material properties listed in Table 2 in order to obtain the reference ply material properties.

Secondly, the reference viscoelastic behaviour of the ply is implemented into the FE model of the multilayer plate in which two both a non-linear modal analysis and a harmonic analysis are performed to calculate the reference damped natural frequencies and the reference harmonic response. The reference damped eigenfrequencies are listed in Table 5.

For modal and harmonic analyses the frequency samples vary between f LB = 100 Hz and f U B = 7500 Hz: for the reference solution n f = 5 damped natural frequencies fall into this interval. The sampling of the considered spectrum is made according to the sequence reported in Table 6.

Nat. freq. Value

The value of δ i is computed according to the following formula:

δ i (f in ) = -4.34 × 10 -12 (f in ) 2 + 2.6 × 10 -7 (f in ) + 6.53 × 10 -4 , (34)
which has been chosen in order to have a value of δ i (f in ) that increases with the frequency. Indeed, by looking at the viscoelastic effect on the amplitude of the FRF

Frequency intervals [Hz]

N. of sampled spectrum points

1 - f ref 1n -f LB f ref 1n f 1n < f < f 1n -δ 1 11 f 1n + δ 1 < f < f 2n -δ 2 11 f 2n + δ 2 < f < f 3n -δ 3 11 f 3n + δ 3 < f < f 4n -δ 4 11 f 4n + δ 4 < f < f 5n -δ 5 11
f 5n + δ 5 < f < 1 - f U B -f ref 5n f ref 5n f 8n 11
f 1n -δ 1 < f < f 1n + δ 1 6 f 2n -δ 2 < f < f 2n + δ 2 6 f 3n -δ 3 < f < f 3n + δ 3 6 f 4n -δ 4 < f < f 4n + δ 4 6 f 5n -δ 5 < f < f 5n + δ 5 6
Table 6: Sampling sequence for FRF calculation.

(evaluated at the fours sampled points highlighted in Figure 8), as shown in Figure 9, it is possible to observe that the damping effect is more pronounced at high frequencies.

The exciting nodal force (as illustrated in Figure 7) has a value |F z | = 1 N and it is not applied at the plate center, in order to be able to excite even and odd modes. Finally, the number of sampling frequencies are N s = 87, for each plate point, whose number is N p = 329.

Results of the inverse problem

The optimisation process has been performed by selecting the main optimisation parameters tuning the behaviour of the GA ERASMUS as a result of a statistic analysis to evaluate their effects on the optimum solutions, according to the best practices discussed in [START_REF] Montemurro | Identification of electromechanical properties of piezoelectric structures through evolutionary optimisation techniques[END_REF]. The parameters governing both the GA and the deterministic algorithm are listed in Tables 7 and8, respectively.

The GA ERASMUS is run with two populations, each one composed of 50 individuals evolving along 100 generations. The exchange of information among populations is realised by using a ring-type operator every 10 generations: the probability of success of the ring-type operator is automatically computed by the GA. As far as the constrainthandling technique is concerned, the Automatic Dynamic Penalisation (ADP) method is used, see [START_REF] Montemurro | The automatic dynamic penalisation method (ADP) for handling constraints with genetic algorithms[END_REF].

The choice of using multiple populations with a small number of individual, is due to the fact that the main goal is to find the global minimum without increasing too much the computational time. In this way, the GA has the possibility to explore the design domain by exchanging information between best individuals belonging to different populations. More details about the use of multiple populations can be found in [START_REF] Montemurro | A contribution to the development of design strategies for the optimisation of lightweight structures[END_REF].

The inverse problem is solved by considering a fibre volume fraction V F = 0.555 [START_REF]Multi-scale modelling of composite material systems. The art of predictive damage modelling[END_REF] and a fibre diameter equal to d f = 6.8 µm [56]. The RVE dimensions are obtained as follows:

a 3 = d f 4 2π V f , a 2 = a 3 , a 1 = a 2 /4. (35)
For each point in the design space, the FE analyses constituted by the union of a numerical homogenisation analysis the solution of the non-linear eigenvalues problem and the computation of the harmonic responses, need about 104 s to be executed (on an Intel R Xeon R 2.70 GHz CPU with two processors and with a RAM of 128 GB). This analysis must be performed for each individual at each iteration, which implies an overall time of about 26.3 days, to get an optimum solution. The optimum solutions obtained from the genetic calculation and the local gradientbased optimisation in terms of microscopic material properties for both fibre and matrix are listed in Table 9, while the relative eigenfrequencies values are summarised in Table 10. As it can be easily inferred from the analysis of these results, the microscopic viscoelastic properties of the optimum solution are in good agreement with the reference data: the absolute percentage difference ranges from 0.10% for E f 1 to 11.94% for E f 2 . Only the material parameters E m 1 and G f 12 by significant percentage errors which are equal to 15.49% and 19.73%, respectively. This is a quite expected result because of the kinematic model at the basis of ANSYS shell elements FSDT. In fact the effect of these material parameters on the dynamic response at the macroscopic scale is negligible and this is also due to the particular stacking sequence (orientation angles and thickness) used for the reference structure: the considered plate is not thick enough to observe a significant influence of G f

Micro

Nat. freq. f ref in [Hz] f in [Hz] GA results f in [Hz] Gradient-
12 and E m 1 on its dynamic response. Nevertheless, both the damped eigenfrequencies and the FRF at each sample point are in excellent agreement with the reference values and the numerical results found at the end of the optimisation perfectly match the reference data with an absolute percentage difference ranging from 3.58 × 10 -2 % (for the first mode) to 2.19 × 10 -1 % (for the fourth mode).

Conclusions and perspectives

In this paper, an extension of the Multi-Scale Identification Strategy (MSIS) (initially presented in [START_REF] Cappelli | Characterisation of composite elastic properties by means of a multi-scale two-level inverse approach[END_REF]) is proposed. The MSIS is here applied to characterise the viscoelastic behaviour of the matrix and the elastic behaviour of the fibres by exploiting the information included into the dynamic response of the composite at the macroscopic scale. The proposed MSIS shows several features that make it a general methodology, that can be easily applied for different classes of materials and structures, e.g. multilayer, fabrics, etc, for identification of material properties without performing destructive tests.

In this study, the multi-scale inverse problem has been solved by means of a "one-shot" hybrid optimisation strategy. The multi-scale inverse problem is stated as an equivalent constrained non-linear programming problem (CNLPP) aiming at minimising the distance between the numerical and reference harmonic responses for the considered multilayer composite plate.

The scales transition is ensured by means of the strain energy homogenisation method for periodic media, which has been generalised to the viscoelastic case. In this way, the ply viscoelastic properties can be computed and used to build the FE model of the multilayer plate. At the microscopic scale the matrix viscoelastic behaviour is described through the Bagley-Torvik model, that requires only four material parameters.

Moreover, the modal and harmonic analyses performed on the multilayer plate at the macroscopic scale are non-linear due to the viscoelastic behaviour of the ply. The main issue is related to the non-linear modal analysis: no dedicated solvers are available in commercial FE software. To this purpose, the Arnoldi's method [START_REF] Voss | An Arnoldi Method for Nonlinear Eigenvalue Problems[END_REF] for non-linear eigenvalue problems has been coded into the MATLAB R environment and interfaced with the ANSYS code.

The effectiveness of the proposed strategy is evaluated through a numerical benchmark in which a composite laminate made of unidirectional carbon/epoxy pre-preg plies T650/F584 is considered as a reference structure.

The results provided by the MSIS are quite satisfactory: all viscoelastic properties are identified with a good level of accuracy, except the in-plane shear modulus of the fibre G f 12 and the viscoelastic matrix parameter E m 1 which are affected by an absolute percentage error of 19% and 15%, respectively. These errors are mainly due to the very low sensitivity of the objective function to these parameters. On the one hand, this low sensitivity is due to the geometry of the considered laminate which is not thick enough to highlight the influence of these properties on its dynamical response. On the other hand, the laminate stacking sequence plays a fundamental role: the stack considered in this work is a standard symmetric balanced stack taken from the literature which has not been designed to maximise the influence of some material properties on the laminate dynamic behaviour.

The proposed strategy constitutes just a "first attempt": the MSIS needs to be generalised to catch the true behaviour of the material of the constitutive phases. In order to achieve this ambitious goal, research is ongoing in order to include into the MSIS the following aspects:

• validation of the effectiveness of the proposed MSIS to characterise the viscoelastic behaviour of composite materials by exploiting the data resulting from experimental harmonic tests;

• design of a suitable stack to maximise the sensitivity of the objective function Φ(x) to the full set of the material properties to be identified;

• extension of the MSIS to the characterisation of the variability related to some parameters like the fibre volume fraction, misalignments of fibres, variation of the plies orientation angles, etc. (as partially done in [START_REF] Cappelli | Multi-scale identification of the elastic properties variability for composite materials through a hybrid optimisation strategy[END_REF]);

• application of the proposed strategy to different macroscopic specimen geometries and microscopic RVE topologies (in terms of constituent phases configurations).

As far as the experimental validation of the MSIS is concerned, two major difficulties must be faced before programming a campaign of harmonic/modal tests. Firstly, the variability of the material properties should be included into the inverse problem formulation. This requires the development of a suitable numerical model to properly describe the variability related to the viscoelastic behaviour of the composite, at each pertinent scale. Even if some models are available in the literature to properly describe the uncertainty of the elastic properties of both fibre and matrix, few research studies focuses on the modelling of the variability of their damping behaviour (to the best of the authors' knowledge). Therefore, a preliminary numerical/theoretical study should be conducted in order to find/develop a pertinent model to describe the variability characterising the viscoelastic behaviour, at each pertinent scale of the composite. Secondly, experimental results are unavoidably affected by noise. In the literature, one can find several methods/techniques to take into account the influence of noise on the characterisation of the elastic properties of the composite (very often at the mesoscopic scale). However, to the best of the authors' knowledge, the influence of noise on the identification of the parameters governing the viscoelastic behaviour of the microscopic constituents of the composite has not fully investigated yet. This aspect is also of paramount importance and should be addressed before starting an experimental campaign which aims at validating the proposed MSIS.

Due to its versatility, the MSIS can be used to characterise the geometrical parameters of the composite material RVE. The variables defining the shape of the inclusion or its volume fraction can be easily integrated into the vector of optimisation variables, without altering the overall architecture of the MSIS. Furthermore, laminate parameters can be included among the unknowns to be identified, e.g. the orientation angle and the thickness of each ply. Research is ongoing on these aspects as well.

Figure 1 :

 1 Figure 1: Geometrical parameters of the reference composite plate (sizes in mm).

Figure 2 :

 2 Figure 2: Flow-chart of the one-shot MSIS.

Figure 3 :

 3 Figure 3: (a) The reference Representative Volume element (RVE) and (b) details of the RVE mesh.

Figure 4 :

 4 Figure 4: Matrix Young's modulus vs. frequency (amplitude, real and imaginary parts).

Figure 5 :

 5 Figure 5: Frequency-dependent elastic properties (amplitude) of the lamina resulting from the homogenisation.

Figure 6 :

 6 Figure 6: Frequency-dependent elastic properties (phase) of the lamina resulting from the homogenisation.

Figure 7 :

 7 Figure 7: FE model of the multilayer plate.

Figure 8 :

 8 Figure 8: Location of the sample points over the plate for harmonic displacements evaluation.

Figure 9 :

 9 Figure 9: Amplitude and phase of the FRF for the reference solution at the four sampling points highlighted in Figure 8.

Table 1 :

 1 Models, techniques and methods used in this study.

	/ Technique / Method Role
	The Bagley-Torvik model	Used to define the viscoelastic behaviour
		of the matrix
	Strain energy homogenisation	Used	to	perform	the	micro-
	technique	scopic/mesoscopic scale transition
	Finite Element Method	Used to model both the RVE of the com-
		posite and the multilayer plate at mi-
		croscopic and macroscopic scales, respec-
		tively
	Non-Linear Arnoldi's method Used to solve the non-linear eigenvalue
		problem at the macroscopic scale
	GA ERASMUS	Used to perform the global search of op-
		timal solutions for the multi-scale identi-
		fication problem
	Automatic Dynamic Penalisa-	Used to handle optimisation constraints
	tion Method	during the global search stage
	fmincon tool and active-set	Used to perform the local search of opti-
	algorithm	mal solutions for the multi-scale identifi-
		cation problem

Table 3 :

 3 Lower and upper bounds of the design variables for the multi-scale inverse problem.

Table 4 :

 4 Material properties of the fibre and the matrix used to illustrate the effectiveness of the homogenisation procedure.

Table 5 :

 5 Reference natural frequencies.

Table 7 :

 7 Optimisation parameters for the GA ERASMUS.

	Parameters	Value
	Solver algorithm	active-set
	Max function evaluation	10000
	Tolerance on the objective function	10 -15
	Tolerance on the gradient norm	10 -15

Table 8 :

 8 Optimisation parameters for the gradient-based algorithm.

Table 9 :

 9 Optimum solution of the inverse problem: results provided by the GA and the active-set algorithm. The percentage difference is indicated in parentheses.

Table 10 :

 10 The damped eigenfrequencies for the optimum solution of the inverse problem: results provided by the GA and the active-set algorithm. The percentage difference (with respect to the reference values) is indicated in parentheses.

Acknowledgements

This research work has been carried out within the project FULLCOMP (FULLy analysis, design, manufacturing, and health monitoring of COMPosite structures), funded by the European Union Horizon 2020 Research and Innovation program under the Marie Skodowska-Curie grant agreement No. 642121.