Effect of Cohesive Properties on Low-Velocity Impact Simulations of Woven Composite Shells

Citation: Ferreira, L.M.; Coelho, C.A.C.P.; Reis, P.N.B. Effect of Cohesive Properties on Low-Velocity Impact Simulations of Woven Composite Shells. Appl. Sci. 2023,13, 6948. https://doi.org/10.3390/ app13126948 Academic Editor: Giangiacomo Minak Received: 24 May 2023 Revised: 6 June 2023 Accepted: 7 June 2023 Published: 8 June 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). applied sciences Article Effect of Cohesive Properties on Low-Velocity Impact Simulations of Woven Composite Shells Luis M. Ferreira 1,2,* , Carlos A. C. P. Coelho 3and Paulo N. B. Reis 4,* 1Grupo de Elasticidad y Resistencia de Materiales, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino Descubrimientos, 41092 Sevilla, Spain 2Escuela Politécnica Superior, Universidad de Sevilla, C/Virgen de África, 7, 41011 Sevilla, Spain 3Unidade Departamental de Engenharias, Escola Superior de Tecnologia de Abrantes, Instituto Politécnico de Tomar, Rua 17 de Agosto de 1808, 2200-370 Abrantes, Portugal; [email protected] 4University of Coimbra, Centre for Mechanical Engineering, Materials and Processes (CEMMPRE), Advanced Production and Intelligent Systems (ARISE), Department of Mechanical Engineering, 3030-194 Coimbra, Portugal *Correspondence: [email protected] (L.M.F.); [email protected] (P.N.B.R.) Abstract: The effect of the interface stiffness and interface strength on the low-velocity impact response of woven-fabric semicylindrical composite shells is studied using finite element (FE) models generated with continuum shell elements and cohesive surfaces. The intralaminar damage is accounted for using the constitutive model provided within the ABAQUS software, while the interlaminar is addressed utilising cohesive surfaces. The results show that the interface stiffness has a negligible effect on the force and energy histories for values between 10 1 N/mm 3 and 2.43 × 10 6 N/mm 3 . However, it has a significant impact on the delamination predictions. It is observed that only the normal interface strength affects the maximum impact force and the delamination predictions. Increasing its value from 15 MPa to 30 MPa resulted in an 8% growth in the maximum force, and a substantial reduction in the delaminated area. The obtained results serve as guidelines for the accurate and efficient computation of delamination. The successful validation of the FE models establishes a solid foundation for further numerical investigations and offers the potential to significantly reduce the time and expenses associated with experimental testing. Keywords: low-velocity impact; cohesive behaviour; interlaminar properties; finite element method; composite shells 1. Introduction Low-velocity impact tests are commonly employed for assessing the damage resistance and tolerance of fibre-reinforced composite structures. The extensive research in this area focuses on experimental analysis, aiming to comprehend the intricate progressive degradation of the material, considering inter-related failure mechanisms [1,2]. The intrinsic characteristics of the structure and of the impactor, and the prevailing environmental conditions collectively influence the structural response to low-velocity impact events. Visual inspection of the damaged area, measurement of dent depth, and a nondestructive evaluation, such as C-Scan, may be used to assess the damage extension in a structure that has experienced low-velocity impact(s) [3]. Typically, matrix cracking is the initial failure observed at the onset of the impact force. Subsequently, delaminations are induced by the matrix cracking and propagate further as the applied force continues to increase. A substantial decrease in force can frequently be used to detect the emergence of notable delaminations [ 4 ]. Depending on the amount of impact energy, there may be noticeable fibre breakage after the delamination has grown. Due to the flexibility of the matrix and the disorganisation of the broken fibres, persistent indentation may appear after the impact [5]. Appl. Sci. 2023,13, 6948. https://doi.org/10.3390/app13126948 https://www.mdpi.com/journal/applsci Appl. Sci. 2023,13, 6948 2 of 19 Understanding the phenomenon of impact damage in laminated composite structures presents a significant analytical challenge due to its intricate and complex nature. The numerical models built using the finite element (FE) method, on the other hand, are a powerful tool for predicting physical processes. Many different impact configurations and structures can be studied using FE models, especially in situations when experimental testing would be prohibitively expensive or complex [ 6 – 9 ]. The computing requirements increase as the model’s complexity increases. Hence, there is a need to dedicate efforts towards developing time-efficient FE models. Simulating low-velocity impact tests on composite structures poses a challenging problem due to the occurrence of contact and progressive degradation of the material. Consequently, it is crucial to accurately model the interlaminar and intralaminar failure mechanisms within an FE analysis. Additionally, defining several important parameters that can impact the precision of numerical predictions, such as material properties, element types, FE mesh discretization, assembly of the parts, and imposed boundary conditions, becomes imperative [10]. There is a scarcity of numerical studies focusing on the impact dynamics of cylindrical composite shells. Among the notable studies, Kim et al. [ 11 ] conducted research that revealed a relationship between the curvature of composite shell and the contact force. They observed an increase in contact force with a higher shell curvature, accompanied by a decrease in deflection and the contact duration. Furthermore, they found that the velocity of the impactor had a more significant influence on the contact force compared to the mass of the impactor. Another noteworthy study by Choi [ 12 ] demonstrated a similar contact force and deflection. Zhao and Cho investigated the onset and evolution of impact damage in composite shells [13], while Kistler and Wass established scaling relationships between impact energy, momentum, mass, and velocity [ 14 ]. In their investigation, Kumar et al. [ 15 ] employed an FE analysis to examine the impact behaviour of cylindrical composite shells. This paper aims to study the effect of the cohesive properties on the numerical predictions of low-velocity impact events on woven-fabric composite shells. This work is based on the FE modelling strategy presented in [ 16 ] and the experimental testing performed in [ 2 ]. The ABAQUS/Explicit [ 17 ] VUMAT subroutine specifically designed for fabric-reinforced composites was used to capture the intralaminar damage [ 18 ]. Simultaneously, a surface-based cohesive model (S-BCM) was implemented to account for the interlaminar damage. The document is structured as follows: Section 2provides a comprehensive overview of the FE modelling approach employed, encompassing the description of intra- and interlaminar damage models, FE mesh discretisation, and boundary conditions and contact definitions. In Section 3, the accuracy of the numerical predictions regarding delamination, force, and energy histories is thoroughly assessed by comparing the results with experimental evidence, while considering the influence of interface stiffness and strength. Finally, Section 4presents the key findings and main conclusions derived from this study. 2. Finite Element Models In this section, the approaches employed for generating the FE models of semicylindrical composite shells are presented. It encompasses the constitutive models utilised to simulate both intralaminar and interlaminar material impact induced damage on the semicylindrical composite shells. Additionally, it provides a concise description of the element types and sizes, the boundary conditions, and the contact definitions imposed to replicate the experimental setup. 2.1. Finite Element Modelling Approach Nine plain weaves of glass fibres at 0 and 90 degrees with 98 ± 4% g/m 2 in each direction were combined with an AROPOL FS 1962 polyester resin and a MEKP-50 hardener to produce the semicylindrical composite shells shown in Figure 1. This figure also shows a schematic representation of the shell layup generated with continuum shell elements. These Appl. Sci. 2023,13, 6948 3 of 19 elements provide an efficient and accurate way to simulate the behaviour of shell structures under various loading conditions, making them suitable to study composite shells structures. Moreover, they are computationally efficient when compared to solid elements. Figure 1. FE modelling approach based on continuum shell elements and cohesive surfaces: ( a ) Representation of the complete FE model. ( b ) Schematic representation of the shell in the through-thickness direction. The parameter trefers to the thickness of a layer. To model the interlaminar behaviour, cohesive surfaces were added to the interfaces of adjacent layers. The cohesive behaviour was applied to the nodes that were initially in contact, meaning that any new contact that occurred during the analysis did not experience cohesive constrain forces. The cohesive surfaces are directly related to the shell elements and share the same FE mesh size. In this way, they provide a simple way to model the interlaminar connection in negligible small interface thicknesses. The use of this strategy for modelling the interlaminar damage implies that the S-BCM must be defined as a surface interaction property and not as a material property. The main advantage of this approach is the quick and easy modelling of the cohesive connection. Moreover, it often allows for higher stable time increments than cohesive elements and does not add any mass to the FE model. The simulations were conducted using a semiautomatic mass-scaling technique that was uniformly applied to the entire FE model. A target time increment of 1 × 10 −7 was employed during the simulations. This approach resulted in an overall mass increase of 1.8%. However, it also provided a notable benefit in terms of computational efficiency, reducing the computational cost of the solution by approximately 50%. 2.2. Intralaminar Damage Modelling To capture the progression of damage at the intralaminar level, a constitutive model specifically designed for fabric-reinforced composites, was used in ABAQUS/Explicit [ 17 ]. This constitutive model, initially developed by Johnson et al. [ 19 ] and based on the work of Ladeveze and Ledantec [ 20 ], offers comprehensive capabilities for capturing the intralaminar damage phenomena. It can be seamlessly integrated into the analysis workflow as a VUMAT user subroutine, with the user-defined material identified by the text “ABQ_PLY_FABRIC” [ 18 ]. This integration enables effective modelling and simulation of the intralaminar damage process in fabric-reinforced composites. In this damage model, each woven fabric-reinforced layer is modelled as a homogeneous orthotropic elastic material with the capability of sustaining progressive stiffness degradation due to fibre failure and matrix cracking, and plastic deformation under shear loading [16,18]. Appl. Sci. 2023,13, 6948 4 of 19 Notice that in the absence of specific experimental data for the tested material, reference values from similar studies/materials were employed as a starting point. However, recognising the need for a more refined and comprehensive analysis, a parametric study was conducted to thoroughly investigate and evaluate the influence of fracture toughness on the numerical-experimental correlation. The intralaminar damage model utilised in this study adopts the maximum stress failure criterion to determine the initiation of damage in the fibres. Additionally, a damage evolution model based on fracture energies was employed to govern the progressive reduction in stiffness. To describe the mechanical behaviour of the degraded orthotropic material, the Hooke’s law was applied, and its expression is as follows: ε=    S1 (1−d1)S12 0 S21 S2 (1−d2)0 0 0 S6 (1−d12)     σ(1) where Sij represents the components of the compliance matrix for the undamaged orthotropic material. The nominal Cauchy stress tensor is denoted as σ , while ε represents the elastic strain tensor. The damage coefficients, d1 , d2 , and d12 , have specific associations: d1 and d2 relate to fibre fracture along directions 1 and 2, respectively, whereas d12 corresponds to matrix microcracking caused by shear deformation. The elastic domain, at any given time, is calculated based on the damage activation functions, Fα, expressed as: Fα=˜ σα Xα −rα≤0 with α=1±, 2±, 12 (2) where ˜ σα and Xα represent the effective stresses and strengths, respectively, and rα denotes the damage thresholds, initially set to 1. Once the damage onset is reached, indicated by ˜ σα Xα= 1, the evolution of the damage coefficients d1 and d2 (associated with tensile and compressive loading) is determined using Equations (3) and (4) . Similarly, the damage coefficient d12 (related to shear loading) is obtained from Equation (5). d1,2 =1−1 r1,2 e−A1,2(r1,2−1)(3) A1,2 =2g1,2 0Lch G1,2 f−g1,2 0Lch with g1,2 0=X2 1,2 2E1,2 (4) d12 =min[α12 ln (r12),dmax 12 ](5) where r1,2 and r12 represent the damage thresholds for tensile/compressive and shear loading, respectively, G1,2 f refers to the fracture energies per unit area, g1,2 0 corresponds to the elastic energy density at the damage onset, and Lch represents the characteristic length of the element. The shear damage parameter α12 is determined through a calibration procedure outlined in [18]. At the intralaminar level, the shear damage response is driven by the nonlinear behaviour of the matrix due to matrix microcracking. This response involves both stiffness loss and plasticity. The VUMAT subroutine employs yield and hardening functions, as shown in Equations (6) and (7), to define the plasticity response of the matrix: Fpl =|˜ σ12| − ˜ σ0(¯ εpl)≤0 (6) ˜ σ0(¯ εpl) = ˜ σy0+C(¯ εpl)p(7) Appl. Sci. 2023,13, 6948 5 of 19 where ˜ σ12 represents the effective shear stress, ˜ σy0 is the initial effective shear stress, and ¯ εpl represents the plastic strain due to shear deformation. The coefficient and power term in the hardening function are denoted by Cand the superscript p, respectively. The two-step homogenisation methodology described by Liu et al. [ 21 ] was used to estimate the stiffness properties of the layers while the remaining properties were taken from the literature [22–27]. The intralaminar properties employed in this study are listed in Table 1. Table 1. Intralaminar properties defined for each layer of the composite shell [16]. Property Symbol Units Value Mass density ρkg/m31900 Stiffness properties E+,− 1=E+,− 2GPa 21.9 E3GPa 8.6 G12 GPa 3.4 G13 GPa 2.4 ν12 - 0.14 Strength properties X+ 1=X+ 2MPa 250 X− 1=X− 2MPa 200 X12 MPa 40 Fracture toughness G1,2 fN/mm 4.5 Shear plasticity dmax 12 - 1 ˜ σy0MPa 25 C- 800 p- 0.552 2.3. Interlaminar Damage Modelling Among the various damage modes observed in composite laminates subjected to low-velocity impact loads, the debonding of adjacent layers stands out as one of the most significant. In this way, a S-BCM was incorporated in the FE models to account for delamination. This model has features that are comparable to cohesive elements and it is particularly designed for extremely thin interfaces. Moreover, the cohesive behaviour is classified as a surface interaction attribute and is controlled by a bilinear traction–separation (τ−δ) constitutive model, which is illustrated in Figure 2for the pure normal and shear responses, and for the expected mixed-mode response. Notice that the subscripts n , s , and t correspond to the normal, first, and second shear directions, respectively, and that the superscripts 0 and f to the onset and fully delaminated stages, respectively. The interface is defined with an initial elastic response with a slope defined by the normal kn and shear ks,t cohesive stiffnesses, and which is maintained until the failure criterion condition is met. The stress-based quadratic failure criterion shown in Equation (1) is used to determine the damage onset of the interface material: hτni τ0 n2 +τs τ0 s2 +τt τ0 t2 =1 (8) where, in the context of the present study, the variables τn , τs and τt symbolise the interface normal and shear tractions, while τ0 n , τ0 s , and τ0 t represent the respective peak contact stress values associated with these tractions. The Macaulay brackets h i indicate that the compressive stress does not contribute to damage. When the quadratic interaction function involving the nominal stress ratios is equal to one, interlaminar damage is presumed to start occurring. Equation (9) represents the softening response of the cohesive surface where dis the damage coefficient. τi=(1−d)kiδiwith i =n,s,t(9) Appl. Sci. 2023,13, 6948 6 of 19 Figure 2. Bilinear cohesive model for mixed modes (normal and shear modes). From this point on, the energy Gc that was dissipated during the damage process and that corresponds to the region below the τ−δ curves controls the damage progression. Its value is obtained implementing the Benzeggagh and Kenane (B-K) damage evolution criterion under mixed-mode loading [28]: Gc=GIc +(GIIc +GIc)GII +GIII GI+GII +GIII η (10) where the parameters GI , GII , and GIII signify the strain energy release rates associated with mode I, II, and III, respectively. On the other hand, GIc and GIIc denote the critical strain energy release rates, and η represents a material coefficient of the B-K equation. Notice that with this constitutive model, damage evolution depicts the deterioration of the cohesive stiffness. In contrast, damage evolution for cohesive elements describes the deterioration of the material stiffness. The values used to study the effect of the interlaminar properties listed in Table 2were taken from the authors previous work [ 16 ]. The fracture toughness parameters remained unchanged for all simulations, and the effect of the variation in the strength property values was studied assuming kn=ks=kt=106[29]. Table 2. Interlaminar properties defined for the composite shell interfaces. Property Symbol Units Value Stiffness properties kn=ks=ktN/mm310, 103, 105, 106 Strength properties τ0 n(τ0 s=τ0 t)MPa 15 (30, 45, 60) 30 (30, 60) Fracture toughness GIc N/mm 0.3 GIIc =GIIIc N/mm 0.6 η- 1.45 2.4. Finite Element Mesh Eight-node continuum shell elements (SC8R) with a characteristiclength of le=0.3 mm , reduced integration, and stiffness hourglass formulation were used to mesh the semicylindrical composite shell layers [ 17 ]. The equations proposed by Turon et al. [ 30 ] were used to Appl. Sci. 2023,13, 6948 7 of 19 calculate the element‘s characteristic length le , and its effect on the numerical predictions was studied by the authors in [ 16 ]. The elements’ coordinate system was set up considering the material directions in the different composite layers. All the nine layers were modelled with the same thickness t= 0.178 mm, resulting in a total thickness of the composite shell of 1.6 mm. As mentioned before, the cohesive behaviour between the adjacent layers of the laminate was modelled using cohesive surfaces, therefore no elements were required. On the other hand, the impactor was meshed with discrete rigid elements (R3D4) with a characteristic length of le= 0.3 mm. Finally, the supports were modelled as analytical rigid bodies, therefore no element or mesh definitions were necessary. During the simulation process in ABAQUS/Explicit [ 17 ], the analysis automatically halts if severe distortion is detected in deteriorated elements. It is a common practice to remove degraded elements from the FE mesh once a damage variable, often associated with the fibre, reaches a predefined threshold close to one. However, it is crucial to define these thresholds carefully as they can significantly impact the shape of the constitutive law. In the case of the low-velocity impact simulations in this study, no puncture was observed in the experimentally tested specimens [ 2 ]. Therefore, element deletion was not implemented, ensuring that the simulations captured the behaviour without compromising the structural integrity of the model. 2.5. Boundary Conditions and Contact Definitions Identically to the previous numerical study presented by the authors [ 16 ], the generated FE models replicate the experimental testing conditions used in [ 1 , 2 ]. In this lowvelocity impact test, the composite shell specimens were stroked by an impactor with a 10 mm diameter hemispherical head and a mass of 2.826 kg. The energy of 5J was used to promote visible damage without full perforation of the specimens (further details about the experimental testing setup can be found in [ 2 ]). The amplitude of the impactor’s velocity throughout the step time was also defined according to the experimental evidence. All the rotations of the impactor were constrained, and only the displacements along the y-direction were allowed. The fixed rigid body supports represented in Figure 1a simply support the bottom straight faces of the composite shell. Notice that, as presented in detail in [ 16 ], only one-fourth of the semicylindrical composite shell was generated to reduce the computational cost of the simulations. Despite this, due to the significant scale of the 9-layer FE model, i.e., consisting of 403,753 elements and 814,462 nodes, the computational complexity necessitated the utilisation of a computer cluster for the numerical simulations. To expedite the computations, 20 CPUs and 40 GB of RAM were allocated. Despite these computational resources, the solution times ranged from 7 days to 20 days. The surface-to-surface contact interactions between the supports, impactor, and the composite shell were carefully defined using the penalty enforcement contact method available in ABAQUS/Explicit [ 17 ]. To accurately capture the behaviour of completely delaminated interfaces, this contact formulation was implemented across all interfaces of the composite shell, ensuring frictional effects were properly accounted for. Specific friction coefficient values were assigned to the tangential behaviour between different surfaces. For the contact between the metal hemispherical head of the impactor and the top surface of the composite shell, a friction coefficient of µ= 0.3 was specified [ 31 , 32 ]. Additionally, a value of µ= 0.7 was assigned to the contact between the composite laminate surfaces and the metal surfaces of the supports. For all other interfaces within the composite shell, a friction coefficient of µ=0.5 was considered. 3. Results Correctly defining the interface stiffness and interface strength is of utmost importance in ensuring the intended functionality of FE models that incorporate cohesive behaviour. These interlaminar properties can significantly influence both the accuracy, in terms of the correlation with experimental evidence, and the computational cost of the solutions. Appl. Sci. 2023,13, 6948 8 of 19 Therefore, the subsequent subsections focus on analysing the impact of these properties on low-velocity impact simulations. 3.1. Effect of the Interface Stiffness To assess the influence of the interface stiffness on the impact response of the semicylindrical woven composite shells, simulations were carried out using cohesive stiffness values ranging from 10 1 N/mm 3 to 2.43 × 10 6 N/mm 3 ( kn=ks=kt ), while maintaining the normal and shear interface strengths constant ( τ0 n= 15 MPa and τ0 s=τ0 t=30 MPa ). Notice that the highest value considered for the cohesive stiffness, that is, 2.43 × 10 6 N/mm 3 , was obtained by employing the equation proposed by Turon et al. [30]: k=αE3 t(11) where t is the thickness of each ply, E3 is the through-thickness elastic modulus, and α is a parameter larger than 1, which, in this study, was set with a value of 50, as suggested in [ 30 ]. The cohesive stiffness proposed by Camanho et al. [ 29 ] (10 6 N/mm 3 ) was also analysed in this study. The numerical and experimental [ 2 ] force–time, energy–time, and force–displacement curves are shown in Figures 3–5, respectively. It is possible to appreciate that the cohesive stiffness has no impact on the force and energy histories for values between 10 3 N/mm 3 and 2.43 × 10 6 N/mm 3 , and that there is a satisfactory numerical–experimental correlation within this range. The force–time and force–displacement response curves are significantly affected by employing smaller cohesive stiffness values, such as 10 1 N/mm 3 , because the stiff connection between the laminate layers is not ensured. It was noted, however, that the energy–time response curve was unaffected because the numerical output results for the impact energy were calculated based solely on the mass and velocity of the impactor along the established step time, which were defined on the FE models in accordance with the experimental evidence. On the other hand, high cohesive stiffness values caused convergence issues. For instance, adopting a cohesive stiffness of 10 9 N/mm 3 resulted in immediate convergence failure. It should be noted that these findings are in good agreement with those presented by Song et al. in [33]. Figure 3. Effect of the interface stiffness on the force–time results with τ0 n= 15 MPa and τ0 s=τ0 t=30 MPa. Appl. Sci. 2023,13, 6948 9 of 19 Figure 4. Effect of the interface stiffness on the energy–time results with τ0 n= 15 MPa and τ0 s=τ0 t=30 MPa. Figure 5. Effect of the interface stiffness on the force–displacement results with τ0 n= 15 MPa and τ0 s=τ0 t=30 MPa. If the cohesive stiffness is maintained between 10 3 N/mm 3 and 2.43 × 10 6 N/mm 3 , the impact bending stiffness (IBS), which is determined by the slope of the ascending branch of the force–displacement curves represented in Figure 5[ 34 ], is also hardly affected. Once more, for 10 1 N/mm 3 , the IBS significantly declines and completely breaks apart from the experimental evidence. It was observed that, in terms of computing efficiency, an increase in cohesive stiffness is accompanied by a significant rise in the number of increments required to compute the solution. As an example, it was observed that increasing the cohesive stiffness from 105N/mm3to 106N/mm3raised the computational cost of the solution by 138%. The cohesive stiffness was observed to have a limited impact on the response curves when maintained within the range of 10 3 –2.43 × 10 6 N/mm 3 . However, it was identified as a critical factor affecting the delamination predictions. To understand and quantify this impact, the effect of the cohesive stiffness on the delamination initiation and progression Appl. Sci. 2023,13, 6948 16 of 19 Figure 16. Effect of the interface stiffness with τ0 n= 30 MPa and kn=ks=kt= 10 6 N/mm 3 : (a) CSQUADSCRT output variable. (b) CSDMG output variable. Figure 17. Effect of interface strength on the delamination predictions with kn=ks=kt= 10 6 N/mm 3 , represented throughout the FE model. 4. Conclusions In this study, the influence of the interface stiffness and strength on the low-velocity impact response of woven fabric-reinforced composite shells was investigated. For that purpose, FE models that incorporated continuum shell elements and cohesive surfaces to account for the intra- and interlaminar damage, respectively, were generated in ABAQUS/Explicit. Appl. Sci. 2023,13, 6948 17 of 19 The simulations were performed using interface stiffness values ranging from 10 1 N/mm 3 to 2.43 × 10 6 N/mm 3 and interface strength values ranging from 15 MPa to 60 MPa. To be able to analyse the results, the force–time, energy–time, and force–displacement curves were presented. Additionally, in order to assess the extent of cohesive degradation and delamination within the FE models, the variation in the percentage of nodes initiating cohesive degradation and the percentage of nodes that become fully delaminated were plotted. It was observed that, for values between 10 3 N/mm 3 and 2.43 × 10 6 N/mm 3 , the interface stiffness has a negligible effect on the force and energy histories, and, within this range, there is a satisfactory numerical–experimental correlation. The use of stiffness values outside this range resulted in convergence issues (high values of interface stiffness) or no numerical–experimental correlation (low values of interface stiffness). Although the interface stiffness has a minimal impact on the response curves’ accuracy, it has a significant impact on the delamination predictions. Only for an interface stiffness value of 10 5 N/mm 3 do completely delaminated nodes become clearly visible. Furthermore, it is worth noting that a significant increase in computational cost was observed with the increase in interface stiffness. The results show that the shear interface strength values (peak contact stresses) have a negligible impact on the force and energy histories. On the other hand, the normal interface strength affects the numerical predictions of the maximum impact force. Increasing its value from 15 MPa to 30 MPa resulted in an 8% increase in the maximum impact force, and thus a deviation from the experimental data. Only the variation in the normal interface strength was found to have an effect on the delamination predictions. It was observed that the percentage of fully delaminated nodes decreases with the increase of the normal interface strength; however, the percentage of nodes that initiate cohesive degradation remains almost unchanged. In addition, it was observed that varying the interface strength had minimal impact on the computational cost of the solutions. The results obtained from this study provide valuable insights and recommendations for accurately and efficiently simulating delamination phenomena. By successfully validating the FE models used in this research, a strong basis is established for conducting future numerical investigations. It also presents an opportunity to reduce the reliance on costly and time-consuming experimental testing. Furthermore, the utilisation of numerical simulations enables researchers to investigate delamination in a more detailed and comprehensive manner. The FE models can capture intricate delamination mechanisms, predict failure modes, and provide detailed information about stress distribution and energy absorption during the process. This information can be utilised to enhance design guidelines, develop mitigation strategies, and optimise the performance and reliability of laminated shell structures. Author Contributions: Conceptualisation, L.M.F. and P.N.B.R.; methodology, L.M.F., P.N.B.R. and C.A.C.P.C.; software, L.M.F.; validation, L.M.F., P.N.B.R. and C.A.C.P.C.; formal analysis, L.M.F. and P.N.B.R.; investigation, L.M.F. and P.N.B.R. and C.A.C.P.C.; writing—original draft preparation, L.M.F. and P.N.B.R.; writing—review and editing, L.M.F. and P.N.B.R. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The data presented in this study are available on request from the corresponding author. Acknowledgments: This research was sponsored by national funds through FCT-Fundação para a Ciência e a Tecnologia, under the project UIB/00285/2020 and LA/P/0112/2020. Conflicts of Interest: The authors declare no conflict of interest. Appl. Sci. 2023,13, 6948 18 of 19 References 1. Ferreira, L.; Coelho, C.; Reis, P. Impact Response of Semi-Cylindrical Composite Laminate Shells Under Repeated Low-Velocity Impacts. In Proceedings of the 2022 Advances in Science and Engineering Technology International Conferences (ASET), Dubai, United Arab Emirates, 21–24 February 2022; pp. 1–5. [CrossRef] 2. Reis, P.; Sousa, P.; Ferreira, L.; Coelho, C. Multi-impact response of semicylindrical composite laminated shells with different thicknesses. Compos. Struct. 2023,310, 116771. [CrossRef] 3. Amaro, A.M.; Reis, P.N.B.; de Moura, M.F.S.F.; Santos, J.B. Damage detection on laminated composite materials using several NDT techniques. Insight 2012,54, 14–20. [CrossRef] 4. ASTM D7136/D7136M-05; Standard Test Method for Measuring the Damage Resistance of a Fiber Reinforced Polymer Matrix Composite to a Drop Weight Impact Event. American Society for Testing and Materials. ASTM: Philadelphia, PA, USA, 2005. 5. Richardson, M.; Wisheart, M. Review of low-velocity impact properties of composite materials. Compos. Part A Appl. Sci. Manuf. 1996,27, 1123–1131. [CrossRef] 6. Ferreira, L.; Coelho, C. Modelling Progressive Damage in NCF Composites using the Continuum Damage Mechanics Method. In Proceedings of the 2022 Advances in Science and Engineering Technology International Conferences (ASET), Dubai, United Arab Emirates, 21–24 February 2022; pp. 1–4. [CrossRef] 7. Ferreira, L.M.; Graciani, E.; París, F. Modelling the waviness of the fibres in non-crimp fabric composites using 3D finite element models with straight tows. Compos. Struct. 2014,107, 79–87. [CrossRef] 8. Ferreira, L.M.; Graciani, E.; París, F. Predicting failure load of a non-crimp fabric composite by means of a 3D finite element model including progressive damage. Compos. Struct. 2019,225, 111115. [CrossRef] 9. Ferreira, L.M.; Graciani, E.; París, F. Three dimensional finite element study of the behaviour and failure mechanism of non-crimp fabric composites under in-plane compression. Compos. Struct. 2016,149, 106–113. [CrossRef] 10. De Moura, M.; Campilho, R.; Amaro, A.; Reis, P. Interlaminar and intralaminar fracture characterization of composites under mode I loading. Compos. Struct. 2010,92, 144–149. [CrossRef] 11. Kim, S.J.; Goo, N.S.; Kim, T.W. The effect of curvature on the dynamic response and impact-induced damage in composite laminates. Compos. Sci. Technol. 1997,57, 763–773. [CrossRef] 12. Choi, I.H. Finite element analysis of low-velocity impact response of convex and concave composite laminated shells. Compos. Struct. 2018,186, 210–220. [CrossRef] 13. Zhao, G.; Cho, C. On Impact Damage of Composite Shells by a Low-Velocity Projectile. J. Compos. Mater. 2004 ,38, 1231–1254. [CrossRef] 14. Kistler, L.S.; Waas, A.M. Experiment and Analysison the Response of Curved Laminated Composite Panels Subjected to Low Velocity Impact. Int. J. Impact Eng. 1998,21, 711–736. [CrossRef] 15. Kumar, S.; Nageswara Rao, B.; Pradhan, B. Effect of Impactor Parameters and Laminate Characteristics on Impact Response and Damage in Curved Composite Laminates. J. Reinf. Plast. Compos. 2007,26, 1273–1290. [CrossRef] 16. Ferreira, L.M.; Coelho, C.A.C.P.; Reis, P.N.B. Numerical Simulations of the Low-Velocity Impact Response of Semicylindrical Woven Composite Shells. Materials 2023,16, 3442. [CrossRef] [PubMed] 17. ABAQUS Analysis User’s Manual, Version 6.14; Dassault Systèmes Simulia Corp: Providence, RI, USA, 2014 . 18. ABAQUS/Explicit. VUMAT for Fabric Reinforced Composites; Dassault Systèmes Simulia Corp, Providence, RI, USA, 2008. 19. Johnson, A.; Pickett, A.; Rozycki, P. Computational methods for predicting impact damage in composite structures. Compos. Sci. Technol. 2001,61, 2183–2192. [CrossRef] 20. Ladeveze, P.; LeDantec, E. Damage modelling of the elementary ply for laminated composites. Compos. Sci. Technol. 1992 , 43, 257–267. [CrossRef] 21. Liu, X.; Rouf, K.; Peng, B.; Yu, W. Two-step homogenization of textile composites using mechanics of structure genome. Compos. Struct. 2017,171, 252–262. . [CrossRef] 22. Sridharan, S.; Pankow, M. Performance evaluation of two progressive damage models for composite laminates under various speed impact loading. Int. J. Impact Eng. 2020,143, 103615. [CrossRef] 23. Bodepati, V.; Mogulanna, K.; Rao, G.S.; Vemuri, M. Numerical Simulation and Experimental Validation of E-Glass/epoxy Composite Material under Ballistic Impact of 9 mm Soft Projectile. Plast. Impact Mech. 2017,173, 740–746. [CrossRef] 24. Esnaola, A.; Elguezabal, B.; Aurrekoetxea, J.; Gallego, I.; Ulacia, I. Optimization of the semi-hexagonal geometry of a composite crush structure by finite element analysis. Compos. Part B Eng. 2016,93, 56–66. [CrossRef] 25. Kostopoulos, V.; Markopoulos, Y.; Giannopoulos, G.; Vlachos, D. Finite element analysis of impact damage response of composite motorcycle safety helmets. Compos. Part B Eng. 2002,33, 99–107. [CrossRef] 26. Hoo Fatt, M.S.; Lin, C. Perforation of clamped, woven E-glass/polyester panels. Compos. Part B Eng. 2004 ,35, 359–378. [CrossRef] 27. Alonso, L.; Martínez-Hergueta, F.; Garcia-Gonzalez, D.; Navarro, C.; García-Castillo, S.; Teixeira-Dias, F. A finite element approach to model high-velocity impact on thin woven GFRP plates. Int. J. Impact Eng. 2020,142, 103593. [CrossRef] 28. Benzeggagh, M.; Kenane, M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos. Sci. Technol. 1996,56, 439–449. [CrossRef] 29. Camanho, P.P.; Davila, C.G.; de Moura, M.F. Numerical Simulation of Mixed-Mode Progressive Delamination in Composite Materials. J. Compos. Mater. 2003,37, 1415–1438. [CrossRef] Appl. Sci. 2023,13, 6948 19 of 19 30. Turon, A.; Dávila, C.; Camanho, P.; Costa, J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng. Fract. Mech. 2007,74, 1665–1682. [CrossRef] 31. Schön, J. Coefficient of friction of composite delamination surfaces. Wear 2000,237, 77–89. [CrossRef] 32. Bresciani, L.; Manes, A.; Ruggiero, A.; Iannitti, G.; Giglio, M. Experimental tests and numerical modelling of ballistic impacts against Kevlar 29 plain-woven fabrics with an epoxy matrix: Macro-homogeneous and Meso-heterogeneous approaches. Compos. Part B Eng. 2016,88, 114–130. [CrossRef] 33. Song, K.; Dávila, C.G.; Rose, C.A. Guidelines and parameter selection for the simulation of progressive delamination. In Proceedings of the 2008 Abaqus Users’ Conference, Newport, RI, USA, 19–22 May 2008; Volume 41, pp. 43–44. 34. David-West, O.; Nash, D.; Banks, W. An experimental study of damage accumulation in balanced CFRP laminates due to repeated impact. Compos. Struct. 2008,83, 247–258. [CrossRef] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.