Prediction of Load-Bearing Capacity of Composite Parts with Low-Velocity Impact Damage: Identification of Intra- and Inter-Ply Constitutive Models

Assessments of residual load-carrying capacity are often conducted for composite structural components that have received impact damage. The availability of a verified simulation methodology can provide significant cost savings when such assessments are required. To support the development of a reliable and accurate simulation methodology, this study investigated the predictive capabilities of a stacked solid-shell finite element model of a cylindrical composite component with a damage mechanics-based description of the intra-ply material response and a cohesive contact model used for simulation of the inter-ply behavior. Identification of material properties for the model was conducted through mechanical characterization. Special attention was paid to understanding the influence of non-physical parameters of the intra- and inter-ply material models on predicting compressive failure load of damaged composite cylinders. Calibration of the model conducted using the response surface methodology allowed for identifying rational values of the non-physical parameters. The results of simulations with the identified and calibrated finite element model showed reasonable correlation with experimental data in terms of the predicted failure loads and post-impact and post-failure damage modes. The investigated modeling technique can be recommended for evaluating the residual load-bearing capacity of flat and curved composite parts with impact damage working under the action of compressive loads.


Introduction
Impact damage induced by hailstone impact, tool, or equipment dropping can lead to severe reductions in composite structures' load-carrying capacity. Aerospace companies and manufacturers of other products, in which composite materials are extensively used, spend considerable resources to determine the level of degradation of composite parts' load-bearing capacity that have received impact damage during operation or during assembly, as well as the permissible degree of damage at which the replacement of an expensive structural member is unnecessary. Usually, such assessments are based on the integrated application of experimental destructive and non-destructive methods, which, in turn, also requires considerable financial investments and time. Understandably, the availability of a verified simulation approach capable of predicting the residual load-carrying capacity of composite parts with impact damage would provide significant cost savings and accelerate the decision making when such assessments are required.
Previous studies considered modeling the impact damage and post-impact load-bearing capacity of damaged composite parts. Multiple publications are focused on the first aspect of the problem-predicting damage induced by the impact on composites. Among them, Mouli et al. [1] compared macro-and mesoscale damage models and reported better correlation with experiments for the latter (11-20% overprediction of max impact force with the mesoscale model as compared with 29-35% underprediction with the macroscale model based on Hashin criterion). Moncayo et al. [2] used the commercial software LS-DYNA R10 (Livermore Software Technology Corporation, Livermore, CA, USA) to predict damage induced by low-speed impact in carbon/epoxy laminates. It was found that, to achieve correlation with experimental data in terms of the extent of damage, an artificial adjustment of parameters of the inter-ply tiebreak contact-based damage model was needed. Aymerich et al. [3] used interface elements with a bilinear cohesive failure model in Abaqus/Explicit to predict delamination in cross-ply laminates. After calibrating the cohesive model, fair agreement with experiments could be achieved in terms of shape, orientation, and size of the delaminated zone. Good agreement with experiments in terms of the measured and predicted damage area (within 16%) was reported in [4]. Long et al. [5] compared prediction of shapes of delamination in different interlayers of a laminate (obtained using the bilinear traction-separation cohesive contact) with experimental data obtained using the ultrasound C-Scan inspection. The importance of having multiple through-the-thickness cohesive interfaces for accuracy of numerical predictions was emphasized. In a recent numerical study, Zhang et al. [6] identified that over 30% of low-speed impact energy can be absorbed through the inter-ply failure (delamination), whereas absorption of the rest of the energy was associated with the intra-ply damage. A comprehensive comparison of several intra-ply material models widely used in impact simulations of composites can be found in [7,8]. Several studies have been focused on improvement of the computational efficiency of lowspeed impact simulations (see, for example, [9,10]), as well as controllability of impact experiments [11]. For the post-impact load-bearing capacity, most of the reported studies were focused on modeling the compression after impact (CAI) behavior of flat composite panels. Examples of this work include CAI simulations of pre-loaded composite plates [12], variable angle tow laminates [13], and flat sandwich panels [14,15].
Unlike previous studies, where mainly simple plate-like geometries were considered, this study is focused on developing a simulation model capable of accurately predicting the load-bearing capacity (failure load) of impact damaged composite structural members. In particular, curved structural components in the form of glass fabric-reinforced composite cylinders were considered. Both intra-and inter-ply properties of the cylinders' material were measured and reported, as described in Section 3.1. Subsequently, two-phase experiments, a schematic of which is shown in Figure 1, were conducted with the composite cylinders. In phase 1, the cylinders were subjected to low-velocity impacts by a 20 mm hemispherical-tip rigid striker with different impact energies (5.15 J and 7.65 J). In the second phase, damaged composite cylinders were tested in quasi-static compression until failure. Both phases of the physical experiments were then simulated in LS-DYNA and results of numerical predictions were compared with the tests outcomes. The influence of nonphysical parameters of the damage mechanics-based intra-ply and cohesive zone inter-ply models was studied and identification of their rational values was performed.

Materials and Manufacturing
A composite material containing approximately 50% of epoxy resin by volume and reinforced by T-26 structural glass fiber fabric was used. The woven fabric had the nominal areal weight of 285 ± 12 g/m 2 , and 11-12 and 6-8 yarns per cm in the warp and fill directions, correspondingly. Samples for material characterization were cut from flat panels manufactured using out-of-autoclave vacuum bagging. Similarly, vacuum bagging was employed in fabricating the composite cylinders. The main stages of the latter process are illustrated in Figure 2. After curing, composite parts were trimmed to have the overall length of 160 mm. The warp direction of all fabric layers was oriented at 90 degrees with respect to the cylinder axis.

Hand lamination
Bleeder/breather plies Vacuum bagging Manufactured part

Experimental Methods
All quasi-static experiments, including the material characterization tests and post-impact compression of the composite cylinders, were conducted using an MTS Systems Corporation test frame equipped with a 100 kN load cell. Crosshead displacement of 2 mm/min was used.

Material Characterization
For material characterization, both intra-and inter-ply properties of the woven glass fabric/epoxy matrix vacuum-bagged composite were experimentally determined. All intra-ply properties tests were conducted according to the corresponding ASTM standards ( Figure 3). Particularly, tensile properties of the composite material along warp and fill directions were determined using the procedure described in the ASTM D 3039 standard [16]. For the shear properties, specimens with a ±45 degree layup were tested in tension according to the ASTM D 3518 test procedure [17]. In these tests, a biaxial extensometer MTS 632.85F-14 was used for strain measurements along and perpendicular to the loading direction. The glass fiber reinforced plastic (GFRP) composite's compressive properties were obtained according to the ASTM D 3410 standard using the loading and alignment fixtures described in the standard [18]. ASTM   The composite material's inter-ply properties were characterized using the double-cantilevered beam (DCB) and end-notched flexure (ENF) tests, which allowed for measuring the critical strain energy release rate in mode I (GIC) and mode II (GIIC), correspondingly. The former test was conducted following the method described in the ASTM D 5528-01 standard [19], according to which the ends of DCB specimens were subjected to controlled opening displacements while the load and delamination length were recorded ( Figure 4). For the latter test, the method based on subjecting endnotched specimens to three-point bending was employed. A description of this test procedure can be found in [20].

Testing of Composite Cylinders
A simple drop-weight apparatus, as shown in Figure 5a, was built to induce controlled-energy impact damage to the composite cylinders. As part of this apparatus, a simple mechanism preventing secondary impacts to the cylinders was implemented. This mechanism fixed the stainless-steel striker after bouncing from the surface of the cylinder for the first time. The striker in the fixed position above the composite specimen right after impact is shown in Figure 5b. In all impact tests, strikers were released from a 1.5 m height. To induce different levels of damage, strikers of 350 g and 520 g were used. In order to provide support during the impact tests, metallic tubes of the same external diameter as the internal diameter of the composite cylinders were inserted in the test specimens from each end to a depth of 10 mm. Following the impact experiment, the damaged composite cylinders were tested until failure under unidirectional compression using the 100 kN MTS test frame. To achieve uniform distribution of the compressive force across the loaded sides of the specimens, thin rubber pads were inserted between the compressive platens of the test frame and the composite cylinders, as can be seen in Figure 5c.

Modeling Techniques
An explicit finite element model of the two-stage process, involving impact damaging of the composite cylinders and compression after impact, was developed in LS-DYNA. The model's major features are illustrated in Figure 6. The steel striker in the model was represented by a rigid body with the assigned translational mass (350 g or 520 g, depending on the modeled impact conditions) and the vertical speed just before impact of 5.425 m/s applied using the *INITIAL_VELOCITY_RIGID_BODY keyword. The *CONTACT_ERODING_ SURFACE_TO_SURFACE contact algorithm was used to simulate the striker and the composite cylinder's interactions. The walls of the cylinders were modeled using three layers of stacked TSHELL elements, each representing two "physical" plies of the woven glass fabric oriented at 90 degrees, that is, with the warp direction being perpendicular to the cylinder's axis. To represent delamination, the two contact interfaces between the three TSHELL layers were modeled using *CONTACT_AUTOMATIC_ONE_WAY_SURFACE_ TO_SURFACE_TIEBREAK with OPTION 9. This contact algorithm is equivalent to using zero-thickness cohesive zone elements and is based on the fracture model with bilinear traction-separation law, mixed-mode delamination criterion, and damage formulation [21]. This delamination model's main parameters include normal and shear failure stresses at the interface between the adjacent layers (NFLS and SFLS parameters), mode I and mode II critical strain energy release rates (GIc and GIIc), and the normal (CN) and tangential (CT = CT2CN × CN, where CT2CN is a coefficient between 0 and 1) stiffness of the material in the interlaminar region. The bilinear law corresponding to a simple case of crack opening (mode I) is exemplified in Figure 7. A detailed description of the mixed-mode loading treatment in this model can be found in [21]. After satisfaction of the failure criterion (δ > δult), the failed interface (master segment-slave node pair) can only resist compressive forces. Tabs of the composite cylinders were represented by an additional layer of TSHELL elements and connected to the body of the cylinder using the simple *CONTACT_TIED_SURFACE_TO_SURFACE method. Overall evaluation of the contacts between the stacked TSHELL layers was conducted by running implicit eigenvalue analysis and observing the deformation modes corresponding to different eigenvalues. As can be seen in Figure 8, modeled composite cylinders deform as a single entity, which confirms that all contacts were appropriately engaged. MAT058 or *MAT_LAMINATED_COMPOSITE_FABRIC-a damage mechanics-based model, which accounts for both pre-and post-peak softening of composite plies-was used to model the intra-ply behavior of the woven glass fabric-reinforced composite cylinders. Details of this model's implementation can be found in [22], and its brief summary is provided below (see Equations (1)-(5)). Depending on the type of failure surface, this model may be used to represent composite materials with unidirectional layers, complete laminates, and woven fabrics. In this study, the following set of failure criteria was used to simulate the behavior of the fabric-reinforced composite (a.k.a failure surface FS = −1): Here, XT,C and YT,C represent the material's strength in the warp (X) and fill (Y) directions in tension ("T"; used if σij > 0) or compression ("C"; used if σij < 0), and S is the fabric shear strength. It should be noted that the effective stresses ( ) in the above expressions are related to the nominal stresses through the damage parameters , also known as area loss parameters, such that, where damage evolution with straining is assumed as = 1 − exp − ⋅ ; and m, , and are the parameters controlling the shape of the stress-strain response, strain, and strain at maximum directional stress, respectively. Thus, the components of the constitutive tensor can be represented as functions of the damage parameters and the properties of the undamaged layer: Such continuous damage mechanics-based (CDM) formulation provides a smooth increase of damage and, on failure initiation, prevents the instantaneous drop of stresses in the failing element. It should also be noted that the two damage parameters and assume different values for tension ( ) and compression ( ) . Additional non-physical parameters associated with MAT058, as well as reasons for their choice, are described in Table 1. Chosen to be significantly higher than any directional strain at failure initiation.

SLIMT1
Factor to determine the minimum stress limit after stress maximum (fiber tension).
-0.10 A recommended value [21] SLIMC1 Factor to determine the minimum stress limit after stress maximum (fiber compression).

SLIMT2
Factor to determine the minimum stress limit after stress maximum (matrix tension).
-0.10 A recommended value [21] SLIMC2 Factor to determine the minimum stress limit after stress maximum (matrix compression).

SLIMS
Factor to determine the minimum stress limit after stress maximum (shear). -1.00 A recommended value [21] The post-failure response of composite in MAT058 is governed by an array of stress limit factors (SLIM_), which represent the amount of residual strength the composite retains after the element's complete failure. For example, even a completely crushed composite can usually retain some resistance to compressive loading, which can be conveniently represented by stress limit factors in compression along warp (SLIMC1) and fill (SLIMC2) directions. As a result of the compressive nature of loading experienced by the cylinders, it was expected that SLIMC1 and SLIMC2 would be the most influential parameters among all stress limit factors. Correspondingly, values of SLIMT1, SLIMT2, and SLIMS were chosen based on the software developer's recommendations provided in [15], whereas the values of SLIMC1 and SLIMC2 were obtained from the range of 0.1-1.0 via calibration with experimental data, as will be discussed in Section 5.2. The upper bound of this range (SLIMC1 = SLIMC2 = 1.0) corresponds to the developer's recommendation. However, the author's experience of using MAT058 suggests that the "optimal" values of these parameters in compression-dominant loading cases may be significantly different from those recommended by the developer (see, for example, the work of [7], in which the value of 0.375 for SLIMC1 and SLIMC2 was obtained through calibration for the case of axial crushing of CFRP energy absorbers).
The idea of stress limit factors is schematically illustrated in Figure 9, where E11C corresponds to the strain at compressive strength XC and GMS is the strain at in-plane shear strength S. With failure surface FS = −1, MAT058 provides an option for defining the shear stress-shear strain diagram as a bi-linear curve, which allows for the improved approximation of nonlinear shear behavior of fabric-reinforced materials, which is usually observed experimentally. This model is schematically illustrated in Figure 9 (right). Parameters TAU1 and GAMMA1 represent shear stress and shear strain, correspondingly, at which behavior of a fabric-reinforced composite becomes nonlinear.
The parameters of the delamination model used in this study, as well as the rationale for their choice, are summarized in Table 2. It should be noted that this delamination model uses such a measure of fracture toughness as the strain energy release rate (Gc) as the input parameter. In addition, the elastic stiffness (slope) and the peak stress are required for complete definition of the bilinear law (see Figure 7) for each mode of fracture (normal or shear), while the ultimate separation displacement (δult) is calculated based on the provided values of initial stiffness, peak stress, and fracture toughness. Numerical studies have shown that the fracture toughness has to be accurate, while the initial stiffness and the peak stress can be varied without affecting the overall results. For example, in [23,24], a constant value for initial stiffness has been used, while the peak stress has been adjusted accordingly to keep the fracture toughness (the area under the triangle) correct. Similarly, a constant value of initial normal stiffness (CN parameter in Table 2) was used in this study, while an appropriate value of normal peak stress (NFLS parameter in Table 2) was obtained from the range provided in Table 2 via calibration with experimental data, as will be discussed in Section 5. Initial tangential stiffness and shear peak stress were related to the corresponding values of the normal stiffness and normal peak stress (see the rationale for choosing CT2CN and SFLS in Table 2).  ), which would be a reasonable estimate in the case of interlaminar failure by adhesive mechanism (cracks formed at the interface between the epoxy in the interlaminar resin-rich region and fibers in the layer adjacent to it). upper bound-the ultimate strength of bulk epoxy resin (~60 MPa), which would be a reasonable estimate in the case of interlaminar failure by cohesive mechanism. In addition, a scaling factor of 0.30 was used to account for the mesh dependency observed for this delamination model (see the recommendation provided in [25] for meshes with element sizes between 2 and 3 mm). A particular value from the specified range was chosen via calibration with experimental data, as will be discussed in Section 5.2.

Assumed as SFLS = NFLS/ 3 (von Mises criterion)
G_Ic, kJ/m2 0.24 Measured experimentally, see Table 3 in Section 5.1 G_IIc, kJ/m2 1.96 Measured experimentally, see Table 3  In addition, the condition for CN > CNmin must be ensured (see [15]), where CNmin = (1/2) × (NFLS 2 )/(G_Ic).This condition is satisfied for the listed set of parameters of the delamination model. As can be seen in Figure 6 (see Section 4 in which the model was described), SPC constraints were applied on one side of the cylinder, while rigid 1D elements connected at a master node were used to define the boundary conditions on the other side of the cylinder. A displacement imitating the test frame crosshead movement was applied to the master node. The corresponding loading rate profile is shown in Figure 10. As can be deduced from the figure, no compressive loading was applied in the first 10 ms when the cylinders were subjected to drop-weight impact. Then, the speed of the master node was gradually increased to 20 mm/s within the next 10 ms and remained at that level for the rest of the simulation. This approach allowed simulating both loading cases (impact and compression after impact) within the same analysis, thus ensuring the seamless transition from impact to compressive loading. It should be noted that, compared with the experiment, somewhat higher compressive loading speed was used in simulations to reduce the computational time. Overall, the developed model contained 43,600 nodes and 54,000 elements. A typical simulation with 8 SMP threads on a desktop computer with 16 Gb RAM and intel i7 processor took approximately 24 h.

Experimental Results
Material properties obtained during the material characterization campaign described in Section 3.1, including the fracture mechanics properties, are summarized in Table 3. The glass fabric composite exhibited a highly nonlinear shear response, as can be seen in Figure 11. Considering this, the MAT058 failure surface FS = −1 was used with the bilinear approximation for the shear stressshear strain diagram (see Figure 11).  Compressive force-displacement diagrams for the two experimentally tested pre-damaged composite cylinders are shown in Figure 12. The initial nonlinearity of the compressive response visible in the figure should be attributed to the deformation of the rubber pads inserted between the compressive platens and the cylinders to provide a uniform load distribution (see Section 3.2). It should be noted that, because of the presence of the rubber pads in experiments and the lack of them in simulations, no direct comparison of the predicted and experimentally obtained forcedisplacement curves was possible. Therefore, only maximal (peak) forces at cylinders' failure were used to compare the simulations and the experiments. Maximal compressive force endured by the composite cylinder impacted by a 350 g striker corresponded to approximately 30 kN, while the cylinder impacted by a 520 g mass could support the force of approximately 28.5 kN before failure. This corresponds to approximately 70% and 67% of the failure load of undamaged composite cylinders with the given layup and dimensions (42.78 kN, as measured in an additional experiment).

Load-Bearing Capacity: Calibration of SLIMC and NFLS Parameters
SLIMC parameters of MAT058 (see Table 1) directly affect intra-ply compressive behavior of the composite represented by this material model, while NFLS and SFLS parameters of the delamination model (see Table 2) may affect the extent of delamination and inter-ply damage, thus reducing the buckling strength of the compressed cylinders. It should be noted that, in each simulation, the value of the SFLS parameter was determined as a function of NFLS parameter's value, as described in Section 4 (see Table 2). Also, it was assumed that SLIMC factors in the fill and warp directions are of the same magnitude, that is, SLIMC1 = SLIMC2, which seems to be a reasonable assumption for a fabric-based material. As a result, rational values of the following two parameters required identification through calibration with experimental data: SLIMC and NFLS.
To conduct the calibration, a response surface for the predicted max compressive force endured by the damaged composite cylinder (520 g striker impact), as a function of SLIMC and NFLS parameters, was built using a genetic aggregation algorithm. To build the response surface, sampling points were chosen according to the value ranges provided for SLIMC and NFLS parameters in Tables 1 and 2, correspondingly, and using the face-centered central composite design DOE (Design Of Experiments) plan. This resulted in, overall, nine combinations of SLIMC and NFLS parameters listed in Figure 13, for which numerical analyses were performed using the LS-DYNA simulation model described in Section 4. The resulting response chart is shown in Figure 14, in which the predicted max force is normalized by the experimentally measured value (the ratio of unity means the exact equality of predicted and experimentally measured failure loads). Goodness of fit for the obtained response surface was evaluated and was found to be adequate (see Figure 15). The scatter chart in Figure 15 is intended to show if the generated response surface correctly fits the simulation data based on which it was generated; the closer the points are to the diagonal trend line, the better the response surface fits the points.  As can be deduced from Figure 14, the model is sensitive to the compressive stress limit factor; depending on the choice of SLIMC parameter, it can either significantly (±40%) under-or overpredict the ultimate compressive force of the damaged composite cylinder, as compared with the experimentally measured value. At the same time, the effect of the NFLS parameter is insignificant with low-to-medium values of the SLIMC factor (less than 0.6) and becomes more noticeable only when SLIMC is large (0.6 < SLIMC < 1.0). Such behavior of the model is somewhat intuitive and can be explained by the fact that the predicted extent of delamination is relatively low at the given impact conditions (see Section 5.3), and failure of damaged composite cylinders under compressive loading is mainly governed by the intra-ply damage. However, as the residual compressive strength of composites after failure increases with an increase of SLIMC, inter-ply damage becomes the factor that triggers compressed cylinder failure.
According to the response chart presented in Figure 14, the following three combinations of SLIMC and NFLS parameters can provide good correlation between the predicted and experimentally measured ultimate load-bearing capacities of the composite cylinder damaged by a 520 g striker: (1) SLIMC = 0.60, NFLS = 9.0 MPa; (2) SLIMC = 0.57, NFLS = 13.5 MPa; and (3) SLIMC = 0.56, NFLS = 18.0 MPa. To verify these predictions of the response surface, LS-DYNA simulations were conducted with all three identified parameter sets. These simulations, in addition to 520 g striker impacts, considered compression of the cylinders damaged by 350 g strikers. The results of these simulations are presented in terms of force-displacement diagrams in Figure 16a-c. In this figure, the peak forces determined experimentally (30 kN for the cylinder damaged by a 350 g striker and 28.5 kN for the 520 g striker-damaged cylinder) are shown by the horizontal level lines. The red and green curves correspond to compressive tests of cylinders damaged by 350 g and 520 g strikers, correspondingly.
As can be seen in these figures, all three combinations provide similar estimates of the loadbearing capacities of the composite cylinders. However, the simulation with the first set of parameters for the case of the 520 g impact (green curve in Figure 16a) slightly overpredicts (by approximately 4%) the experimentally observed peak force. It was, therefore, judged that the models with higher values of NFLS and lower SLIMC (cases b and c in Figure 16) provide more accurate predictions of the load-bearing capacity. It should also be noted that the compressive forces in all simulations quickly drop to the level of approximately 10 kN after cylinder failure, which correlates well with the physical experiments ( Figure 12).

Predicted and Experimentally Observed Damage
Predicted impact damage (simulations with SLIMC = 0.57, NFLS = 13.5 MPa), as well as the impact damage visible on the tested specimens, is illustrated in Figures 17 and 18. In these figures, damage is shown for the outer ply only, while delamination is only shown for the interface between the outer and middle plies. As can be deduced from Figures 17 and 18, visible damage (elliptical area marked by red arrows in Figures 17 and 18) at the considered impact energies mainly corresponds to matrix cracking produced by shear stress during impact, while the predicted extent of delamination is relatively small. Traces of damage resembling the predicted transverse damage pattern can also be noticed on the experimental images (marked by green arrows in Figures 17 and 18; also noticeable in Figure 19). It should be noted that precise correlation between the extent of visible and predicted damage was not expected, as not all areas in which impact-induced microcracks are present may be clearly seen by the naked eye. This can explain the somewhat larger extent of predicted damage, as compared with the experimental images. In our future studies, an ultrasound inspection (C-scan) will be employed to address this visualization issue.

Conclusions
This study investigated the predictive capabilities of a stacked finite element model of a cylindrical composite component with a CDM-based description of the intra-ply material response (LS-DYNA MAT_LAMINATED_COMPOSITE_FABRIC model, MAT058) and a cohesive contact model used for simulation of the inter-ply behavior (SURFACE_TO_SURFACE_TIEBREAK contact, OPTION 9). Identification of material properties for the model was conducted through mechanical characterization (ASTM standards D3039, D3410, D5528-01 tests; ENF test). Special attention was paid to understanding the influence of non-physical parameters of the intra-and inter-ply material models on predicting compressive failure load of damaged composite cylinders. Calibration of the model conducted using the response surface methodology allowed for identifying rational values of the non-physical parameters.
The conducted investigation suggests that, with the widely used MAT058 model and the delamination model based on SURFACE_TO_SURFACE_TIEBREAK contact (OPTION 9), material properties obtained from conventional mechanical tests do not provide all the required information for accurately predicting the load-carrying capacity of composite structural components experiencing post-impact compressive loading. In addition to the tests, identification of non-physical parameters of the intra-and inter-ply models is required to achieve good correlation with experimental data. In this study, such an identification was conducted for the compressive stress limit factors (SLIMC) of MAT058 and the peak stress parameter (NFLS) of the delamination model. A response surface, relating the predicted cylinders' load-bearing capacity with the SLIMC and NFLS parameters, was built and used to determine their rational values. As a result, it was found that the models are sensitive to the compressive stress limit factor; depending on the choice of SLIMC, they can either significantly under-or overpredict the ultimate compressive force of the damaged composite cylinder, as compared with the experimentally measured value. The effect of the NFLS parameter was found to be insignificant in combination with low-to-medium values of the SLIMC factor (less than 0.6) and more noticeable when SLIMC was large (0.6 < SLIMC < 1.0). The numerical model identified using this method was able to accurately predict the failure loads observed in the physical experiments, including the test used in the parameters' identification process (compression of the composite cylinder pre-damaged by a 520 g striker), and the test that was not used for the parameters' calibration (compression of the composite cylinder pre-damaged by a 350 g striker). In addition, fair correlation was observed between the experimental and the predicted post-impact and post-failure damage modes of the composite cylinders. Therefore, the investigated modeling technique can be recommended for evaluating the residual load-bearing capacity of composite parts with impact damage working under the action of compressive loads.
Future work will be focused on extending the analysis and verification studies to prepreg-based carbon fiber composite systems, as well as other post-impact loading conditions (bending and torsion). It is also planned to employ non-destructive ultrasound inspection for more precise measurements of the impact damage.

Funding:
The authors would like to acknowledge financial support of this study provided by Samara University and the University of Windsor.

Conflicts of Interest:
The authors declare no conflict of interest.