Numerical–Experimental Correlation of Impact-Induced Damages in CFRP Laminates

Andrea Sellitto 1,* , Salvatore Saputo 1, Francesco Di Caprio 2, Aniello Riccio 1 , Angela Russo 1 and Valerio Acanfora 1 1 Department of Engineering, University of Campania Luigi Vanvitelli, via Roma 29, 81031 Aversa, Italy; salvatore.saputo@unicampania.it (S.S.); aniello.riccio@unicampania.it (A.R.); angela.russo@unicampania.it (A.R.); valerio.acanfora@unicampania.it (V.A.) 2 C.I.R.A.—Italian Aerospace Research Centre, via Maiorise snc, 81043 Capua, Italy; f.dicaprio@cira.it * Correspondence: andrea.sellitto@unicampania.it; Tel.: +39-081-5010-407


Introduction
Composite materials are characterized by high mechanical properties values if compared to classic metallic materials. However, their damage behavior is difficult to predict and control. As a matter of fact, impact events on a composite structure are still a "hot topic" for research, nowadays. The different impact phenomena can be divided into three main groups, based on the velocity of the impactor: low velocity impacts, ballistic impacts, and hyper-velocity impacts. In particular, Low Velocity Impacts (LVIs) often result in Barely Visible Impact Damages (BVIDs), which can be hardly detected by the naked eye, but may relevantly decrease the material strength [1][2][3]. Damages occurring as a consequence of low velocity impacts include both intra-laminar (fiber and/or matrix failure) and inter-laminar damages (delaminations). Among the others, delaminations, which can growth under service loads, can seriously reduce the load carrying capability of the overall structure [4,5]. For this reason, a large number of studies, focused on the prediction of the damages (especially delaminations) induced by LVI, can be found in the literature. In these studies, different approaches, both analytical-numerical and between the impacted area and the rest of the specimen. The results, in terms of intra-laminar and inter-laminar damages, have been compared with the experimental ones, obtained by means of Ultrasonic C-Scan tests. Further numerical-experimental correlations on the impactor displacement and velocity, the internal energy, and the force exerted during the impact have been performed as well.
Numerical results, obtained considering different Abaqus element types, were compared to each other. The Finite Element Model has been discretized by means of SC8R Continuum shell elements and C3D8R Solid elements [36]. Inter-laminar damages were modeled by a Cohesive Zone Model based approach, while intra-laminar damage was modeled by means of Hashin's failure criteria and gradual material properties degradation rules. However, the intra-laminar damage model was only provided for Continuum shell elements in Abaqus. Since no intra-laminar damage model was provided for Solid elements, a user-defined material model, able to take into account the intra-laminar damage, was implemented in a user subroutine VUMAT.
In Section 2, the theoretical background on the analytical model and on the intra-laminar damage model implemented in the VUMAT are presented. In Section 3, the test cases are introduced, while in Section 4 the results are presented and discussed.

Classical Laminated Plate Theory
A multi-degree of freedom analytical model [1,2,37] was introduced in order to study the impact event. A simply supported plate was analytically solved by using the Classical Laminated Plate Theory. The solutions for the analytical approach were obtained in terms of contact force, maximum displacement, impactor velocity, and kinetic energy. Then, a system of ordinary differential equations (ODE), representing the simply supported plate, coupled with a differential equation representing the impactor motion, was obtained.
A laminated composite plate, with in-plane dimensions a and b, respectively, along the x-and y-directions, was considered. The plate was considered simply supported along the four edges, and a concentrated force was applied at the center (x 0 = a/2, y 0 = b/2). Under these assumptions, the equation of motion along the transverse direction can be expressed as the following system of ordinary differential equations [1]: In Equation (1), F is a concentrated applied force, while ω mn , which is the natural frequency of the system for the strain modes m,n, is given by: ω mn = π 4 a 4 I 1 (D 11 m 4 + 2(D 12 + 2D 66 )m 2 n 2 r 2 − 4D 16 m 3 nr − 4D 26 mn 3 r 3 + D 22 n 4 r 4 ), (2) where D ij are the terms of the bending stiffness matrix, I 1 is the mass per unit length, and r = a/b. As reported in [1], when both values of m and n are negative, the effects of deflection in the plate become negligible. The choice of high positive values for both values is related to the phenomenon being analyzed: as the values of m and n increase, the rotary inertia and shear deformation effects becomes significant.
Moreover, the following Equation (3), governing the impactor dynamic, has been considered: where M and .. z(t) are, respectively, the impactor mass and acceleration, F c (t) is the contact force, and t is the time variable.
Appl. Sci. 2019, 9,2372 4 of 21 A proper contact law has been chosen to couple Equations (3) and (1). In particular, the contact force F c (t) has been related to the indentation α = z(t) − w 0 (x 0 ,y 0 ,t), expressed as the difference between the impactor displacement z(t) and the laminate central point displacement w 0 (x 0 ,y 0 ,t). In this work, a linearized form of the Hertz Contact Law was adopted [37], therefore the contact force history, exerted during the impact phenomenon, can be expressed as: where k y represents the linearized contact stiffness. The boundary conditions of the system of ODEs are: α mn (0) = 0, where V is the initial velocity of the impactor.

Composites Damage Criteria
In order to account for the damage behavior of the laminate, a User Material was developed in the ABAQUS FEM environment capable to consider the fiber and matrix damage mechanisms in tension and compression [4,19,29]. Each failure mode consisted of two different phases. Initially, a linear mechanical behavior is considered up to the damage initiation threshold. Then, the damage evolution up to the complete damage is evaluated.
The criteria used for evaluate the damage initiation are based on Hashin's (fiber and matrix tension and fiber compression) and Puck-Shurmann [38] (matrix compression) failure criteria, as reported in Equations (7)- (10).
where X T and X C are the fiber tensile and compressive strengths, Y T is the matrix tensile strength, and S 12 and S 23 are respectively the longitudinal and transverse shear strengths. Moreover, the parameters introduced in Equation (10) have been defined with respect to the potential fracture plane; in particular, S A 23 is the transverse shear strength in the potential fracture plane, µ nt and µ nl are respectively the friction coefficients in the transverse and longitudinal directions, σ nn is the normal stress, and τ nt and τ nl are respectively the shear stresses in the transverse and longitudinal directions. Indeed, experiments have demonstrated that unidirectional laminates experience shear fail [39] when subjected to transverse compressive loads, with a fracture plane oriented with an angle θ f = 53 • ± 2 • [21]. Hence, to correctly evaluate the damage status of the laminate, the stress/strain components in the θ f oriented fracture plane have been considered, instead of the nominal 0 • oriented nominal plane. Therefore, the stress components in the general fracture plane L-N-T oriented with a fracture angle θ f , as in Figure 1, can be obtained as a function of the stress components defined in the lamina coordinate system 1-2-3: Hence, the parameters related to the fracture angle introduced in Equation (10) can be defined [40]: where YT is the matrix compressive strength. Each mode is characterized by the failure behavior reported in Figure 2. In particular, a linear elastic behavior with an initial stiffness K can be observed up to location A, which identifies the damage onset. Then, each location B on the segment AC identifies a partial damage condition characterized by a degraded stiffness Kd, up to location C, where the complete failure can be observed. Hence, the parameters related to the fracture angle introduced in Equation (10) can be defined [40]: where Y T is the matrix compressive strength. Each mode is characterized by the failure behavior reported in Figure 2. In particular, a linear elastic behavior with an initial stiffness K can be observed up to location A, which identifies the damage onset. Then, each location B on the segment AC identifies a partial damage condition characterized by a degraded stiffness K d , up to location C, where the complete failure can be observed. Appl. Sci. 2019, 9, x 6 of 22 To take into account the material stiffness degradations in locations B, the damage coefficients di [36] are introduced for each failure mode: , , , Indeed, the damaged stiffness matrix CD can be expressed as a function of the damage coefficients as: , and df, dm, and ds are respectively the fiber, matrix, and shear damage coefficients: However, this methodology has been found strongly dependent on the domain discretization. Hence, to reduce mesh dependencies, the element characteristic length LC has been introduced. In this work, the solution proposed by Bazant and Oh [39] and reported in Equation (19) has been adopted: where Aip is the Area of the element corresponding to the ip-th integration point and θ is the fracture angle. Then, the equivalent displacements was computed as a function of the element characteristic length LC: Constitutive relation adopted for fiber and matrix failure modes in tension and compression.
To take into account the material stiffness degradations in locations B, the damage coefficients d i [36] are introduced for each failure mode: Indeed, the damaged stiffness matrix C D can be expressed as a function of the damage coefficients as: where , and d f , d m , and d s are respectively the fiber, matrix, and shear damage coefficients: However, this methodology has been found strongly dependent on the domain discretization. Hence, to reduce mesh dependencies, the element characteristic length L C has been introduced. In this work, the solution proposed by Bazant and Oh [39] and reported in Equation (19) has been adopted: where Aip is the Area of the element corresponding to the ip-th integration point and θ is the fracture angle. Then, the equivalent displacements was computed as a function of the element characteristic length L C : where G T IC and G C IC are the tensile and compressive Mode I fracture toughness, G T IIC and G C IIC are the tensile and compressive Mode II fracture toughness, and:

Analytical Set-Up
The analyzed test case is a simply supported laminated plate impacted by a 3.58 kg mass steel sphere. A 10 J impact energy is considered, leading to an impactor initial velocity equal to 2.364 m/s. The plate is made by a unidirectional carbon-fiber laminate, with in-plane dimensions in the x-and y-directions, respectively, a = 150 mm and b = 100 mm. The stacking sequence of the laminate is [45/-45/0/0/90/0] 2S , while the mechanical properties of the adopted material system and the ply thickness are reported in Table 1. A spherical impactor with an 8 mm in radius was considered. The linear contact stiffness was set equal to k y = 7.8755 × 10 3 N/mm, according to the equations reported in [34]. In these analyses, aimed to preliminary assess the elastic behavior of the structure under impact loading conditions, only the elastic behavior has been taken into account. Therefore, no damage has been considered.

Experimental Set-Up
The experimental impact test was performed according to the ASTM D7136 standard [34], as shown in Figure 3. A spherical impactor with an 8 mm in radius was considered. The linear contact stiffness was set equal to ky = 7.8755 ×•10 3 N/mm, according to the equations reported in [34]. In these analyses, aimed to preliminary assess the elastic behavior of the structure under impact loading conditions, only the elastic behavior has been taken into account. Therefore, no damage has been considered.

Experimental Set-Up
The experimental impact test was performed according to the ASTM D7136 standard [34], as shown in Figure 3. The three tested specimens were cut from a single CFRP plate by using the waterjet technique. The CFRP plate is composed of prepreg high modulus carbon fiber and thermoset epoxy resin, cured in autoclave. The test was performed with the CEAST Fractovis plus drop tower test machine. The plates were normally impacted at the center. The impactor was a 16 mm diameter hemisphere made of hardened steel, with a mass of 3.58 kg. A 10 J impact energy was considered. Ultrasonic inspections were performed on the specimen prior to the impact test to assure the lack of manufacturing defects, and after the experimental test to evaluate the damaged area.

Preliminary Analyses
The analytical approach results expressed in terms of impactor displacement, contact force, impactor velocity, and kinetic energy versus time are shown respectively in Figures 4-7. The results were obtained by considering an increasing number of modes. The three tested specimens were cut from a single CFRP plate by using the waterjet technique. The CFRP plate is composed of prepreg high modulus carbon fiber and thermoset epoxy resin, cured in autoclave. The test was performed with the CEAST Fractovis plus drop tower test machine. The plates were normally impacted at the center. The impactor was a 16 mm diameter hemisphere made of hardened steel, with a mass of 3.58 kg. A 10 J impact energy was considered. Ultrasonic inspections were performed on the specimen prior to the impact test to assure the lack of manufacturing defects, and after the experimental test to evaluate the damaged area.

Preliminary Analyses
The analytical approach results expressed in terms of impactor displacement, contact force, impactor velocity, and kinetic energy versus time are shown respectively in             Figure 6 shows the trends of the impact velocity considering different m and n values. Although the velocity had the same initial and final values, different m and n led to different velocity trends. Moreover, as a direct consequence of the variation of the impact velocity, the kinetic energy, shown in Figure 7, showed different behavior when low values of m and n were considered. However, as m and n increase, the graphs shown in Figures 4-7 appear to be perfectly overlapped.
Then, the numerical results were compared with the analytical ones. The m,n = 10 analytical results were chosen as the reference analytical solution.
In order to define the numerical model which best-fit the analytical results, several parameters were investigated, such as the element type (reduced integration continuum shell elements SC8R and reduced integration solid elements C3D8R), the in-plane element size (1 mm, 2 mm, and 4 mm), and  Figure 6 shows the trends of the impact velocity considering different m and n values. Although the velocity had the same initial and final values, different m and n led to different velocity trends. Moreover, as a direct consequence of the variation of the impact velocity, the kinetic energy, shown in Figure 7, showed different behavior when low values of m and n were considered. However, as m and n increase, the graphs shown in Figures 4-7 appear to be perfectly overlapped.
Then, the numerical results were compared with the analytical ones. The m,n = 10 analytical results were chosen as the reference analytical solution.
In order to define the numerical model which best-fit the analytical results, several parameters were investigated, such as the element type (reduced integration continuum shell elements SC8R and reduced integration solid elements C3D8R), the in-plane element size (1 mm, 2 mm, and 4 mm), and the number of elements along the thickness direction. In particular, Continuum Shell elements were discretized with 1, 12, and 24 elements in the thickness direction, while Solid elements were discretized with 24 elements in the thickness direction (one element per ply). Hence, more than one orientation is associated to the continuum shell elements belonging to the models characterized by less than 24 elements in the thickness direction.
In the numerical analyses, the plate was considered simply supported on the edges, as reported in Figure 8. In the framework of preliminary analyses, aimed to compare the analytical model with the numerical ones to determine the numerical parameters which best represent the mechanical behavior of the plate subjected to low velocity impact, both the inter-laminar and intra-laminar damages were neglected.
Although the velocity had the same initial and final values, different m and n led to different velocity trends. Moreover, as a direct consequence of the variation of the impact velocity, the kinetic energy, shown in Figure 7, showed different behavior when low values of m and n were considered. However, as m and n increase, the graphs shown in Figures 4-7 appear to be perfectly overlapped.
Then, the numerical results were compared with the analytical ones. The m,n = 10 analytical results were chosen as the reference analytical solution.
In order to define the numerical model which best-fit the analytical results, several parameters were investigated, such as the element type (reduced integration continuum shell elements SC8R and reduced integration solid elements C3D8R), the in-plane element size (1 mm, 2 mm, and 4 mm), and the number of elements along the thickness direction. In particular, Continuum Shell elements were discretized with 1, 12, and 24 elements in the thickness direction, while Solid elements were discretized with 24 elements in the thickness direction (one element per ply). Hence, more than one orientation is associated to the continuum shell elements belonging to the models characterized by less than 24 elements in the thickness direction.
In the numerical analyses, the plate was considered simply supported on the edges, as reported in Figure 8. In the framework of preliminary analyses, aimed to compare the analytical model with the numerical ones to determine the numerical parameters which best represent the mechanical behavior of the plate subjected to low velocity impact, both the inter-laminar and intra-laminar damages were neglected.  In Figures 9-13, the numerical and analytical results, in terms of displacement of the impactor and contact force, are compared. In particular, the results shown in Figure 9 were obtained by using a shell element formulation, while the results in Figures 10-12 were obtained by using a continuum shell formulation considering, respectively, 1, 12, and 24 elements in the thickness direction. Finally, in Figure 13, the results obtained by using a solid element formulation, with 24 elements in the thickness direction, are introduced.
Appl. Sci. 2019, 9, x 11 of 22 In Figures 9-13, the numerical and analytical results, in terms of displacement of the impactor and contact force, are compared. In particular, the results shown in Figure 9 were obtained by using a shell element formulation, while the results in Figures 10-12 were obtained by using a continuum shell formulation considering, respectively, 1, 12, and 24 elements in the thickness direction. Finally, in Figure 13, the results obtained by using a solid element formulation, with 24 elements in the thickness direction, are introduced.  According to the Figures 9-13, Shell and Solid elements' models were not particularly influenced by the in-plane element size. On the other hand, Continuum shell elements' models were more sensitive to the in-plane element size as the number of elements in the thickness direction increased. In particular, the numerical results of the Continuum shell model with 12 and 24 elements in the thickness direction associated to coarser mesh (size 2 mm or 4 mm) differed noticeably from the analytical results. This is mainly due to hourglass problems, which can be relevant for the reduced integration scheme elements. Unfortunately, the reduced integration scheme is mandatory in Abaqus/Explicit when Continuum shell elements were adopted. To alleviate hourglass deformation problems that could arise in the numerical models, the default hourglass controls available in Abaqus /Explicit [36] was adopted.
In general, the results highlight that the best results were obtained by the Solid element model discretized with one element per ply and an in-plane element size equal to 1 mm; however, a good numerical-analytical correlation was obtained also by the Continuum shell mode, discretized by 12 elements in the thickness direction and an in-plane element size equal to 1 mm. Indeed, this last configuration was chosen for the analysis presented hereafter due to its excellent compromise between low computational cost and agreement with the analytical solution.
Despite their simplifications, the analytical models available in the literature well describe the behavior of a CFRP plate subjected to impact loading conditions. Hence, the numerical model was validated by comparisons with the analytical one, and the used modeling approach, neglecting in this preliminary phase the damages, has been assessed. Once the accuracy of the impact numerical model is verified, some considerations should be made about the relevant difference between the According to the Figures 9-13, Shell and Solid elements' models were not particularly influenced by the in-plane element size. On the other hand, Continuum shell elements' models were more sensitive to the in-plane element size as the number of elements in the thickness direction increased. In particular, the numerical results of the Continuum shell model with 12 and 24 elements in the thickness direction associated to coarser mesh (size 2 mm or 4 mm) differed noticeably from the analytical results. This is mainly due to hourglass problems, which can be relevant for the reduced integration scheme elements. Unfortunately, the reduced integration scheme is mandatory in Abaqus/Explicit when Continuum shell elements were adopted. To alleviate hourglass deformation problems that could arise in the numerical models, the default hourglass controls available in Abaqus/Explicit [36] was adopted.
In general, the results highlight that the best results were obtained by the Solid element model discretized with one element per ply and an in-plane element size equal to 1 mm; however, a good numerical-analytical correlation was obtained also by the Continuum shell mode, discretized by 12 elements in the thickness direction and an in-plane element size equal to 1 mm. Indeed, this last configuration was chosen for the analysis presented hereafter due to its excellent compromise between low computational cost and agreement with the analytical solution.
Despite their simplifications, the analytical models available in the literature well describe the behavior of a CFRP plate subjected to impact loading conditions. Hence, the numerical model was validated by comparisons with the analytical one, and the used modeling approach, neglecting in this preliminary phase the damages, has been assessed. Once the accuracy of the impact numerical model is verified, some considerations should be made about the relevant difference between the boundary conditions implemented in the analytical model and the ones implemented in the experimental test. Actually, simply supported conditions have been considered in the analytical and preliminary numerical models; however, the boundary conditions of the experimental test, suggested by the ASTM standards and reported in Figure 14, are substantially different. These boundary conditions are a consequence of a support fixture whose degrees of freedom are restrained and four rigid clamps, also modeled as rigid bodies, which hold the specimen still during impact and avoid any rebound. boundary conditions implemented in the analytical model and the ones implemented in the experimental test. Actually, simply supported conditions have been considered in the analytical and preliminary numerical models; however, the boundary conditions of the experimental test, suggested by the ASTM standards and reported in Figure 14, are substantially different. These boundary conditions are a consequence of a support fixture whose degrees of freedom are restrained and four rigid clamps, also modeled as rigid bodies, which hold the specimen still during impact and avoid any rebound. A further numerical model obtained with the implementation of the Experimental boundary condition (BC-E) was introduced to understand the effects of boundary conditions on the numerical results. In Figure 15, the analytical, experimental and numerical results of the shell model obtained by considering both analytical (BC-A) and experimental (BC-E) boundary conditions were compared. A further numerical model obtained with the implementation of the Experimental boundary condition (BC-E) was introduced to understand the effects of boundary conditions on the numerical results. In Figure 15, the analytical, experimental and numerical results of the shell model obtained by considering both analytical (BC-A) and experimental (BC-E) boundary conditions were compared. In Figure 15, a comparison between experimental, analytical, and numerical results, in terms of force as a function of the time, is reported. The comparison shows that in the first part of the curve both numerical and analytical models were in agreement with the experimental results. However, after a first phase, the analytical model and the FEM BC-A numerical model greatly differed from the experimental result. On the other hand, the FEM BC-E numerical model was characterized by the same mechanical response in terms of stiffness of the laminate, up to the maximum force peak. Indeed, neglecting the damages in the numerical model resulted in grater force peak, if compared with the experiment. Once the accuracy of the numerical model, in terms of elastic response, was established, the numerical models which provided the best results compared to the analytical solution was considered for the final impact test: final Continuum shell model and Final 3D solid model. Since Abaqus does not provide any intra-laminar damage models for Solid elements, a VUMAT was introduced. On the other side, Abaqus built-in Hashin's criteria were used for model discretized by means of Continuum shell elements. For both the final continuum shell model and the final solid FEM model, gradual material properties degradation rules have been adopted for intralaminar damage progression and inter-laminar damage were implemented through cohesive elements layers placed at plies interfaces. In Figure 15, a comparison between experimental, analytical, and numerical results, in terms of force as a function of the time, is reported. The comparison shows that in the first part of the curve both numerical and analytical models were in agreement with the experimental results. However, after a first phase, the analytical model and the FEM BC-A numerical model greatly differed from the experimental result. On the other hand, the FEM BC-E numerical model was characterized by the same mechanical response in terms of stiffness of the laminate, up to the maximum force peak. Indeed, neglecting the damages in the numerical model resulted in grater force peak, if compared with the experiment. Once the accuracy of the numerical model, in terms of elastic response, was established, the numerical models which provided the best results compared to the analytical solution was considered for the final impact test: final Continuum shell model and Final 3D solid model. Since Abaqus does not provide any intra-laminar damage models for Solid elements, a VUMAT was introduced. On the other side, Abaqus built-in Hashin's criteria were used for model discretized by means of Continuum shell elements. For both the final continuum shell model and the final solid FEM model, gradual material properties degradation rules have been adopted for intra-laminar damage progression and inter-laminar damage were implemented through cohesive elements layers placed at plies interfaces.

Final Solid Model
This Final Solid model was created by introducing two different meshed areas connect by a global-local conditions. Indeed, in the specimen impacted area, where most of the damage is supposed to develop (Figure 16) a ply-by-ply refined discretization was considered with each ply represented by a single layer of 3D solid elements C3D8R with a reduced integration scheme. On the contrary, the external global domain has been modeled with a coarser continuum shell mesh. With this discretization scheme, an accurate study can be performed by reducing, at the same time, the computational effort [40].
global-local conditions. Indeed, in the specimen impacted area, where most of the damage is supposed to develop (Figure 16) a ply-by-ply refined discretization was considered with each ply represented by a single layer of 3D solid elements C3D8R with a reduced integration scheme. On the contrary, the external global domain has been modeled with a coarser continuum shell mesh. With this discretization scheme, an accurate study can be performed by reducing, at the same time, the computational effort [40]. For 3D solid elements, the intra-laminar damage model introduced in Section 2.2 was applied by means of a VUMAT. Moreover, since the inter-laminar damages are likely to occur only between plies with different fiber orientation, zero-thickness cohesive layers were placed only in these locations. As highlighted in Figure 17, cohesive layers were connected to the adjacent layer by means of tie-constraint. Table 2 reports the mechanical properties of the cohesive material system, which were experimentally evaluated by means of three-point bending (En, Et, Es), double cantilever beam (Nmax, GIc), and end notched flexure tests (Tmax and Smax, GIIc and GIIIc).   For 3D solid elements, the intra-laminar damage model introduced in Section 2.2 was applied by means of a VUMAT. Moreover, since the inter-laminar damages are likely to occur only between plies with different fiber orientation, zero-thickness cohesive layers were placed only in these locations. As highlighted in Figure 17, cohesive layers were connected to the adjacent layer by means of tie-constraint. Table 2 reports the mechanical properties of the cohesive material system, which were experimentally evaluated by means of three-point bending (E n , E t , E s ), double cantilever beam (N max , G Ic ), and end notched flexure tests (T max and S max , G IIc and G IIIc ).
supposed to develop (Figure 16) a ply-by-ply refined discretization was considered with each ply represented by a single layer of 3D solid elements C3D8R with a reduced integration scheme. On the contrary, the external global domain has been modeled with a coarser continuum shell mesh. With this discretization scheme, an accurate study can be performed by reducing, at the same time, the computational effort [40]. For 3D solid elements, the intra-laminar damage model introduced in Section 2.2 was applied by means of a VUMAT. Moreover, since the inter-laminar damages are likely to occur only between plies with different fiber orientation, zero-thickness cohesive layers were placed only in these locations. As highlighted in Figure 17, cohesive layers were connected to the adjacent layer by means of tie-constraint. Table 2 reports the mechanical properties of the cohesive material system, which were experimentally evaluated by means of three-point bending (En, Et, Es), double cantilever beam (Nmax, GIc), and end notched flexure tests (Tmax and Smax, GIIc and GIIIc).    Both solid and cohesive elements have been discretized by 1 × 1 mm 2 elements. The external global domain, which has been modeled by one layer of Continuum shell SC8R elements along thickness direction, has been discretized with an element in-plane size of 4 mm. The specimen is constrained by a fixture and four clamps, modeled with R3D4 rigid elements. The impactor was modeled as a sphere of 16 mm diameter modeled with 3D Solid C3D8R elements. A surface-to-surface contact was used between the impactor and the specimen. During the analysis, the completely damaged elements were removed from the model to increase the computational efficiency. In particular, when cohesive elements were removed, plies can touch each other at the interfaces. Hence, a general contact on the whole model was defined (except for the already defined impact surfaces) to take into account this behavior. A friction coefficient with a value of 0.5 was introduced for both impactor-to-ply contact and ply-to-ply contact. Enhanced hourglass control was set to decrease hourglass problems related to the reduced integrated elements.

Final Continuum Shell Model
In the Final Continuum shell model, 12 continuum shell elements along thickness direction were used to model the specimen (Figure 18). For each element, a composite layup composed of two plies was defined.
thickness direction, has been discretized with an element in-plane size of 4 mm. The specimen is constrained by a fixture and four clamps, modeled with R3D4 rigid elements. The impactor was modeled as a sphere of 16 mm diameter modeled with 3D Solid C3D8R elements. A surface-to-surface contact was used between the impactor and the specimen. During the analysis, the completely damaged elements were removed from the model to increase the computational efficiency. In particular, when cohesive elements were removed, plies can touch each other at the interfaces. Hence, a general contact on the whole model was defined (except for the already defined impact surfaces) to take into account this behavior. A friction coefficient with a value of 0.5 was introduced for both impactor-to-ply contact and ply-to-ply contact. Enhanced hourglass control was set to decrease hourglass problems related to the reduced integrated elements.

Final Continuum Shell Model
In the Final Continuum shell model, 12 continuum shell elements along thickness direction were used to model the specimen (Figure 18). For each element, a composite layup composed of two plies was defined. The intra-laminar damages have been accounted by means of Hashin's failure criteria, while cohesive layers placed at plies interfaces and introduced in Table 2 were used to model inter-laminar damage. The contact laws introduced in Section 4.2 were adopted.

Numerical-Experimental Correlation
In this section, a comparison between the experimental results and the numerical ones for the both investigated configurations (continuum shell element and solid elements) is performed considering both intra-laminar and inter-laminar damages. Three 10 J impact experimental tests were performed. The damaged area was acquired by means of the ultrasonic C-Scan. Figure 19 shows the C-scan of one representative specimen. The intra-laminar damages have been accounted by means of Hashin's failure criteria, while cohesive layers placed at plies interfaces and introduced in Table 2 were used to model inter-laminar damage. The contact laws introduced in Section 4.2 were adopted.

Numerical-Experimental Correlation
In this section, a comparison between the experimental results and the numerical ones for the both investigated configurations (continuum shell element and solid elements) is performed considering both intra-laminar and inter-laminar damages. Three 10 J impact experimental tests were performed. The damaged area was acquired by means of the ultrasonic C-Scan. Figure 19 shows the C-scan of one representative specimen. A linear pattern obtained by means of Abaqus/Viewer was used to evaluate the experimental damaged area. The tested specimens were comparable in terms of force, displacement, and kinetic energy as a function of the time, and inter-laminar damages. In particular, the experimental damaged area was about 250 mm 2 , while the numerical damaged area, evaluated by using different numerical models, was in the range 166-609 mm 2 . The continuum shell model uses a lower number of interfaces between the plies, if compared with the solid element discretization model, and two plies with different orientations are discretized with only one element through the thickness. This assumption could influence the plate mechanical response, leading to different results in terms of delaminated area. A linear pattern obtained by means of Abaqus/Viewer was used to evaluate the experimental damaged area. The tested specimens were comparable in terms of force, displacement, and kinetic energy as a function of the time, and inter-laminar damages. In particular, the experimental damaged area was about 250 mm 2 , while the numerical damaged area, evaluated by using different numerical models, was in the range 166-609 mm 2 . The continuum shell model uses a lower number of interfaces between the plies, if compared with the solid element discretization model, and two plies with different orientations are discretized with only one element through the thickness. This assumption could influence the plate mechanical response, leading to different results in terms of delaminated area.
In Figures 20 and 21, the numerical results of the final Continuum shell model in terms of inter-laminar and intra-laminar damages are respectively reported. A linear pattern obtained by means of Abaqus/Viewer was used to evaluate the experimental damaged area. The tested specimens were comparable in terms of force, displacement, and kinetic energy as a function of the time, and inter-laminar damages. In particular, the experimental damaged area was about 250 mm 2 , while the numerical damaged area, evaluated by using different numerical models, was in the range 166-609 mm 2 .
The continuum shell model uses a lower number of interfaces between the plies, if compared with the solid element discretization model, and two plies with different orientations are discretized with only one element through the thickness. This assumption could influence the plate mechanical response, leading to different results in terms of delaminated area.
In Figures 20 and 21, the numerical results of the final Continuum shell model in terms of interlaminar and intra-laminar damages are respectively reported.  According to Figure 20, the simulation over predicts the size of the delamination. This overprediction can be related to the shell formulation, which is not able to correctly predict the forces acting in the out-of-plane direction, and to the discretization in the thickness direction. Indeed,  A linear pattern obtained by means of Abaqus/Viewer was used to evaluate the experimental damaged area. The tested specimens were comparable in terms of force, displacement, and kinetic energy as a function of the time, and inter-laminar damages. In particular, the experimental damaged area was about 250 mm 2 , while the numerical damaged area, evaluated by using different numerical models, was in the range 166-609 mm 2 .
The continuum shell model uses a lower number of interfaces between the plies, if compared with the solid element discretization model, and two plies with different orientations are discretized with only one element through the thickness. This assumption could influence the plate mechanical response, leading to different results in terms of delaminated area.
In Figures 20 and 21, the numerical results of the final Continuum shell model in terms of interlaminar and intra-laminar damages are respectively reported.  According to Figure 20, the simulation over predicts the size of the delamination. This overprediction can be related to the shell formulation, which is not able to correctly predict the forces acting in the out-of-plane direction, and to the discretization in the thickness direction. Indeed, According to Figure 20, the simulation over predicts the size of the delamination. This over-prediction can be related to the shell formulation, which is not able to correctly predict the forces acting in the out-of-plane direction, and to the discretization in the thickness direction. Indeed, assuming 12 elements in the thickness direction reduces the number of interfaces where the delamination can occur, leading to larger inter-laminar damages in the remaining interfaces. Figures 22 and 23 show the predicted and measured contact force, impactor displacement, internal energy, and impactor velocity as a function of the time. According to Figure 22a, the numerical curve deviates from the experimental one at about 0.3 ms due to anticipated numerical prediction of the onset of delaminations, resulting in a longer contact duration. This effect can be due to the lower number of interfaces. Moreover, the numerically overestimated impactor displacement (6.77% as reported in Figure 22b) and the underestimated residual internal energy of the numerical model (−28.52% as reported in Figure 23a), confirm the overprediction of damages observed in the simulation. Hence, in the adopted numerical model, significant differences can be observed between the numerical and experimental inter-laminar and intra-laminar damages, as shown in Figures 20 and 21. In particular, the stiffness of the continuum shell element is underestimated if compared with experimental results. Indeed, the maximum displacement of the numerical result is lower than the experimental result in conjunction with a lower peak force.
The numerical results of the solid element model in terms of inter-laminar and intra-laminar damages are reported in Figures 24-27. According to Figure 22a, the numerical curve deviates from the experimental one at about 0.3 ms due to anticipated numerical prediction of the onset of delaminations, resulting in a longer contact duration. This effect can be due to the lower number of interfaces. Moreover, the numerically overestimated impactor displacement (6.77% as reported in Figure 22b) and the underestimated residual internal energy of the numerical model (−28.52% as reported in Figure 23a), confirm the over-prediction of damages observed in the simulation. Hence, in the adopted numerical model, significant differences can be observed between the numerical and experimental inter-laminar and intra-laminar damages, as shown in Figures 20 and 21. In particular, the stiffness of the continuum shell element is underestimated if compared with experimental results. Indeed, the maximum displacement of the numerical result is lower than the experimental result in conjunction with a lower peak force.
The numerical results of the solid element model in terms of inter-laminar and intra-laminar damages are reported in Figures 24-27.        Figure 24 shows a good agreement between the predicted numerical delaminated area (267 mm 2 ) and the experimental one. Moreover, the delaminations through the specimen sections are reported in Figure 27b. A good agreement between numerical and experimental results has been reached in terms of displacement-time, energy-time, and impact velocity-time as highlighted in Figure 28a,b and Figure 29b. The numerical model with solid element approach is able to predict with good  Figure 24 shows a good agreement between the predicted numerical delaminated area (267 mm 2 ) and the experimental one. Moreover, the delaminations through the specimen sections are reported in Figure 27b. A good agreement between numerical and experimental results has been reached in terms of displacement-time, energy-time, and impact velocity-time as highlighted in Figure 28a,b and Figure 29b. The numerical model with solid element approach is able to predict with good approximation the maximum peak force, with a deviation of −8.90% with respect to the experiment, although it is not able to accurately predict the entire load history. Moreover, the proposed numerical model better predicts the maximum displacement (deviation of 0.37%) and reproduces the overall stiffness of the whole plate. According with Figure 29a, a good correlation in terms of internal energy vs time has been reported. In particular, the solid model compared to the shell one reproduces, in a good way, the amount of energy dissipated in the form of friction, viscosity, and especially due to breakages, is better computed by the solid model (−15.41%) respect to the shell one (−28.52%).
Appl. Sci. 2019, 9, x 19 of 22 approximation the maximum peak force, with a deviation of −8.90% with respect to the experiment, although it is not able to accurately predict the entire load history. Moreover, the proposed numerical model better predicts the maximum displacement (deviation of 0.37%) and reproduces the overall stiffness of the whole plate. According with Figure 29a, a good correlation in terms of internal energy vs time has been reported. In particular, the solid model compared to the shell one reproduces, in a good way, the amount of energy dissipated in the form of friction, viscosity, and especially due to breakages, is better computed by the solid model (−15.41%) respect to the shell one (−28.52%).

Conclusions
In this work, the influence of modeling parameters on the accuracy of results in simulating the impact response of a composites plate subjected to low velocity impact was investigated. Indeed, a numerical procedure has been introduced, able to assess the relevance of modeling parameters for the simulation of intra-laminar and inter-laminar damage onset and progression. In the first stage, an analytical procedure for the prediction of the impact event has been introduced. The results have been used to calibrate the numerical models in terms of element types, element size, and element discretization along the thickness. Continuum shell and solid elements formulations have been analyzed. After this first screening based on analytical-numerical comparisons, continuum shell element model (12 elements through the thickness) and a solid element model (24 elements through the thickness, therefore one element for ply) were selected as the most suitable combination to represent the impact behavior of the specimen. The built-in Abaqus intra-laminar damage model has been used for Continuum Shell elements, while a user defined material model has been developed to simulate the intra-laminar damages in Solid elements, by means of Hashin's failure criteria and gradual material properties degradation rules. Cohesive zone model elements were adopted to simulate the inter-laminar damage for both the FEM models. The results obtained highlight some

Conclusions
In this work, the influence of modeling parameters on the accuracy of results in simulating the impact response of a composites plate subjected to low velocity impact was investigated. Indeed, a numerical procedure has been introduced, able to assess the relevance of modeling parameters for the simulation of intra-laminar and inter-laminar damage onset and progression. In the first stage, an analytical procedure for the prediction of the impact event has been introduced. The results have been used to calibrate the numerical models in terms of element types, element size, and element discretization along the thickness. Continuum shell and solid elements formulations have been analyzed. After this first screening based on analytical-numerical comparisons, continuum shell element model (12 elements through the thickness) and a solid element model (24 elements through the thickness, therefore one element for ply) were selected as the most suitable combination to represent the impact behavior of the specimen. The built-in Abaqus intra-laminar damage model has been used for Continuum Shell elements, while a user defined material model has been developed to simulate the intra-laminar damages in Solid elements, by means of Hashin's failure criteria and gradual material properties degradation rules. Cohesive zone model elements were adopted to simulate the inter-laminar damage for both the FEM models. The results obtained highlight some important aspects. The firstly introduced analytical method provides a very fast approach in terms of computational costs to preliminary assess the impact behavior of the composite plate. In the secondly introduced numerical model, which neglects the damages and uses a plain strain formulation, the best ratio between accuracy of the results and computational costs is represented by a reduced number of elements in the thickness direction. Then, Abaqus built-in failure criteria for intra-laminar and inter-laminar damage onset and propagation have been used to preliminary assess the damage behavior of the impacted panel. Finally, the solid element formulation, considering both intra-laminar and inter-laminar damages, has been introduced. In particular, user-defined intra-laminar failure criteria (Hashin's for fiber and tensile matrix, Puck-Shurmann for compressive matrix, and a non-linear behavior of the matrix shear) have been considered, allowing a more accurate prediction of the impact response of the plate. All the presented models provide an increasing degree of accuracy, starting from the analytical to the numerical solid model, being characterized, however, by increasing computational cost.
The results demonstrate the superior capabilities of the solid element model in conjunction with the user material model to simulate the low velocity impact phenomenon on the composite plate, with respect to the shell element model which substantially anticipate the damage onset and overestimates the damages. In conclusion, the solid elements model gives a more realistic representation of the impact behavior of the panel. Indeed, the shell elements formulation does not allow to accurately evaluate the stresses in the normal direction of the plane (which is the predominant direction of deformation during an impact phenomenon) and neglects the effects of the shear.