M Lezgy-Nazargah

P Vidal

O Polit

A 1D nonlinear finite element model for analysis of composite foaminsulated concrete sandwich panels

Introduction

Due to specific properties such as high thermo-acoustic efficiency, high blast resistance, high workability and low self-weight characteristics, foam-insulated concrete sandwich panels are nowadays applied extensively in constructions of various civil engineering structures. As shown in Fig. 1, the FICSP structures are made of insolation foam, concrete, shear connectors, longitudinal rebar and welded wire reinforcements (WWR). The connectivity between the upper and lower concrete layers, which are referred as wythes, is provided using shear connectors. Different materials like polystyrene, fiberglass and rock wool are used for insulation layer. The expanded polystyrene (EPS) and extruded expanded polystyrene (XEPS) are the most common materials which are used for the insolation foams [START_REF] Hee | Composite behavior of a novel insulated concrete sandwich wall panel reinforced with GFRP shear grids: effects of insulation types[END_REF]. Dayton superior delta tie, Thermomass rod, standard steel C-clip and universal Teplo tie are some of various types of shear connectors [START_REF] Naito | Precast/ prestressed Concrete Experiments Performance on Non-load Bearing Sandwich Wall Panels[END_REF]. Without shear connectors, there is no connection between concrete wythes and insulation layer. In this case, the FICSP has non-composite action and slip occurs at the interfaces between concrete wythes and insulation layer when it is bent. When a FICSP has sufficient number of shear connectors, its constituent layers will act together as a single unit to resist bending. In this case, it will have a fully composite action and no slip occurs at the interface between concrete wythes and insulation layer. It has been observed that generally the composite action of FICSPs is somewhere between the fully composite action and the non-composite action, and is called the partial composite action.

The mechanical behavior of FICSP structures is to some extent similar to traditional sandwich structures which are fabricated by fastening two thin but stiff face-sheets to a thick but lightweight core. In both of them, there is a large difference between the mechanical properties of core and upper/lower face-sheets. Different laminated theories are presented until now for the analysis of traditional sandwich structures. A finite element based on the variables separation method was presented by Vidal et al. [START_REF] Vidal | Thermo-mechanical analysis of laminated composite and sandwich beams based on a variables separation[END_REF] for the thermo-mechanical analysis of sandwich beams. Vidal and Polit [START_REF]A refined sinus plate finite element for laminated and sandwich structures under mechanical and thermomechanical loads[END_REF] proposed a 6-node triangular finite element for analysis of multilayered composite and sandwich plates subjected to mechanical and thermomechanical loads. The formulation of this reference is based on the sinus model with layer refinement and it includes the transverse normal effects. An efficient finite element method based on global-local theory was introduced by Lezgy-Nazargah and colleagues [START_REF] Lezgy-Nazargah | Efficient coupled refined finite element for dynamic analysis of sandwich beams containing embedded shear-mode piezoelectric layers[END_REF][START_REF] Lezgy-Nazargah | A refined high-order global-local theory for finite element bending and vibration analyses of the laminated composite beams[END_REF][START_REF] Lezgy-Nazargah | Coupled refined layerwise theory for dynamic free and forced responses of piezoelectric laminated composite and sandwich beams[END_REF] for bending and vibration analyses of laminated composite and sandwich beams. A large variety of plate theories including classical, higher order, zigzag, layerwise, and mixed theories were described and assessed by Carrera and Brischetto [START_REF] Carrera | A survey with numerical assessment of classical and refined theories for the analysis of sandwich plates[END_REF] to evaluate the bending and vibration of sandwich structures. By adding a zig-zag function to displacement fields of known theories, Brischetto et al. [START_REF] Brischetto | Improved bending analysis of sandwich plates using a zig-zag function[END_REF] developed higher order theories for bending analysis of sandwich flat panels. A very flexible variable kinematics modeling T technique based on sublaminate Generalized Unified Formulation was proposed by D'Ottavio [START_REF] D'ottavio | A sublaminate generalized unified formulation for the analysis of composite structures[END_REF] for the analysis of laminated composite and sandwich structures. Choe et al. [START_REF] Choe | An efficient approach to investigate the postbuckling behaviors of sandwich structures[END_REF] investigated the post-buckling behavior of sandwich structures. They used a novel one-dimensional layer-wise model using Euler-Bernoulli beam theory in the skins and higher-order kinematics in the core. Shen et al. [START_REF] Shen | Accurate predictions of bending deflections for soft-core sandwich beams subject to concentrated loads[END_REF] studied bending deflections of soft-core sandwich beams subject to concentrated loads. Two correction factors to the classical deflection formula were proposed by them based on the higher-order theory approach. By expanding the displacement field into Fourier series, Huang et al. [START_REF] Huang | A fourier-related double scale analysis on the instability phenomena of sandwich plates[END_REF] studied the instability phenomena of sandwich plates. Saoud and Grognec [START_REF] Saoud | Post-buckling analysis of elastoplastic sandwich columns by means of an enriched 1D finite element model[END_REF] introduced an advantageous 1D finite element model to analyze the post-buckling behavior of sandwich beams. They modeled the thin skin layers as Timoshenko-Reissner beams. In this reference, the complex behavior of core layer is represented by specific kinematics involving hyperbolic functions. The buckling behavior of two-layer shear-deformable beams with partial interaction was studied by Grognec et al. [START_REF] Grognec | Exact buckling solution for two-layer Timoshenko beams with interlayer slip[END_REF]. They considered Timoshenko kinematic hypotheses for both layers. For representing the shear connection, they assumed a continuous relationship between the interface shear flow and the corresponding slip. Challamel et al. [START_REF] Challamel | Out-of-plane behavior of partially composite or sandwich beams by exact and finite element methods[END_REF] studied the out-of-plane vibrations of composite and sandwich beams with interlayer slip for general boundary conditions by exact and finite element methods. For more details about various theories and numerical models on sandwich beams structures, the interested readers can refer to [START_REF] Hu | Review and assessment of various theories for modeling sandwich structures[END_REF][START_REF] Sayyad | Bending, buckling and free vibration of laminated[END_REF][START_REF] Mf | A review on plate and shell theories for laminated and sandwich structures highlighting the finite element method[END_REF][START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF]. However, most of these available laminate theories are not appropriate for the structural analysis of FICSP systems. In contrast to traditional sandwich structures which are made from metallic face sheets, FICSPs consist of concrete wythes. Failure mode in traditional sandwich structures is local buckling of face-sheets while the concrete cracking/crushing is the dominant failure modes in FICSPs. Traditional sandwich structures have a relatively high composite action and can be considered in some extent as a "fully composite" structure. In FICSP structures, stiffness and spacing of shear tie connectors will cause partial composite action to occur between the foam and concrete layers.

Benayoune et al. [START_REF] Benayoune | Structural behavior of eccentrically loaded precast sandwich panels[END_REF] studied the structural behavior of FICSPs made of steel truss shear connectors experimentally. Six full-scale FICSP with variable slenderness ratio were tested under eccentric loads till failure. Woltman et al. [START_REF] Woltman | Investigation of various GFRP shear connectors for insulated precast concrete sandwich wall panels[END_REF] investigated the behavior of FICSP made of Glass fiber-reinforced polymer shear connectors via experimental tests. Naito et al. [START_REF] Naito | Precast/ prestressed Concrete Experiments Performance on Non-load Bearing Sandwich Wall Panels[END_REF] studied experimentally the blast resistance of FICSP structures. Flexural behavior of FICSPs with FRP shear connectors was studied numerically and experimentally by Chen et al. [START_REF] Chen | Experimental investigation and finite element analysis of flexural behavior of insulated concrete sandwich panels with FRP plate shear connectors[END_REF]. Carbonari et al. [START_REF] Carbonari | Flexural behaviour of lightweight sandwich panels composed by concrete and EPS[END_REF] studied experimentally the flexural behavior of FICSPs containing inclined shear connectors under different boundary conditions. Using 2D finite element model and considering both geometrical and material nonlinearities, Benayoune et al. [START_REF] Benayoune | Response of pre-cast reinforced composite sandwich panels to axial loading[END_REF] studied the axial strength capacity of FICSP structures. They also carried out full scale tests in order to validate their finite element model. Bai and Davidson [START_REF] Bai | Analysis of partially composite foam insulated concrete sandwich structures[END_REF] introduced an analytical solution for the analysis of simply supported FICSP structures in elastic regime. They considered the effects of both partial composite action and shear deformation of insulation layer in formulations. With considering the nonlinear behavior of foam, concrete and rebar, Kang [START_REF] Kang | Composite and non-composite behaviors of foam-insulated concrete sandwich panels[END_REF] used a 3D nonlinear finite element model (ABAQUS) for simulating the structural behavior of FICSPs. In [START_REF] Kang | Composite and non-composite behaviors of foam-insulated concrete sandwich panels[END_REF], the multi-point constraints method is used for modeling the shear connectors.

The review of open literature shows that most of research works on FISCP structures are experimental studies. The 2D and 3D finite element models introduced for the structural analysis of FISCPs are computationally expensive methods and cannot be used for the practical design purposes. To fill this gap, a new displacement-based materially nonlinear 1D finite element model with low number of degrees of freedom (DOFs) is introduced for the flexural analysis of FISCP structures. Thanks to introducing a Sine function in the assumed displacement fields, the presented finite element model takes into account the transverse shear deformations induced in concrete wythe and middle insolation layers. The partial composite action of FISCP has been modeled based on spring layer model concept. Shear correction factor is not needed in the presented nonlinear formulation and the zero conditions of the transverse shear stress on the external plane of concrete wythes are fulfilled exactly. Max stress failure criterion is used for the detection of failure in insulation layer. Elasto-plastic model is used for modeling the material nonlinearity of steel rebar. The material acting manner of concrete after tensile crack is modeled via smeared crack method while the elasto-plastic model is employed for capturing the compressive behavior of concrete. A standard 1D displacement-based nonlinear finite element model is derived for solving the governing differential equations. The developed element has three nodes and 22 DOFs. The full-scale experimental data, the 3D finite element results and other analytical approaches have been used to validate the present formulation.

Coordinate system and kinematic

As shown in Fig. 2a, the considered FISCP structure is made of three main constituent elements: top concrete wythe, insulation layer and bottom concrete wythe. In this figure, the chosen coordinate systems for constituent elements are also shown. The coordinate systems are taken so that the coordinate plane y-z lies on the cross-section symmetry plane of the FISCP structure. Each constituent element of FISCP is divided into a series of sub-layers. Although the material property and strain at each constituent element are independent on the number of sub-layers, the sub-layer discretization is required for the purpose of materially nonlinear analysis. As shown in Fig. 2b, longitudinal steel rebar embedded in top and bottom concrete wythes are treated as being steel layers with equivalent thickness t A b / s s = (where A s is the total cross-section area of steel rebar in each wythe). These equivalent steel layers are placed at the centerline of the longitudinal steel rebar.

The following expressions are defined for representing the axial (U i ) and transverse displacement (W i ) in insolation layer and two concrete wythes: where subscript i = t, m, b denote the top concrete wythe, insolation layer and bottom concrete wythe, respectively. Superscript k represents the sub-layer number. u i 0 is the axial displacement at mid-plane of each constituent element. h i denotes the total thickness of insulation/wythe layers. The shear-bending rotation of constituent elements around the y axis is represented with x ( ) i . w denotes the transverse displacement of FICSPs whose variations are assumed to be constant along the thickness direction of FICSP system. From Eq. ( 1), it is seen that the assumed kinematic has only seven global unknown variables:

U x y z u x z w x x h z h x w x x ( , , ) ( ) ( ) sin ( ) ( ) i k i i i i i i i ( ) 0 = + + (1.a)
u b 0 , w, b , u m 0 , m , u t 0 , t .
Using the well-known strain-displacement relations, axial and shear strains at the kth sub-layer of FICSP can be determined from the following relations:

x y z u x x z w x x h z h x x w x x ( ) ( , , ) ( ) ( ) ( )sin ( ) ( ) xx k i i i i i i i i ( ) 0 2 2 2 2 = + + (2.a) x y z z h x w x x ( ) ( , , ) cos ( ) ( ) xz k i i i i i ( ) = + (2.b)
It is easy to find from Eq. ( 2) that the transverse shear stress xz becomes zero at z h /2

t t = and z h /2 b b = .
Thus, the considered kinematic does not need the shear correction factor. The geometrically nonlinear effects are not considered in the present formulation since FICSP structures are not generally designed to operate in this regime.

Due to the partial composite action of shear ties, relative sliding will appear at the interface between the insolation layer and concrete wythes. This phenomenon has been modeled using distributed elastic shear springs which are placed at the interfaces between insulation layer and concrete wythes. The stiffness of these distributed springs may be approximated based on the real tangential stiffness of shear ties and distance between them. The slip at the interface between insulation layer and top concrete wythe can be determined as follow:

s U x z h U x z h u x u x h h h h w x x h x h ( , /2) ( , /2) ( ) ( ) ( /2 /2 / / ) ( ) ( / ) ( ) ( / ) tm t t t m m m t m m t t m t t m m 0 0 = = = = + + (3)
Similarly, the relative sliding at the interface between the bottom concrete wythe and insulation layer can be approximated as:

s U x z h U x z h u x u x h h h h w x x h x h ( , /2) ( , /2) ( ) ( ) ( /2 /2 / / ) ( ) ( / ) ( ) ( / ) mb m m m b b b m b b m m b m m b b 0 0 = = = = + + (4)

Constitutive relations

Insulation layer

Before the occurrence of failure, the following linear elastic behavior is assumed for the kth sub-layer of insulation layer:

E G 0 0 xx k xz k xx m xz m xx k xz k ( ) ( ) ( ) ( ) = (5.a) or D m k m k m k ( ) ( ) ( ) = (5.b)
where G xz m and E xx m are shear and Young's modulus of insulation layer, respectively. The Maximum stress criterion is used for detecting the occurrence of failure in the insulation layer. Based on this criterion, the failure occurs at a point in the insulation layer when the stress components satisfy one of the below conditions [START_REF] Lezgy-Nazargah | Assessment of refined high-order global-local theory for progressive failure analysis of laminated composite beams[END_REF]:

X ( 0) xx t xx > (6.a) X ( 0) xx c xx < (6.b) S | | xz xz (6.c)
where X c , X t and S xz denote uniaxial compressive strength, uniaxial tensile strength and shear strength in x-z plane, respectively. After failure, the elastic constants of insulation layer degrade to zero [START_REF] Lezgy-Nazargah | Assessment of refined high-order global-local theory for progressive failure analysis of laminated composite beams[END_REF].

Steel rebar

It is assumed that no slip occurs at the interface between steel rebar and concrete wythes. The constitutive relations for the kth equivalent steel layer can be expressed as follow:

E 0 0 xx k xz k s E v xx k xz k ( ) ( ) 2 (1 ) ( ) (
)

s s = + ( 7
)
s and E s are Poisson's ratio and Young's modulus of the steel layers, respectively. The above relation can be written as the following matrix form:

D s k s k s k ( ) ( ) ( ) = (8)
Elasto-plastic model is considered for the material behavior of steel layers after yielding. In this condition, the incremental constitutive equations revert to

D d d s k t s k s k ( ) ( ) ( ) = ( 9
)
where

{ }{ } { } { } D D D D D t s k s k s k T s k T s k ( ) ( ) ( ) f( ) f( ) ( ) f( ) ( ) f( ) s k s k s k s k s k s k s k s k ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) =
and f ( ) is the von Mises yield function. The expression of f ( ) is given in Appendix A.

Concrete wythes

Before cracking/crushing, the following linear elastic behavior is assumed for the kth concrete layer: where c and E c denote Poisson's ratio and Young's modulus of concrete.

E 0 0 xx k xz k c E v xx k xz k ( ) ( ) 2 (1 )
In this study, smeared crack modeling is adopted for representing the behavior of concrete after tensile cracking. According to this constitutive model, the incremental stress-strain relations at a point after concrete cracking revert to:

d d d d 0 0 0 k k E v k k 11 ( ) 13 ( ) 2(1 ) 11 ( ) 13 (
)

c c = + (11.a) or D d d c k t c k c k ( ) ( ) ( ) = (11.b)
where denotes the shear retention factor whose range is between 0 and 1. The subscripts 1 and 3 are direction perpendicular and parallel to the tensile crack, respectively. In this study, the constant value 0.001 = is assumed for the shear retention factor. Interested readers can refer to [START_REF] Lezgy-Nazargah | An efficient materially nonlinear finite element model for reinforced concrete beams based on layered global-local kinematics[END_REF] for more details about the proper selection of shear retention factor.

Concerning the material behavior of concrete after compressive crushing, elasto-plastic model is employed.

Finite element formulation

The standard finite element method is used for solving the governing differential equations of the FICSPs. By adopting suitable interpolation schemes for the unknown displacement field variables, the finite element formulation is derived from the following virtual work principle: A i , P ¯i and P ¯iV denote the cross-section area, surface and body force vectors, respectively. i = t, m, b and k ¯is the equivalent stiffness of the shear ties. l e and V e denote the length and volume of elements, re- spectively. e is the total number of elements.

Since the unknown field variables of the FICPS are only function of x coordinate, a 1D finite element formulation is adopted. The variations of transverse displacement w(x) within each element are approximated using Hermite shape functions. Other displacement field variables are approximated using quadratic Lagrangian shape functions. Thus, one can write the vector of displacement variables with respect to the nodal displacement vector as below:

u Nu e = ( 13
)
where It should be pointed here that the linear or cubic Lagrangian shape functions may be used for interpolating of u i 0 and i (i = t, m, b). However, previous published results [START_REF] Vidal | A family of sinus finite elements for the analysis of rectangular laminated beams[END_REF][START_REF] Lezgy-Nazargah | NURBS-based isogeometric analysis of laminated composite beams using refined sinus model[END_REF][START_REF] Beheshti-Aval | A refined sinus finite element model for the analysis of piezoelectric-laminated beams[END_REF] show that the selection of the same order of interpolation for w x , and i in the expression of transverse shear strain ( ) i 13 , prevents from the appearance of shear locking phenomenon.

u u w u u { } T b b m m t t 0 0 0 = , N N N N [ ]
The kinematic relations considered for the constituent elements of the FICSP (Eqs. ( 1)) can be written as the following matrix form: By substituting Eq. ( 13) into Eq. ( 14), the vector of global displacements U i can be written with respect to the nodal displacement vector u e :

U A u t t = (14.a) U A u m m = (14.b) U A u b b = (14.c) where { } U U W i i k i k T ( ) ( ) = , and ( ) z
U A Nu u e e i i i = = ( 15
)
The strain components of FICSP can be rewritten as the below matrix form: Based on Eqs. ( 3) and ( 4), the slips appeared at the interface between the concrete wythe and insulation layer can be written as the following matrix form:

C u t t = ( 16
( ) By introducing the interpolation functions (Eq. ( 13)) in the virtual work principle (Eq. ( 12)), the element load vector is obtained as

A u s 0 0 1 1 tm h h h h d dx h h tm spring 2 2 m t t m m t = + = (18.a) ( ) A u s 1 1 0 0 mb h h h h d dx h h mb
F P P dz dx d d dx d d ¯¯i e nc z z i T i V i T 1 1 1 1 k k 1 = + + + + ( 20
)
After introducing the constitutive relations (Eqs. (5)(11)) as well as the strain-nodal variables relations (Eq. ( 17)) in the virtual work principle, the following expression is obtained for the element stiffness matrix: n t , n m , n b , n sb and n st are the number of sub-layers of top concrete wythe, insulation layer, bottom concrete wythe, top equivalent steel layer, and bottom equivalent steel layer, respectively. is the local coordinate of FICSP elements. Summing the element contributions, the following final nonlinear system of equations are obtained from the application of the virtual work principle:

K D D D D D dz dx d d dz dx d d dz dx d d dz dx d d dz dx d d k dx d d k dx d d ē n z z t T c k t n z z m T m k m n z z b T c k b n z z t T s k t n z z b T s k b tm T tm mb T mb 1 1 ( ) 1 1 ( ) 1 1 ( ) 1 1 ( ) 1 1 ( ) 1 1
K u F [ ]{ } { } T T = ( 22
)
where

K K F F u u [ ] , { } , { } T T e e
e e e e

= = =

Since the nonlinear behavior of insolation layer, steel and concrete is considered in the formulation, the above system of equations are nonlinear and should be solved using the usual incremental-iterative methods. By assuming that the estimate solution at the ith iteration is u { } i , the solution at the next step can be estimated as below using the first-order Taylor series:

u u u { } { } { } i i m i 1 ( ) 1 = + + + (23)
In the above equation,

u { } m i ( ) 1 +
is the mth estimation to the incremental vector displacement u { } i 1

+ . It can be obtained from the resolution of the linear system:

K u F [ ] { } { } T T m i m i i ( 1) 1 ( ) 1 1 = + + + ( 24
)
where

K D D D D D dz dx d d dz dx d d dz dx d d dz dx d d dz dx d d k dx d d k dx d d [ ] T m i e n z z t T c k t n z z m T m k m n z z b T c k b n z z t T s k t n z z b T s k b tm T tm mb T mb ( 1) 1 1 1 ( ) 1 1 ( ) 1 1 ( ) 1 1 ( ) 1 1 ( ) 1 1

Numerical results

In this section, some FICSPs were analyzed until failure stage by means of the presented 1D nonlinear finite element model. In order to assess the accuracy of the present numerical model in predicting the structural responses of FICSP systems, the obtained numerical data were compared with the full-scale experimental data, results of the 3D finite element (ABAQUS) analysis and the analytical results reported in open literature.

Example 1

A simply supported FICSP with span 3.05 m was analyzed by means of the present nonlinear finite element model. It is submitted to a uniform pressure. The geometrical characteristics and cross-section properties of the considered FICSP are shown in Fig. 4. The insulation layer is made of XEPS material with E MPa 40.9

=

. Table 1shows the mechanical properties assumed for the longitudinal steel rebar and concrete wythes. Three different values k ¯= 952 kN/m 3 , 8490 kN/m 3 , and 374,678 kN/m 3 are assumed for the stiffness of shear ties.

Due to the structural symmetry of the considered FICSP system, only the half of it was modeled. Based on the convergence mesh study, the FICSP was axially discretized into 20 elements with equal length. The applied boundary conditions at the support of FICSP model is w 0 = . At the middle of FICSP system, the considered boundary conditions are

w u u u 0 x b m t b m t , 0 0 0 = = = = = = = .
The cross-section of FICSP was also divided into 20 sub-layers. Seven sub-layers were assumed for each concrete wythe while the insolation layer was divided into six sublayers. As shown in Fig. 5, the selected number of sub-layers is adequate for a nonlinear finite element analysis. Although the kinematics relation given in ( 1) is written individually for each constituent layer (and not for sub-layers), the sub-layer discretization is also needed due to the consideration of nonlinear behavior of material. In contrast to concrete layer, the foam insulation layer does not experience significant stress and it will remain in the elastic regime. To this reason, the lower number of sub-layer division is assumed for it. For different values of shear ties' stiffness, the mid-span deflection of FICSP under different pressure state is depicted in Fig. 6. The pressure-slip diagram at the end of the FICSP is also shown in Fig. 7. The depicted graphs of Figs. 6 and7 indicate that the stiffness of shear ties has a remarkable effect on the structural responses of FICSP systems. It should be pointed here that the FICSP system of this example has previously been investigated by Naito et al. [START_REF] Naito | Precast/ prestressed Concrete Experiments Performance on Non-load Bearing Sandwich Wall Panels[END_REF] via experimental tests. For each FICSP system with a distinct k ¯, two similar specimens were fabricated and tested by Naito and colleagues in their laboratory. The results corresponding to these two specimens are shown in depicted graphs of Figs. 6 and7. The considered FICSPs of the present example have also been analyzed by Bai and Davidson [START_REF] Chen | Experimental investigation and finite element analysis of flexural behavior of insulated concrete sandwich panels with FRP plate shear connectors[END_REF] using an analytical elastic model. In Figs. 6 and7, the results of these researchers are also shown and compared with the present results. It can be seen that the present finite element model gives results which are in good agreement with reported experimental data. In comparison to the analytical results of Bai and Davidson, the results of present model are more close to the experimental results. In contrast to the analytical model of Bai and Davidson which assumes a linear elastic behavior for materials, the present finite element model considers the nonlinear behavior of all constituent elements of the FICSPs.

Example 2

The structural response of a FICSP system till collapse stage is evaluated in this section using the proposed nonlinear finite element model. The geometrical and structural details of considered FICSP system are shown in Fig. 8. The FICSP is simply supported and is under the action of a uniform pressure. Table 2 gives the mechanical properties corresponding to longitudinal steel rebar and concrete wythes. The insulation layer is made of XEPS material. The distributed equivalent stiffness of shear ties is approximated as k ¯= 28,310 kN/m 3 [START_REF] Naito | Evaluation of shear tie connectors for use in industrial concrete sandwich panels[END_REF].

The pressure-deflection graph of the FICSP system predicted by the nonlinear finite element model of present study is depicted in Fig. 9. In this figure, the static test results performed at the University of Missouri [START_REF] Newberry | Simulation of prestressed concrete sandwich panels subjected to blast loads[END_REF], and the analytical results of Kang [START_REF] Kang | Composite and non-composite behaviors of foam-insulated concrete sandwich panels[END_REF] based on the ACI standard design code are also shown for more comparison. It is seen that the results obtained from the nonlinear finite element model of present study are in good agreement with the reported experimental data. It can be also observed that the analytical methods based on ACI standard design code can only predict the fully composite and non-composite behavior of FICSPs systems. Thus, such methods are not able to model the partial composite behavior of FICSPs. As seen in Fig. 9, the fully composite and non-composite behavior make, respectively, the upper and lower boundaries of the capacity curve of the FICSPs.

In the following, the accuracy of presented finite element formulation in predicting the detailed (local) structural behavior of FICSP's components is evaluated. To this aim, variations of mid-span axial stress at upper and lower concrete wythes against the mid-span deflection are shown in Fig. 10. Fig. 11 shows distributions of axial (inplane) stress in upper steel rebar at the ultimate pressure stage. Finally, axial stress-strain history on insulation layer at the middle of span is also shown in Fig. 12.

It can be seen from Figs. 10-12 that there exist a good correlation between the experimental and present results in predicting the local structural behavior of FICSP systems. Concerning the numerical results corresponding to the axial stress in upper steel rebar (Fig. 11), it can be seen that the slope in the middle part of graph is not the same as the experiment. This discrepancy may be due to neglecting the structural effect of WWR in the 1D finite element formulation of the present model. Indeed, that part of wire mesh which is transverse to the longitudinal rebar cannot be modeled in the framework of a 1D finite element model. Assuming a perfect bonding between the longitudinal rebar and surrounding concrete may be other reason which leads to the difference between the experimental and present results.

The above obtained numerical results confirm that the introduced finite element formulation of the present study is an appropriate tool for global-local structural analysis of FICSP systems.

Example 3

As a final example, a cantilever FICSP system till failure stage was analyzed via the nonlinear finite element formulation introduced in Section 4. The length of the considered FICSP is 4 m and it is under the action of a distributed pressure. The cross-section details and the mechanical properties of the constituent elements of the considered FICSP system of the present example are the same as example 1. The distributed equivalent stiffness of shear ties which provides the connectivity between the upper and lower concrete wythes is assumed to be k ¯= 0 kN/m 3 (non-composite action), 123 10 3 × kN/m 3 (partially composite action), and 123 10 6 × kN/m 3 (fully composite action). In Fig. 13, variations of the reaction force against the tip deflection of FICSP are depicted for three different interface conditions (i.e. partially composite, non-composite, and fully composite actions). It is seen again that the ultimate capacity of FICSP system is extremely dependent on the stiffness of shear ties. In Fig. 13, the present results were compared with nonlinear 3D finite element results. ABAQUS software was employed for carrying out a 3D finite element analysis. Concrete whythes and insulation layer are modeled using 8-node brick elements (C3D8R). Longitudinal steel rebar as well as the WWRs are modeled × × brick elements. Longitudinal steel rebar and WWR were divided into 50 and 10 elements, respectively. Since longitudinal steel rebar and the WWRs are embedded in the concrete layers, a perfect bond is assumed at the interface between them. Due to the occurrence of slip at the interfaces between the concrete and insulation layers, assumption of a perfect bond is not correct for describing the interfacial behavior between them. In order to capture the partial composite action in the 3D model of considered FICSP, traction-separation laws (cohesive behavior) are used for describing the interfacial behavior between the concrete wythes and insulation layer. Tractionseparation laws describe the interaction between two surfaces by defining a relative displacement at each contact point. For more details about interfacial traction-separation law, interested readers can refer to [START_REF] Lezgy-Nazargah | The effects of different FRP/concrete bond-slip laws on the 3D FE modeling of retrofitted RC beams -A sensitivity analysis[END_REF]. The ideal elasto-plastic and concrete damage plasticity (CDP) models are respectively employed for representing the material nonlinearity of steel and concrete. CDP uses concepts of isotropic damaged elasticity in combination with isotropic tensile and compressive plasticity to represent the inelastic behavior of concrete. The material parameters used for the CDP model are given in Table 3. The expressions corresponding to yield function of the CDP model are given in Appendix A. The present finite element model estimates the maximum load capacity of the cantilever FICSP in non-composite, partially composite and fully composite conditions as 7.74 kN, 14.80 kN, and 21.92 kN, respectively. In case of 3D finite element analysis, the aforementioned values are 7.63 kN, 14.15 kN and 22.05 kN. The percent error between the ABAQUS and present results is less than 4.5. It is worthy to be pointed here that the converged mesh of ABAQUS model has about 109,500 DOFs while this value for the 1D finite element formulation of the present study is only 288 DOFs. This comparison proves the 1D finite element formulation introduced in the present study is an efficient tool for nonlinear structural analysis of FICSP systems.

Through-the-thickness distributions of stress and displacement components of FICSP at a section with distance 1 m from the fixed end are depicted in Figs. 15-17. Note that the transverse shear stress was calculated at the post-processing level by integrating from equilibrium equations. For comparison purposes, the numerical results are also given in Table 4. It can be observed from depicted graphs of Figs. 15-17 that the present finite element formulation has a good accuracy in predicting the local structural behavior of FICSP structures. The small discrepancies between the present and 3D finite element results is mainly due to difference in employed constitutive law used for modeling the nonlinear demeanor of concrete materials. In can be observed from Figs. 15-17 that the interface properties between concrete and foam (i.e. k ¯) has a significant influence on the local structural responses of FICSP systems. For low values of k ¯, significant slip occurs at the interface between concrete wythe and foam layer. This slip reduces with increasing of shear ties' stiffness so that for k ¯123 10 6 = × kN/m 3 Fig. 14. The refined mesh used for 3D finite element model (ABAQUS).

Table 3

The values of parameters used for the CDP model. no slip occurs at the interface between layers. When k ¯123 10 6 = × kN/ m 3 , the maximum transverse shear stress situates at the foam layer. With reducing of k ¯, the value of transverse shear stress in the foam layer reduces. In case of non-composite action, no shear stress appears in the foam layer and maximum shear stress occurs at the center of concrete wythes. Figs. [START_REF] Grognec | Exact buckling solution for two-layer Timoshenko beams with interlayer slip[END_REF][START_REF] Challamel | Out-of-plane behavior of partially composite or sandwich beams by exact and finite element methods[END_REF][START_REF] Hu | Review and assessment of various theories for modeling sandwich structures[END_REF] show that the through-the-thickness distribution of axial stress is extremely affected by the interface characteristic between foam and concrete. In case of fully composite action, the considered FICSP system is bent as an integral cantilever beam. Top concrete wythe experiences tensile stresses while the lower wythe undergoes compressive stresses. With reducing of k ¯, every concrete layer is inclined to bend individually.

The tensile crack patterns developed at the concrete wythes are shown in Fig. 18 for non-composite interface condition (k ¯= 0). In this figure, the tensile crack patterns corresponding to 3D finite element analysis are also shown. It is seen that the patterns of tensile cracks predicted via the present 1D finite element are similar to those predicted by ABAQUS.

Conclusions

A simple and efficient displacement-based 1D finite element formulation was introduced for predicting the nonlinear structural behavior of FICSP structures. Non-polynomial expressions are considered for representing the axial displacement of concrete wythes and insolation layer. The assumed displacement field takes into account the transverse shear deformations induced in the constituent components of FICSP structures. Without using any shear correction factor, the presented formulation guarantees the zero conditions of transverse shear stress on exterior plane of concrete wythes. The partial composite action of shear connectors is represented using the concept of spring layer model. In comparison to available 2D and 3D finite element models, the proposed approach is computationally low cost.

CDP model

The yield function of the CDP model can be written as the following form [START_REF] Lezgy-Nazargah | The effects of different FRP/concrete bond-slip laws on the 3D FE modeling of retrofitted RC beams -A sensitivity analysis[END_REF]: where b0 and c are the biaxial and uniaxial compressive strength of the concrete, respectively. The function ) pl is given as where K c is the ratio of the tensile to the compressive meridian.

f P ( ¯, ) 1
( ) ¯( ) ¯( ) (1 ) (1

Fig. 1 .

 1 Fig. 1. Main constituent elements of the typical FISCP structures.

Fig. 2 .

 2 Fig. 2. Coordinate system and typical geometry of FICSPs: (a) side view, (b) cross-section.

  for N1 , N3 , N2 , u e 1 , u e 3 and u e 2are given in Appendix A. The considered three-nodded element and corresponding nodal variables are shown in Fig.3.

Fig. 3 .

 3 Fig. 3. The three-node FICSP element and related degrees of freedom.

  of Eq. (13) into the above equations gives:

+

  Based on the history of stress at each sub-layer element,

Fig. 4 .

 4 Fig. 4. Geometrical characteristics of considered FICSP of example 1: (a) loads and boundary conditions (b) cross-section details-unit (mm).

Fig. 5 .

 5 Fig. 5. Effect of the number of cross-section sub-layers on the mid-span deflection of FICSP (k ¯= 8490 kN/m 3 ).

Fig. 6 .

 6 Fig. 6. Mid-span deflection of FICSP versus the applied pressure: (a) k ¯= 952 kN/m 3 , (b) 8490 kN/m 3 , (c) 374,678 kN/m 3.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Pressure-slip history at the end of the FICSP: (a) k ¯= 952 kN/m 3 , (b) 8490 kN/m 3 . Fig. 8. Geometrical characteristics and cross-section details of FICSP: unit (mm).

Fig. 10 .

 10 Fig. 10. Stress-deflection history at the mid-span of FICSP: (a) upper concrete wythe, (b) lower concrete wythe.

Fig. 11 .

 11 Fig. 11. Variations of stress in upper steel rebar at the ultimate load stage.

Fig. 12 .

 12 Fig. 12. Axial stress-strain history on insulation layer at mid-span.

Fig. 13 .

 13 Fig. 13. Tip deflection of cantilever FICSP versus the reaction force.

using 2 -

 2 node truss elements. The refined mesh used for 3D finite element model is shown in Fig.14. The number of elements used for meshing each concrete whyth is 6 30 50 × × . The mesh of insulation layer includes5 30 50

Fig. 15 .

 15 Fig.15. Through-the-thickness distributions of displacement and stress components for cantilever FICSP: fully composite case.

Fig. 16 .

 16 Fig. 16. Through-the-thickness distributions of displacement and stress components at cantilever FICSP: partially composite case.

Fig. 17 .

 17 Fig. 17. Through-the-thickness distributions of displacement and stress components at cantilever FICSP: non-composite case.

  It has only 22 DOFs per element. The presented finite element can accurately model the fully composite, non-composite and partially composite actions of FICSP structures. The obtained numerical results show that the formulation introduced in this study can estimate both global (deflection, flexural capacity) and local (through-the-thickness distributions of stress components) structural behaviors of FICSP structures at different loading stages with enough accuracy. Future studies are toward the development of efficient plate/shell elements for predicting the structural responses of FICSP structures till failure stage.

Fig. 18 .

 18 Fig. 18. Patterns of tensile crack at cantilever FICSP: (a) present finite element model, (b) 3D finite element model-(k ¯= 0).

P

  ¯is the hydrostatic pressure stress; Q ¯denotes the Mises equivalent effective stress; . is the Macauley bracket; pl and t are equivalent plastic strains in compression and tension. Parameter is defined as the follow:

  )

  ¯c and ¯t are, respectively, the effective compressive and tensile cohesion stresses. The parameter can be calculated as the below:

Table 1

 1 Mechanical characteristics of concrete wythes and steel rebar -Example 1.

	Mechanical property	Concrete wythes	Steel rebar
	Modulus of elasticity (GPa)	36	200
	Yield strength (MPa)	-	480
	Compressive strength (MPa)	49	-
	Tensile strength (MPa)	2.8	-
	Poisson's ratio	0.18	0.3

Table 2

 2 Mechanical characteristics of concrete wythes and steel rebar -Example 2.

	Mechanical property	Concrete wythes	Steel rebar
	Modulus of elasticity (GPa)	24.8	200
	Yield strength (MPa)	-	480
	Compressive strength (MPa)	34	-
	Tensile strength (MPa)	2.8	-
	Poisson's ratio	0.18	0.3

Fig. 9. history of FICSP-example 2.

Table 4

Displacement (in mm) and stress (in MPa) components of the cantilever FICSP-x = 1000 mm. (A.1) y is the uniaxial yield stress.

Matrix formulation

where


Philippe Vidal

Laurent Gallimard

Olivier Polit

come

Multiresolution strategies for the modeling of composite shell structures based on the variable separation method

Multiresolution strategies for the modeling of composite shell structures based on the variable separation method P. Vidal, L. Gallimard, O. Polit

INTRODUCTION

Nowadays, composites are widely used in the industrial field (eg, aerospace, automotive, marine, and medical industries) due to their excellent mechanical properties, especially their high specific stiffness and strength. They are used for the manufacturing of primary structures that have to ensure the integrity of those in service. Moreover, complex geometries such as curved structures can be involved. Thus, accurate knowledge of displacements and stresses is required to design reliable systems. One way consists in considering three-dimensional (3D) modeling. However, due to the complexity of such numerical simulations, it is suitable to represent the problem as a two-dimensional (2D) model, leading to the construction of shell theories. There are two ways to define the approximation of the displacement field. A "pure shell model" can be considered, in which the displacement is associated with the local curvilinear vectors and strain and stress are deduced using differential geometry. [START_REF] Dau | An efficient C 1 finite element with continuity requirements for multilayered/sandwich shell structures[END_REF] Alternatively, the shell-like solid approach [START_REF] Zienkiewicz | The Finite Element Method: The Basis[END_REF] to obtain a shell finite element (FE) is widely used in commercial software, as it is more simple. In this case, the displacement vector is defined in the global Cartesian frame, and Jacobian matrix transformation is used to express strain and stress with respect to the reference frame defined on the middle surface in order to introduce the constitutive law. In this approach, differentiation is simplified, and the curvatures are not directly calculated. [START_REF] Vidal | An efficient finite shell element for the static response of piezoelectric laminates[END_REF] Moreover, the designer has to solve some specific problems: (i) uncertainties in the manufacturing process of such materials and (ii) the optimization of the structure (material and geometry) to improve the performance. These call for the resolution of numerous problems. Hence, the development of efficient computational models for the analysis of shells appears to be of major interest.

According to published research, various theories based on the FE method for composite shells have been developed. In the following, most of the mentioned works refer to the pure shell model. Thus, it is nowadays admitted that two families of models 4 can be identified.

• Equivalent single-layer models (ESL). Here, the classical shell theory (CST/Koiter) and first-order shear deformation theory (Nagdhi) models are commonly used. Assumptions on the strain to derive different shell models can be found in the work of Leissa. [START_REF] Leissa | Vibration of shells[END_REF] The CST suffers from inaccurate results when applied to composites because both transverse shear and normal strains are neglected. The reader can refer to the works of Rao, [START_REF] Rao | Rectangular laminated anisotropic shallow thin shell finite element[END_REF] Jeyachandrabose and Kirkhope, [START_REF] Jeyachandrabose | Explicit formulation of two anisotropic, triangular, thin, shallow shell elements[END_REF] and Qatu and Leissa [START_REF] Qatu | Bending analysis of laminated plates and shells by different methods[END_REF] for discussions on shallow laminated shells. The first-order shear deformation theory is the most popular model due to the possibility of using a C 0 FE; however, it needs shear correction factors, and the transverse normal strain is always neglected (cf the works of Reddy, [START_REF] Reddy | Bending of laminated anisotropic shells by a shear deformable finite element[END_REF] Chakravorty et al, [START_REF] Chakravorty | Finite element free vibration analysis of doubly curved laminated composite shells[END_REF] Hossain et al, [START_REF] Hossain | A finite element formulation for the analysis of laminated composite shells[END_REF] and Asadi et al [START_REF] Asadi | Static and vibration analyses of thick deep laminated cylindrical shells using 3D and various shear deformation theories[END_REF] ). Thus, higher-order shear deformation theories have been developed to overcome these drawbacks. Different kinematics including five, [START_REF] Reddy | A higher-order shear deformation theory of laminated elastic shells[END_REF] seven, [START_REF] Balah | Finite element formulation of a third order laminated finite rotation shell element[END_REF][START_REF] Sgambitterra | A robust three-node shell element for laminated composites with matrix damage[END_REF] or nine parameters [START_REF] Kant | Estimation of interlaminar stresses in fibre reinforced composite cylindrical shells[END_REF] are considered. Note that various theories based on Carrera's unified formulation (CUF) are addressed in the works of Carrera et al. [START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF][START_REF] Carrera | Historical review of zig-zag theories for multilayered plates and shells[END_REF][START_REF] Cinefra | Shell finite elements with different through-the-thickness kinematics for the linear analysis of cylindrical multilayered structures[END_REF] In the equivalent single-layer model context, a simple way to improve the estimation of the mechanical quantities consists in adding one zig-zag function (Murakami) in the expression of the displacement to introduce the slope discontinuity at the interface between two adjacent layers. It allows to describe the so-called zig-zag effect. This has been carried out in the works of Bhaskar and Varadan, [START_REF] Bhaskar | A higher-order theory for bending analysis of laminated shells of revolution[END_REF] Jing and Tzeng, [START_REF] Jing | Refined shear deformation theory of laminated shells[END_REF] Brank, [START_REF] Brank | On composite shell models with a piecewise linear warping function[END_REF] and Ganapathi et al. [START_REF] Ganapathi | Vibration analysis of laminated cross-ply oval cylindrical shells[END_REF] • Layer-wise models (LW). Here, the expression of the mechanical quantities is written over each layer. Some of these works are based on a linear distribution of the in-plane displacements through each layer, without taking into account the transverse normal stress. The transverse displacement is assumed constant across the thickness, such as in the works of Botello et al [START_REF] Botello | A layer-wise triangle for analysis of laminated composite plates and shells[END_REF] and Zinno and Barbero, [START_REF] Zinno | A three-dimensional layer-wise constant shear element for general anisotropic shell-type structures[END_REF] or in each layer, as in the work of Seide and Chaudhuri. [START_REF] Seide | Triangular finite element for analysis of thick laminated shells[END_REF] However, this type of approach fails to predict accurate transverse stresses, unless when using special post-treatments (see, eg, the work of Zinno and Barbero 25 ). Thus, higher-order approaches, taking into account the transverse normal effect, have been developed. The 3D behavior law is used. The second-, third-, and fourth-order expansions are derived in the works of Reddy, [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells: Theory and Analysis[END_REF] Basar and Ding, [START_REF] Basar | Interlaminar stress analysis of composites: layer-wise shell finite elements including transverse strains[END_REF] and Grigolyuk and Kulikov. [START_REF] Grigolyuk | General direction of development of the theory of multilayered shells[END_REF] In this framework, Kulikov and Plotnikova [START_REF] Kulikov | Advanced formulation for laminated composite shells: 3D stress analysis and rigid-body motions[END_REF] have developed the sampling surfaces method (see also the previously mentioned work [START_REF] Cinefra | Shell finite elements with different through-the-thickness kinematics for the linear analysis of cylindrical multilayered structures[END_REF] ). In all the aforementioned works, it should be noted that the number of unknowns depends on the number of layers, affecting the performance in terms of computational cost.

As an alternative, refined models have been developed in order to improve the accuracy of ESL models, avoiding the additional computational cost of LW approaches. On the basis of physical considerations and after some algebraic transformations, the number of unknowns becomes independent of the number of layers. The deduced model can be derived from the CST [START_REF] Soldatos | A unified formulation of laminated composite, shear deformable, five-degrees-of-freedom cylindrical shell theories[END_REF] and the higher-order shear deformation theory with a third-order theory [START_REF] Yasin | An efficient layerwise finite element for shallow composite and sandwich shells[END_REF][START_REF] Cho | Efficient higher-order shell theory for laminated composites[END_REF][START_REF] Shariyat | Non-linear dynamic thermo-mechanical buckling analysis of the imperfect laminated and sandwich cylindrical shells based on a global-local theory inherently suitable for non-linear analyses[END_REF][START_REF] Yasin | An efficient layerwise finite element for shallow composite and sandwich shells[END_REF][START_REF] Shu | Refined theory of laminated shells with higher order transverse shear deformation[END_REF] or the sinus model. [START_REF] Dau | An efficient C 1 finite element with continuity requirements for multilayered/sandwich shell structures[END_REF] The number of parameters varies from 5 to 15. Thus, for all these studies, the computational cost decreases, but the transverse normal effect is introduced in none of them. This is achieved in the subsequent works of Zhen and Wanji [START_REF] Zhen | A global-local higher order theory for multilayered shells and the analysis of laminated cylindrical shell panels[END_REF] and Icardi and Ferrero. [START_REF] Icardi | Multilayered shell model with variable representation of displacements across the thickness[END_REF] Note also the possibility of managing different kinematics through the thickness in the framework of the CUF approach. [START_REF] Carrera | Shell elements with through-the-thickness variable kinematics for the analysis of laminated composite and sandwich structures[END_REF] The aforementioned literature deals with only some aspects of the broad research activity about composite shells. An extensive assessment of the different approaches for various theories and/or FE applications can be found in other works. [START_REF] Kapania | A review on the analysis of laminated shells[END_REF][START_REF] Noor | Assessment of computational models for multilayered composite shells[END_REF][START_REF] Gilewski | A survey of finite element models for the analysis of moderately thick shells[END_REF][START_REF] Yang | Survey of recent shell finite element[END_REF][START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF][START_REF] Reddy | Shear deformation plate and shell theories: from Stavsky to present[END_REF][START_REF] Hohe | Advances in the structural modeling of elastic sandwich panels[END_REF][START_REF] Qatu | Review of recent literature on static analyses of composite shells: 2000-2010[END_REF][START_REF] Caliri | A review on plate and shell theories for laminated and sandwich structures highlighting the finite element method[END_REF] For the present topics, the fundamental subject on the shear and membrane locking of the shell is not addressed here.

Over the past years, an a priori method called proper generalized decomposition (PGD) [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF] has shown interesting features in the reduction model framework. It has been successfully used in the context of the separation of coordinate variables for different types of composite structures in other works. [START_REF] Savoia | A variational approach to three-dimensional elasticity solutions of laminated composite plates[END_REF][START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF][START_REF] Vidal | Assessment of a composite beam finite element based on the proper generalized decomposition[END_REF][START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF][START_REF] Vidal | Shell finite element based on the proper generalized decomposition for the modeling of cylindrical composite structures[END_REF] Note that other reduction model techniques are available in open literature. A posteriori methods that consist in computing a reduced basis from snapshots using proper orthogonal decomposition or singular value decomposition are not addressed hereafter. The main goal of this work consists in assessing a multiresolution strategy based on PGD to model cylindrical composite shell structures. Classically, the present approach is based on the separation representation of the unknown displacements. In this way, the three displacements are written under the form of a sum of products of 2D polynomials of (𝜉 1 , 𝜉 2 ) and one-dimensional (1D) polynomials of z. A piecewise fourth-order Lagrange polynomial is chosen for the three unknown z-functions. As far as the variation of the in-plane coordinates is concerned, a 2D eight-node quadrilateral FE is used for the three unknown (𝜉 1 , 𝜉 2 ) functions. It should be noted that the LW unknown functions of z are global for the whole shell. Finally, the deduced nonlinear problem implies the resolution of two linear problems alternatively. This process yields 2D and 1D problems in which the number of unknowns is smaller than that in a classical layer-wise approach. The interesting feature of this approach lies in the possibility of having a higher-order z-expansion and refining the description of the mechanical quantities through the thickness without increasing the computational cost. This is particularly suitable for modeling composite structures.

This approach can be advantageously used to perform numerous computations in which some geometric or material parameters change, with a low computational cost. It can be achieved in an optimization or a reliability process for instance. Once a first computation is carried out for a fixed set of parameters, it is possible to reuse the 2D in-plane (𝜉 1 , 𝜉 2 ) functions previously built. Thus, only new 1D problems have to be solved for a new set of parameters. This method is an extension of the so-called preliminary stage introduced in the works of Bussy et al [START_REF] Bussy | The large time increment method for numerical simulation of metal forming processes[END_REF] and Boisse et al [START_REF] Boisse | A new approach in non-linear mechanics: the large time increment method[END_REF] to solve nonlinear problems with a time-space decomposition of the solution. If needed, only few 2D functions are computed to obtain an accurate solution.

We now outline the remainder of this paper. First, the shell definition, the differential geometry, and the mechanical formulation are described. Then, the principles of PGD are given. The particular assumption on the displacements yields a nonlinear problem that is solved using a fixed-point method. The classical greedy algorithm is recalled, and a multiresolution version (1D problem) is given. FE discretization is also described. Finally, numerical tests are performed on L-angle specimens that can be considered as representative structures. The variations of the stacking sequences, the geometry of the shell, and the number of layers are addressed in the multiresolution process. The efficiency of the approach is proved, and the accuracy of the results is assessed by a comparison with the layer-wise LD4 reference solution.

SHELL DEFINITIONS AND DIFFERENTIAL GEOMETRY

A shell  with a middle surface  and a constant thickness e (see Figure 1) is defined by 58

 = { M ∈  3 ∶ ⃗ OM(𝜉, 𝜉 3 = z) = ⃗ Φ(𝜉) + z ⃗ a 3 ; 𝜉 ∈ Ω; - 1 2 e ≤ z ≤ 1 2 e } ,
where the middle surface can be described by a map ⃗ Φ from a parametric bidimensional domain Ω as

⃗ Φ ∶ Ω ⊂  2 -→  ⊂  3 𝜉 = (𝜉 1 , 𝜉 2 )  -→ ⃗ Φ(𝜉). (1)
In Figure 1, the map ⃗ Φ describing the shell middle surface (in gray) and the local basis vectors are presented. The basis vectors ⃗ a i are defined for a point on , and the basis vectors ⃗ g i are defined for a generic point of the shell. For a point on the shell middle surface, the covariant basis vectors defining the tangent plane to the middle surface are usually obtained as follows:

⃗ a 𝛼 = ⃗ Φ(𝜉 1 , 𝜉 2 ) ,𝛼 ; ⃗ a 3 = ⃗ a 1 × ⃗ a 2 ||⃗ a 1 × ⃗ a 2 || , ( 2
)
where ⃗ a 3 is the unit normal vector to the surface  (see Figure 1). In Equation ( 2) and further on, Latin indices i, j, … take their values in the set {1, 2, 3}, whereas Greek indices 𝛼, 𝛽, … take their values in the set {1, 2}.

a 𝛼𝛽 = ⃗ a 𝛼 .⃗ a 𝛽 a 𝛼𝛽 = ⃗ a 𝛼 .⃗ a 𝛽 b 𝛼𝛽 = ⃗ a 𝛼,𝛽 .⃗ a 3 b 𝛽 𝛼 = ⃗ a 𝛽 .⃗ a 3,𝛼 . (3)
For a generic point of the shell, covariant basis vectors must be defined, and we have

⃗ g 𝛼 = ⃗ OM(𝜉, z) ,𝛼 = ( 𝛿 𝛽 𝛼 -zb 𝛽 𝛼 ) ⃗ a 𝛽 = 𝜇 𝛽 𝛼 (z) ⃗ a 𝛽 and ⃗ g 3 = ⃗ a 3 , (4)
where 𝛿 𝛽 𝛼 is the Kronecker symbol and b 𝛽 𝛼 is the mixed form of the second fundamental form. This basis ⃗ g i , illustrated in Figure 1, must be used to define quantities for any point of the shell. The form 𝜇 𝛽 𝛼 (z) introduced in Equation ( 4) defines the transport from the shell middle surface to any point of the shell and is associated with the curvature variation along the thickness direction z of the shell. The inverse tensor of the mixed tensor 𝜇 𝛽 𝛼 is denoted m 𝛽 𝛼 and is defined as

m 𝛽 𝛼 = (𝜇 -1 ) 𝛽 𝛼 = 1 𝜇 { 𝛿 𝛽 𝛼 + z ( b 𝛽 𝛼 -2H𝛿 𝛽 𝛼 )} , ( 5
)
where we have introduced the determinant of the mixed tensor 𝜇 = det (𝜇 𝛽 𝛼 ) = 1 -2H𝜉 3 + (𝜉 3 ) 2 K and the invariants of the second fundamental form H = 1 2 tr(b 𝛽 𝛼 ) and K = det (b 𝛽 𝛼 ). Finally, the surface element d and the volume element d are given by

d = √ a d𝜉 1 d𝜉 2 d = 𝜇 ddz, ( 6
)
where a is the determinant of the first fundamental form a 𝛼𝛽 . The geometry of a shell can also be defined using a contravariant or mixed forms. Furthermore, covariant and contravariant differentiation involving Christoffel symbols is not detailed here, and readers can refer to the book of Bernadou. 58

REFERENCE PROBLEM DESCRIPTION

The definition of the strain field

For geometrically linear elastic analysis, the components of the strain tensor 𝜀 ij expressed in the contravariant basis ⃗ a i are obtained through covariant differentiation, denoted |, as follows:

𝜀 = 𝜀 i𝑗 (⃗ a i ⊗ ⃗ a 𝑗 ) with 2 𝜀 𝛾𝜆 = m 𝛽 𝜆 (u 𝛾 | 𝛽 -b 𝛾𝛽 u 3 ) + m 𝛼 𝛾 (u 𝜆 | 𝛼 -b 𝜆𝛼 u 3 ) 2 𝜀 𝛾3 = u 𝛾 | 3 + m 𝛼 𝛾 ( u 3,𝛼 + b 𝜆 𝛼 u 𝜆 ) 𝜀 33 = u 3,3 . (7)
The mixed tensor m 𝛽 𝜆 carries out the transport from any point of the shell to the shell middle surface, that is, from ⃗ g i to ⃗ a i .

Constitutive relation

The stress tensor is obtained from the strain tensor using the constitutive equations. For this purpose, all these tensors must be referred to the covariant and contravariant basis vectors, ⃗ a i and ⃗ a i , respectively, associated with the middle surface of the shell. In case of laminated shells composed of orthotropic plies, this reference frame ensures to consistently take into account the different material orientations of the layers. The tensorial relation is

𝜎 i𝑗 = C i𝑗kl 𝜀 kl with 𝜎 i𝑗 ( ⃗ a i ⊗ ⃗ a 𝑗 ) , C i𝑗kl ( ⃗ a i ⊗ ⃗ a 𝑗 ⊗ ⃗ a k ⊗ ⃗ a l ) , 𝜀 kl ( ⃗ a k ⊗ ⃗ a l ) . ( 8
)
It should be noted that the stress tensor is defined in the covariant reference frame, whereas the strain components are defined in the contravariant frame. If the frame is assumed to be orthonormal, then the covariant and contravariant components are equal, that is, the superscript and subscript are interchangeable. With this assumption and for a shell made of NC perfectly bonded orthotropic layers, the 3D constitutive law of the kth layer is expressed in matrix notation by

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝜎 11 (k) (𝜉, z) 𝜎 22 (k) (𝜉, z) 𝜎 33 (k) (𝜉, z) 𝜎 23 (k) (𝜉, z) 𝜎 13 (k) (𝜉, z) 𝜎 12 (k) (𝜉, z) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ C 11 (k) C 12 (k) C 13 (k) 0 0 0 C 22 (k) C 23 (k) 0 0 0 C 33 (k) 0 0 0 C 44 (k) 0 0 sym C 55 (k) 0 C 66 (k) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝜀 11 (k) (𝜉, z) 𝜀 22 (k) (𝜉, z) 𝜀 33 (k) (𝜉, z) 𝛾 23 (k) (𝜉, z) 𝛾 13 (k) (𝜉, z) 𝛾 12 (k) (𝜉, z) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (9)
ie, 𝝈 (k) = C (k) 𝜺 (k) ,

where we denote C (k) i𝑗 as the 3D stiffness coefficients of layer k, stress vector 𝝈, and strain vector 𝜺.

The weak form of the boundary value problem

The shell is submitted to a surface force density t defined over a subset 𝜕 F of the boundary. We assume that a prescribed displacement u = u d is imposed on Γ D = 𝜕 -𝜕 F . The mechanical problem is based on the variational principle, which is given as follows: find u(M) ∈ U (space of admissible displacements) such that

a(u, 𝛿u) = b(𝛿u), ∀𝛿u ∈ 𝛿U (10) with a(u, 𝛿u) = ∫ ∫ ∫  𝛿𝜺 T 𝝈 d b(𝛿u) = ∫ ∫ 𝜕 F 𝛿u • t d𝜕,
where U is the space of admissible displacements, ie, U = {u ∈ (H 1 ()) 3 ∕u = u d on Γ D }. We have also

𝛿U = {u ∈ (H 1 ()) 3 ∕u = 0 on Γ D }.
For the present work, 𝜕 F is considered as the top or bottom surface of the shell, that is, z = z F with z F = ±e∕2. For the sake of clarity, the body forces are neglected in the developments.

APPLICATION OF PGD TO A CYLINDRICAL SHELL

In this section, the application of PGD is developed in order to carry out numerous shell analyses. This work is an extension of the previous studies on beam, plate, and shell structures. [START_REF] Vidal | Assessment of a composite beam finite element based on the proper generalized decomposition[END_REF][START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF][START_REF] Vidal | Shell finite element based on the proper generalized decomposition for the modeling of cylindrical composite structures[END_REF]

Cylindrical geometry

A cylindrical shell (see Figure 2) is described using the following map:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ x 1 (𝜉 1 , 𝜉 2 ) = R cos ( 𝜉 1 R ) x 2 (𝜉 1 , 𝜉 2 ) = R sin ( 𝜉 1 R ) x 3 (𝜉 1 , 𝜉 2 ) = 𝜉 2 . ( 11
)
Following Equation (3), the nonzero terms for the covariant and mixed forms are

a 11 = a 22 = 1 b 11 = b 1 1 = - 1 R 𝜇 1 1 = 1 + z R m 1 1 = ( 1 + z R ) -1 (12)
and

𝜇 = 1 + z R .

FIGURE 2

Middle surface of the cylindrical shell

The displacement

Let us denote (u

1 (𝜉 1 , 𝜉 2 , 𝜉 3 = z), u 2 (𝜉 1 , 𝜉 2 , 𝜉 3 = z), u 3 (𝜉 1 , 𝜉 2 , 𝜉 3 = z))
as the curvilinear components of the displacement field associated with the contravariant basis vectors ⃗ a i . Let Ω and Ω z be the bidimensional domain associated with the mid-surface of the shell (see Equation ( 1)) and the unidimensional domain associated with the normal fiber, respectively. This displacement solution is constructed as the sum of N products of functions of in-plane coordinates and the transverse coordinate (N ∈ N is the order of the representation), ie,

u = ⎡ ⎢ ⎢ ⎢ ⎣ u 1 (𝜉, 𝜉 3 = z) u 2 (𝜉, 𝜉 3 = z) u 3 (𝜉, 𝜉 3 = z) ⎤ ⎥ ⎥ ⎥ ⎦ = N ∑ i=1 ⎡ ⎢ ⎢ ⎢ ⎣ 𝑓 i 1 (z) v i 1 (𝜉) 𝑓 i 2 (z) v i 2 (𝜉) 𝑓 i 3 (z) v i 3 (𝜉) ⎤ ⎥ ⎥ ⎥ ⎦ = N ∑ i=1 ⎡ ⎢ ⎢ ⎢ ⎣ 𝑓 i 1 (z) 𝑓 i 2 (z) 𝑓 i 3 (z) ⎤ ⎥ ⎥ ⎥ ⎦ • ⎡ ⎢ ⎢ ⎢ ⎣ v i 1 (𝜉) v i 2 (𝜉) v i 3 (𝜉) ⎤ ⎥ ⎥ ⎥ ⎦ , ( 13
)
where

(𝑓 𝑗 1 , 𝑓 𝑗 2 , 𝑓 𝑗 3 ) are defined in Ω z and (v 𝑗 1 , v 𝑗 2 , v 𝑗 3 )
are defined in Ω. The "•" operator is Hadamard's element-wise product. In this paper, a classical eight-node FE approximation is used in Ω, and an LW description is chosen in Ω z as it is particularly suitable for modeling a composite structure.

The strain field for the cylindrical composite structure

The strain field in Equation ( 7) is free of any approximated shell kinematics. These strain components are simplified using Equation ( 12), and we recover the following relations:

𝜀 11 = ( 1 + z R ) -1 ( u 1,1 + 1 R u 3
)

𝜀 22 = u 2,2 𝜀 33 = u 3,3 𝛾 23 = u 2,3 + u 3,2 𝛾 13 = u 1,3 + ( 1 + z R ) -1 ( u 3,1 - 1 R u 1
)

𝛾 12 = u 1,2 + ( 1 + z R ) -1 u 2,1 , (14)
where the covariant derivative becomes a classical derivative for the case of a cylindrical shell.

Equation (13) must be introduced at this level in order to obtain the compatible strain expansion along the normal coordinate z of the shell. Hence, we obtain

𝜺(u) = N ∑ 𝑗=1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝜇 -1 ( 𝑓 𝑗 1 v 𝑗 1,1 + 1 R 𝑓 𝑗 3 v 𝑗 3 ) 𝑓 𝑗 2 v 𝑗 2,2 ( 𝑓 𝑗 3 ) ′ v 𝑗 3 ( 𝑓 𝑗 2 ) ′ v 𝑗 2 + 𝑓 𝑗 3 v 𝑗 3,2 ( 𝑓 𝑗 1 ) ′ v 𝑗 1 + 𝜇 -1 ( 𝑓 𝑗 3 v 𝑗 3,1 -1 R 𝑓 𝑗 1 v 𝑗 1 ) 𝑓 𝑗 1 v 𝑗 1,2 + 𝜇 -1 𝑓 𝑗 2 v 𝑗 2,1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 15
)
where the prime stands for the classical derivative

( 𝑓 ′ i = d𝑓 i dz
) , () ,𝛼 stands for the partial derivative, and 𝜇 = 1 + z R .

The problem to be solved

The resolution of Equation ( 10) is based on a greedy algorithm. If we assume that the first m functions have been already computed, the trial function for the iteration m + 1 is written as

u m+1 = u m + ⎡ ⎢ ⎢ ⎣ 𝑓 1 v 1 𝑓 2 v 2 𝑓 3 v 3 ⎤ ⎥ ⎥ ⎦ = u m + f • v, (16)
where

(v 1 , v 2 , v 3 ) and ( f 1 , f 2 , f 3 )
are the functions to be computed and u m is the associated known set at iteration m defined by

u m = m ∑ i=1 ⎡ ⎢ ⎢ ⎢ ⎣ 𝑓 i 1 v i 1 𝑓 i 2 v i 2 𝑓 i 3 v i 3 ⎤ ⎥ ⎥ ⎥ ⎦ . ( 17
)
The test functions are

𝛿 (f • v) = 𝛿f • v + 𝛿v • f (18) with v = ⎡ ⎢ ⎢ ⎣ v 1 v 2 v 3 ⎤ ⎥ ⎥ ⎦ f = ⎡ ⎢ ⎢ ⎣ 𝑓 1 𝑓 2 𝑓 3 ⎤ ⎥ ⎥ ⎦ . ( 19
)
Introducing the test function defined by Equation ( 18) and the trial function defined by Equation ( 16) into the weak form given by Equation (10), the two following equations to be solved can be deduced:

• for the test function 𝛿f, we have

a(v • f, v • 𝛿f) = b(v • 𝛿f) -a(u m , v • 𝛿f) ∀𝛿f; (20)
• for the test function 𝛿v, we have

a(f • v, f • 𝛿v) = b(f • 𝛿v) -a(u m , f • 𝛿v) ∀𝛿v. (21)
This coupled nonlinear problem is solved using a classical strategy based on a fixed-point method. For each problem, only an unknown 1D or 2D function has to be found. Hence, the approach leads to the process given in Algorithm 1. The fixed-point algorithm is stopped when the distance between two consecutive terms is sufficiently small (cf Vidal et al [START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF] ).

Multiresolution algorithm with a preliminary stage

In this section, a multiresolution algorithm is given, using the advantages of the present variable separation method. The approach described here can be considered as an extension of the so-called preliminary stage introduced in the works of Bussy et al [START_REF] Bussy | The large time increment method for numerical simulation of metal forming processes[END_REF] and Boisse et al [START_REF] Boisse | A new approach in non-linear mechanics: the large time increment method[END_REF] to improve the convergence rate of the LArge Time Increment (LATIN) method. Once some couples (at least two) are built, the process consists in updating only all the 1D z-functions 𝑓 𝑗 i (z), with the basis functions v 𝑗 i (𝜉) being known. Thus, assuming that N 1 couples are built, the problem to be solved can be written as a

( N 1 ∑ i=1 v i • f i , N 1 ∑ i=1 v i • 𝛿f i ) = b ( N 1 ∑ i=1 v i • 𝛿f i ) ∀𝛿f i . ( 22
)
The associated multiresolution algorithm is given in Algorithm 2. The first step consists in building N 1 couples for an initial set of parameters. Then, these 𝜉-functions v p are reused in the following calculations involving new sets of parameters. Thus, only the z-functions are updated. The associated problem is 1D. If the residual error remains high, only few new couples are computed to achieve the convergence in the so-called enrichment stage.

Remark 1. The convergence rate of the classical Algorithm 1 can be improved by introducing the preliminary stage at each iteration of the algorithm to update the 1D z-functions. The new associated algorithm is given in Algorithm 3.

FE discretization

To build the shell FE approximation, a discrete representation of the functions (v, f ) must be introduced. We use a classical FE approximation in Ω and a polynomial expansion in Ω z . The elementary vector of degrees of freedom (dofs) associated with one element Ω e of the mesh in Ω is denoted q v e . The vector of dofs associated with the polynomial expansion in Ω z is denoted q f . The displacement fields and the strain fields are determined from the values of q v e and q 𝑓 e by

v e = N 𝜉 q v e ,  e v = B 𝜉 q v e , f e = N z q 𝑓 e ,  e 𝑓 = B z q 𝑓 e , (23)
where

 e T v = [v 1 v 1,1 v 1,2 v 2 v 2,1 v 2,2 v 3 v 3,1 v 3,2 ] and  e T 𝑓 = [ 𝑓 1 𝑓 ′ 1 𝑓 2 𝑓 ′ 2 𝑓 3 𝑓 ′ 3 ] .
The matrices N 𝜉 , B 𝜉 , N z , and B z contain the interpolation functions, their derivatives, and the Jacobian components.

FE problem to be solved on Ω

For the sake of simplicity, the function f(k) , which is assumed to be known, will be denoted f, and the function ṽ(k) to be computed will be denoted v. The strain in Equation ( 21) is defined in matrix notations as

𝜺( f • v) = Σ z ( f ) v ( 24
)
with

Σ 𝐳 ( f ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 f1 ∕𝜇 0 0 0 0 f3 ∕(𝜇R) 0 0 0 0 0 0 0 f2 0 0 0 0 0 0 0 0 0 f ′ 3 0 0 0 0 0 f ′ 2 0 0 0 0 f3 f ′ 1 -f1 ∕(𝜇R) 0 0 0 0 0 0 f3 ∕𝜇 0 0 0 f1 0 f2 ∕𝜇 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 25
)
The variational problem defined on Ω from Equations ( 6) and ( 21) is

∫ Ω 𝛿 𝐯 T k z ( f ) 𝐯 √ a dΩ ⏟⏟ ⏟ d𝜉 1 d𝜉 2 = ∫ Ω 𝛿v T t z ( f ) √ adΩ -∫ Ω 𝛿 𝐯 T 𝝈 z ( f , u m ) √ adΩ ( 26
)
with

k z ( f ) = ∫ Ω z Σ z ( f ) T CΣ z ( f )𝜇dz (27)
t z ( f ) = f • t 𝜇 | z=z F (28) 𝝈 z ( f , u m ) = ∫ Ω z Σ z ( f ) T C𝜺(u m )𝜇dz. (29)
Note that the calculation of Equations ( 27) and ( 29) is performed using an analytical integration. The introduction of the FE approximation given by Equation ( 23) in the variational Equation ( 26) leads to the linear system

K z ( f ) q v =  v ( f , u m ), (30)
where

• q v is the vector of the nodal displacements in the contravariant basis,

• K z ( f ) is the stiffness matrix obtained by summing up the elements' stiffness matrices

K e z ( f ) = ∫ Ω e B T 𝜉 k z ( f )B 𝜉
√ adΩ e , and •  v ( f , u m ) is the equilibrium residual obtained by summing up the elements' residual load vectors

 e v ( f , u m ) = ∫ Ω e N T 𝜉 t z ( f ) √ adΩ e -∫ Ω e B T 𝜉 𝝈 z ( f , u m ) √ adΩ e .

FE problem to be solved on Ω z

For the sake of simplicity, the function ṽ(k-1) , which is assumed to be known, will be denoted ṽ, and the function f(k) to be computed will be denoted f. The strain in Equation ( 20) is defined in matrix notations as

𝜺(ṽ • 𝑓 ) = Σ 𝛏 (ṽ) f (31)
with

Σ 𝛏 (ṽ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ṽ1,1 ∕𝜇 0 0 0 ṽ3 ∕(𝜇R) 0 0 0 ṽ2,2 0 0 0 0 0 0 0 0 ṽ3 0 0 0 ṽ2 ṽ3,2 0 -ṽ 1 ∕(𝜇R) ṽ1 0 0 ṽ3,1 ∕𝜇 0 ṽ1,2 0 ṽ2,1 ∕𝜇 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 32
)
The variational problem defined on Ω z from Equation ( 20) is

∫ Ω z 𝛿 T f k 𝜉 (ṽ) f 𝜇dz = 𝛿f T t 𝜉 (ṽ) 𝜇 | z=z F -∫ Ω z 𝛿 T f 𝝈 𝜉 (ṽ, u m )𝜇dz (33)
with

k 𝜉 (ṽ) = ∫ Ω Σ 𝛏 (ṽ) T CΣ 𝛏 (ṽ) √ adΩ (34) t 𝜉 (ṽ) = ∫ Ω ṽ • t √ adΩ ( 35
)
𝝈 𝜉 (ṽ, u m ) = ∫ Ω Σ 𝛏 (ṽ) T C𝜺(u m ) √ adΩ. ( 36
)
The introduction of the FE discretization given by Equation ( 23) in the variational Equation (33) leads to the linear system K 𝜉 (ṽ) q 𝑓 =  𝑓 (ṽ, u m ),

where q f is the vector of dofs associated with the polynomial expansion in Ω z , K 𝜉 (ṽ) is a stiffness matrix defined by Equation (38), and  𝑓 (ṽ, u m ) is an equilibrium residual defined by Equation (39), ie,

K 𝜉 (ṽ) = ∫ Ω z B T z k 𝜉 (ṽ)B z 𝜇dz (38)  𝑓 (ṽ, u m ) = N T z t 𝜉 (ṽ)𝜇 | z=z F -∫ Ω z B T z 𝝈 𝜉 (ṽ, u m )𝜇dz. ( 39
)
Remark 2. For the preliminary stage, the same type of computation has to be performed. The system to be solved can be written under the following form:

K PS (v 1 , … , v m-1 ) q 𝑓 1m-1 =  𝑓 (v 1 , … , v m-1 ), (40)
where

• K PS (v 1 , … , v m-1 ) is the assembly of matrices K 𝜉 (v i , v 𝑗 ) = ∫ Ω z B T z k 𝜉 (v i , v 𝑗 )B z 𝜇dz with k 𝜉 (v i , v 𝑗 ) = ∫ Ω Σ 𝛏 (v i ) T CΣ 𝛏 (v 𝑗 ) √ adΩ, •  𝑓 (v 1 , … , v m-1
) is written with the submatrices N T z t 𝜉 (v i )𝜇 | z=z F , and

• q 𝑓 1m-1 T = [q 𝑓 1 T • • • q 𝑓 m-1 T ].
In fact, the computational cost of this stage is mainly due to the calculation of the integrals. The order of magnitude is proportional to N z × N x × N y × N couple , where N x , N y are the number of elements in the 𝜉 1 , 𝜉 2 directions, respectively. N z is the number of numerical layers (number of elements in Ω z ), and N couple is the number of couples built in the PGD process. Nevertheless, as shown in the work of Vidal et al, [START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF] the estimation of the total computational time is proportional to N 3

x × N 𝑦 × N couple (for N x > N y , and a high value of N x × N y , N z ). Thus, the additional cost of the preliminary stage remains low.

NUMERICAL RESULTS

In this section, an eight-node quadrilateral FE based on the Serendipity interpolation functions is used for the unknowns depending on the in-plane coordinates. The geometry of the shell is approximated by this classical FE in the parametric space. The geometrical transformation is based on an explicit map ⃗ Φ. A Gaussian numerical integration with 3 × 3 points is used to calculate the elementary matrices.

In the framework of multiparametric processes, it is required to reduce significantly the cost of the analyses, which can be numerous. Hereafter, as explained in Section 4.5, we take advantage of the possibility offered by the present method to reuse the basis of space 𝜉 functions built during the first resolution. Indeed, once a first analysis has been carried out for a set of geometric/material parameters chosen a priori, the following are based on the built space basis, and only the z-functions are updated. If the residual error remains too high, only few new couples are computed to achieve the convergence. Note that this idea has already been used in other works [START_REF] Boucard | A multiple solution method for non-linear structural mechanics[END_REF][START_REF] Allix | A new multi-solution approach suitable for structural identification problems[END_REF][START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF][START_REF] Néron | Time-space PGD for the rapid solution of 3D nonlinear parametrized problems in the many-query context[END_REF] in the framework of nonlinear problems with a time-space decomposition to perform multiparametric analysis. Here, it is extended to the space decomposition in the framework of the composite shell structure.

Several static tests are presented, validating our approach (denoted VS-LD4) and evaluating its efficiency. First, the dedicated approach is assessed by varying the stacking sequences of the laminates. Usual stacking sequences used in aeronautics are addressed. Then, the multiresolution process is carried out for different geometries of the shell. Finally, the possibilities of the approach are illustrated by considering a variation of both the geometry and the number of layers. We pay special attention to the assessment of the out-of-plane stresses that are of relevant importance in the damage phenomenon of such structures. Note that the accuracy of the results is assessed by a comparison with a fourth-order layer-wise model (denoted LD4), referring to the systematic work of Carrera and Carrera's unified formulation. [START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices[END_REF][START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 2: numerical implementations[END_REF] All of the following numerical assessments are based on the configuration proposed in the work of Charrier et al. [START_REF] Charrier | Determination of the out-of-plane tensile strength using four-point bending tests on laminated L-angle specimens with different stacking sequences and total thicknesses[END_REF] It concerns L-angle specimens that are usually used in the junction between different perpendicular panels of the aeronautical structures. It is described as follows.

Geometry: a very thick composite cylindrical shell with 𝜙 = 𝜋∕2. The number of layers (NC) with the stacking sequences and the radius R are precise in each section. All layers have the same thickness. Boundary conditions: a shell clamped along one straight edge and free on the other sides, subjected to a constant global pressure: q(𝜉) = q 0 = 1 MPa. Material properties: the T700GC/M21 laminates are considered with the following characteristics: E L = 115 GPa , E T = 8.5 GPa , G LT = 4.5 GPa , G TT = 3.0357 GPa , 𝜈 LT = 0.32, 𝜈 TT = 0.4, where L refers to the fiber direction, T refers to the transverse direction. Mesh: The mesh consists of (N x = 16) × (N y = 16) elements in the 𝜉 1 and 𝜉 2 directions, respectively. It is considered as a suitable mesh for the present test cases. Number of dofs: N do f xy = 3(N x + 1)(N y + 1) and N do f z = 12 × 𝛼NC + 3 are the number of dofs of the two problems associated with v i 𝑗 and 𝑓 i 𝑗 , respectively. 𝛼 is the number of numerical layers per physical layer. Hence, the total number of dofs is N do f xy + N do f z . Reference values: LD4 model. Different stacking sequences (reported in the work of Charrier et al [START_REF] Charrier | Determination of the out-of-plane tensile strength using four-point bending tests on laminated L-angle specimens with different stacking sequences and total thicknesses[END_REF] ), which can be considered as a large class of lay-up types, are addressed. They are described as follows:

• LU1: any laminate [0 • ∕45 • ∕90 • ∕ -45 • ∕0 • ∕45 • ∕ -45 • ∕90 • ], it has no symmetry; • LU2: a highly oriented laminate [0 • 3 ∕45 • ∕90 • 2 ∕ -45 • ∕0 •
] ns , usually used to be subjected to high uniaxial tensile loading;

• LU3: a quasi-isotropic laminate [(45 • ∕90 • ∕ -45 • ∕0 • ) 2 ] ns ,
usually used in order to be subjected to complex multiaxial loadings;

• LU4: a highly disoriented laminate [45 • 2 ∕0 • ∕ -45 • 3 ∕90 • ∕45 • ]
ns designed mainly to be subjected to in-plane shear loading, including 75% of ±45 • plies.

The different geometries of the shell involved in the subsequent examples are summarized in Table 1. In particular, n is given. Each configuration will be denoted LUI-SPJ, where I ∈ {1, 2, 3, 4} and J ∈ {1, 2, 3}.

Initial solution

First, the LU2-SP2 specimen with 32 layers is considered as the initial set of parameters. Eight couples are built to obtain the solution. The distributions of the in-plane functions v 𝑗 i (𝜉) and the thickness functions 𝑓 𝑗 i (z) are given in Figures 3 and4. In Figure 3B, v 𝑗 1 (𝜉) is associated to the displacement along the radius in the curvilinear coordinate system. (A) (B) In Figure 4B, v 𝑗 3 (𝜉) is associated to the transverse displacement. This function is plotted with respect to the undeformed shape of the structure (red cross) in the plane 𝜉 2 = L y ∕2, and the amplification ratio 𝛼 i is given for all couples, allowing us to compare them against each other. 𝑓 𝑗 i (z) are normalized such that max |𝑓 𝑗 i (z)| = 1. The two first functions represent the main contributions to the displacement solution. They can be considered as global modes. On the contrary, the following functions involve local corrections that are localized close to the clamped or free edges. For the latter, the distribution of the associated z-functions through the thickness is rather complex. Thus, highly refined models are required to capture the variations of the mechanical quantities.

FIGURE 3 Couples for LU2-SP2. A, 𝑓 𝑗 1 (z); B, v 𝑗 1 (𝜉 1 , 𝜉 2 = L 𝑦 ∕2) [

Multiresolution: different lay-up configurations

Once the eight couples are built for the LU2-SP2 specimen, they are now used for two subsequent configurations, namely, the LU3 and LU4-SP2 shells. The new z-functions associated to the LU3-SP2 shell are given in Figure 5. The change of stacks requires the correction of the thickness functions. It should be noted that the two first 𝑓 𝑗 3 (z) functions remain quasi-constant and linear. To assess the accuracy of the results, the error rate of the local stresses 𝜎 11 , 𝜎 13 , and 𝜎 33 at 𝜉 1 = L 𝜉 ∕4, 𝜉 2 = L y ∕2 is first evaluated, as it is located in an area with the highest stresses. The results are summarized in Table 2. The calculation of these stresses with the initial functions drives an important error rate. This can be reduced significantly by updating the eight z-functions 𝑓 𝑗 i (z). The maximum error rate becomes 3.7%. The accuracy can be improved by enriching the basis with only three functions (error of less than 1.8%). For further assessment, the distributions of 𝜎 11 and 𝜎 13 through the thickness are compared for the initial basis and the updated version with three new couples in Figure 6. Significant corrections are provided to the initial solution. The accuracy of the results becomes very good with respect to the LD4 model. The complex behavior of the out-of-plane stress is also well described. The continuity of this stress is fulfilled. For this test case, the of magnitude of the transverse shear and normal stresses is the same in each configuration. Note that the maximum values are located at about z = -e∕4. To illustrate the capability of the present method to provide quasi-3D results, the distribution of stresses over the plane 𝜉 2 = L y ∕2 is also shown in Figures 7,8, and 9. The results appear in very good agreement with the LD4 model. Even the local effects, such as the clamped-edge ones, are rather similar. Ignoring this effect, the area where solicitation is most significant is also very well described (around 𝜉 1 = L 𝜉 ∕4 and z = -e∕4). Finally, despite the strong gap between the initial solution (only the initial functions) and the final solutions (with three new couples) illustrated in Figure 6, a very good accuracy level is achieved with only few new couples. The method can be carried out for a large variety of stacking sequences.

Multiresolution: different geometries

The same procedure is tested, considering a variation of the geometry. In this study, the involved parameters are the thickness of each layer e ply , the radius of the shell R, and the opening angle 𝜙. The initial basis consists of 13 couples for the LU4-SP1 configuration (16 layers). A simple mapping of the 1D domain (for e ply and R) or the 2D domain (for 𝜙) is achieved for the following computations. The results are given in Tables 3,4, and 5. As previously, a local error on the stresses is provided. Wide domains of variation of the geometry have been considered. For a variation of 10% of the parameters, no additional couple is required to obtain an accurate solution; only the z-functions have to be updated. Beyond this value, only one or two couples are built despite the high geometric modification. Note that the basis 𝜉-functions seem to be sufficiently rich to represent the changes of the opening angle as the update of the z-functions allows us to decrease drastically the error rate (see Table 5).

Multiresolution: different geometries with different numbers of layers

In this section, the multiresolution process is performed with different numbers of layers and different total thicknesses. We take advantage of the separation of the solution between the in-plane and out-of-plane functions. Thus, using the 13 in-plane functions computed for the LU4-SP1 shell (16 plies), the solution can be built for the LU2-SP2 (32 plies; e = 8.384 mm) and LU1-SP3 (eight plies; e = 2.096 mm) configurations. Table 6 shows that two new couples allow us to have an error rate of less than 1.1%. The comparison of the initial solution (LU4-SP1) and the final solution (LU1-SP3 [eight plies]) is given in Figure 13. The solutions are deeply different. This figure shows also that the distribution of the in-plane and out-of-plane stresses through the thickness is very accurate when compared with the LD4 model. This test case allows us to show the possibilities of the method to address many shell configurations.

CONCLUSION

In this paper, PGD is advantageously performed in the framework of the multiresolution process for modeling composite shell structures. Indeed, the in-plane functions built for one fixed set of parameters can be reused for the subsequent configurations, reducing the computational cost to a 1D problem. If needed, only few new couples are computed. The performance of the method is illustrated for a wide type of shell structures. The changes of stacking sequences, geometry (radius, thickness, and open angle), and number of layers have been addressed. The good accuracy of the results is shown by a comparison with a quasi-3D layer-wise reference solution. Thus, the potential of the approach seems to be very interesting in the framework of multiparametric modeling of composite structures.

ORCID

P. Vidal https://orcid.org/0000-0002-9634-4443

FIGURE 4 11 FIGURE 5

 4115 FIGURE 4 Couples for LU2-SP2. A, 𝑓 𝑗 3 (z); B, v 𝑗 3 (𝜉 1 , 𝜉 2 = L 𝑦 ∕2) [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 6

 6 FIGURE 6 Distribution of 𝜎 11 (𝜉 1 = L 𝜉 ∕2, 𝜉 2 = Ly∕2, z) and 𝜎 13 (𝜉 1 = L 𝜉 ∕2, 𝜉 2 = Ly∕2, z) along the thickness -multiresolution process from LU2-SP2 (𝜉) basis functions -32 layers. A, LU3; B, LU4 [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 7

 7 FIGURE 7 𝜎 11 (𝜉 1 , 𝜉 2 = L y ∕2, z) -LU3-SP2 -32 layers. A, VS-LD4; B, LD4 [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 8 FIGURE 9

 89 FIGURE 8 𝜎 13 (𝜉 1 , 𝜉 2 = L y ∕2, z) -LU3-SP2 -32 layers. A, VS-LD4; B, LD4 [Colour figure can be viewed at wileyonlinelibrary.com]

10FIGURE 11

 11 FIGURE 11 𝜎 11 (𝜉 1 , 𝜉 2 = L y ∕2, z) -LU4-SP1 -16 layers. A, VS-LD4; B, LD4 [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 12 𝜎

 12 FIGURE 12 𝜎 13 (𝜉 1 , 𝜉 2 = L y ∕2, z) -LU4-SP1 -16 layers. A, VS-LD4; B, LD4 [Colour figure can be viewed at wileyonlinelibrary.com]For further comparison, the distribution of the out-of-plane stresses through the thickness is presented in Figure10for e ply = 0.1 mm and R = 18 mm. It can be shown that the results are in excellent agreement with the reference LD4 solution. The top/bottom conditions are fulfilled. Then, the distribution of 𝜎 11 and 𝜎 13 over the plane 𝜉 2 = L y ∕2 for 𝜙 = 𝜋 4 is given in Figures11 and 12. It can be also inferred from these figures that the accuracy of the results is very satisfactory. The localization of the in-plane stresses is well captured near the clamped edge.

2 FIGURE 13

 213 FIGURE 13 Distribution of 𝜎 11 (𝜉 1 = L 𝜉 ∕4, 𝜉 2 = Ly∕2, z), 𝜎 13 (𝜉 1 = L 𝜉 ∕4, 𝜉 2 = Ly∕2, z), and 𝜎 33 (𝜉 1 = L 𝜉 ∕2, 𝜉 2 = Ly∕2, z) along the thicknessmultiresolution process from LU4-SP1 𝜉 basis functions -geometry -16 layers [Colour figure can be viewed at wileyonlinelibrary.com]

  indices is used, and the partial derivative is denoted by () ,𝛼 . A shell is characterized by the first fundamental form a 𝛼𝛽 and the second fundamental form b 𝛼𝛽 . Their covariant, contravariant, and mixed-form definitions are given by

	convention on repeated
				The summation
		g 3
	a 3	a 2	g 2	g 1
		a 1
	e
	FIGURE 1 The map ⃗ Φ and the local basis vectors ⃗ a i and ⃗ g i for a shell panel

TABLE 1

 1 Geometry of the shell

	n No. of Layers Layer Thickness, mm	𝜙	R, mm L y , mm
	SP1 1	16	0.262	𝜋∕2	6.095	20
	SP2 2	32	0.262	𝜋∕2	12.19	40
	SP3 -	8	0.262	𝜋∕2	6.095	20

TABLE 2

 2 Multiresolution process from LU2-SP2 𝜉 basis functions -32 layers

	Lay-Up		LU3									LU4
	f	Initial Basis i Not Updated f j j i Updated Only New	( f	j i , v	i	j	)	f	Initial Basis i Not Updated f j j i Updated Only New	(	f	j i , v	i	j	)
	𝜎 max 11 (𝜉 1 = L 𝜉 ∕4, 𝜉 2 = L 𝑦 ∕2, z)	-16.21 (49%)	-31.81 (0.4%)	-32.04 (0.3%)		-18.91 (43%)	-33.12 (0.3%)	-32.86 (1%)
	𝜎 max 13 (𝜉 1 = L 𝜉 ∕4, 𝜉 2 = L 𝑦 ∕2, z)	-3.22 (10%)	-2.89 (1.1%)	-2.87 (1.8%)		-3.18 (11%)	-2.78 (3%)	-2.85 (0.5%)
	𝜎 max 33 (𝜉 1 = L 𝜉 ∕4, 𝜉 2 = L 𝑦 ∕2, z)	-2.32 (9%)	-2.50 (2.3%)	-2.57 (0.6%)		-2.32 (5.8%)	-2.37 (3.7%)	-2.46 (0.2%)
	Number of new couples	0	0		3						0	0		3

TABLE 3

 3 Error rate on stresses at (𝜉 1

	-multiresolution process from LU4-SP1 𝜉 basis functions -parameter e ply -16 layers	e ply = 0.2 mm (-24%) e ply = 0.1
	z)
	y ∕2,
	= L 𝜉 ∕4, 𝜉 2 = L	0.2358 mm (-10%)
		e ply =

mm (-62%) Initial Basis Initial Basis Initial Basis f j i Not Updated f j i Updated New ( f j i , v j i ) f j i Not Updated f j i Updated New

	)
	j
	i
	j i Updated New ( f j i , v	20% 0.8%	36% 1.2%	16% 0.9%	0 2
	j i Not Updated f	83%	35%	42%	0
	f
	)
	j
	i
	( j f i , v	0.6%	0.4%	0.5%	1
		4.0%	5.6%	2.5%	0
		38%	31%	10%	0
		-	-	-	-
		1.4%	1%	0.6%	0
		17%	15%	4%	0
		𝜎 max 11	𝜎 max 13	𝜎 max 33	New couples

TABLE 4

 4 Error rate on stresses at (𝜉 1

			)
			j
			i
	= L 𝜉 ∕4, 𝜉 2 = L y ∕2, z) -multiresolution process from LU4-SP1 𝜉 basis functions -parameter R -16 layers	R = 6.7 mm (+10%) R = 8 mm (+31%) R = 18 mm (+195%)	j i Not Updated f j i Updated New ( Initial Basis f j i , v j i ) f j i Not Updated f Initial Basis j i Updated New ( f j i , v j i ) f j i Not Updated f Initial Basis j ( f j i , v i Updated New	24% 1.8% -55% 6.6% 1.4% 96% 0.8% 47%	22% 0.8% -48% 6.1% 0.6% 55% 2.3% 26%	6% 0.8% -18% 4.1% 0.06% 55% 1.5% 37%	0 0 -0 0 1 0 0
			f
				𝜎 max 11	𝜎 max 13	𝜎 max 33	New couples

TABLE 5

 5 Error rate on stresses at (𝜉 1 = L 𝜉 ∕4, 𝜉 2 = L y ∕2, z) -multiresolution process from LU4-SP1 𝜉 basis functions -parameter 𝜙 -16 layers

		𝝓 = 𝟑 𝟒	𝝅 𝟐 (-25%)							𝝓 = 𝟏 𝟐	𝝅 𝟐 (-50%)
	f	Initial Basis i Not Updated f j j i Updated Only New	(	f	j i , v	j i	)	f	Initial Basis i Not Updated f j j i Updated Only New	(	f	j i , v	i j	)
	𝜎 max 11	122%	0.6%	-					880%	0.8%	0.7%
	𝜎 max 13	97%	0.9%	-					447%	2.7%	1.7%
	𝜎 max 33	34%	0.3%	-					30%	0.1%	0.9%
	New couples	0	0	-					0	0		1

TABLE 6

 6 Error rate on stresses at (𝜉 1 = L 𝜉 ∕4, 𝜉 2 = L y ∕2, z) -multiresolution process from LU4-SP1 𝜉 basis functions -different lay-ups

	f	32 Layers LU2-SP2 i Updated New j ( f j i , v	j i	)	f	8 Layers LU1-SP3 i Updated New j ( f j i , v	i j	)
	𝜎 max 11	1.2%	0.7%				15%	0.9%
	𝜎 max 13	5.2%	0.04%				19%	1.1%
	𝜎 max 33	3.7%	0.6%				9%	0.1%
	New couples	0	2				0	2