:

The present paper presents an investigation and analysis study of the effect of the head shapes of the impactor on the damages observed during low-velocity impact on T700/Epoxy composite laminate. Two types of impactors were investigated: hemispherical and flat-face. A new criterion based on the LARC05 damage model was linked as a LARC_VUMAT subroutine to simulate the impact and explore the effects of the head form shape through a three-dimensional finite element model. To properly analyze the problem, the history time of the mechanical responses, such as impact forces, bending, principal, and residual stresses, are highlighted and assessed. Additionally, a comparison with the experimental data found in the literature was performed to check the validity and accuracy of the considered finite element model. The damage occurring in the T700/Epoxy plates is illustrated for each impactor head shape. The mechanical response curves and all kinds of damage of the presented simulations are in perfect agreement with the experiments. The proposed VUMAT is efficient in the prediction of fiber kinking, matrix cracking, fiber splitting, and fiber tension of a laminate, and more importantly, it is easy to implement for other types of materials and the reproducibility of the analysis is assured.


INTRODUCTION
Composite structures are increasingly used in the manufacturing of aircraft and automobile structural components, mainly due to the balanced ratio of high strength to low weight [1]. However, they remain vulnerable to defects induced by the non-miscible and inhomogeneous composition of the material. Moreover, the thin thickness of composite structures makes them susceptible to damage, especially under impact loading [2]. For such a loading case, the damage may barely be visible, but it might nevertheless considerably reduce strength. Accordingly, a better understanding of the behavior under impact loading is required for composite structure [3]. To meet this need, considerable efforts have been made to study and optimize the impact performance of composite laminates [4,5]. Many researchers use experimentation to give a clear explanation of the impact behavior of composite laminates, although relatively few comparative studies have been reported on the impactor head shape under low-velocity impact conditions. Impacts commonly occur accidentally during manufacturing, maintenance operations, etc. Up to now, research has fallen beneath three categories: impact response of composite laminates, damage characterization, and contact mechanisms in impact-loaded composite structures [6].
Authors in [7,8] investigated and experimentally emphasized the impact of impactor head-shape on the impact response of thin woven carbon/epoxy laminates. They report that different impactor shapes have a large effect on the damage mechanism. Authors in [9] reported that different impactor head-shapes produce different degrees of different types of damage, such as fiber breakage, matrix cracking, and delamination, which can affect the residual properties of the composite. Authors in [10] were interested in the effect of impactor diameter. They showed that this parameter significantly affects the impact and post-impact compression response of composite materials. Authors in [11], from impact tests on woven glass and carbon fiber composites at low impact velocities using a drop tower having a hemispherical impactor, found that impact damage was primarily a function of impact energy rather than velocity or mass separately. Furthermore, regardless of the mass of the impactor used, the measured peak forces reveal a correlation with the impact energy, indicating that the low-velocity impact of these samples is a quasi-static process [12]. However, according to the results of their tests, the relationship between peak load and the size of the damage was independent of the diameter of the impactor. In most of the mentioned studies, impact tests or simulations were performed using hemispherical impactors according to ASTM D7136 [13]. Since different shaped objects might collide with the composite structures, the head shape of impactors must be considered in the experiments [14].
Recently, it has been pointed out that the contact area between the impactor and the composite test sample is a crucial factor in the response of the impacted composites. In addition, impactors with a larger contact area yielded a higher initial peak load, highest maximum load, shortest contact time for the impact event timescale, and more extensive delamination [15,16]. Despite this, the overall relationship between internal and external damage to impacted composite laminates is difficult to describe in a simple way [17,18], especially at various impact conditions. Developing of analytical and numerical models is another alternative for investigating the impact response of composite laminates. These techniques provide a better understanding of the physics of impact effects and are more efficient computationally since they are comparatively processable with simple geometries and boundary constraints [19]. For example, to simulate the behavior of composite laminates in high-mass, low-velocity impacts, a simplified analytical model for predicting the onset and propagation of impact damage in composites was proposed in [20,21]. Authors in [22] distinguished the relation between analytical and numerical models for finite element simulation. As a result, they report that analytical models are mainly used to model the overall response to an impact. In [23], two analytical design methods are distinguished. One based on energy equilibrium, where the behavior remains quasi-static at low energy and low velocity [24] and another that uses the notion of semi-infinite space which only takes into account the plate deflection. The effectiveness of this model is established when the strain front does not reach the plate boundaries. Otherwise, the model is no longer suitable, as the dynamics of the plate must be considered [25].
The literature information on the effects of impactor head shapes upon the impact response of laminated composite plates is rather poor. The effect of the impactor head shape is an essential parameter to be studied in order to determine the response of laminates under different impactor shapes as well as the induced damage. So, the main objective of the current study is to numerically determine the influence of impactor head shapes on the low-velocity impact performance of composite laminate plates and compare them to the results of an experimental study found in the literature [6]. For this purpose, a three-dimensional non-linear finite element model was developed to simulate the effect of the impactor head shapes on the impact behavior of a composite laminate plate associated to a user subroutine based on LARC05 criterion and was implemented under the commercial software ABAQUS/Explicit. In terms of impact forces-time, deflectionstime, predicted stresses, deformation plots, and failure damages ascertainment in the composite plate, the finite element analysis findings agree well with the experimental results of [6]. In the predicted stress results, the stresses in the three main directions (S 11 , S 22 , and S 33 ) are illustrated as a replacement for the commonly used von Mises stresses, to show the stress state in each layer according to the fiber direction as well as in the strain illustration case.

II. NUMERICAL MODELLING TECHNIQUES
In order to maintain the same experimental protocol conducted by [1], all the simulated composite specimens possessed the same lay-up of +45 /0 /−45 /0 . The simulated plates are made of T700/Epoxy pre-pregs. The dimensions of the composite specimens were 150×100×3mm 3 according to the ASTM D7136 standard [2]. Two impactor head shapes were taken into consideration: An impactor with a Hemispherical head end (HemI) and a Flat-Face end Impactor (FFI). The two impactors were made of steel. All mechanical data concerning the plate and the impactors are shown in Table  I. The modeling techniques for the composite failure and the impactor head shapes were investigated and ABAQUS/Explicit was used for the impact simulation. The composite plate was modeled using 337296 C3D8I elements and for the impactor R3D4 elements were used as follows: 3953 elements for the HemI and 3638 for the FFI. An investigation of mesh convergence was conducted and the optimal size of 1mm was considered. All details are illustrated in Figure 1(a).
The boundary conditions for the constrains applied during the simulations are similar to those used in the experiments of [6]. More details are presented in Figure 1(b) and Table II    The z-direction of the supports of the plate is pinned. The z-direction of the impactors is unconstrained. Ensure impactor-to-plate contact Interaction conditions were attributed using standard contact algorithm

Loading
A predefined field in terms of velocity was applied into the impact: In addition, for the HemI and FFI, a point inertia mass of 5.20kg and 5.25Kg was assigned.

Simulation time
Time of impact is calculated by:

Contact time in experiments number of increments
Furthermore, considering the damages, a VUMAT user subroutine modeling based on the LARC05 criterion [26] was implemented for the explicit scheme. This failure criterion allows the distinction of 4 different damage initiation mechanisms: matrix cracking, fiber kinking, fiber splitting, and fiber tension, as follows: ) 5" cos . + ) 5 sin .
The angle α represents the direction of the critical plane perpendicular to the local 23-plane. The α value is numerically found such that it maximizes 2) Fiber Tension 1) 55 7 03

3) Fiber Splitting
Terms " , 5" , and ) "" are the traction components of the fracture plane. where I Y is the direction of thickness under transverse compressive load, which is the experimentally measurable material property and its value is usually 53°.
Through the VUMAT, the signs of σ 11 for all the above fiber failure mechanisms must be firstly checked. If σ 11 is positive, then the fiber tension criterion 8 9:;< is calculated. Else, for both fiber splitting and kinking, the stresses are rotated to the kinking plane of the fibers, and subsequently these resultant stresses are rotated to the misalignment plane of fibers. Through this process, the angles ψ and φ are varied together until a peak value of 8 <ABC9 or 8 C; is reached. More details are available in Appendix C.

A. Impact Forces and Bending
The typical profile of the forces and bending curves of the simulated composite samples obtained from the two shapes were compared to the experiments of [1] and are shown in Figure. 2. For the same energy level of impact of 15J, it was found that:  On average, a peak impact force of 6.92kN was recorded for the flat head shape impactor upon contacting the plate. Compared to the HemI, where the average maximum impact load was 4.86kN, this is roughly 70% higher, see Figure 2(a).
 Regarding contact time, the time duration of the impact event changes for each head shape of impactor. Accordingly, the FFI struck the plate with an average contact time of 6.4ms, a reduction of 0.8ms (or approximately 12.50%) compared to the HemI impactor, where the average contact time of 7.2 was registered.
 Observing the records of bending shown on Figure 2(b), the highest values are recorded at the mid-contact time and are 5.26 and 4.48mm for the HemI and FFI respectively. The reduction in bending between the two shapes results from the load distribution over the contact area. For the HemI, the motion process of the impacts is associated with the formation of a spherical cavity in the medium. However, for the FFI this process is eliminated, and a large contact area is put at stake.

B. Main Stresses
Compared to experiments that require stress gauges to measure stresses, the numerical model presented in this study provides an easy way to pick up stresses. Because the presented composite plate is balanced, only the stress levels of the upper layers of the composite are presented. Figure 3 shows the principal stress distribution in the three plate directions for each layer. The main observations are:  A detailed overview of the stress status of each layer is available within the model.
 The stresses have a quasi-symmetric distribution between the upper and lower layers of the composite.
 All the impact resistance is absorbed by the two upper layers. The same findings will be observed but conversely in layers [ 55 and [ 5" .
 As the impacts reaches their maximum, both inner layers L 11 and L 12 undergo traction, while the others are subjected to compression in both X and Y-directions.         The peak value of stresses for ) 55 and ) "" is registered through the first layer (L 1 /+45°) in contact with impactors. However, for ) the peak value exists in the inner-layer of the composite plate (L 12 /+0°).

www.etasr.com Rabouh et al.: Bridging the Effect of the Impactor Head Shape to the Induced Damage during Impact at …
 On average, the peak values are respectively for ) 55 , ) "" and ) equal to 750MPa, 62.35MPa, and 4.04MPa obtained from the simulated impact of the FFI impactor.

www.etasr.com Rabouh et al.: Bridging the Effect of the Impactor Head Shape to the Induced Damage during Impact at …
For the HemI, the recorded values are 1147MPa, 172MPa, and 96.5MPa.
 The values of ) 55 , ) "" are doubling when the contact surface changes from flat to curved, which is confirmed by the high increase in the value of ) in the case of the HemI simulation.
 In each stress versus time curve, the maximum stresses are smaller in the order of FFI and HemI, and the contact time is longer in the same order.
 An interaction between the main stresses is shown in Figure  4.

C. Residual Stresses
After the separation of the impactor from the plate, residual stresses 1) :< 3 remain. In Figure 5, the residual stress contours are presented according to the FE simulation on the T700/Epoxy +45 /0 /−45 /0 composite plate after the impact using both proposed impactor head shapes, through which it can be clearly observed that:  The impactor with the larger contact area provides very high stress waves across the whole area. Nevertheless, similar stress wave patterns are concentrated in the area around the contact area and the supports for both HemI and FFI.
 The highest values of ) :< recorded for the FFI and HemI are respectively 59.76MPa and 40.96MPa along the largest dimension of the plate (X-direction, see Figure 5(a)). This reflects the fact that the bending stresses are more dominant and important at this direction. On the other hand, through the depth of the plate, the values of ) :< decreased as a result of the indentation mechanism, as shown on Figure 7. In fact, the occurrence of this phenomenon can be explained by the storage of a portion of the impact energy by strain mechanisms in the matrix (indentation) and fiber breakage dissipation.

D. Expected Strain
In Figure 6, the results for the strain distribution contours across thickness ε 33 are illustrated. From the simulations at the peak displacement along the focused mid-zone we notice:  For both HemI and FFI, it is found that there is a good correlation between the strain distribution maps ε 33 and the extent and location of intralaminar matrix damage in the specimens, as shown in Figures 5-7.
 In the case of the HemI, the strain, which was initially localized and located at the center region under the impactor, was extended to the peripheral regions during the interaction between the impactor and the plate, while the FFI has initiated the damage in the composite plate at the region in contact with the ring edge of this impactor, as shown in Figure 5(c).
 In addition, it is noticeable that the curvature is less in the central region under the FFI compared to the HemI.

E. Damage Analysis
For the sake of brevity, only the damages caused by the HemI are presented in this section. However, in Appendix B, the reader can find detailed information on the FFI damages. All observations about this impactor will be highlighted during the damage results discussion.
In Figure 7, the typical damages generated using the implemented VUMAT_LARC subroutine of the LARC05 damage criterion in Abaqus/Explicit for both the energy level and the impactor head shape under consideration are illustrated. Based on the outcomes, it is found that:  The model displays the most important mechanisms of damage that can be observed when a composite plate is impacted. Moreover, compared to the existing damage criteria such as those of [4,27,28], the implemented VUMAT_LARC gives a large insight of the kinking and the slitting of fibers in the zone under impactors (see Figures 7(a),(d)) as well as the extend of the matrix cracking (see Figure 7(b)), which cannot be identified using the previously mentioned criteria.

www.etasr.com Rabouh et al.: Bridging the Effect of the Impactor Head Shape to the Induced Damage during Impact at …
Case 1: HemI Case 2: FFI (a) Residual X-stresses after impact.
(b) Residual X-stresses after impact.
(c) Residual Z-stresses after impact.  During the analysis of the impact results, damages are observed on both the front and back surfaces. This observation is highly observed for fiber traction in Figure  7(b), which supports the experimental findings.
 The visible damages on the frontal face impacted of the plate are only impact dents of different sizes. The damage differences in the frontal surface mainly occur due to the difference in contact sizing under different impact conditions (impactor head shapes, interactions, BCs, etc.).

www.etasr.com Rabouh et al.: Bridging the Effect of the Impactor Head Shape to the Induced Damage during Impact at …
 Indeed, in the case of the HemI, damages like fiber kinking and matrix cracking, are initially localized at the center region under the impactor. Then they are extended to the peripheral regions during the interaction between the impactor and the plate. Moreover, it is observed that the local damages are growing trough the plate depth, as shown in Figures 7(a),(c). On the other hand, the FFI damages are initiated at the contact region evolving in a ring shape on the edge of this impactor, as shown in Figure 8. Accordingly, it is found that the FFI damages are less pronounced than the HemI damages.
 It is also noticed that the weakest damage mechanism is the fiber splitting, compared to the other damage mechanisms which are, in decreasing order, matrix cracking, fiber kinking, and fiber tension.
As a final note, some key results from the proposed modeling predictions are summarized in Table III, i.e. maximum impact forces Z A , maximum bending, duration of contact with the impactor, and the CAREAs. The findings clearly demonstrate the capabilities of the model and the power of implementation via VUMAT_LAR to accurately and efficiently capture the impact responses and the damaged area among the composite specimens.

IV. CONCLUSION
In an effort to better understand the effects of the impactor head geometry upon low velocity impact, a three-dimensional finite element model was developed and linked to a VUMAT subroutine to simulate the impact response of the T700/Epoxy composite laminates when struck by the two suggested impactors. Based on the results of these studies, the following conclusions can be made [29]:  The integration of the LARC damage model for material response has shown its effectiveness in determining very accurately the predicted key parameters of the impact event, such as (a) peak impact forces and bending, (b) duration of the impact event, (c) the stresses state, (d) the different mechanisms of damage, and (e) contact area shape (CAREA).
 The composite plate impacted by the Flat-Face Ιmpactor (FFI) showed a higher impact force peak during impact, whilst suffering a lower out-of-plane bending than the Hemispherical Impactor (HemI).
 The HemI's contact time is considerable and may affect the extent of the damages.
 The analysis of the stress state allows clarification of the extent of the damage and how it is caused. Furthermore, the analysis reveals that the material has accumulated a proportion of the impact energy after the impact, which is highlighted through the existence of the residual stresses.
 The extent of the damaged area was slightly greater when using a HemI. Indeed, as the contact area increases, the depth of impact decreases, which implies the shape of the impactor head. For a small contact area impactor, the visual impact damages are worse on the frontal face than on a large face in the case of an equienergetic impact, as well as on the rear face.
 The determined impact responses and damage characteristics are an important factor to assess the effect of the impactor head shape on the low-velocity impact and how they correlate.

www.etasr.com Rabouh et al.: Bridging the Effect of the Impactor Head Shape to the Induced Damage during Impact at …
 Finally, it can be said that the implemented model strongly agrees with the experimental findings of [1].
Our investigations have enhanced the understanding of the behavior of laminates under low velocity impact and also enabled the elaboration and validation of a three-dimensional finite element model based on the LARC criterion to predict the response of the T700/Epoxy laminate [30]. The proposed model will be important for future research into the predictability and the enhancement of the impact behavior of composite laminates. Young's modulus in X-direction MPa ε22 Young's modulus in Y-direction MPa ε33 Young's modulus in Z-direction MPa G12 Shear modulus in XY-plane direction MPa G13 Shear modulus in XZ-plane direction MPa G23 Shear modulus in YZ-plane direction MPa X T Tensile strengths in X-direction MPa Y T Tensile strengths in Y-direction MPa X C Compressive strengths in X-direction MPa Y C Compressive strengths in Y-direction MPa Transverse to fiber direction shear strengths (S12) MPa Tensile strains in X-direction (%) (%) Tensile strains in Y-direction (%) (%) Compressive strains in X-direction (%) (%) Compressive strains in Y-direction (%) (%) Transverse to fibre direction shear strains (%) b 5" Poisson's ratio modulus in XY-plane direction b 5 Poisson's ratio modulus in XZ-plane direction b " Poisson's ratio modulus in YZ-plane direction c 8 Volumetric fraction of fibers (%) d Density g/mm 3 S11, S22, S33 Main stresses in directions X, Y, and Z MPa APPENDIX B The FFI damages for the presented criterion mentioned in Section V are presented below.