N Mathieu
email: norman.mathieu@univ-lorraine.fr

C Czarnota

H Obeid

S Mercier

Prediction of flatness defects and of the stable configuration of thin multilayer assemblies due to chemical shrinkage

Keywords: Effective shrinkage strain, Flatness defects, Variational principle, Analytical modeling, FE simulation

The manufacturing process of multimaterial and multilayer assemblies involving pre-impregnated laminates consist of heating the composite structure at high temperature, typically of the order of 200 • C, at which polymerization occurs. During the curing, a permanent deformation, called chemical shrinkage strain, is generated and may strongly influence the future flatness of the assembly. Moreover, the cooling generates additional thermal deformations, which also participate into the manifestation of flatness defects at room temperature. To predict warpage or flatness defects, the chemical shrinkage strain needs to be precisely determined.

This work proposes an analytical approach dedicated to flatness prediction of multilayer composites taking into account shrinkage strain generated during processing. Our contribution also aims at predicting and analyzing stable and unstable solutions of flatness defects. The proposed analytical model, developed for any multilayer composite, is obtained from an extension of the classical laminate theory (CLT). Geometrical nonlinearities are also accounted for. The analytical approach relies on trial fields for strain and displacements, and on total potential energy minimization. The theory is applied to a bilayer laminate consisting of a cured layer made of epoxy/glass fiber composite and of a pre-impregnated one of the same material. Results obtained from the analytical modeling are validated by numerical simulations. Influence of material parameters is also analyzed for this configuration. Finally, from experimental measurements of curvatures on bilayer composite samples, an inverse method and a minimization procedure, the proposed analytical development provides an estimate of the effective shrinkage strain, which is responsible for the flatness defect. Illustration of this strategy is exemplified by considering a bilayer composite manufactured for this work.

Introduction

Because of differences in thermo-mechanical properties between layers, a multimaterial and multilayer assembly tends to bend in order to equilibrate residual stresses induced for instance during the manufacturing process. When the multilayer composite is made of some pre-impregnated laminates (i.e. partially cured matrix), the heating at sufficiently high temperature induces a permanent deformation (chemical 1 and thermal shrinkages) which needs to be accomodated by the elastic strains. A flatness defect thus emerges in the composite structure. The goal for manufacturing industries is to minimize the warpage development after processing but also in service operation. This optimization is strongly related to the chemical and thermal shrinkages. The combined effect of both contributions will be examined in this paper by referring to the effective shrinkage strain as the permanent strain resulting from the curing step. This feature has to be precisely determined.

As an application of our analysis, attention is paid in this paper on structures with dimensions and constituent materials used in the Printed Circuit Board (PCB) industry. For microelectronics, predicting flatness defects is an important issue. A PCB consists of an assembly of conductive thin copper tracks, embedded in insulating layers (typically made of epoxy resin or polyimide composite reinforced with glass fibers). For some specific applications such as aeronautics, aerospace industry or military sector, Printed Circuit Boards are subjected to severe conditions during their manufacturing process and their service life.

The manufacturing process of PCBs has received particular attention since flatness defects usually emerge during this stage and pertain during subsequent thermal loadings. It may contain more than 70 steps, depending on the complexity of the multilayer structure (see [START_REF] Schuerink | Modeling printed circuit board curvature in relation to manufacturing process steps[END_REF]. One of the key step, which is illustrated in Fig. 1, consists in the assembly of materials of different natures under thermo-mechanical loading following a four-stage sequence: • 1 heating: a pressure P 0 (around 2MPa) is imposed and the sample is heated from room temperature (T room = 20 • C) to T max (around 200 • C);

• 2 curing: the temperature is maintained at T max to cure the resin totally. This step leads to an effective shrinkage strain (permanent deformation at the level of the ply) that is induced by the curing of the resin and the mechanical interaction with glassfibers;

• 3 cooling: the sample is cooled from T max to room temperature (leading to thermal shrinkage);

• 4 demoulding: the pressure is removed at room temperature.

Many authors used the Classical Laminate Theory (CLT ) to determine the final flatness of a multilayer and multimaterial plate (see [START_REF] Schuerink | Modeling printed circuit board curvature in relation to manufacturing process steps[END_REF]. [START_REF] Nakagawa | Optimum design of printed circuit board to reduce deformation in reflow process by a global optimization method[END_REF] combined the CLT with an optimization algorithm in order to find the optimum design of a PCB that minimizes its deformation after reflow process. To go further, [START_REF] Abouhamzeh | Closed form expression for residual stresses and warpage during cure of composite laminates[END_REF] proposed a modified CLT model to predict the warpage and residual stresses in fiber metal laminates composed of aluminum, epoxy resin and prepreg layers. In this work, a deformation due to chemical shrinkage was added in the CLT and the epoxy resin stiffness modulus was assumed to increase during curing. [START_REF] Kravchenko | Chemical and thermal shrinkage in thermosetting prepreg[END_REF] and [START_REF] Kravchenko | Cure history dependence of residual deformation in a thermosetting laminate[END_REF] distinguished two kinds of shrinkage, as also reported in the forming process of Fig. 1: chemical shrinkage which consists of volume reduction as the result of crosslinking and depends on the degree of the resin cure; thermal shrinkage which results from cooling. Nawab et al. (2013b) conducted a literature survey about all techniques used for the characterization of resin effective shrinkage. They found that the most important issue is to identify the chemical shrinkage of resin without being influenced by thermal expansion. Among experimental approaches, we may cite [START_REF] Kravchenko | Digital image correlation measurement of resin chemical and thermal shrinkage after gelation[END_REF] who used Digital Image Correlation to record in situ stress-free strain fields in a thermosetting polymer submitted to curing process. This technique allowed to obtain relations between the chemical shrinkage strain and the dependence of the coefficients of thermal expansion to the degree of cure. Another recent method proposed by [START_REF] Tsai | Cure strain monitoring in composite laminates with distributed optical sensor[END_REF] also relied on in situ measurements using Distributed Optical Sensors (DOS).

The aim was to measure the strain history in pre-impregnated thermoset composites during their curing process. By coupling experimental data with the laminate theory, the authors identified effective shrinkage values (-2500µε in the transverse direction against ∼ 0µε in the fiber direction). Note that 1D long fiber composite was considered. However, only partial determination was obtained since, owing to the DOS technique, the effective shrinkage (e.g. the combined chemical and thermal shrinkage strains) could only be captured when the sensor bonds to the epoxy resin (i.e. when the temperature overcomes the glass transition T g ).

It is also important to stress that for large deformation of the composite, the linear modeling issued from the CLT is no longer sufficient to evaluate the structure flatness after manufacturing process (see [START_REF] Cantera | Modelling of [0/90] laminates subject to thermal effects considering mechanical curvature and through-the-thickness strain[END_REF]. [START_REF] Gigliotti | Loss of bifurcation and multiple shapes of thin [0/90] unsymmetric composite plates subject to thermal stress[END_REF] proved, from a comparison between finite element calculations and experimental observations, that geometrical nonlinearities have to be accounted for to properly predict the deformed shape of plates. A so-called Extended Classical Lamination Theory (ECLT), based on the total potential energy minimization, was then developed to find the nonlinear deformation of a mutlilayer plate submitted to thermal loading (see the pioneering work of [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF]. Several authors followed this way for estimating the final plate curvature and discuss the stability of solutions for manufactured samples having different stacking sequences (e.g. see [START_REF] Eckstein | Morphing high-temperature composite plates utilizing thermal gradients[END_REF]Nawab et al., 2013a). Furthermore, multistability was found to occur for particular laminate designs. For instance, thermally actuated morphing structures were constructed by [START_REF] Eckstein | Multi-mode morphing using initially curved composite plates[END_REF]. They were considering initially curved composite laminates with temperature dependent properties. In their approach, thermal gradient through thickness was accounted for. They showed a complex multistability behavior which could not have been obtained with initially flat laminates.

A second route of flatness investigation has recoursed to Finite Element analyses to model the behavior of laminates. [START_REF] Li | Bistable hybrid symmetric laminates[END_REF] predicted unstable shapes for hybrid symmetric lay-up laminates. To predict stable shape, a buckling step followed by a simulation with initial imperfection (deduced from buckling modes) were implemented by the authors. [START_REF] Lee | Effect of anisotropic thermo-elastic properties of woven-fabric laminates on diagonal warpage of thin package substrates[END_REF] highlighted the importance of considering both orthotropic elastic constitutive law and anisotropic thermal expansion (with a non-zero shear coefficient) to predict thin package substrate deformation and diagonal warpage which usually emerges as an unstable solution in CLT approaches. [START_REF] Macurova | Stress and deflection development during die embedding into printed circuit boards[END_REF] predicted stresses and deflection of PCBs after embedding process. [START_REF] Groh | Extreme mechanics in laminated shells: New insights[END_REF] included Finite Element discretization into numerical continuation algorithms to study the multistability of [0 • /90 • ] laminate during its cool-down. Their approach was termed extreme mechanics in the sense that the use of arc-length techniques provides bifurcation detection and several branch-switchings showing possible snapping behavior between five different modes. [START_REF] Kuang | A computational framework for multi-stability analysis of laminated shells[END_REF] recently investigated the stability of cylindrical laminated shells. They were subjected to a central point load, with free and hinged edges. The stability of various configurations (symmetric and antisymmetric structures) was captured by an asymptotic numerical method.

In this paper, we restrict our attention to flatness defects generated during the manufacturing process described in Fig. 1. Layers made of the same material, only differ in the initial resin-state. Some layers are cured prior to the manufacturing process while others are pre-impregnated materials. For those plies, the polymerization at elevated temperature leads to permanent deformation. The observed effective shrinkage strain is the result of the interaction between the B-stage resin cure and the elastic fiberglass. Flatness defects that may develop in composite laminates due to heterogeneity in material properties (i.e. different coefficients of thermal expansion and/or elastic properties), are thus not considered in the present work but could be addressed in a future work. The paper is organized as follows. In section 2, the theory based on Extended Classical Lamination Theory is developed for thin multilayer composites submitted to thermal loading. Section 3 presents a parametric study where the effects of orthotropic elastic stiffness and of effective shrinkage strain tensors are examined in terms of curvatures and maximum displacement. The occurrence of stable and unstable solutions is observed. In section 4, the theoretical approach is used to identify effective shrinkage strains from experimental tests conducted on bilayer samples. The capability of the modeling to predict final stable shapes is discussed. Path-following algorithms were implemented to analyse the complex multistability of the bilayer laminate. One of the main outcome of the proposed analytical approach is a new straightforward experimental method to identify effective shrinkage strains without having recourse to sophisticated devices. Finally, the analysis developed in this paper furnishes a fundamental tool to study flatness defects and to investigate stable or unstable solutions. The proposed work may help PCB designers to anticipate defect development.

The proposed methodology can be transposed to any industrial domains where large composite structures are frequently used and when the warpage control is crucial.

Analytical modeling

Potential energy minimization

An analytical modeling is proposed to predict the shape of a thin multilayer composite structure under general thermal loading. Illustration will be conducted with the one of Fig. 1. Cured and pre-impregnated layers are present in the stack-up sequence. The approach is based on the minimization of the total potential energy developed in [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF] (see also [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF][START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF][START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF], and is extended by accounting for the effect of permanent effective shrinkage strain.

Consider a thin rectangular laminate (length L, width l and thickness h, with h/L < h/l << 1) in which the plane z = 0 coincides with the midplane of the plate. The cartesian coordinates of any material point in the plate are (x, y, z) with

-L 2 ≤ x ≤ L 2 , -l 2 ≤ y ≤ l 2 , -h 2 ≤ z ≤ h 2 .
The structure is composed of n layers. The reduced stiffness tensor, which varies from one layer to another, is denoted by Q and the following vector notation is introduced (see [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells: Theory and Analysis[END_REF], for notation details):

ε =    ε xx ε yy γ xy    , ε th =    ε th xx ε th yy γ th xy    and Q =    Q 11 Q 12 Q 16 Q 12 Q 22 Q 26 Q 16 Q 26 Q 66    (1)
with ε is the total strain. Due to cure of the resin in some plies and thermal loading, the effective thermal strain in our analysis ε th is decomposed in two parts:

ε th = ∆T α + ε sh (2)
where the first term on the right hand side is induced by the temperature variation ∆T = T -T 0 , T 0 being a reference temperature. The second one mimics the consequence of the volumetric contraction of the resin due to cure and the associated mechanical interaction with glassfibers (denoted as effective shrinkage strain ε sh ). α is the reduced tensor of thermal expansion coefficients. When the assembly is composed of layers that have all been cured prior to the forming process, ε sh vanishes and the modeling reduces to approaches found in the literature, see [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF]. Note also that if some layers have not been cured, the effective thermal strain tensor may vary from one layer to another, even if all layers exhibit identical coefficients of thermal expansion. The total potential energy of the composite structure under plane stress condition subjected to a thermal loading is written as:

Π = L/2 -L/2 l/2 -l/2 h/2 -h/2 1 2 ε • Q • ε -ε • Q • ε th dzdydx, (3)
with • being the contracted product.

We assume that strains are linear through the thickness and geometrical nonlinearities are accounted for so that, using vector notation, the total strain is written as:

ε (x, y, z) = ε 0 (x, y) + zκ (x, y) (4) with ε 0 =    ε 0xx ε 0yy γ 0xy    =      ∂u 0 ∂x + 1 2 ∂w 0 ∂x 2 ∂v 0 ∂y + 1 2 ∂w 0 ∂y 2 ∂u 0 ∂y + ∂v 0 ∂x + ∂w 0 ∂x ∂w 0 ∂y      and κ =    κ xx κ yy κ xy    =     -∂ 2 w 0 ∂x 2 -∂ 2 w 0 ∂y 2 -2 ∂ 2 w 0 ∂x∂y     . (5)
In Eq. ( 5), u 0 , v 0 , w 0 denote the midplane displacements, ε 0xx , ε 0yy , γ 0xy the membrane strains and κ xx , κ yy , κ xy are the curvatures. For ease of reading the dependance in x and y will be intentionally omitted in the following.

Combining Eqs (3-5), the effective potential energy is expressed in terms of membrane strains ε 0 and curvatures κ:

Π = L/2 -L/2 l/2 -l/2 h/2 -h/2 1 2 ε 0 • Q • ε 0 + 2zε 0 • Q • κ + z 2 κ • Q • κ - ε 0 • Q • ε th + zκ • Q • ε th dzdydx. ( 6
)
Each layer is assumed homogeneous. So the through-thickness integration is replaced by a summation over the number of layers (see elsewhere in any handbook for the CLT ). The current label of each individual layer is k (1 ≤ k ≤ n). The global extensional stiffness tensor A, the bending-extensional coupling stiffness tensor B and the bending stiffness tensor D are classically defined as:

A = n k=1 Q (k) (z k+1 -z k ) , B = n k=1 Q (k) z k+1 2 -z k 2 2 , D = n k=1 Q (k) z k+1 3 -z k 3 3 (7)
with Q (k) standing for the reduced stiffness tensor of the k th layer. z k and z k+1 are the positions of the bottom and top surfaces of layer k.

The thermal forces N th and thermal moments M th per unit surface are given by:

N th = n k=1 Q (k) • ε th (k) (z k+1 -z k ) and M th = n k=1 Q (k) • ε th (k) z k+1 2 -z k 2 2 (8)
with ε th (k) being the global thermal strain related to the k th layer, defined in Eq. (2).

Combining Eqs (6-8), the total potential energy can be written in the condensed form:

Π = L/2 -L/2 l/2 -l/2   1 2 ε 0 κ T A B B D ε 0 κ - ε 0 κ T N th M th   dydx. ( 9
)
For a given composite structure and material parameters, the unknown vectors ε 0 and κ are determined from the minimization of the total potential energy. Since the solution cannot be obtained in a closedform expression, a Rayleigh-Ritz method is adopted in the present work, see also [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF] or [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF]. More specifically, the dependency of ε 0 and κ upon the position (x, y) is approximated by a truncated series expansion and the unknown quantities ε 0 and κ are expressed in terms of m independent parameters p i (i = 1, . . . , m): ε 0 = ε 0 (x, y, p 1 , p 2 , ..., p m ) and κ = κ(x, y, p 1 , p 2 , ..., p m ).

(10)

Note that the order of the truncated series and the number of independent parameters are not linked.

Adopting the form (10) for the six unknown components (ε 0xx , ε 0yy , γ 0xy of ε 0 and κ xx , κ yy , κ xy of κ), the minimization of Π leads to the following necessary conditions:

∂Π ∂p i = 0 for i = 1, . . . , m (11)
which bring a system of m nonlinear equations that need to be satisfied simultaneously. Solving (11) provides the identification of parameters p i (i = 1, . . . , m) which in turns allows to determine each component of ε 0 and κ from Eq. ( 10). Note that Eq. ( 11) can define a critical point which may lead to a stable (minimum) or unstable (maximum) solution. A solution is said to be stable if the following condition on the Hessian matrix is fulfilled:

∂ 2 Π ∂p i ∂p j 1≤i,j≤m is definite positive. ( 12
)
The modeling developed here is similar to some approaches of the literature where minimization of the potential energy is involved and Eqs (9,11-12) can be found for instance in [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF]. However, to the best of our knowkledge, none of them has considered the effect of effective shrinkage strain (permanent deformation due to curing) on the occurrence of flatness defects and on the loss of stability.

Polynomial approximation of displacements and strains

Various approaches can be found in the literature where polynomial functions were used to expand either midplane displacements, or membrane strains and curvatures, or combined displacements and strains.

• In the approach of [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF], the midplane displacements u 0 , v 0 , w 0 were expressed in terms of four independent parameters:

u 0 = cx - a 2 x 3 6 - ab 4 xy 2 , v 0 = dy - b 2 y 3 6 - ab 4 x 2 y , w 0 = 1 2 ax 2 + by 2 . ( 13
)
Using Eq. ( 5), the membrane strains and curvatures are derived:

ε 0xx = c - ab 4 y 2 , ε 0yy = d - ab 4 x 2 , γ 0xy = 0 , κ xx = -a , κ yy = -b , κ xy = 0. ( 14
)
In [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF], four parameters are representing the six components of the mid-plane strain and of the curvature. With this simple model, the shear strain γ 0xy and the shear curvature κ xy are vanished.

• In order to introduce non-zero shear strain in the x -y plane, [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF] proposed a six parameter formulation for the midplane displacements where:

u 0 = cx - a 2 x 3 6 -exy 2 , v 0 = dy - b 2 y 3 6 -f x 2 y , w 0 = 1 2 ax 2 + by 2 . ( 15
)
Using Eq. ( 5), one obtains:

ε 0xx = c -ey 2 , ε 0yy = d -f x 2 , γ 0xy = -2 e + f - ab 2 xy , κ xx = -a , κ yy = -b , κ xy = 0. ( 16
)
The approach of [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF] reduces to the modeling proposed by [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF] when e = f = ab/4. Note also that with the proposed trial displacement fields, the predicted curvature κ xy is still zero. Clearly, some configuration of unsymmetric composites present a non vanishing shear curvature. Such extensions were proposed in the literature.

• From another point of view, [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF] proposed a quadratic expansion of the membrane strains ε 0xx , ε 0yy and of the vertical displacement w 0 :

ε 0xx = ε x 00 + ε x 20 x 2 + ε x 02 y 2 + ε x 11 xy , ε 0yy = ε y 00 + ε y 20 x 2 + ε y 02 y 2 + ε y 11 xy, ( 17
)
w 0 = w 20 x 2 + w 02 y 2 + w 11 xy. ( 18
)
Based on the trial fields (18) and the definition of ε 0 in Eq. ( 5), the in-plane displacements are found: [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF] considered that the in-plane displacements can be well reproduced by a polynomial expression of order three. The shear strain is obtained from Eq. (5-19-20). This model involves 14 independent parameters since the condition u 01 = v 10 is imposed to avoid rigid-body rotation. From the proposed kinematic description of the composite plate, the curvatures are constant and given by:

u 0 =   ε 0xx - 1 2 ∂w 0 ∂x 2   dx + u 01 y + u 03 y 3 , ( 19
)
v 0 =   ε 0yy - 1 2 ∂w 0 ∂y 2   dy + v 10 x + v 30 x 3 , (20)
κ xx = -2w 20 , κ yy = -2w 02 , κ xy = -2w 11 . ( 21
)
By contrast with the modeling of [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF], the modification introduced by [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF] leads to a shear curvature κ xy which is no longer null.

• [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF] complemented the work of [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF] by a further enrichment of the formulation of the out-of-plane displacement for the midplane surface w 0 :

w 0 = w 00 + w 10 x + w 01 y + w 20 x 2 + w 02 y 2 + w 11 xy + w 21 x 2 y + w 12 xy 2 + w 22 x 2 y 2 . ( 22
)
From a mathematical point of view, w 0 is now searched as a product of two parabolas with axes lying along the principal directions. The in-plane membrane strains are still given by Eq. ( 17). This approach involves 20 independent parameters. In practice, w 00 is set to 0 to reflect that the center of the plate is kept fixed during the loading. In this way, rigid-body translation no more exists. With such a formulation, the curvatures are varying in the plane according to:

κ xx = -2 w 20 + w 21 y + w 22 y 2 , κ yy = -2 w 02 + w 12 x + w 22 x 2 , κ xy = -2 (w 11 + 2w 12 y + 2w 21 x + 4w 22 xy) . ( 23
)
To compare curvature predictions given by [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF] and the three previous models adopted [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF], [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF] and [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF]), it is of interest to evaluate the average curvatures κ

= 1 L×l L/2 -L/2 l/2 -l/2 κ (x, y) dydx. From relationship (23), one obtains: κxx = -2w 20 - w 22 6 l 2 , κyy = -2w 02 - w 22 6 L 2 , κxy = -2w 11 . ( 24
)
The parameter labels in Eqs (13-24) are the original labels provided by the corresponding authors [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF], [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF], [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF] and [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF]. In the following, to be consistent with the notation of section 2.1, they will be identified as p i .

Note the progressive enrichment in the description of the deformation from [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF] to [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF]. While in [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF] the shear strain is set to zero, it is bilinear in [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF] and of a quadratic form in [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF] and [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF]. In addition, curvatures κ xx and κ yy are homogeneous in [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF], [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF] and [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF], whereas it varies in the plane when the modeling of [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF] is adopted. In particular, regarding curvature κ xy , the approach of [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF] allows to account for the development of inhomogeneous curvature that may occur for some specific cases (e.g. a moment applied to a plate edge generates local variation of curvature). It provides also the ability to consider general boundary conditions (not only free surfaces) or compound plates i.e. different stacking sequences joined side by side.

Obviously, the degree of complexity of the solving procedure increases as the number of parameters raises. If one can easily solve Eq. ( 11) up to 14 independent parameters, using tools available in numerical softwares (e.g. Matlab here or Mathematica in [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF], the 20 parameter model needs advanced numerical strategies to capture the solution. The latter is obtained through a procedure combining arc-length method and deflation technique. Appendix A provides details of various procedures we have implemented to identify parameters p i (i = 1, . . . , m) with m = 4, 6, 14 or 20. In addition, a validation of the theoretical approach and solving procedure is proposed in Appendix B where results calculated with the model of [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF] are compared against finite element simulations obtained using Abaqus/Standard.

In our work, the four sets of parameters involved in the polynomial interpolations have been considered to predict the plate deflection due to effective thermal strain and more specifically the role of shrinkage strain. The following section compares results obtained from the above cited approaches. They will be referred to as 4p-model, 6p-model, 14p-model and 20p-model for the model of [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF] (Eqs 13-14), [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF] (Eqs 15-16), [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF] (Eqs 17-18) and Mattioni et al. (2009) (Eqs 17,22-23) respectively.

Effect of the enrichment in the polynomial kinematic trial fields

The influence of shrinkage process on the final flatness of manufactured thin multilayer plates is explored in this section. The adopted configuration consists in a bilayer rectangular plate of initial length L = 470mm, width l = 300mm and thickness h = 0.6mm, as displayed in Fig. 2. The bottom and top layers, made of the same material, differ from their initial state, the first one being already cured prior to the manufacturing process.

The two layers have the same thickness (t = 0.3mm). The constituents of the bilayer structure have an orthotropic elastic behavior. The stacking of the assembly is [0 • /0 • ] (no misorientation between the two layers). As we consider the manufacturing process described in Fig. 1, the sole permanent deformation originates from the curing at T max . As the two layers have the same in-plane coefficients of thermal expansion, it appears that the four-stage process does not need to be accounted for and the analysis can be restricted to a single step of virtual cure where instantaneous shrinkage is admitted to occur at room temperature. This assumption has been verified with numerical simulations (results not shown here) where the simplified virtual cure approach has led to the same final state in terms of stresses, strains and shape as the one obtained considering the complete four-stage process. Since virtual cure will be considered, only values of Young's moduli, Poisson's ratios and shear moduli at room temperature are relevant for this study. Table 1 provides material parameters adopted in this work, which are representative of the configuration tested in section 5. Material characterization was performed on an electrodynamic tensile machine using a contactless strain measurement technique(DIC), see [START_REF] Girard | Experimental and numerical characterization of thin woven composites used in printed circuit boards for high frequency applications[END_REF] for more details.

The out-of-plane properties are not needed in the present CLT based approach but have been used in Finite Element simulations. They were obtained following the numerical homogenization procedure developed in [START_REF] Girard | Experimental and numerical characterization of thin woven composites used in printed circuit boards for high frequency applications[END_REF].

Tab. 1 -In-plane elastic properties of the fiberglass reinforced epoxy laminate at room temperature (20

• C). E x [GPa] E y [GPa] ν xy G xy [GP a] 20.20 22.71 0.160 4.75
The bilayer plate is thus submitted to a virtual cure at T = T 0 = T room so ∆T = 0 • C in Eq. ( 2) and the effective thermal strain reduces to the effective shrinkage strain. Since only the top layer is uncured prior to the loading, the effective thermal strain is expressed as:

ε th = ε sh =         0 0 0    for - h 2 ≤ z ≤ 0 ;    ε sh xx ε sh yy γ sh xy    for 0 ≤ z ≤ h 2      . ( 25
)
In this section, calculations are performed with ε sh xx = ε sh yy = ε sh , γ sh xy = 0 and ε sh is varying in the range [-10 -4 , 0]. Depending on the adopted modeling and prescribed value of ε sh , the bilayer composite can deform into one of the five shapes illustrated in Fig. 3 (the notations S1 to S5 are introduced).

Figs 3(a-b-c) display the deformed shape with κ xy = 0. The vertical displacement w 0 is maximum at the four corners. More specifically, Fig. 3(a d-e) reflect cases where the maximal vertical displacement is located at two opposite corners with κ xy > 0 in Fig. 3(d), κ xy < 0 in Fig. 3(e). Note that shapes S4 and S5 cannot be predicted by the 4p-model and the 6p-model for which κ xy = 0 by construction. Fig. 4 shows the evolution of curvatures κ xx , κ yy and κ xy as a function of |ε sh | calculated from the four approaches presented in section 2.2. For the 20p-model, average values given in Eq.( 24) are considered. Stable solutions (i.e. minimum of total potential energy) are represented with solid lines whereas unstable solutions (i.e. maximum) are displayed with dashed lines. It first appears from Fig. 4(a) to Fig. 4(g) that results are not depending on the adopted modeling route for values of |ε sh | up to 10 -5 . In that range, there is only one (stable) solution with κ xx /κ yy ≃ 1 so the deformed configuration corresponds to shape S1 as

L 2 x 0 -L 2 -l 2 0 y l 2 w 0 (a) S1 L 2 x 0 -L 2 -l 2 0 y l 2 w 0 (b) S2 L 2 x 0 -L 2 -l 2 0 y l 2 w 0 (c) S3 L 2 x 0 -L 2 -l 2 0 y l 2 w 0 (d) S4 L 2 x 0 -L 2 -l 2 0 y l 2 w 0 (e) S5
ε sh xx = ε sh yy = ε sh = -10 -4 , γ sh xy = 0. Material parameters are provided in Table 1. The plate dimension is L = 470mm, l = 300mm, h = 0.6mm. The 20p-model is used.
depicted in Fig. 3(a). When |ε sh | is raised, κ xx is shown to increase in magnitude (κ xx being negative), while κ xy = 0. As |ε sh | is increased, other shapes can emerge depending on the adopted modeling strategy. Let us first concentrate on the 4p-model, see Figs 4(a) and 4(b). Following the continuous evolution starting from the lowest value of |ε sh |, the curvature κ xx is negative and tends to increase in magnitude until |ε sh | reaches approximately 2 10 -5 . For larger values of |ε sh |, the curvature decreases in magnitude to asymptotically reach a zero value for large shrinkage strain. This reflects that the plate tends to adopt an alternative shape once |ε sh | is larger than a critical value. In fact, one observes that three different shapes are obtained for a value of |ε sh | around 3 10 -5 . Two shapes (S2 and S3) are stable solutions for which Eq. ( 12) is satisfied, while the shape (S1) is unstable. Note that this last solution was stable for smaller values of |ε sh |. With the 20p-model, more information on the curvature can be obtained as the displacement description includes higher order terms. The variation of curvature along the plate length and width is displayed in Fig. 5. Evaluation is provided by the 14p and 20p-models for the configuration of Fig. 3 with |ε sh | = 10 -4 . As previously shown in Fig. 4, stable solutions S4 and S5 obtained from the two models coincide when the average curvatures κxx , κyy and κxy are considered, see Eq. ( 24). This founds explanation in Figs 5(a) and 5(b) where the local curvatures for the stable solutions predicted by the 20p-model are shown to slightly vary within the plate. By contrast, unstable solutions S1, S2 and S3 exhibit heterogeneous curvatures, with a particular marked effect for the shape S3. In that case, κ xx ≃ -10 -4 mm -1 for the 14p-model, close to the value κ xx (y = 0) predicted by the 20p-model at the plate center. For the same case, the amplitude of the curvature is twice larger at the edges with κ xx (y = ± l 2 ) ≃ -2 10 -4 mm -1 . Similar observations are made for κ yy in Fig. 5(b) with still a strong heterogeneity (even slightly larger). Recall that the hypothesis of constant curvatures has been removed in the 20p-model developed by [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF] to account for other boundary conditions than free edges. Therefore, if the focus is only made on the stable solution for the deformed plate, 20p-model does not bring a real added value for the resolution of our problem where edges are stress free. To be convinced, Figs 5(c) and 5(d) present a contour plot of the shear curvature κ xy through the plate for shapes S4 and S5. Recall that these shapes are symmetrical with w S4 0 (x, y) = w S5 0 (-x, y) which implies that κ S4 xy (x, y) = -κ S5 xy (-x, y). Similarly, w S4 0 (x, y) = w S5 0 (x, -y) implies that κ S4 xy (x, y) = -κ S5 xy (x, -y). It appears that the average shear curvature κxy (±3.52 10 -4 mm -1 ) is a good approximation of the local curvature at the center (saddle point) but also over a large area of the plate. Iso-curvature contours, shown as solid lines, confirm that the maximum deviation to the mean is observed in the corners of the plate. Nevertheless, the deviation remains limited with less than 10% difference between the value at the corners and the average value κxy . Predictions of the four models are compared. For all possible solutions S1 to S5, the maximal value is located at the plate corners as illustrated in Fig. 3. Configurations S4 and S5 which stand for symmetric solutions, are superimposed in Fig. 6 and display the highest value of w max 0 (e.g. 3 -4 times greater than the one of shape S3). It also appears from Fig. 6 that, for a given shape S2 to S5, w max 0 is less dependent on the chosen approach when |ε sh | becomes large. Note that this observation does not hold for the shape S1 which presents enhanced deflection when the 20p-model is adopted, see Fig. 6(c). Fig. 6 confirms that the 14p-model is sufficient to describe the stable solution in terms of maximum deflection of the plate, at least with our adopted geometrical configuration.

S1 S1 ✁ ✁ ☛ ✁ ✁ ☛ S1 S2 S3 S1 S2 S3 S1 S1 ✁ ✁ ☛ ✁ ✁ ☛ S1 S2 S3 S4 S5 S1 S2 S3 S4 S5 S1 S1 S1 ✁ ✁ ☛ ✁ ✁ ☛ ✁ ✁ ☛ S1 S2 S3 S4 S5 S1 S2 S3 S4 S5 S1 S2 S3 S4 S5
Finally, this section is complemented by results obtained with the 20p-model up to a large shrinkage strain of |ε sh | = 5 10 -4 . Fig. 7 shows that no additional deformed configuration appears at least for |ε sh | < 5 10 -4 . In addition, for the stable solution, the three curvatures and the maximum displacement evolve linearly with the magnitude of |ε sh | in the range [5 10 -5 ; 5 10 -4 ]. This result is of great interest and is one of the main outcome of the paper. Indeed, for this stable configuration, the measurement of the curvature can lead to a direct estimate of the effective shrinkage strain.

Influence of material parameters

In the previous section, analytical results were obtained considering the geometrical configuration of Fig. 2 with layer orientations [0 • /0 • ] and given material parameters: elastic moduli of Table 1 and effective shrinkage strains ε sh xx = ε sh yy = ε sh , γ sh xy = 0. We have seen previously that shapes S4 and S5 present a non-vanishing shear curvature together with principal curvatures of similar amplitude. We propose to investigate first the role of elastic orthotropy on curvature predictions. In a second part, curvatures are evaluated when the curing process induces different directional shrinkage strains. Only the 20p-model will be considered in this section.

Influence of elastic stiffness ratio (E x /E y )

The stiffness ratio E x /E y is varied in Fig. 8 between 0.5 and 2.0 to observe the influence of elastic orthotropy on the plate deflection (while holding E x = 20.2GP a). Results of section 2.3 (E x /E y = 0.89) are reported on graphs using magenta dots. In this section, we keep displacement w max 0 in terms of the stiffness ratio E x /E y . Clearly, one observes that the number of solutions evolves when E x /E y is varied. Four domains can be identified. In domain 1 of Fig. 8, when E x /E y < 0.63, the unstable shape S3 (with κxy = 0) and stable shapes S4 (κ xy > 0) and S5 (κ xy < 0) are present. Note that S4 and S5 configurations have same κxx , κyy and w max 0 so only one blue solid line is visible in figs 8 a), b) and d). When E x /E y = 0.63 (i.e. at the frontier between two domains), the additional unstable solution S2 appears which in turns gives rise to the unstable shape S1 as E x /E y becomes larger than 0.63. In region 2 , where 0.63 < E x /E y < 1.70, five solutions exist: two of them are stable with κxy = 0 and three are unstable with κxy = 0. These various stable and unstable solutions have been illustrated in the previous section for the case of E x /E y = 0.89. Note that the unstable solution S1 finally merges with S3 when E x /E y = 1.70. The third domain 3 (E x /E y > 1.70) exhibits three shapes: unstable shape S2 with κxy = 0 and two stable solutions S4 and S5 with κxy = 0. If the elastic stiffness ratio keeps increasing, the three remaining solutions finally merge when E x /E y = 1.89, the lower limit of domain 4 . As an important outcome from the analysis is that a unique stable solution exists for large values of E x /E y : shape S2 or coil set (preponderant curvature κ xx , with κ xy = 0). It is important to emphasize that curvatures are gradually evolving with the stiffness ratio so that, in general, mixed shapes (e.g. shape S3 -S4 with κxx << κyy and κxy = 0) are likely to occur. Regarding the maximum displacement, one can see from Fig. 8(d) that the highest value of w max 0 is always obtained for stable solutions. In addition, as observed for curvatures, the evolution of the maximal vertical displacement for unstable solutions shows two nose points.

Recall that with the 20p-model, all curvatures are depending on the position within the plate, see Eq. ( 23). In order to examine the level of heterogeneity developed for the stable solutions, thin black lines standing for the minimum and maximum values of each curvature within the plate are displayed on respective graphs. It appears from Fig. 8 that stable solutions provide relatively homogeneous curvatures for κ xx (y) and κ yy (x), and a very limited heterogeneity for κ xy (x, y). In fact, it is observed that the 14pmodel provides very similar evolutions for the stable solutions over the total range of stiffness ratio tested here (see Fig. S2 in supplementary material). This result confirms that the 14p-model, which relies on a constant value of curvatures, is sufficient to describe stable solutions in terms of curvatures and maximum displacement, at least for the studied configurations.

In order to have a clear view on the evolution of stable solutions over the range of E x /E y studied here, a movie showing the evolution of the shape in terms of E x /E y is provided in supplementary material. The investigation over a larger range of stiffness ratio has revealed that no additional solutions appeared. Note also that for very small ratio E x /E y (domain 1 ) and for very large value of E x /E y (domain 4 ), the curvatures and maximal displacement reach plateau values. This is illustrated in Fig. S1 of supplementary material. For small and large values of the stiffness ratio, the configuration is representative somehow of laminates with longitudinal fibers. So the coil set configuration is expected. Nevertheless, it is expected for real materials that the effective shrinkage strain would be smaller in the direction of large Young's modulus.

Influence of the effective shrinkage ratio (ε sh

xx /ε sh yy )

In most laminates used in industry (as in the PCB industry), the size and shape of the fill and warp yarns are usually different, see for instance [START_REF] Girard | Experimental and numerical characterization of thin woven composites used in printed circuit boards for high frequency applications[END_REF]. Therefore, the elastic behavior is orthotropic. Thus, it is mostly that the effective shrinkage strains in the fill or warp direction may differ. We propose in this section to illustrate the role of the difference in shrinkage strains on the deformed configuration of the bilayer materials. For that, the ratio ε sh xx /ε sh yy is varied with still ε sh xy = 0. The effective shrinkage strain in the longitudinal direction is fixed, ε sh xx = -10 -4 , and the ratio ε sh xx /ε sh yy is varied in the range [0.2, 5]. The elastic properties of the layers are still given in Table 1, so the stiffness ratio is E x /E y = 0.89. Fig. 9 shows the evolution of the three curvatures and the maximal vertical displacement as a function of the effective shrinkage ratio. A semi-log plot is proposed.

Five domains for the number of possible solutions can be identified. For effective shrinkage ratio between 0.79 and 1.26 (domain 3 ), five solutions exist for the plate deflection problem, two being stable with non-zero shear curvature. This domain includes the case of ε sh xx /ε sh yy = 1 explored in section 2.3 (magenta points). Outside this interval (domains 2 and 4 ), two unstable shapes disappear. If the shrinkage ratio is lower than 0.79 and keeps decreasing, then the three remaining solutions (two stables and one unstable) will tend to merge into the crossbow shape S3. This unique configuration exists in domain 1 when ε sh xx /ε sh yy < 0.52. On the other hand, if the shrinkage ratio is greater than 1.26 and further increases, a coil set (shape S2) will emerge and becomes the unique possible solution in the domain 5 when ε sh xx /ε sh yy > 2.11. Of particular note from Fig. 9 is the high sensitivity of κxx and κxy to a slight variation in the shrinkage ratio near the frontier between domains 1 and 2 . This implies that the final shape of two almost identical plates, i.e. made of materials exhibiting a small contrast in terms of effective shrinkage ratio, may reveal relatively marked differences. Recall also that mixed shapes may exist in a region close to boundaries of various domains where the switching from one shape to another one is occurring. The dependence of the number of unstable solutions on the effective shrinkage ratio is also reflected in Fig. 9(d) showing the evolution of w max 0 . Again, the largest maximal vertical displacement is found for the stable solution. It has been checked that the level of curvature heterogeneity for the stable solutions is still consistent with trends observed in Fig. 8. Stable solutions have relatively homogeneous curvatures for κ xx (y) and κ yy (x), and a very limited fluctuation for κ xy (x, y) (results not presented here).

Identification of effective shrinkage strains

As a practical application, the theoretical modeling is used to identify in-plane effective shrinkage strains for a bilayer laminate made of materials used in the PCB industry. Such information is missing in the literature and so the prediction of the PCB warpage could be improved with the proposed work. In this section, experimental data are first presented. A minimization procedure is next introduced.

Bilayer samples and experimental deflection

Four bilayer composite structures following the configuration of Fig. 2 were manufactured according to the process described in Fig. 1. The two layers are made of epoxy resin reinforced with woven glass fibers. One layer is already cured and the second one is still in a pre-impregnated stage (the resin being in B-stage). As reported by [START_REF] Wijskamp | Residual stresses in non-symmetrical carbon-epoxy laminates[END_REF], this elemental configuration is found to be a good testing structure to accurately estimate the effective shrinkage caused by the polymerisation of the resin at elevated temperature.

After processing, all plates exhibited at room temperature similar deformation illustrated in Fig. 10. More specifically, two specimens, labeled N o 1 and N o 3 in the following, deformed into the shape S4 (κ xy > 0) and two others, N o 2 and N o 4, were deflecting according to the shape S5 (κ xy < 0). Experimental configurations, where the maximal vertical displacement is located at two opposite corners, correspond to stable solutions predicted by the analytical approach (with the 14p-model and 20p-model). As a matter of fact, 4p-and 6p-models, for which the shear curvature is zero, are not able to predict experimental observations. The deflection of each specimen was assessed with a Vic3D measurement system (digital image correlation using a pair of stereo-mounted cameras) after preparation of a speckle pattern made of a spray paint. Samples were observed in an upright position and the spatial position of each material point was measured with respect to a reference point (RP) located at the plate center, see Fig. 11.

Experimental data for each sample reduce to a point cloud which is approximated by the second order polynomial function in x and y of the form proposed by [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF], see Eq. ( 22). Fig. 12 illustrates the accuracy of the approximations for samples N o 1 and N o 2 with corresponding coefficients of determination R 2 close to one (see Table 2).

Table 2 contains, for all specimens, average curvatures deduced from polynomial approximations of the deformed shape. One notices that the four samples present almost same curvatures along x and y with κxx being negligible with respect to the others and κyy ≃ -10 -3 mm -1 . Table 2 also reveals a discrepancy in the shear curvature κxy with an amplitude almost three times higher for specimen N o 4 when compared to specimen N o 3. One would expect that samples would reveal closer similarities. However, one should also note that even if all samples are identical in their composition and have followed the same manufacturing process, some small differences may exist at various scales, which in turn may lead to different sets of effective shrinkage strains ε sh xx and ε sh yy . These differences are believed to be at the origin of the measurement mismatch between shear curvatures. Indeed, in Fig. 9(c), the shear curvature κxy was shown to be highly sensitive to a slight variation in the effective shrinkage ratio. This trend was observed for some specific values of ε sh xx /ε sh yy (close to 0.5 when the value of ε sh xx is taken equal to -10 -4 ).

Determination of effective shrinkage strains

The determination of the effective shrinkage strains of the composite material is carried out considering all samples. We assume that the effective shrinkage strain ε sh xy is playing a minor role and a vanishing value is assumed for this parameter. Indeed, as the fill and warp yarns are perpendicular for materials used in the PCB industry (see for instance [START_REF] Girard | Experimental and numerical characterization of thin woven composites used in printed circuit boards for high frequency applications[END_REF] for microtomography views of the internal structure of similar materials), it is expected that the curing of the resin will generate only effective shrinkage strains Fig. 11 -Configuration for the deflection measurement. The spatial position of each material point within the deformed plate is measured, using DIC technique, with respect to a reference point located at the plate center. From all positions, the surface profile is approximated by a polynomial expansion as proposed in [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF].

Two offset cameras

✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✣ ❏ ❏ ❏ ❏ ❏ ❏ ❏ ❏ ❏ )
Vertical position of the plate with speckle pattern

✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✮ RP ε sh
xx and ε sh yy . Thus only these two components are left to be identified. An optimization procedure relying on the minimization of a dimensionless cost function f cost which is written in a least squares sense as:

f cost ε sh xx , ε sh yy = r i=1 β i z exp i (x i , y i ) -w 0 x i , y i , ε sh xx , ε sh yy 2 , ( 26
)
is adopted to evaluate the two unknown shrinkage strains ε sh xx and ε sh yy . z exp i is the value of the polynomial approximation at a position (x i , y i ) of a regular grid attached to the L × l plate (15 × 15 points are used in this work) and w 0 is the vertical displacement as defined in the 20p-model. r is the number of grid points (r = 225 for the adopted 15x15 grid) and β i is a weighting factor. Various expressions could have been adopted for the family (β i ). As all points of the grid seem to have similar importance for the evaluation of the deformed shape, we adopt identical values for β i :

β i = 1 r j=1 [z exp j (x j , y j )] 2 (27)
The minimization is realized for each sample using a built-in Matlab function to identify shrinkage strains (lsqnonlin function with Levenberg-Marquardt algorithm). Identified values of ε sh xx and ε sh yy as well as amplitudes of the cost function at the end of the optimization procedure, are listed in Table 3. It appears that parameters ε sh xx , ε sh yy identified for each sample are hardly dispersed, which means that the determination of the shrinkage strains is satisfactory. Note that the ratio ε sh xx /ε sh yy obtained by the optimization procedure is close to 0.5. This may explain the large fluctuation in κ xy measured experimentally. Indeed recall that Fig. 9(c) shows that a slight change in the effective shrinkage ratio (in the vicinity of the frontier between domains 1 and 2 ) may generate a large change in κ xy . Final estimations of the two effective shrinkage strains are obtained using the mean and related standard errors (adjusted using a Student's correction factor of 1.2): ε sh xx identified = -1.90 10 -4 ± 2.7 10 -6 ; ε sh yy identified = -3.43 10 -4 ± 1.7 10 -5 .

(28)

To confirm the validity of the identification process, identified shrinkage strains of Table 3 are used in the 20p-model. Analytical predictions obtained for samples N o 1 to N o 4 are compared with experimental curvatures in Fig. 13 and in Table 4. 3.

Tab. 4 -Comparison of curvatures from measured points and from predictions using 20p-model and shrinkage values of Table 3. Results with mean values of Eq. ( 28) are also reported for comparison.

Sample

Method κxx [mm -1 ] κyy [mm -1 ] κxy [mm -1 ] N o 1
Measure 8.65 10 -6 -9.93 10 -4 2.23 10 -4 N o 1 Prediction -1.46 10 -5 -9.74 10 -4 2.15 10 -4 N o 2 Measure -4.30 10 -5 -9.00 10 -4 -5.40 10 -4 N o 2 Prediction -7.08 10 -5 -8.37 10 -4 -4.77 10 -4 N o 3 Measure 1.78 10 -5 -1.03 10 -3 2.11 10 -4 N o 3 Prediction -1.32 10 -5 -9.93 10 -4 2.05 10 -4 N o 4 Measure -8.85 10 -5 -8.50 10 -4 -6.42 10 -4 N o 4 Prediction -1.21 10 -4 -7.94 10 -4 -6.13 10 -4 Prediction (using Eq. 28) -5.35 10 -5 -9.01 10 -4 ± -4.28 10 -4

Even if some discrepancy can be noticed, the negligible amplitude of the curvature κxx is restituted, while theoretical values of κyy and κxy are in a close agreement with experiments. It appears that the theoretical prediction of each curvature calculated using the effective shrinkage strains of Eq. 28 falls close to the corresponding mean value calculated from all specimens. More specifically, the average value κxx and κyy calculated from analytical results of each specimen is -5.49 10 -5 mm -1 and -8.99 10 -4 mm -1 respectively. These values fall close to the theoretical predictions -5.35 10 -5 mm -1 and -9.01 10 -4 mm -1 deduced from the mean effective shrinkage strains of Eq. 28, see Table 4.

To the best of authors' knowledge, effective shrinkage strains for the material tested in this paper are not available in the literature. With the present method, an efficient methodology is proposed, which opens the possibility to investigate precisely the development of warpage in thin complex multimaterial composites.

Conclusion

An analytical modeling, extending CLT approaches, has been developed to analyze flatness defects occuring during the assemby of multilayer composites. The warpage originates from difference in coefficients of thermal expansion and in effective chemical shrinkage strains between layers. The theory was applied to a symmetric bilayer structure under virtual cure, but is general enough to be used for any multilayer composites submitted to more general plane stress thermo-mechanical loadings. Validation, conducted from a comparison with finite element calculations, was proposed to check the predictive capabilities of our analytical approach. Even if some discrepancies were noticed in the level of heterogeneity of curvatures obtained in numerical simulations, the plate deflection was well captured by the modeling.

To detect all possible stable and unstable configurations of the laminates by analytical means, trial fields for the displacement and strain available in the literature were adopted. It was seen that for limited effective shrinkage, the stable solution is representative of a coil shape and all trial fields are salient for this case. It has been shown as the magnitude of the effective shrinkage strain is enlarged, only the modeling route with enriched expressions of the trial fields can capture configurations with varying curvature in the plane. As expected, models with a limited number of parameters cannot give access to all possible configurations. In the case of trial fields with 14 or 20 parameters, all configurations S1 to S5 are found. It could be of interest to define additional enrichments of the trial fields so that another configurations could be obtained. It could be also of interest to account for non-vanishing shear curvatures which seem to be responsible of additional stable solutions. As the number of parameters becomes large, the solution is obtained through an in-house numerical solving procedure combining arc-length method and deflation technique.

The combined effects of stiffness and effective shrinkage ratio have been examined by performing a parametric analysis. The relationship between material properties and curvatures of the plate, intimately related to flatness defects, have been uncovered for a bilayer laminate. The existence and stability of solutions are depending on the material properties. Interestingly, asymptotic behavior was observed for large mismatch in the stiffness ratio or in the effective shrinkage strain ratio.

The analysis developed in this paper furnishes a fundamental tool to study flatness defects, investigate stable or unstable solutions. In line with this concern, an inverse identification method was implemented to identify effective shrinkage strains from a series of manufactured samples. Values identified show only a slight variation between samples. This result also participates to the validation of our modeling.

In the future, the model could be extended to account for the effect of the ply misorientation and stacking of a general laminate. In addition, of prime interest for PCB manufacturers is to analyze the thermo-mechanical response of multilayer structures containing holes. This can be envisioned with the development of finite element calculations of a more complex assembly, for which effective shrinkage strains identified in this work would be used as input parameters to model the manufacturing step.

A Numerical resolution of the total potential energy minimization problem

This appendix details how the system of nonlinear equations ( 11) was numerically treated to obtain the solution of the total potential energy minimization problem presented in section 2. Matlab software was used in this work but another software could have been used. Let us first mention that for [START_REF] Hyer | Calculations of the room-temperature shapes of unsymmetric laminates[END_REF], [START_REF] Jun | Effect of residual shear strain on the cured shape of unsymmetric cross-ply thin laminates[END_REF] and [START_REF] Dano | Snap-through of unsymmetric fiber-reinforced composite laminates[END_REF] models (respectively 4, 6 and 14 parameters), one can use vpasolve function to solve the symbolic system (11) for a given shrinkage strain tensor ε sh . This procedure brings directly all the solutions which ensure the gradient nullity. As seen in section 2.3, there can be one, three or five solutions for the deformed bilayer plate, depending on the adopted modeling and the value of the components of ε sh . For illustrative purpose, results in the Appendix are obtained with ∆T = 0 • C, γ sh xy = 0 and ε sh is varying in the range [-10 -4 , 0].

A.1 Arc-length method

For the [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF] approach (20 parameters), an alternate approach was needed due to the large size of the system (20 unknowns). Therefore, we have developed our own procedure. As results exhibit critical points (bifurcation and nose points) and disconnected branches, we have chosen an arc-length method (continuation method originally developed by [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF][START_REF] Riks | The application of newton's method to the problem of elastic stability[END_REF] which is recommended for the analysis of load-deflection structure equilibrium. Continuation methods are also widely used for the computation of bifurcation points in fluid mechanics, see [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF] for instance. A load coefficient λ is introduced in Eq. ( 2) which, under the condition ∆T = 0 considered in the results section, is written as follows:

ε th (k) = λ ε sh (k) with λ ∈ [0; 1] . (29)
This leads to a slightly modified total potential energy expression:

Π = L/2 -L/2 l/2 -l/2   1 2 ε 0 κ T A B B D ε 0 κ -λ ε 0 κ T N th M th   dydx = Π p, λ . ( 30
)
p standing for the vector of parameters p i , i = 1, . . . , m. In above relationships, ε sh (k) is the prescribed shrinkage. N th and M th are thermal forces and thermal moments per unit surface associated to ε sh (k) , see Eq. ( 8).

The minimization of Π is now expressed as: In Newton-Raphson algorithms, only p is unknown and λ is linearly evolving from 0 to 1 to reach the solution. By contrast, in continuation methods both p and λ are unknowns and simultaneously varied.

For the path-following technique, an arc-length parameter s is then introduced to control the evolution of p and λ: p = p (s) and λ = λ (s) . ( 32)

Assume now that the point p k , λ k = p (s = s k ) , λ (s = s k ) , with k ∈ IN, satisfies the system of nonlinear equations (31). Continuation method consists of finding p k+1 , λ k+1 such as:

   f p k+1 , λ k+1 = 0 p k+1 -p k T • ṗk + (λ k+1 -λ k ) λk = ∆s (33)
where (•) = ∂(•) ∂s . The second equation (arc-length equation) complements the system which is now composed of (m + 1) equations for (m + 1) unknowns.

The resolution of system (33) requires three steps: tangent vector determination, prediction and correction.

i. Considering that the initial tangent vector (i.e. for k = 0) is defined as:

   f p 0 ṗ0 + f λ 0 λ0 = 0 ṗ0 2 + λ2 0 = 1 (34)
the tangent vector ṗk , λk , k = 0, is first determined from the following incremental relationship: ii. Then an increment of arc-length parameter is imposed as ∆s k = s k+1 -s k (adaptative increment size) and the tangent predictor is built as follows :

   f p k ṗk + f λ k λk = 0 ṗk-1 T • ṗk + λk-1 λk = 1 such as ṗk 2 + λ2 k = 1 (normalization), (35)
p k+1 (0) = p k + ∆s k ṗk and λ k+1 (0) = λ k + ∆s k λk ; (36)
iii. Finally, an incremental correction step is carried out, based on the Newton's method:

p k+1 (j+1) = p k+1 (j) + ∆p k+1 (j) and λ k+1 (j+1) = λ k+1 (j) + ∆λ k+1 (j) , (37)
where ∆p k+1 (j) and ∆λ k+1 (j) are obtained from solving the following system:

  f p k+1 (j) f λ k+1 (j) ṗk T λk   ∆p k+1 (j) ∆λ k+1 (j) = -   f p k+1 (j) , λ k+1 (j) ṗk T • p k+1 (j) -p k + λk λ k+1 (j) -λ k -∆s k   . ( 38
)
Iterations end when ∆p k+1 (j) 2 + ∆λ k+1 (j) 2 ≤ η 2 , η being a given precision threshold.

In practice, a conjugate gradient method is used to solve system (34) whereas the bordering algorithm is employed for resolution of systems (35,38) (see [START_REF] Keller | The bordering algorithm and path following near singular points of higher nullity[END_REF], for details about this algorithm). In our calculations, the value adopted for the precision threshold in the correction step was η = 10 -10 . This procedure was applied from the starting point p 0 = 0, λ 0 = 0 to get the first solution branch of system (31) (see Fig. 14). Note that this first branch could have been also determined with a classical Newton-Raphson algorithm.

A.2 Disconnected branches

The analytical model shows that beyond certain shrinkage values, several different shapes appear and some solutions become unstable. These multiple solutions are classically obtained by bifurcation point detection (nullity of a test function at this point) and branch switching performance (with continuation on each detected branch). But, in some cases, there is no bifurcation point and the solution branches are disconnected (as seen in Fig. 4 for instance). It is therefore impossible to detect others branches using the approach described in section A.1. To overcome this difficulty, deflation technique is considered (see [START_REF] Farrell | The computation of disconnected bifurcation diagrams[END_REF], for details on the so-called deflated continuation). This method can find multiple solutions from only one single initial point, even when this point is not a bifurcation point. This approach thus is well adapted to our problem. Its principle is based on the modification of the system of nonlinear equations. We assume that a solution p 1 , λ 1 is already known (end of the first continuation branch for instance). The initial problem f p, λ 1 = 0 is replaced by:

g p = M p, p 1 f p, λ 1 = 0 ( 39
)
where M is a m × m diagonal matrix defined by:

M p, p 1 = 1 p -p 1 q + σ 1, (40)
1 being the second-order identity tensor. σ is a shift parameter (used to improve convergence) and || • || is the 2-norm of a vector. The Newton-Raphson method is chosen to find another solution p 2 of g p = 0. Hence, from Eq. (39), p 2 is also solution of Eq. (31) (since f p 2 , λ 1 = 0). Arc-length continuation is then conducted to solve system (31) from the new solution p 2 , λ 1 , in the range λ ∈ [0; 1]. This procedure (deflated continuation) is repeated until no other solution can be found. Fig. 15 illustrates the sequence of alternating continuation and deflation steps after determination of the first continuation branch using the arc-length method of section A.1. From the end point of the continuation branch, represented in Fig. 15(a) by a plain square, the deflation technique is applied and a second solution is found (plain cercle in Fig. 15(a)). The arc-length method is next used to build the second continuation branch (from the plain circle to the triangle). One notes that this new branch presents a nose point so that classical Newton-Raphson algorithms are not able to capture this path entirely. An additional solution is obtained from the deflation technique, using Eq. ( 39) and the third solution (triangle) as the initial point. This solution is represented by a plain diamond in Fig. 15(b) and becomes the starting point of a new continuation branch which needs to be determined using an arc-length method. Indeed, the third point leads to two distinct solutions (plain and hollow diamonds) in the shear curvature κxy and an examination of the diagram inserted in Fig. 15(b) reveals a bifurcation point (hollow cercle). Finally, a last deflation step is run with the hollow diamond as a starting point. No additional solution is achieved and the numerical resolution of system (31) is stopped.

B Numerical simulation and validation of the proposed model

To verify the proposed analytical model of plate deformation due to shrinkage, finite element calculations have been carried out with the commercial software Abaqus/Standard (6.13-2 version). First, the meshing strategy is discussed in section B.1 and comparison is exemplified in section B.2.

B.1 Finite Element modeling

Tab. 5 -Out-of-plane elastic properties of the fiberglass reinforced epoxy laminate at room temperature (20

• C). E z [GPa] ν xz ν yz G xz [GP a] G yz [GP a] 16.34 0.451 0.451 2.71 2.71
The bilayer thin structure [0 • /0 • ] of Fig. 2 is adopted for the analysis. Both layers have orthotropic elastic behavior with moduli listed in Table 1 andin Table 5. The structure is loaded during one step of virtual cure for which the theoretical approach predicts flatness defects due to shrinkage. Geometrical nonlinearities are included in the calculation. A Fortran user-subroutine UEXPAN has been written to reproduce the shrinkage occurring in the upper layer (see Kpobie et al., 2016, for details). Boundary conditions consist in considering that the center of the plate is clamped (no translation and rotation). Because of aspect ratio of the studied structure (h << L, l) and of bending phenomenon that occurs, the plate is discretized using composite shell finite elements. Simpson method is chosen for numerical integration through thickness (with 15 points in each layer). Triangular and quadrilateral elements with linear and quadratic interpolation were tested. Following the designation used in Abaqus, these elements are S3R (three noded triangular element), S4R (four noded quadrilateral element), STRI65 (six noded triangular element) and S8R (eight noded quadrilateral element). The laminate is meshed with 96x60, 48x32, 48x30 and 24x16 elements (element aspect ratio close to one) when having adopted element type S3R, STRI65, S4R, and S8R. Note that mesh size independence has been checked. Fig. 16 shows the deformed plate for ε sh xx = ε sh yy = -10 -4 and ε sh xy = 0 (results are obtained using the four meshes specified above). Note that from our analytical modeling, the stable configuration is a S5 shape, see Fig. 3. From Fig. 16, it is observed that results are not depending on the degree of interpolation but differ from triangular to quadrangular elements. Fig. 16 reveals that a structure meshed with triangular elements S3R or STRI65 deforms into shape S5 (maximal vertical displacement at two opposite corners and κ xy < 0), while a coil set deformation (shape S2) is obtained when quadrangular elements S4R or S8R are used. In fact, the shape reached with quadrangular elements may be an artifact of the mesh structuring which, in this configuration, is aligned with principal directions of curvatures. This is confirmed from additional calculations conducted with a bias introduced in the mesh made of quadrangular elements which yield identical results as that obtained using triangular elements, see Fig. 16. As a consequence, it is believed that results obtained with the structured mesh using quadrangular elements suffer from a lack of physical meaning and this configuration is rejected for the comparative analysis. In addition, it has been proved that numerical simulations of the plate considered in this work do not benefit from a quadratic interpolation. Based on the above analysis, linear triangular elements S3R are selected to conduct a comparison between analytical and numerical results.

B.2 Comparison between analytical and numerical results

Fig. 17 shows a comparison of the average curvatures and maximum vertical displacement between numerical simulation (FEM ) and analytical modeling (stable and unstable solutions) obtained with the 20p-model. Curvatures, based on the FEM are obtained from second order derivative of the plate deflection and further approximated by the same second-order polynomial as in [START_REF] Mattioni | Multistable composite plates with piecewise variation of lay-up in the planform[END_REF]. A very good agreement is observed for κxx , while differences remain for κyy , κxy and w max 0 . In fact, the critical value of ε sh corresponding to the occurrence of a new solution (different shape), is larger in FEM which creates a minor offset in the evolution of average curvatures. The difference in the critical value of ε sh may be explained by the complex form of local curvatures within the plate. Fig. 18 depicts the variation of curvature κ yy (respectively κ xx ) along the plate length (respectively the plate width), when ε sh = -10 -4 . As it was already highlighted in Figs. 5(a) and 5(b), the 20pmodel leads to quasi homogeneous curvatures for the stable solution, while the numerical simulation is shown to predict curvatures that significantly change with the position in the plate. In addition, in the analytical 20p-model, κ xx (respectively κ yy ) is depending on y (respectively x) only, and is approximated using a quadratic form. By contrast, curvatures obtained from finite element calculations behave in a more complex manner. Moreover, one observes that results from FEM are not necessarily symmetric with respect to the x-axis and the y-axis while by construction, results are symmetric in the 20p-model. It thus appears that the analytical approach is too restrictive in terms of curvature approximation and so trial fields need to be enriched to reflect higher order dependencies. Finally, results are quite similar to those obtained from FEM which illustrates the good predictive capabilities of the proposed modeling.

Fig. 19 displays a comparison between the predicted deflections w 0 (x, y). It is seen that the difference ∆w 0 (x, y) = w analytical 0 (x, y) -w F EM 0 (x, y) between the two predictions of the plate deflection is limited and a maximum of 1mm difference is found in the vicinity of the four corners. A 6.6% accuracy is thus predicted for the deflection of the four corners.

Fig. 1 -

 1 Fig. 1 -Schematic time evolution of the imposed temperature (left axis) and pressure (right axis) during manufacturing of PCB.

Fig. 2 -

 2 Fig.2-Thin bilayer composite. The geometrical configuration adopted in the present work is : L = 470mm, l = 300mm, h = 0.6mm.

  ) depicts situations where |κ xx | ≃ |κ yy | while Fig. 3(b) (respectively Fig. 3(c)) is related to cases where | κxx κyy | ≫ 1 (respectively | κxx κyy | ≪ 1). To some extent, S2 of Fig. 3(b) and S3 of Fig. 3(c) can be seen as 'orthogonal' shapes. Figs 3(

Fig. 3 -

 3 Fig.3-Deformed shapes of the bilayer plate at the end of the virtual cure process when ε sh xx = ε sh yy = ε sh = -10 -4 , γ sh xy = 0. Material parameters are provided in Table1. The plate dimension is L = 470mm, l = 300mm, h = 0.6mm. The 20p-model is used.

  Figs 4(a) and 4(b) reveal that the 6p-model provides similar solutions in terms of numbers, shapes and stability. The main differences lie in (i) the critical value of ε sh where ∂κ xx /∂|ε sh | = 0 (|ε sh | ≃ 2 10 -5 for the 4p-model and |ε sh | = 3 10 -5 for the 6p-model) and (ii) the value of ε sh from which several solutions emerge (|ε sh | = 3 10 -5 for the 4p-model and |ε sh | = 5 10 -5 for the 6p-model). An interesting feature induced by a further enrichment in the modeling is illustrated in figures 4(c)-4(d) by comparing the predictions of the 6p-model and the 14p-model. Due to non vanishing shear curvature for the 14p-model, two additional shapes are found when |ε sh | is close to the critical value |ε sh | = 3 10 -5 . These shapes, corresponding to S4 and S5 of Fig. 3, are defined by an identical value of κ xx , a ratio κ xx /κ yy close to unity but two opposite values for κ xy , see Fig 4(g). Obviously, these solutions do not exist in 6p-model where κ xy = 0 is imposed by construction of the model. More important is that these 'two symmetric' solutions (in the sense that w S4 0 (x, y) = w S5 0 (-x, y) and w S4 0 (x, y) = w S5 0 (x, -y)) are the unique stable solutions among the various possible shapes predicted by the approach. As a consequence, the 14p-model predicts a unique family of stable solutions (if one considers S4 and S5 as similar shapes) for the whole range of shrinkage values tested in the present work. Figs 4(e-f-g) compare results obtained with 10

  Variation of the longitudinal curvature along the plate width for the five deformed shapes. Variation of the transversal curvature along the plate length for the five deformed shapes.

  Shear curvature κ xy for shape S4.

Fig. 5 -

 5 Fig. 5 -Variation of curvatures in the plate (obtained by 20p-model). Stable solutions are represented in figs.(a) and (b) with solid lines whereas unstable solutions are displayed with dashed lines. The effective shrinkage strain is set to ε sh = -10 -4 . The plate dimensions are those of Fig. 3.

Fig. 6 Fig. 6 -

 66 Fig. 6 -Evolution of maximal vertical displacement as a function of shrinkage strain |ε sh |. Stable solutions are represented with solid lines whereas unstable solutions are displayed with dashed lines.

Fig. 7 -

 7 Fig. 7 -Evolution of average curvatures a) κxx , b) κyy , c) κxy (see Eq. 24) and d) maximal displacement for the 20p-model as a function of shrinkage strain. The range of variation of |ε sh | is 0; 5 10 -4 . Stable solutions are represented with solid lines whereas unstable solutions are displayed with dashed lines.

Fig. 8 -

 8 Fig. 8 -Evolution of curvatures a) κxx , b) κyy , c) κxy and d) maximal vertical displacement as a function of the elastic orthotropy ratio E x /E y for |ε sh xx | = |ε sh yy | = 10 -4 . In Figs a-c), dashed lines correspond to the average value of the curvature for an unstable solution. Stable solution is visualized with a solid line. Average curvatures are represented by a thick blue line while the minimal and maximal curvature values are displayed with thin black lines. Magenta dots are identifying the case E x /E y = 0.89 studied in section 2.3. In all calculations, E x = 20.20GP a.

Fig. 9 -

 9 Fig. 9 -Evolution of curvatures a) κxx , b) κyy , c) κxy and d) maximal vertical displacement as a function of effective shrinkage strain ratio ε sh xx /ε sh yy for |ε sh xx | = 10 -4 (semi-log scale). Results are obtained considering the 20p-model for which the average values of curvatures are in blue line. Magenta points are identifying the case ε sh xx = ε sh yy studied in section 2.3.

Fig. 10 -

 10 Fig. 10 -View of a deformed bilayer sample (specimen N o 1). The maximum deflection (around 15mm for this specimen) is located at two opposite corners. The specific sample presents a S4 shape.

Fig. 12 -

 12 Fig. 12 -Comparison between measured points and polynomial approximation (meshgrid) Tab. 2 -Curvatures of the bilayer samples calculated from polynomial approximation of the measured point cloud Sample Coefficient of κxx [mm -1 ] κyy [mm -1 ] κxy [mm -1 ] determination R 2 N o 1 0.9953 8.65 10 -6 -9.93 10 -4 2.23 10 -4 N o 2 0.9974 -4.30 10 -5 -9.00 10 -4 -5.40 10 -4 N o 3 0.9973 1.78 10 -5 -1.03 10 -3 2.11 10 -4 N o 4 0.9988 -8.85 10 -5 -8.50 10 -4 -6.42 10 -4

Tab. 3 -

 3 Identified effective shrinkage strains and cost function values at the end of the minimization procedure.

Fig. 13 -

 13 Fig. 13 -Comparison between experimental (blue points) and theoretical (surface) plate deflections. Theoretical results have been obtained using the 20p-model and shrinkage values of Table3.

  i = 1, . . . , m rewritten as f p, λ = 0 where f : IR m+1 → IR m . (31)

  k , λ k and f λ k = ∂f ∂λ p k , λ k ;

Fig. 14 -

 14 Fig. 14 -First continuation branch obtained by numerical resolution of system (31), a) average curvature κxx , b) average curvature κyy versus shrinkage strain |ε sh xx |

Fig. 15 -

 15 Fig. 15 -Second (a) and third (b) continuation branches obtained with the deflation technique.

Fig. 16 -

 16 Fig. 16 -Plate deformation obtained from numerical simulations using various types of finite elements a) shape S5, b) shape S2. U 3 corresponds to the vertical displacement in millimeters.

Fig. 17 -Fig. 18 -

 1718 Fig. 17 -Evolution of curvatures a) κxx , b) κyy , c) κxy and d) maximal vertical displacement as a function of shrinkage strain. Comparison of the predictions obtained from analytical modeling and finite element calculations.

Fig

  Fig. 19 -a) Comparison of the two plate deflections captured by the analytical approach and the simulations. b) Map of the difference between the two predictions ∆w 0 (x, y) = w analytical 0

Acknowledgments

The authors would like to gratefully thank the Région Grand-Est, the Département de la Moselle, the Metz Métropole, CIMULEC, SYSTRONIC and CSI Sud Ouest through the NIT foundation for their financial support. The authors would also like to acknowledge the support of the French government through the National Research Agency ANR and the LabCom programm (ANR-14-LAB7-0003-01) for the creation of LabCom LEMCI.

Data availability

The raw and processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.


Nanying Liang

Laurent Bougrain

Decoding Finger Flexion using Amplitude Modulation from Band-specific ECoG

EEG-BCIs have been well studied in the past decades and implemented into several famous applications, like P300 speller and wheelchair controller. However, these interfaces are indirect due to low spatial resolution of EEG. Recently, direct ECoG-BCIs attract intensive attention because ECoG provides a higher spatial resolution and signal quality. This makes possible localization of the source of neural signals with respect to certain brain functions. In this article, we present a realization of ECoG-BCIs for finger flexion prediction provided by BCI competition IV. Methods for finger flexion prediction including feature extraction and selection are provided in this article. Results show that the predicted finger movement is highly correlated with the true movement when we use band-specific amplitude modulation.

Introduction

A brain-computer interface (BCI) is a new communication and control option for those with severe motor disabilities [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. A BCI system translates brain signals into commands for a computer or other devices. Electroencephalography (EEG) based BCIs are the most studied non-invasive interfaces mainly due to the fine temporal resolution, ease of use, portability and low set-up cost of EEG recordings [START_REF] Hughes | EEG in Clinical Practice[END_REF]. Unfortunately, non-invasive implants produce noisy signals because of the damping effect introduced by the skull. Further, its spatial resolution is low due to the size of the electrode. Thus, it is difficult to use EEG signals for directly decoding the task events.

Electrocorticography (ECoG) has recently emerged as a promising recording technique for use in brain-computer interfaces (BCI) and has been extensively investigated [START_REF] Schalk | Decoding two-dimensional movement trajectories using electrocorticographic signals in humans[END_REF][START_REF] Sanchez | Extraction and localization of mesoscopic motor control signals for human ECoG neuroprosthetics[END_REF][START_REF] Chin | Identification of arm movements using electrocorticograhpic signals[END_REF][START_REF] Pistol | Prediction of arm movement trajectories from ECoG-recordings in humans[END_REF]. ECoG electrodes arrays were initially implanted (under the skull but over the surface of the cortex) for severe epileptic patients for presurgical planning, in order to identify the sources generating epileptic seizures [START_REF] Asano | Origin and propagation of epileptic spasms delineated on electrocorticography[END_REF]. The spatial resolution of ECoG signals is higher than EEG with a better signal-to-noise ratio (SNR). Typically, the diameter of one ECoG electrode is of 4 mm with 1 cm inter-electrode distance. Therefore ECoG can provide a spatial resolution of approximately 1 cm [START_REF] Asano | Origin and propagation of epileptic spasms delineated on electrocorticography[END_REF]. Spatial resolution plays an important role in BCI [START_REF] Sanchez | Extraction and localization of mesoscopic motor control signals for human ECoG neuroprosthetics[END_REF]. The fine spatial resolution of ECoG provides an opportunity for decoding brain functions directly and therefore possibly leads to the implementation of direct neural interfaces, which is difficult to be accomplished through EEG-based BCIs.

In order to study the application of ECoG in BCIs, several research groups have recorded ECoG signals when the participants performed certain kind of tasks related to the brain functional areas that the implanted electrode arrays had covered. The testing tasks include center-out reaching or pointing task [START_REF] Sanchez | Extraction and localization of mesoscopic motor control signals for human ECoG neuroprosthetics[END_REF], finger flexion [START_REF] Schalk | Decoding two-dimensional movement trajectories using electrocorticographic signals in humans[END_REF] and cursor trajectory [START_REF] Pistol | Prediction of arm movement trajectories from ECoG-recordings in humans[END_REF]. Some of these groups are interested in the use of frequency features as the inputs of decoder [START_REF] Pistol | Prediction of arm movement trajectories from ECoG-recordings in humans[END_REF], and others use amplitude modulation [START_REF] Sanchez | Extraction and localization of mesoscopic motor control signals for human ECoG neuroprosthetics[END_REF], or both [START_REF] Schalk | Decoding two-dimensional movement trajectories using electrocorticographic signals in humans[END_REF]. For features extracted from frequency domain, it was found that different frequency bands have different correlation coefficients to the event tasks [START_REF] Schalk | Decoding two-dimensional movement trajectories using electrocorticographic signals in humans[END_REF][START_REF] Chin | Identification of arm movements using electrocorticograhpic signals[END_REF]. Further, the implanted ECoG electrode array is usually a grid with a size from 6*6 to 8*8. It covers a cortical zone involving several different functions of the brain. Thus, frequency band and electrode selection are necessary for identifying good descriptors.

In this paper, we first present the material used, the ECoG data set adopted from BCI competition IV in section 2. The methods for feature extraction and selection are presented in 3.1. The construction of a linear regression model for decoding is done in 3.2. Simulation results on this competition data set are reported in section 4. We draw conclusions with discussions on further work in section 5.

Material

Recently several paradigms of ECoG-BCIs have been realized focusing on motor tasks. Here, we just adopt the ECoG-BCI data set provided by BCI competition IV [START_REF] Schalk | Decoding two-dimensional movement trajectories using electrocorticographic signals in humans[END_REF][START_REF] Miller | Prediction of finger felxion: 4 th brain-computer interface data competition[END_REF]. Three subjects participated in the recording were epileptic patients. Each subject was implanted a subdural electrode array with the arrangement of 8*6 or 8*8 electrodes. The ECoG signals were recorded when the subjects performed a finger movement task. The subjects were instructed to move one certain finger by the corresponding word displayed on a computer screen. The execution of finger movement lasted 2 seconds and it was followed by a 2-second resting period. There were 30 movement stimulus for each finger resulting in 600-second recording for each subject. The ECoG signals were recorded through the general-purpose BCI system BCI2000 [START_REF] Schalk | BCI2000: A general-purpose brain-computer interface (BCI) system[END_REF], bandpass filtered between 0.15 to 200 Hz and sampled at 1000 Hz. The finger movements were recorded using a dataglove sampled at 25 Hz. Figure 1 shows an example of the visualization of the ECoG signals and the finger movement time course from subject 1. Due to space limitation, only a subset of electrodes is displayed.

Decoding Workflow

Feature Extraction and Selection

The evidence of sensorimotor ECoG dynamics has been reported in several specific frequency bands including sub-bands (1-60 Hz), gamma band (60-100 Hz), fast gamma band (100-300 Hz) and ensemble depolarization (300-6k Hz) [START_REF] Sanchez | Extraction and localization of mesoscopic motor control signals for human ECoG neuroprosthetics[END_REF]. Being inspired by the firing rate coding scheme used in spike train decoding, Sanchez proposed the band-specific amplitude modulation (AM) as the descriptor for ECoG signal decoding, which is defined as the sum of the power of the voltage of the band-specific ECoG signals in a time bin ∆t:

x(t n ) = ∆t t=0 v 2 (t n + t) (1)
where t n+1 = t n + ∆t. We simply let ∆t = 40ms so the resulting bandspecific AM feature inputs have the same sampling rate (i.e. 25Hz) as that of the dataglove position measurements. It seems that the sensitivity profile of ECoG for each frequency band is specific to a given task [START_REF] Sanchez | Extraction and localization of mesoscopic motor control signals for human ECoG neuroprosthetics[END_REF]. For our study, we also found that each finger movement was correlated to two or three specific electrodes. Therefore, feature selection including frequency band and electrode site selection is necessary for removing the unrelated input descriptors. For each finger and each subject, we use a stepwise feature selection procedure based on the method of train and validation (i.e. 3/5 of available data are used for training and 2/5 for validation). In the future, we will refer to this method as FD-ECoG AM for frequency decomposition with feature extraction scheme.

Next, we want to compare it with other methods for estimating AM features: raw-ECoG AM feature and PCA-ECoG AM feature. For computing the raw-ECoG AM feature, we repeat the above procedure but exclude the frequency decomposition operation. For PCA-ECoG AM feature, we first apply the principal component analysis (PCA) method to the whole ECoG data set, and then for each principal component, we calculate the corresponding AM feature. After obtaining the descriptor features, we are ready to present a decoder for predicting finger movement in the following section.

Linear Regressor

We have no a priori information about the relationship between the descriptors and the target signals or the interaction between descriptors. Thus, we simply applied a linear model as a decoder for its robustness property. We also noticed that other advanced methods have been used for ECoG signals decoding, for example, the Kalman filter [START_REF] Pistol | Prediction of arm movement trajectories from ECoG-recordings in humans[END_REF] and nearest neighbour classifier (NNC) [START_REF] Chin | Identification of arm movements using electrocorticograhpic signals[END_REF]. But these methods are not suited for the case being studied here because the first method needs a finger model which we do not have and the second one is suitable for classification. Linear models have been applied to ECoG signal analysis by other research groups with successful outcome [START_REF] Schalk | Decoding two-dimensional movement trajectories using electrocorticographic signals in humans[END_REF][START_REF] Sanchez | Extraction and localization of mesoscopic motor control signals for human ECoG neuroprosthetics[END_REF]. Here, we adopt the linear model with the form as follows:

d(t) = W T x(t) ( 2
)
where d is the finger position as measured by a dataglove. x(t) is the tapdelay AM feature vector. The coefficients W of the model are trained by the Wiener solution in Equation ( 3):

W = E(x T x) -1 E(x T d) (3)
where E is the expected mean. The number of tap-delay is optimized with the value of 25 tap-delays for our case (i.e., using the present AM input and the previous 1-second AM inputs for predicting the present finger flexion). Further, in order to improve the stability for estimating the coefficients of the Wiener model, we replace the inverse operation in Equation ( 3) with the pseudo-inverse.

Results

First, we present the simulation results of the linear decoder combined with different feature extraction methods as summarized in Table 1 1 .

Table 1 shows that the method based on AM features generated from bandspecific ECoG signals outperforms the others. It can be inferred that the sensitivity profile of ECoG is both band-specific and channel-specific to a task and subject. Therefore feature selection is important for the ECoG signal decoding. For example, in the sense of correlation, channel 1 ranks highest for flexion decoding for thumb and index finger, and channel 39 ranks highest for middle, ring and little finger for subject 1. We provide an example for the predicted finger movement for subject 3 based on the method of FD-ECoG AM as shown in solid red curve in Figure 2. For comparison purpose, the corresponding true finger movement time course is plotted in the dash blue curve.

Subj

Conclusion and Further Work

In this paper, we have studied a linear decoding scheme for predicting finger movement from ECoG signals. The approach presented in this article won the BCI competition IV addressed to this problem. The experimental results show that, with the fine spatial resolution of ECoG signals, it is promising to implement a direct neural interface in applications, for example, dedicated robotic arm/hand control and simple gesture language communication such as cued language. Moreover, when comparing the results obtained by the method taking into account the frequency-specific features and those not, it can be inferred that the sensitivity profile of ECoG signals is band-specific and channel-specific to a task and subject.

In the stepwise feature selection procedure, we did not consider the correlation between band-specific ECoG signals. It is suggested that incorporating the feature correlation into feature selection, for example, using correlation feature selection (CFS) method [START_REF] Hall | Correlation-based feature selection for machine learning[END_REF], may produce an optimal compact feature set.

We have no a priori information about the relationship between the descriptors and the target signals or the interaction between descriptors and thus make a simple assumption of linear regression between the descriptors and the target signals. We try to introduce some user-defined interaction terms into the regression model. It may improve the decoding accuracy using less number of regressors. Thus, non-parameter models incorporating the interaction terms into the decoding model, for example, using multivariate adaptive regression splines (MARS) [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] will be investigated in the future.

Fig. 1 :

 1 Fig. 1: Example ECoG signals from channel 46 to 50 (the first 60 seconds in the training data set) from subject 1. The last 5 rows show the corresponding finger movement time course.

Fig. 2 :

 2 Fig.2: Predicted (solid red) and real (dash blue) finger movement time courses for the first 60 seconds in validation data set from subject 3.

Table 1 :

 1 . Performance of the linear decoder combined with different feature extraction methods in terms of correlation coefficient between the predicted and true finger movement for each finger and subject. The last column presents the averaged result for each method and subject.

		AM feature	Thumb Index Middle Ring Little	Av.
		Raw ECoG AM	0.26	0.26	0.25	0.40	0.32	0.30
	1	PCA-ECoG AM	0.22	0.29	0.22	0.33	0.23	0.26
		FD-ECoG AM	0.56	0.76	0.37	0.62	0.48	0.56
		Raw ECoG AM	0.46	0.31	0.38	0.41	0.25	0.36
	2	PCA-ECoG AM	0.30	0.25	0.28	0.27	0.19	0.26
		FD-ECoG AM	0.64	0.40	0.50	0.57	0.45	0.51
		Raw ECoG AM	0.63	0.46	0.54	0.56	0.39	0.52
	3	PCA-ECoG AM	0.31	0.34	0.44	0.42	0.35	0.37
		FD-ECoG AM	0.73	0.68	0.78	0.68	0.64	0.70
	Thumb
	Index
	Middle
	Ring
	Little
	0	10	20	30	40		50	60
				Time (s)

The performance is evaluated only with the validation data set because the finger position measurements because the testing data set are not available yet by the BCI competition IV.