Ahmed Boukar

Stephane Corn
email: stephane.corn@mines-ales.fr

Pierre R L Slangen

Patrick Ienny

Finite element modelling of low velocity impact test applied to biaxial glass fiber reinforced laminate composites

Keywords:

The aim of this paper is to describe the dynamic behavior of biaxial glass fiber reinforced laminate composites under low velocity impact test through finite element modelling. Experimental investigations by impact test performed using an instrumented drop weight testing machine were conducted on three-point bending composite samples in order to assess their impact damage resistance. Moreover, the experimental setup allowed the visualization of real-time damage progression of the impacted laminate composite via high-speed camera Phantom V2512 enabling to capture 83000 frames per second. Dynamic strain fields were extracted by Digital Image Correlation (DIC) method. Based on the experimental results, a numerical study of the impacted specimens was developed as a user subroutine VUMAT integrated to ABAQUS/Explicit in order to precisely capture the progressive dynamic failure of the laminate composite under impact test. In the proposed model, the damage and failure of each ply are accounted by a Hashin 3D damage-based behavior, and a cohesive zone model is employed to capture the onset and progression of inter-laminar delamination. A good experimental-numerical correlation was obtained for peak force and failure modes.

Introduction

Composite materials are being used increasingly in many industrial sectors, such as aerospace, aeronautics, defence, sports and buildings. The benefits of composites are related to their excellent mechanical strength and stiffness, considering their lightness. However, impact damages in service may severely reduce the strength and stability of the laminated composite structure. Even when impact does not lead to visible damage, the impact energy dissipation may induce internal cracking networks, limiting the integrity of the structure [START_REF] Abrate | Impact on Composite Structures[END_REF][START_REF] Sr R G Z | Impact Behaviour of Fibre-Reinforced Composite Materials and Structures[END_REF]. Composite laminates present a brittle behavior and may be subject to internal damage in the form of matrix cracking, fiber fracture and interlaminar delamination when they are subjected to object impacts [START_REF] Tomita | Effect of lay-up sequence on mechanical properties and fracture behaviour of advanced CFRP laminate composite[END_REF]. These internal damages may considerably reduce the capacity of the composite laminate to undergo additional loads [START_REF] Hitchen | The effect of stacking sequence on impact damage in a carbon fibre /epoxy composite[END_REF].

In the literature, we can find several studies on the face-on impact of laminated composites [START_REF] Zhang | Finite element modeling of damage development in cross-ply composite laminates subjected to low velocity impact[END_REF][START_REF] Tuo | Damage and failure mechanism of thin composite laminates under low-velocity impact and compression-after-impact loading conditions[END_REF][START_REF] Lopes | Low-velocity impact damage on dispersed stacking sequence laminates . Part II : Numerical simulations[END_REF][START_REF] Zhou | Numerical investigation on the repeated low-velocity impact behavior of composite laminates[END_REF][START_REF] Li | Assessment of failure criteria and damage evolution methods for composite laminates under low-velocity impact[END_REF][START_REF] Sun | Barely visible impact damage in scaled composite laminates: Experiments and numerical simulations[END_REF][START_REF] Zhang | Simulating low-velocity impact induced delamination in composites by a quasi-static load model with surface-based cohesive contact[END_REF]. Many researches have focused on the response of damage mechanisms and the extent of damage due to impact and its effect on the residual strength of the structure. Studies on edge-on impact have been scarce in comparison to face-on impact. The edge-on impact is most frequently associated to impact on the edge of a stringer on a stiffened panel. The bending test under low velocity impact test is well suited to study the edge-on impact of laminated composite. The edge-on impact measurements presented in this paper were carried out on biaxial glass fiber-reinforced laminated composite to understand the edge-on impact failure modes correctly and develop models to assess the structural integrity during the impact event.

The characterization and the assessment of impact damage is usually done after the impact using visual inspection, ultrasonic C-scan [START_REF] Caminero | Damage resistance of carbon fibre reinforced epoxy laminates subjected to low velocity impact: Effects of laminate thickness and ply-stacking sequence[END_REF][START_REF] Caminero | Internal damage evaluation of composite structures using phased array ultrasonic technique: Impact damage assessment in CFRP and 3D printed reinforced composites[END_REF] or other NDI methods. To improve standard test, high-speed cameras have been used as practical tools to capture the impact events for various types of composites materials and projectiles [START_REF] Kirugulige | Measurement of transient deformations using digital image correlation method and high-speed photography: application to dynamic fracture[END_REF]. The method of digital image correlation (DIC) is employed in conjunction with high-speed imaging to study full-field surface displacements. For example, Puech et al. [START_REF] Puech | Investigating the impact behaviour of short hemp fi bres reinforced polypropylene biocomposites through high speed imaging and fi nite element modelling[END_REF] and Ramakrishnan et al. [START_REF] Ramakrishnan | Experimental assessment of low velocity impact damage in flax fabrics reinforced biocomposites by coupled high-speed imaging and DIC analysis[END_REF] applied coupled high-speed imaging and DIC analysis to study the initiation and propagation of cracks in short natural fibers reinforced biocomposite plates during the low-velocity impact test. Pan et al. [START_REF] Pan | Full-field transient 3D deformation measurement of 3D braided composite panels during ballistic impact using singlecamera high-speed stereo-digital image correlation[END_REF] conducted a full-field transient of 3D deformation measurement of 3D braided composite panels during ballistic impact using a single high-speed camera and stereo-digital image correlation. Some researchers used two high-speed cameras in order to recover all the 3D data [START_REF] Cuynet | Damage characterisation of fl ax fi bre fabric reinforced epoxy composites during low velocity impacts using high-speed imaging and Stereo Image Correlation[END_REF]: Flores et al. [START_REF] Flores | Highspeed 3D digital image correlation of low-velocity impacts on composite plates[END_REF] conducted a 3D digital image correlation using two high-speed digital cameras at 20,000 fps to capture the out-of-plane displacements and in-plane strains.

The finite element simulation to assess the response of composite materials to impact loads has been widely considered in the literature in order to understand the complexity of internal damage and predicted in a relatively short time when composite is subject to impact loading. The numerical modelling of laminated composites can address micro, meso and macro scales. The continuum damage model (CDM) has become the most popular simulation approach, which considers damage initiation and stiffness degradation [START_REF] Zhou | Numerical investigation on the repeated low-velocity impact behavior of composite laminates[END_REF][START_REF] Cui | A critical study on impact damage simulation of IM7/8552 composite laminate plate[END_REF].

Developed initially from classical lamination theory, failure criteria for composites materials have been subject to continuous development and have been organized into two general groups: non-interactive and interactive criteria. Maximum stress and strain criteria, which do not contain shear stress, are usually used to evaluate fiber damage produced by low-velocity impact [START_REF] Liu | Predicting impact induced delamination of FRP laminates[END_REF][START_REF] Liao | Multi-scale modelling of dynamic progressive failure in composite laminates subjected to low velocity impact[END_REF]. As for interactive criteria, the polynomial Tsai-Wu criteria [START_REF] Tsai | A General Theory of Strength for Anisotropic Materials[END_REF] are presently being employed to capture the in-plane matrix cracking of composites structures under impact loads. With regard to damage initiation, the interactive criteria with separate expressions to assess fiber and matrix damage under tensile and compression loading have been widely used, in particular the Hashin [24], Chang-Chang [START_REF] Chang | Laminated Composites Containing Stress Concentrations[END_REF] and Hou [START_REF] Hou | Prediction of impact damage in composite plates[END_REF] criteria. As the unique among these criteria, the 3D Hashin criteria contains the normal stress component σ 33 to correctly capture the through-thickness stress and its distribution during ply failure.

Additionally, an appropriate damage evolution model needs to be defined in order to describe the loss of stiffness after damage initiation (where failure criteria are satisfied). To simulate the process of damage accumulation, several researchers used a predefined constant for the stiffness degradation rule. Obviously, this approach is not applicable for a wide range of situations due to the involvement of mechanical parameters. An exponential damage evolution law is widely used to control the loss of stiffness and damage surface growth for each failure mode [START_REF] Khan | Effect of metal layer placement on the damage and energy absorption mechanisms in aluminium/glass fibre laminates[END_REF][START_REF] Yang | In-plane compression response of woven CFRP composite after low-velocity impact: Modelling and experiment[END_REF][START_REF] Han | Damage evolution analysis of open-hole tensile laminated composites using a progress damage model verified by AE and DIC[END_REF]. The equivalent displacement and stress method in the damage evolution model is now widely employed for predicting the progressive damage behavior of composite materials [START_REF] Li | Assessment of failure criteria and damage evolution methods for composite laminates under low-velocity impact[END_REF][START_REF] Li | Impact behaviour of pultruded GFRP composites under low-velocity impact loading[END_REF]. Among these methods, the fracture energy-based methods that involve equivalent displacement damage evolution will be chosen to model the loss of stiffness during the impact test. Among all damage modes, delamination is considered to be the most critical damage mode since it may propagate invisibly between adjacent layers and lead to unexpected failure of the structure. To model this failure mode, two main methods have been adopted by researchers in literature to predict delamination in laminated composites, such as the virtual crack closure technique VCCT [START_REF] Krueger | Virtual crack closure technique: History, approach, and applications[END_REF] in which the propagation occurs along the predefined surface and cohesive zone model CZM [START_REF] Camanho | Mixed-Mode Decohesion Finite Elements in for the Simulation Composite of Delamination Materials[END_REF] using a bilinear traction-separation law and fracture energy criterion to predict the damage initiation and propagation, which can surmount the limitations of the VCCT method [START_REF] Lopes | Low-velocity impact damage on dispersed stacking sequence laminates . Part II : Numerical simulations[END_REF][START_REF] Sun | Barely visible impact damage in scaled composite laminates: Experiments and numerical simulations[END_REF][START_REF] Yang | In-plane compression response of woven CFRP composite after low-velocity impact: Modelling and experiment[END_REF].

This work is devoted to studying the dynamic behavior of laminated composite subjects to low velocity impact test. A Non-Crimp Fabric (NCF) biaxial glass fiber/polyester laminate was chosen as a material for this study. Impact tests were performed on samples using a drop-weight testing machine in order to characterize the key parameters influencing their impact damage resistance. Moreover, the implementation of a high-speed camera allowed the visualization of real-time damage progression and the measurement of full-field displacement of the specimen during the impact test. Based on these results and on the mechanical tensile behavior of the laminate composites, a three-dimensional finite element model has been implemented in a user-defined material subroutine VUMAT built on ABAQUS/Explicit to predict the impact behavior of these composites.

Materials and methods

Materials and process

The GFRP laminate materials tested in this work were manufactured from a biaxial Non-Crimp Fabric (NCF) glass fiber 0 • /90 • FGE100 with a density of 2.55g/cm 3 manufactured by FORMAX. This fabric is made of two layers having an areal weight of 283g/m 2 , stitched together by polyester threads of 6g/m 2 , leading to a total areal weight of 572g/m 2 for the fabric. The layers are glass fiber bundles with a linear density of 1200tex in direction 0 • and 300tex in direction 90 • . The polyester (Synolite 8488-G-2) provided by DMS Composite Resins was used as a matrix. From relevant datasheet values, this polyester resin has a density of 1.05g/cm 3 , viscosity of 80-90mPa.s, tensile modulus of 3.5GPa and tensile strength of 66MPa.

The laminate was manufactured by vacuum infusion of the resin into the stacked NCF-reinforcement. Across this paper, the fabric composites are denoted [θ] ns , where θ is the angle (in degree) between the warp direction and the plate side edge (direction 1) Fig. 1, s denote that laminates are symmetric about their geometric midplane and n the number of (NCF) plies that is shared by both halves of the laminate. A total of 6 plies were stacked together, resulting in 3mm thickness, with two main lay-up configurations (cross-ply [0/90] 3s, angle-ply [±45] 3s) .

The glass fiber weight fraction was determined by measuring the fibers weight before the infusion process and then the weight of the resulting composite. Values of 558g and 813g, respectively led to a glass fiber weight fraction of 68.6wt%.

Panels were cut into rectangular coupons of 250 × 25 mm 2 for the tensile tests and 80 × 10 mm 2 for the impact tests relative to the specimen thickness.

Microstructure characterization

Observations of polished cross-sections of NCF laminate samples were performed with a FEI Quanta 200 Scanning Electron Microscope (SEM) in order to assess average values of the diameter d f of the glass fibers and their volume fraction V f inside the yarns as well as the shape and dimensions of the yarns cross-section. The cross-section of the NCF laminate Fig. 2(a) exhibits the yarns in the two orthogonal directions and the resin-rich zones in the intra-yarn between two consecutive plies. The section shape of yarn were considered as elliptic with major axes and minor axes of average length 2.8 mm and 0.285 mm respectively, and 0.35 mm average length of resin-rich region present between yarns. SEM images processed using the ImageJ software were applied a thresholding operation consisting of separating fibers and matrix Fig. 2 (b)-(c), leading to V f = 70% and d f = 13.46 μm.

Mechanical testing 2.3.1. Static tensile test

Tensile mechanical properties of the biaxial NCF laminated composite were evaluated with an MTS testing machine (model criterion C45.105) equipped with 100 KN load cell and DIC equipment to measure the strain. Samples were loaded at a crosshead speed of 1mm/min. For each orientation (0 • , 90 • and 45 • ), three specimens were tested. Sample dimensions were 250mm length, 25mm width and around 3mm thickness.

Low-velocity impact test

A schematic of the experimental setup is shown in Fig. 3. It consisted of an Instron CEAST 9340 drop-weight impact machine for delivering low-velocity impact and a Phantom V2512 high-speed camera for capturing speckle pattern images in real-time in order to track the progression of cracks at the surface of samples during impact. The impact kinetic energy was varied by modifying the impactor height and weight.

The setup also included a CEAST Data Acquisition and Analysis System (DAS) that generates a signal pulse to trigger the image capture when the impactor contacts the specimen. Two separate computers were used to record the impact force and the other to control the high-speed camera and store digitized images. The high-speed images were recorded in 10bits at a constant frame rate of 83000 fps and with reduced sensor resolution (640 × 400 pixels). The exposure time was set at 1μs, which is short enough to ensure blur-free images [START_REF] Practices | A Good Practices Guide for Digital Image Correlation[END_REF]. The impact test specimen was painted with a random black pattern on a painted white background to aid in the image correlation.

The drop-weight impact device applies a dynamic three-point bending to an un-notched specimen, all the specimens in this study beeing tested in edgewise (the applied load direction is parallel to the laminate) see Fig. 1(b) according to ISO 179-1. The tup of the impactor is instrumented with a sensor for measuring the dynamic load and thus the acceleration during the impact. Following Newton's law and assuming permanent contact between the impactor and the sample, the force contact history F(t), velocity v(t) and the displacement s(t) of the

Table 1

Predicted elastic properties of the intra-yarn material.

Elastic properties FE homogenization

Young

Table 2

Effective experimental and predicted mechanical properties of laminated NCF composite.

Homogenization Experimental

E 1 = E 2 (GPa) 22.73 23.1 E 3 (GPa) 11.3 - G 12 (GPa) 4.98 5.76 G 13 = G 23 (GPa) 3.65 - v 12
0.141 0.14 v 13 = v 23 0.38 specimen during the impact test is calculated according to the following relationships.

F(t) = m * a(t), V(t) = V 0 - ( 1 m ) ∫ t 0 F(t)dt and s(t) = ∫ t 0 v(t)dt (1)
With m as the mass of impactor, V(t) as the instantaneous velocity of the impactor during the test, and V 0 as the initial impact velocity.

DIC on high-speed camera images

A 2D local DIC technique was used to obtain the full-field displacements and strains at the centre region of the specimen. A rectangular Region Of Interest (ROI) of 320 pixels (10 mm) in the horizontal direction and 256 pixels (8 mm) in the vertical direction was selected, with its centre at the middle of the specimen lateral face, comprising the expected region of cracks occurrence and propagation. The images postprocessing using a in-house DIC software called CinEMA led to the subpixel accurate components of the two-dimensional displacements field.

Results and discussions

Elastic properties assessment

The described experimental procedure provided part of the elastic mechanical properties of the biaxial NCF composite, while a multiscale homogenization approach seemed appropriate for assessing its out-ofplane elastic modulus. Fig. 4 depicts the successive steps of this multiscale homogenization process detailed hereinafter.

Homogenization at the microscopic intra-yarn scale

Each yarn, embedded in NCF fabric composites, consists of many fibers densely packed and bonded together by the matrix. Yarns are commonly modelled as unidirectional composites, whose mechanical behavior is considered to be transversely isotropic. In order to determine the homogenized elastic properties of this intra-yarn material, a periodic FE numerical homogenization based on strain energy equivalence [START_REF] Xia | A unified periodical boundary conditions for representative volume elements of composites and applications[END_REF] is applied to a parallelepiped Representative Elementary Volume (RVE) with a hexagonal pattern (generated using the Python script language of Abaqus Finite Element software), as shown Fig. 5. Based on the outputs of the aforementioned microstructure characterization, the dimensions a 2 and a 3 of the RVE were set such as to comply to 70% volume fraction of fibers of constant diameter d f = 13.46 μm (the not-influent dimension a 1 being set at a 1 = a 2 /4). Glass fibers and the polyester resin were considered isotropic in the FE computation (their elastic properties set according to datasheets). The resulting anisotropic elastic properties of the intra-yarn material predicted by this FE homogenization method are listed in Table 1.

Homogenization at the macroscopic scale

The FE model considered at the laminate macroscopic scale consists of two perpendicular yarns generated using Texgen software in a parallelepiped RVE of dimensions (L1; L2; L3) = (3.6; 3.6; 0.6) Fig. 6. The cross-section of the yarns is elliptical with major and minor axes of lengths 2.8 and 0.285 mm respectively, according to the aforementioned microstructure analysis. The RVE is then meshed, and periodic boundary conditions are applied for the computation.

The computed predicted elastic properties are gathered in Table 2 and compared to their corresponding measured value (if available) as a validation.

As observed from Table 2, the in-plane elastic engineering constant computed using the proposed multiscale numerical homogenization agree well with those measured using the experimental tensile test. Moreover, the multiscale homogenization approach provides the out-of-plane mechanical modulus essential for the 3D finite element damage simulation presented in section 4.

Low-velocity impact test results

Fig. 7. reports some representative load-displacement and energytime responses for [0/90] 3s and [±45] 3s laminates in order to assess and characterize their impact behavior. The load curve response recorded by the instrumented impact test shows various vibrational waves superimposed on the main load-time curve, especially in the rising part of the curve. These oscillations can be interpreted in terms of the main impact impulse, high-frequency oscillation and free vibration [START_REF] Found | Interpretation of signals from dropweight impact tests[END_REF].

The result showed that lay-up architecture affected the impact resistance of (NCF) GFRP laminates. [±45] 3S laminates exhibited lower peak force, longer duration of impact force and higher absorbed energy than [0/90] 3S .

The average absorbed energy measured with the impact test for both materials is shown in Fig. 8, and the failure impact modes are depicted in Fig. 9 and Fig. 10.

Digital image correlation results for the low velocity impact test

Fig. 11 and Fig. 12 show the displacement field of [0/90] 3S and [±45] 3s specimen at labelled events, ①, ②, ③ and ④ described in Fig. 9 and Fig. 10 above. It is observed that the maximum displacement occurs at the top of the sample section travelling vertically down due to impactor displacement, and its distribution is symmetric about the 2axis as expected by the loading conditions.

To underline the presence of the damage process, Fig. 13 presents a comparison of the damage propagation and crack growth at a different time of the impact for both different materials architectures [0/90] 3s and [±45] 3s laminated composite using the in-plane equivalent strain maps obtained from DIC analysis at the surface of the samples. The selected subset size was large enough to ensure that there is a sufficiently distinctive pattern in the ROI. From the current analysis, an investigation of the local stress concentrations was possible in order to detect whether there was any internal damage developed in the ROI or not. For [0/90] specimen, it can be seen that at time of 0.65ms after impact, event ② there is localized damage at the top side of the specimen beneath the point of impact in high compressive impact zone probably due to fiber breakage of the 0 • layers. The main cracks begin to propagate mainly along the 90 • yarn direction until linking up with the pullout of the fiber's procedure by the tensile stress. In the case of [±45] 3s, the damage was primarily due to matrix cracking and delamination. Consequently, fiber/matrix interfacial fracture occurred in the ±45 layers due to high free-edge interlaminar stress in the tensile side (at the bottom of the specimen), propagating through the specimen until it reached the compressive side.

Three-dimensional finite element simulation of the impact test

A user subroutine VUMAT integrated to ABAQUS/Explicit is implemented to precisely capture the dynamic progressive failure of the laminate composite under the impact test (Fig. 18). In the proposed model, each ply is assumed to be homogeneous orthotropic linear elastic material, within damage and failure are accounted with a Hashin 3D damage-based behavior. A bilinear cohesive contact model with quadratic stress failure criteria in ABAQUS is employed between each composite ply with a zero thickness to capture the onset and progression of inter-laminar delamination. A penalty contact algorithm in ABAQUS/ Explicit with a friction coefficient of 0.3 is adopted to describe the tangential behavior for the contact between the impactor and top surface of specimens, while hard contact method is employed to describe the normal behavior. The frictionless contact is adopted to simulate the contact between the specimens and the support. The impactor is modelled with 3D rigid shell elements (R3D4) initially placed at the top center of the specimen with no initial gap. A 3.14Kg mass and 2.9m/s initial velocity are assigned to a reference point (RP) of the impactor. The movement of the impactor is driven according to the experimental velocity response and allowed only in the impact direction (along the 2axis) in Fig. 14. The composite laminate ply was modelled by eight-node solid deformable continuum elements (C3D8R) with three elements through the thickness of each ply. The contact area between the specimen and the impactor has been meshed with 0.25 mm maximum edge length elements which is equal to single element per lamina in the thickness direction. A coarse mesh is considered in other rigions of the finite element model to reduce the computational effort. These mesh sizes were obtained after performing a mesh convergence study, which is illustrated in Fig. 15 for the maximum contact force.

Intra-laminar model

Damage in composites materials can be qualified using stress or strain-based criteria or suitable polynomial criteria like Tsai-Wu or Tsai-Hill [START_REF] Tsai | A General Theory of Strength for Anisotropic Materials[END_REF]. However, the polynomial failure criteria are not applied at the ply level and are only used to predict the failure envelope of the laminate subjected to different multi-axial loading. Hashin [24] proposed failure criteria applied at the ply level to predict four major damage modes in laminates composite.

Based on the continuum damage method (CDM), the intra-laminar damage model presented in this section comprises failure criteria and

Table 3

Equivalent displacement and stress of each failure mode. Compression in direction 2 (σ22≤ 0)

δ2c,eq = lc ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 〈-ε22〉 2 + ε12 2 + ε23 2 √ σ2c,eq = lc[〈 -σ22〉〈 - ε22〉 + ε12σ12 + ε23σ23]/ δt2,eq F t 1 = ( σ 11 X t ) 2 + σ 12 2 S 2 12 + σ 13 2 S 2 13
(2)

Compressive failure direction 1 (σ 11 ≤ 0)

F c 1 = ( σ 11 X c ) 2 (3)
Tensile failure direction 2 (σ 22 + σ 33 ≥ 0) (5)

F t 2 = ( σ 22 + σ 33 Y t
F c 2 = ( ( Y c 2S 23 ) 2 - 1
) ( σ 22 + σ 33 Y c ) + ( σ 22 + σ 33 2S 23
In the above Equations, σ ij (i, j = 1, 2, 3) are the Cauchy stresses tensor components calculated in material coordinates, the tensile and compressive allowable strengths for lamina are denoted by subscripts t and c, respectively. X t , Y t denote the allowable tensile strengths in longitudinal and transverse directions. Similarly, X c , Y c denote the allowable compressive strengths in the longitudinal and transverse direction. Further, S 12 , S 13 and S 23 present the in-plane and out-of-plane shear strengths of laminate composite ply, respectively. F i j = 1 (i = c, t and j = 1, 2) indicate the onset of failure in the corresponding mode.

Damage evolution laws

In order to model the structural stiffness loss due to damage of the material in elements where the index criterion F t 1 , F c 1 , F t 2 and F c 2 are satisfied, a damage evolution law has to be defined. In the present study, the evolution of each damage variable is governed by energy-based damage evolution methods with linear degradation; this crack band model was developed by Bazant and Oh [START_REF] Banzah | Crack band theory for fracture of concrete[END_REF], successfully implemented by Zhou et al. [START_REF] Zhou | Finite element analysis of a modified progressive damage model for composite laminates under low-velocity impact[END_REF] and [START_REF] Zhang | Finite element modeling of damage development in cross-ply composite laminates subjected to low velocity impact[END_REF][START_REF] Li | Assessment of failure criteria and damage evolution methods for composite laminates under low-velocity impact[END_REF][START_REF] Jun | Progressive damage and nonlinear analysis of 3D four-directional braided composites under unidirectional tension[END_REF] in their damage models.

G I = 1 2
σ f eq ε f eq l c [START_REF] Tuo | Damage and failure mechanism of thin composite laminates under low-velocity impact and compression-after-impact loading conditions[END_REF] Where G I is the fracture energy density of failure mode I; σ f eq and ε f eq are respectively the equivalent peak stress and equivalent failure strain; l c is the characteristic length of the element. The characteristic length l c , which equals the cube root of the element volume, is introduced to alleviate mesh dependency during material softening in the finite element model.

The following equation can describe the damage evolution variable for each failure mode. The damage variable may evolve in such a way so that the stress-displacement behaves as shown in Fig. 16.

d I = δ f I,eq ( δ I,eq -δ 0 I,eq ) δ I,eq ( δ f I,eq -δ 0 I,eq ) (d I ∈ [0, 1], I = 1t, 1c, 2t, 2c) (7)
where δ 0 I,eq is the equivalent displacement for damage initiation; δ f I,eq is the equivalent displacement at final failure. The above relation is presented graphically in Fig. 16 (blue line), and they can be computed as fellows:

δ 0 I,eq = δ I,eq ̅̅̅̅ ̅ F I √ (8)
δ f I,eq = 2G I σ 0 I,eq
; σ 0

I,eq = σ I,eq ̅̅̅̅ ̅

F I √ (9)
where F I is the value of damage initiation criteria; σ I,eq and δ I,eq are the equivalent stress and equivalent displacement for different failure

G 1+ f = G 2+ f (N/mm) G 1- f = G 2- f (N/mm)
45 (45-70) 39 (39 -50) [START_REF] Yang | In-plane compression response of woven CFRP composite after low-velocity impact: Modelling and experiment[END_REF][30]

[28] [START_REF] Kinvi-Dossou | A numerical homogenization of E-glass/acrylic woven composite laminates: Application to low velocity impact[END_REF] Interface properties N = S = T (MPa)

G C I (N/mm) G C II = G C III (N/
[C d ] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ (1 -d 1 )C 0 11 (1 -d 1 )(1 -d 2 )C 0 12 (1 -d 1 )(1 -d 2 )C 0 13 0 0 0 (1 -d 1 )(1 -d 2 )C 0 12 (1 -d 1 )(1 -d 2 )C 0 22 (1 -d 1 )(1 -d 2 )C 0 23 0 0 0 (1 -d 1 )(1 -d 2 )C 0 13 (1 -d 1 )(1 -d 2 )C 0 23 (1 -d 1 )(1 -d 2 )C 0 33 0 0 0 0 0 0 (1 -d 3 )C 0 44 0 0 0 0 0 0 (1 -d 3 )C 0 55 0 0 0 0 0 0 (1 -d 3 )C 0 66 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (10)
modes are listed in Table 3, where the symbol < > represents the Macaulay operator 〈x〉 = (x + |x|)/2. Before the initiation of damage, the material behavior is linear elastic, and the stress and strain are related as σ = [C]ε where [C] denotes the elastic matrix. After initiation of damage, the elastic matrix [C d ] of the damaged material is defined as follows: Where C 0 ij are the undamaged stiffness matrix components.

C 0 11 = 1 -v 23 v 32 E 2 E 3 Δ = 1 -v 13 v 13 E 1 E 3 Δ = C 0 22 ; C 0 33 = 1 -v 12 v 21 E 2 E 1 Δ C 0 12 = v 21 -v 31 v 23 E 2 E 3 Δ ; C 0 13 = v 31 -v 21 v 32 E 2 E 3 Δ ; C 0 23 = v 32 -v 12 v 31 E 1 E 3 Δ Δ = 1 -v 12 v 21 -v 23 v 32 -v 13 v 31 -2v 21 v 32 v 13 E 1 E 2 E 3 C 0 44 = 2G 23 = 2G 13 = C 0 55 C 0 66 = 2G 12 (11)
And the expressions d 1 , d 2 and d 3 are respectively the damage factors in longitudinal, transverse and shear direction [START_REF] Zhou | Finite element analysis of a modified progressive damage model for composite laminates under low-velocity impact[END_REF][START_REF] Raimondo | A progressive failure model for mesh-size-independent FE analysis of composite laminates subject to low-velocity impact damage[END_REF]. The damage variable d 1 , d 2 and d 3 are defined as follows.

d 1 = 1 -(1 -d 1t )(1 -d 1c ) d 2 = 1 -(1 -d 2t )(1 -d 2c ) d 3 = 1 -(1 -d 1t )(1 -d 1c )(1 -S t d 2t )(1 -S c d 2c ) (12)
Where d 1t , d 1c , d 2t , d 2c are the damage variable for each failure mode under tensile and compression loads calculated by the damage evolution model using Eq. [START_REF] Lopes | Low-velocity impact damage on dispersed stacking sequence laminates . Part II : Numerical simulations[END_REF].

Additionally, the factors S t and S c in shear damage variable d 3 are introduced to control the loss of shear stiffness due to tensile and compressive matrix failures. In this study, these factors are assumed to be S t = 0.9 and S c = 0.5 [START_REF] Zhou | Finite element analysis of a modified progressive damage model for composite laminates under low-velocity impact[END_REF].

Inter-laminar damage modelling

To predict the delamination between consecutive plies of the lami-

⎨ ⎩ t n t s t t ⎫ ⎬ ⎭ = ⎡ ⎣ k 1 0 0 0 k 2 0 0 0 k 3 ⎤ ⎦ ⎧ ⎨ ⎩ δ 1 δ 2 δ 3 ⎫ ⎬ ⎭ ( 13
)
Where t is the traction stress, δ is the separation displacement, and k represent the initial stiffness of the interface.

Therefore, a quadratic failure criterion and the Benzeggagh-Kenane (BK) criterion are used to predict the initiation and propagation of inter-laminar damage, as shown in Eqs. ( 15)- [START_REF] Ramakrishnan | Experimental assessment of low velocity impact damage in flax fabrics reinforced biocomposites by coupled high-speed imaging and DIC analysis[END_REF].

{ 〈t n 〉 N } 2 + { 〈t s 〉 S } 2 + { 〈t t 〉 T } 2 = 1 (14)
Where t n , t s and t t are the normal and shear tractions; N, S and T denote the interface normal and shear strength.

G c = G c n + ( G c s -G c n ) { G s G T } η ( 15
)
with G c , G c s and G c n are the total, shear and normal critical fracture energy respectively; G s is the dissipated energy in the out-of-plane direction; G T is the total dissipated energy in all three direction;η is the relevant material coefficient in the B-K formula.

The detailed material properties of the biaxial glass fiber reinforced laminate composite used in the present FE simulation are listed in Table 4. Any strain rate dependency of those properties is neglected due to the relatively low deformation speeds occurring in this low-velocity impact study.

Since no reliable source was found for S 23 , the value used was supposed to equal S 12 and S 13 .

Impact tests results and discussion

Global mechanical response

The good performance of the proposed FE model is highlighted in Fig. [START_REF] Flores | Highspeed 3D digital image correlation of low-velocity impacts on composite plates[END_REF]. In which the typical force-time and energy-time responses of [0/ 90] 3s and [±45] 3s laminated composite subject to impact test obtained from experiment results are compared to that predicted by FE. Little governed by compression damage variables at the top side and the tensile damage variable at the bottom side of the specimen during the impact test.

Conclusion

In this paper, an original method combining low velocity impact tests, high-speed imaging, multiscale homogenization and a progressive finite element damage model was developed to characterize and predict the impact damage behaviour of a composite sample under a complex stress state due to the three-point bending impact test. Two main lay-up composite architectures ( [0/90] 3s and [±45] 3S ) have been studied. The results have shown the influence of the lay-up architecture on the impact resistance of (NCF) GFRP laminates. [±45] 3S laminates exhibit lower peak force, longer duration of impact force and higher absorbed energy than [0/90] 3S .

High-speed imaging coupled with the DIC technique was successfully used to measure the deformation of the specimen and to track the evolution of damage during impact. It has been shown that damage in [0/ 90] 3S laminate occurs under the loading contact point, dominated by delamination between adjacent layers and matrix cracking in the compressive side and propagate through the specimen until reaching the tensile side. For [±45] 3S laminates, the damage is primarily due to matrix cracking and delamination, and consequently, fiber/matrix interfacial fracture occurs in the ±45 layers as a result of high free-edge interlaminar stress in the tensile side.

A numerical finite element model accounting for Intra-laminar and Inter-laminar damages has been developed to simulate this test and understand how damage factors evolve inside a specimen. The experimental impact curves have been compared with numerical simulation to evaluate the accuracy of the finite element model. This comparison showed that the proposed method could predict the impact behaviour in terms of peak force, duration of impact contact time, absorbed energy and delamination area.

difference is found among the numerical results of the load-time curve to experiment result. For the energy-time curve, the initial kinetic energy of the impactor begins to transfer to the specimen as soon as the contact occurs. During the rising part of the impact test, the kinetic energy is absorbed by the elastic deformation of the sample; the numerical and experimental energy-time curve shows a good correlation before load achieve the maximum value. After the load peak point, the energy is dissipated by intra-laminar damage and delamination. Compared to the experimental results, the predicted energy curve decreases and shows a difference from the measured value. In the numerical model, the friction coefficients used in the contact between the impactor and the specimen and the properties of the cohesive elements affect absorbed energy and may be the cause of this difference. Moreover, the predicted distribution of delamination shows that the delamination areas, including the complete and partially failed region and the separation between adjacent ply, are comparable to the experimental phenomenology observed Fig. 19e.

Local mechanical response: FE results compared to experimental DIC measurements

Figs. 20 and 21 shows the transverse displacement U2 distribution in the central region for [0/90] 3S and [±45] 3s laminate composite. The magnitude of displacement is high at the center of the laminate in both FE and DIC results, and they are in good correlation.

Finite element damage evolution and failure mechanism process case of [0/90] laminate

It is worth mentioning that once the failure initiates, damage evolution is governed by fracture energy dissipated during failure. That means the damage variable keeps evolving from 0 to 1 irreversibly after damage occurrence, and the stiffness of the damaged element degrades accordingly until the complete failure of the element. Figs. 22 and23 shows deformation and damage development for the two orthogonal directions of the NCF composite under tension and compression loading, respectively. For both transverse tensile damage d 2t and compression damage d 2c , damage occurs initially in the center line at the impact zone Fig. 23, initiated by delamination under compression. The transition from compression to tension passing through specimen width results in shear stress through the sample thickness. The damage progresses with the increase of external load, and the damage areas expand and connect with those of neighbouring zones. As seen in the illustration Fig. 22, the longitudinal tensile damage occurs and propagates on the bottom of the specimen. It may also be produced compressive damage occurring on the top of the specimen produced by the bending moment.

The delamination damage consists of important cracks between plies. During the impact test, the surface under the center loading point exhibits extreme matrix deformation and Mode I of fracture is the most thoroughly characterized fracture mode in impact on edge for laminate composite. Fig. 24 shows the delamination profile between adjacent layers for the sample during the impact test using the cohesive surface damage model CSDMG.

As shown in Fig. 25, all the damage variables factors rise to the peak at the same moment t = 0.7 ms for the [0/90] composite laminate. The high compressive zone in the top side of the specimen was selected to present the time history for transverse tensile and compression damage variables, and the bottom side in the high tensile zone was selected to present the longitudinal tensile and compression damage variables. As shown in Fig. 25(a), unlike the tensile damage variable dt2, the transverse compression damage variable dc2 rises to the max value 1. This is because compressive deformation is more significant than tensile deformation at the top side of the specimen for the impact test. In the meantime, at the bottom side of the specimen, the longitudinal tensile damage variable is easier to generate, and it propagates faster than the compression damage variable. Consequently, the loss of fitness is

Fig. 1 .

 1 Fig. 1. (a) biaxial NCF glass fiber fabric and cutting direction (b) Impact specimen dimension.

Fig. 2 .

 2 Fig. 2. SEM observations. (a) yarns identification and measurements (b) fiber repartition (c) fiber repartition after segmentation (d) fiber diameters distribution.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Experimental setup of impact test and high-speed imaging.

Fig. 5 .

 5 Fig. 5. Yarn model (a) 2D view (b) 3D model dimensions (a 1 = 1.91 μm, a 2 = 7.66 μm and a 3 = 13.267 μm).

Fig. 6 .

 6 Fig. 6. Full view of basic cell with yarn bundles and dimensions in mm.

Fig. 7 .

 7 Fig. 7. Curves of load versus time (a) [0/90], (b) [±45] and energy versus time (c) [0/90], (d) [±45], under impact test.

Fig. 8 .

 8 Fig. 8. Absorbed energy under impact test for [0/90] 3s and [±45] 3s .

Fig. 9 .

 9 Fig. 9. Evolution of failure on [0/90] 3S specimen, under impact test.

Fig. 10 .

 10 Fig. 10. Evolution of failure on [±45] 3s specimen, under impact test.

Fig. 11 .

 11 Fig. 11. DIC transverse displacement field at a different event in [0/90] 3S specimen, under impact test.

Fig. 12 .

 12 Fig. 12. DIC transverse displacement field at a different event in [±45] 3s specimen, under impact test.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Progressive Exx DIC strain of the GFRP under impact test for [0/90] and [±45].

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Maximum contact force versus contact mesh size.

〈ε11〉 2 + ε12 2 + ε13 2 √〈ε22〉 2 + ε12 2 + ε23 2 √

 2222 σ1t,eq = lc[〈σ11〉〈ε11〉 + ε12σ12 + ε13σ13]/δt1,eq Compression in direction 1 (σ11≤ 0) δ1c,eq = lc〈 -ε11〉 σ1c,eq = lc[〈 -σ11〉〈 -σ2t,eq = lc[〈σ22〉〈ε22〉 + ε12σ12 + ε23σ23]/δt2,eq

  2 (σ 22 + σ 33 ≤ 0)

Fig. 17 .

 17 Fig. 17. Traction-separation Law.

Fig. 18 .

 18 Fig. 18. Numerical modelling process.

Fig. 19 .

 19 Fig. 19. Comparison between numerical and experimental impact (a)-(c) force-time curve [0/90] 3s , [±45] 3s (b)-(d) energy-time [0/90] 3s , [±45] 3s (e) delamination area: view from the top side of the sample.

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Comparison of transverse displacement from FE simulation with DIC displacement for [0/90] 3S .

Fig. 22 .

 22 Fig. 22. Damage failure: longitudinal direction (1) in tension and compression.

Fig. 23 .

 23 Fig. 23. Damage failure: transverse direction (2) in tension and compression.

  . The constitutive relationship of the interface is given by ⎧

Fig. 24 .

 24 Fig. 24. Delamination growth during the impact test (top view).

Fig. 25 .

 25 Fig. 25. Damage variable versus time history (obtained from the element at the center of the high compressive zone (a) and element at the center of the high tensile zone at the bottom of sample (b)).

Table 4

 4 Material parameters.

	Property	Value	References
	Density (kg//m 3 )	1538	Measured
	Elastic properties	23.010 ± 0.3	Measured
	E1 = E2(GPa)	11.380	FE homogenization
	E3 (GPa)	5.768 ± 0.6	Measured
	G12 (GPa)	3.680	FE homogenization
	G13 = G23 (GPa)	0.14	Measured
	v12	0.38	FE homogenization
	v13 = v23
	Strength
	σ1t = σ2t (MPa)	450	Measured
	σ1c = σ2c (MPa) S12 (MPa) S13 = S23 (MPa)	280 (260-360) 119 119	[28][27][30] Measured Estimated
	Intra-laminar fracture energies

Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations.

CRediT authorship contribution statement

Ahmed Boukar: Conceptualization, Methodology, Investigation, Writingoriginal draft. Stephane Corn: Conceptualization, Methodology, Investigation, Supervision. Pierre R.L. Slangen: Conceptualization, Methodology, Investigation, Supervision. Patrick Ienny: Conceptualization, Methodology, Investigation, Supervision.

Declaration of Competing Interest

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. We confirm that the manuscript has been read and approved by all named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. We further confirm that the order of authors listed in the manuscript has been approved by all of us. We understand that the Corresponding Author is the sole contact for the Editorial process (including Editorial Manager and direct communications with the office). He is responsible for communicating with the other authors about progress, submissions of revisions and final approval of proofs. We confirm that we have provided a current, correct email address which is accessible by the Corresponding Author and which has been configured to accept email from stephane.corn@minesales.fr

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ijimpeng.2022.104218.