Composite Material Failure Model Updating Approach Leveraging Nondestructive Evaluation Data

A methodology to demonstrate how nondestructive evaluation (NDE) data could be leveraged to update computational modeling approaches related to progressive failure analysis (PFA) of composite materials is presented in this paper. The particular case of a carbon fiber-reinforced polymer composite (CFRP) is chosen as a characteristic case of a material for which a wealth of mechanical behavior, characterization, failure, and NDE data and information is available [1,2]; however, it remains challenging to develop computational models for any component or structure that is built with this material that could be informed by such data in ways that would relate to PFA of the data-driven type [3–5].

The basic reason for such connection between data and PFA using computational models is the fact that when experimental data sets are obtained, the intent in many cases is to comply with existing guidelines for certification and fleet management, as for example prescribed by organizations such as NASA, the U.S. Air Force of the Federal Aviation Administration (FAA) and not to create comprehensive data sets aligned, for example, with the Digital Twin concept [4,6]. Consequently, although data exist on, for example, fracture properties, such data were not necessarily obtained simultaneously with, for example, NDE monitoring or by using a design of experiment approach that would allow its hierarchical use from the specimen to the structural level [4]. In addition, the lack of concurrent development of appropriate physics-based models that could read and be driven by sensing data while being capable to account for variability in loading and material state conditions create challenges in data-driven modeling appropriate for composite structures [5,7].

The fact is that both material characteristics and failure processes in composites are fairly complex and highly variable, depending on manufacturing, processing, geometrical, and testing parameters, which result in several challenges when monitoring or simulating their failure. In general, mechanical and/or thermal loading causes microcracks in the matrix in addition to small-scale delaminations at the layered interfaces. Further increase in such loads leads to the growth of microcracks into mesoscale cracks, fiber breaks, as well as delaminations, which eventually lead to the formation of macroscale cracks and fiber pullouts at the failure state [8–11]. This process indicates that damage initiation and evolution in composite materials are inherently multiscale and progressive in the sense that once initiation occurs there are multiple stages of evolutionary degradation before the final failure occurs [12]. In this context, several different approaches to test, monitor, characterize, and simulate the failure process in composites have been developed, as explained next; however, these vary in several ways related to the length scale targeted, as well as their accuracy and reliability in relating local damage with global behavior, in addition to their potential to be used on actual structures versus only in lab conditions.

In the case of the type of composite material used in this paper, Hallett et al. investigated both un-notched and notched specimens tested in tension [13]. A strong influence of matrix cracking and delamination on the overall damage of the composite was found. Differences both in failure stress and in damage mechanisms were reported in the case of notched specimens due to changes in test configuration (ply and laminate thickness, and hole diameter). Moreover, Nixon-Pearson et al. investigated damage development of open-hole specimens loaded in tension–tension fatigue [14]. Computer tomography images of interrupted tests showed that initially damage starts to propagate out from the edge of the hole in terms of matrix cracks and delamination. Moreover, Schaefer et al. experimentally characterized the failure progression using cross-ply and quasi-isotropic laminate specimens and confirmed matrix-dominated failure modes [15].

In this context, multiple research endeavors have been undertaken with respect to monitoring and assessing the damage evolution process in composite materials using NDE methods [12,16]. Among them and in relation to the methods used in this paper, acoustic emission (AE) has been widely used to evaluate the failure processes [9,17–20]. AE refers to the energy release in the form of transient elastic waves due to a sudden change in a material or structure subjected to external loads (thermal/chemical/mechanical), and therefore, it provides a volumetric assessment of the progressive failure process [21,22]. Moreover, the digital image correlation (DIC) method, which is an optical metrology method that typically employs a pair of digital cameras to feed stereoscopic vision images into computer algorithms that track changes due to deformation, has also been extensively used to monitor strain localizations at tow crossings and boundaries between plies and coupon edges [12,16,23,24]. In post-processing, DIC can also be used to capture general material property information such as elastic constants [16,25]. By nature, DIC provides the surface projection of a much more complex volumetric response. Other NDE methods used include infrared thermography (IR) [26,27], X-ray-based computed tomography [28,29], and electrical resistance methods [30–32]. Accompanying post-processing techniques of NDE data have also been reported to extract information. Such techniques include data mining and pattern recognition methodologies [33–35]. Data mining methodologies are particularly successful when it comes to combatting convolution effects caused by the simultaneous or overlapping activation of multiple damage sources [20,34,36,37]. In this context, several examples of machine learning methods have been reported [38–40]. The performance of such classification algorithms can be divided into three groups: unsupervised and supervised classifications, as well as those that seek relationships to the physical problem [41–43].

Regarding modeling of failure in composites, there are mainly two approaches that could be classified as fracture and damage mechanics; the latter is the approach followed in this paper. The fracture mechanics approach generally seeks to explicitly represent, for example, a crack as an internal discontinuity, while crack equivalent effects, i.e., stiffness and strength reduction, are typically modeled by damage mechanics approaches. Fracture mechanics approaches for composites have been implemented using a variety of numerical methods [44,45]. In the case of damage mechanics approaches, continuum-based failure criteria are widely used to relate stresses and experimental data of material strength to the onset of failure. These include several failure theories such as Hashin [46] and Tsai and Wu [47] which have been proposed to predict the initiation of various damage mechanisms such as fiber breakage, fiber buckling, transverse matrix cracking, and in-plane shear failure.

The use of such failure criteria is not sufficient though to fully describe the quasi-brittle failure of composite laminates that results from the activation of several failure mechanisms. To study the nonlinear response due to the accumulation of damage, continuum damage mechanics (CDM) methods have been reported. In CDM methods, the accumulation of damage is simulated by reducing specific material properties which causes the loss of load capacity in a specific direction [48]. In fact, the macroscopic stress–strain behavior changes gradually during the simulation as a result of a distributed network of microcracks. This idea was first implemented in metals by employing a single scalar damage representing microcrack density over a continuum medium [49]. This concept indicates that several intermediate phases may exist between intact material and final fracture. This theory has been extended by introducing numerous damage variables associated with different failure mechanisms to model damage initiation and evolution in composites [8]. Matzenmiller et al. [50] developed a constitutive model for anisotropic damage to describe the relationship between the damage of the material and the effective elastic properties for stress analysis of structures. Internal variables were employed to describe the evolution of the damage state under loading and as a subsequence the degradation of the material stiffness. In this framework, CDM considers damaged materials as continuous, despite the inherent heterogeneities, micro-cavities, and defects. This nonlinear response to the loading conditions is predicted based on the constitutive relations between macroscopic variables (e.g., stress and strain) and internal variables which model, on a macroscopic scale, the irreversible changes occurring at the microscopic level [51].

In orthotropic or transversely isotropic materials, such as fiber-reinforced composites, the microstructural morphology of the material causes more difficulties in using a CDM approach. Specifically, crack orientations are not only affected by the geometry and boundary conditions. For instance, the interface between fiber and matrix is weaker than the surrounding material, and interfacial decohesion is the first damage mechanism to occur followed by matrix cracking. Hence, damage in composites is a multiple scale phenomenon coupling different scales of damage initiation and progression. In micromechanical models, a representative volume element (RVE) is conducted with subsequent homogenization, to predict evolving material damage behavior [52]. Micromechanical damage approaches can be computationally very expensive, however, since detailed micromechanical analyses need to be conducted in each load step at every integration point of the macroscopic structure, even though they are effective in characterizing damage.

As an alternative, in phenomenological CDM, stiffness degradation is modeled by employing the anisotropic tensorial damage variable using mathematically and thermodynamically consistent models of damage mechanics [48,53]. Such phenomenological CDM models do not explicitly consider microstructural damage. Instead, the damage events evolve at the constituent material scale with some knowledge on the constituent properties. Overall, predictions from computational models, therefore, may encompass a high level of uncertainty, while they might require additional material testing to address the empirical evolution laws.

Based on this introduction, this paper targets the development of a model updating approach that is based on a key difference with respect to the previously reported approaches on composite damage modeling. Specifically, instead of using a number of inputs for material and failure properties based on literature or some experiments [54,55], the presented approach targets the use of NDE data collected during controlled mechanical testing to generate progressive failure inputs to a finite element method (FEM) approach that is first calibrated for its predictions and then uses the NDE data to tune aspects of damage initiation and evolution. Hence, instead of solely using a number of ad hoc experimental inputs to mathematical criteria in damage modeling approaches, the proposed approach uses the NDE data to form specific inputs given the stress–strain curve at the laminate level to tune the activation and increase of a damage parameter which tracks the progressive failure of the composite, while it considers the development of several individual failure modes. With this approach, mechanical and NDE data are leveraged, while the damage model is then used to resolve damage information as a function of the stacking sequence.

A 16-ply carbon fiber-reinforced polymer (CFRP) composed of 8552 epoxy resin reinforced with unidirectional IM7 carbon fiber prepreg sheets (manufactured by Hexcel) was used for this investigation. The 16-ply layup had a [45/0/-45/90]_2s stacking sequence, as shown in Fig. 1(a). The fiber volume fraction for the laminate was nominally 60%. Figure 1(b) shows a schematic of the layup with the 0-deg fiber orientation aligned with the x-axis of the plate and the 90 deg aligned with the y-axis. The composite specimens used for this work were cut using a water jet to produce straight edge (SE) samples per ASTM D3039. All specimens had a final nominal thickness of 5 mm, width of 25 mm, and length of 250 mm in the loading direction (aligned with the y-axis in Fig. 1(b). To quantify the initial condition of the material, a 10 × 25 × 5 mm section was cut from one of the “as- received” specimens for X-ray scanning. A ScanCo µCT100 scanner was used with a 90-kV accelerating voltage resulting in good void contrast and a resolution of 10 μm/voxel. Figure 1(c) shows the results of the scan where the red areas indicate voids/microcracks produced during manufacturing. As the voids do not appear to be uniformly distributed, variation in the initiation and evolution of damage mechanisms and subsequently in the mechanical response during loading was expected from different specimens.

In addition, specific material property sets were used based on literature related to elastic constants and failure properties for the composite lamina (Table 1).

Furthermore, Table 2 provides reported energy release rates that were used in the modeling approach. The physical meaning of the parameters in Tables 1 and 2 is described in Sec. 3.

All specimens were monotonically loaded in tension until failure using an MTS 370.10 Landmark servo-hydraulic load frame equipped with a 100-kN load cell. The load was applied in displacement control at a rate of 2 mm/min based on ASTM 3039, while the specimens were monitored by a combination of AE, DIC, and passive infrared thermography (pIRT). The front surface of the specimen was monitored with DIC to observe the localization and evolution of strain concentrations that occur as a result of the mechanical load. A 5-megapixel stereo camera system (manufactured by GOM) was used with a field of view (FOV) of 135 × 115 mm and a frame rate of 1 Hz so that both in-plane strain and out-of-plane motion could be measured. The required intensity contrast was achieved with a random pattern of black speckles applied on a white background with commercially available spray paint. Strain maps were generated using a subset and step size of 25 and 10 pixels (1.35 and 0.54 mm, respectively), which produced a maximum noise level of 1100μ strain. The average strain in the specimen was determined using a virtual extensometer placed across the monitored gauge length.

The surface temperature was monitored with an infrared camera (FLIR A325sc) having a 320 × 240 resolution and a detectible temperature range between 0 °C and 350 °C. Video files of pIRT data were recorded at 14 Hz, and images were post-processed to observe changes in temperatures between frames. Finally, AE activity was monitored using four sensors (PICO, 150–750 kHz operating frequency) symmetrically placed with respect to the center of the specimen. All AE waveforms were recorded at 10 million samples per second (MSPS) to avoid aliasing of the recorded waveform using a PCI-2 data acquisition board with an analog filter between 100 and 1000 kHz, which represents the closest built-in filter available to the AE sensor range. Peak definition (PDT), hit definition (HDT), and hit lockout time (HLT) were set as 100, 500, and 500 μs, respectively, based on pencil lead break tests and verified through monotonic testing. All waveforms were further filtered in post-processing to the operating range of the sensor using a tenth-order digital Butterworth bandpass filter. The AE equipment was manufactured by Mistras Group.

An investigation of the tested specimens after failure was conducted to identify key damage mechanisms (fiber failure, fiber pull-out, delamination, and matrix cracking) present in relation also to related NDE results in the literature [57–59] to assess the evolution of these mechanisms over the testing process. A side view of a failed specimen in Fig. 2 shows significant amounts of delamination between plies both at the center of the specimen (Fig. 2(a)) and near the grip zone (Fig. 2(b)). The left side of the image in Fig. 2(b) reveals less delamination than the right-hand side which is expected as the gripping occurred just to the left of the image location. Additional mechanisms were observed by examining the specimen in different orientations. Figures 2(c) and 2(d) show the specimen from a top view at the top and bottom portions of the final failure location, respectively.

The individual damaged plies are color-coded by fiber orientation to emphasize the observed damage mechanisms. The red, green, and purple dashed lines indicate the 0-deg, 45-deg, and −45-deg plies, respectively, and show the locations where separation has occurred. The primary damage mode observed is fiber-matrix debonding. The yellow lines indicate fibers aligned with the loading direction (90-deg ply). Many of these fibers were found fractured.

In Eq. (11), E₁₁ is the longitudinal Young's modulus, E₂₂ is the transverse Young's modulus, G₁₂ is the shear modulus, and ν₂₁ and ν₁₂ are Poisson's ratios. Also, Δ = 1 − (1 − d_f)(1 − d_m)ν₁₂ν₂₁. The critical energy release rates (energy dissipation), $Gic,i=XT,XC,YT,YC$⁠, are employed to address how the state damage variables increase (from 0 to 1) after one or more damage initiation criteria are met [64].

It should be noted that only the tension mode failures (longitudinal and transverse) were considered for the maximum element size. Based on the assumed material properties (Tables 1 and 2), the maximum element size was chosen as 0.8 mm in this FEM implementation. Furthermore, to ensure that the explicit code successfully predicts the quasi-static condition, the kinetic energy of the model needs to not exceed 10% of the total internal energy [61]. Also, the quasi-static analysis should be performed without generating waves due to discontinuities in the rate of applied loading. To address this issue [68], a displacement loading with a smooth step amplitude curve was used in which the first and second derivatives are zero at the initial and final data points.

A python code was developed to read the Abaqus output file and create the computational results shown in this paper. The python code reads the required values, namely, the stress and strain in global coordinates for composite plies, longitudinal, and transverse damage variables for composite plies, and damage variable for cohesive plies. Then, average values are calculated for the entire domain. It should be noted here that the summation of the averages of the transverse, longitudinal, and cohesive damages was calculated for each strain step, and then, they were divided by their maximum value to find the normalized overall damage.

The experimental information used for the model updating approach in this investigation consists of nondestructive evaluation data which is correlated with mechanical behavior data. Specifically, representative mechanical testing data obtained by testing the specimens presented in Fig. 1(b) are shown in Fig. 4. The monotonic behavior of this material exhibits initially a linear behavior (Fig. 4(a)) during which strains increase uniformly across the gauge section, as shown between steps (i) and (ii) in the corresponding DIC results (Fig. 4(b)). Once the load value reaches point (ii), an edge-type localization was observed on the top ply which is oriented at 45 deg. After this point, multiple similar strain localizations at a 45-deg inclination were observed. These strain concentrations are directly linked to the separation of fibers from the matrix, and they could be considered a damage precursor in the evolving failure process.

In addition to the surface information obtained from the DIC results, analysis of the thermal evolution during loading revealed additional data trends which corroborate with the DIC results in Fig. 4(a). Specifically, although the minimum and maximum temperature trends do not show any significant change, the standard deviation shows a downward trend after 0.33% strain marked by the black dashed line. This trend does not correspond to any temperature hotspot; however, it could be indicative of subsurface damage which might not be capable to be depicted in the surface maps in Figs. 4(a) and 4(b). In this context, the first thermal localization is observed in Fig. 4(d) at 0.7% strain indicated by the dashed green line labeled as (1), which agrees with the corresponding strain localization observed in Fig. 4(b), while it also appears to have the same 45-deg inclination and to be at the same location as the strain hotspot. Further thermal hotspots are observed in lines labeled (2)–(3), while just prior to final failure (4) an area with a significant rise of temperate (∼ 10 deg with respect to the initial average temperature in the measurement field of view) was found, which corresponds well with the area where the final fracture occurred.

Beyond the DIC and pIRT data, representative results from the application of the AE method are shown in Fig. 5. The AE activity is represented (Figs. 5(a) and 5(b)) by two of its features, namely, the amplitude and the peak frequency of the recorded waveforms, corresponding to the specimen data in Fig. 4. Additional AE features including the cumulative absolute energy, hits, and counts are shown in Fig. 5(c). It is interesting to note that at the strain increment where AE amplitude and activity appear to increase (at 0.5% strain marked in Fig. 5(a)), the cumulative AE features reveal an inflection point that coincides with the initiation of the cumulative data trend in Fig. 5(c). In addition, the corresponding strain accumulation observed in Fig. 4(a) corresponds well with the second inflection point observed in Fig. 5(c). Hence, this collaboration between these three different NDE methods increases the confidence regarding the capability of the NDE methods chosen to monitor failure, as also previously shown by the authors and others in the case of composites.

Based on the NDE data trends presented in Figs. 4 and 5 and given the prior work by the authors [71–74], a data post-processing procedure were followed using outlier analysis to develop a more objective description of the damage initiation and evolution process instead of solely plotting NDE features as a function of loading. To achieve this, 31 AE features from both the frequency and time domains were extracted and normalized to describe the real-time recorded waveforms. The feature space was reduced based on a combination of correlation and principal component analyses. These features were then used in a Mahalanobis squared distance (MSD) outlier analysis. Specifically, the cumulative sum of Mahalanobis squared distance (CMSD) values was used as shown in Fig. 6(a), and the normalized values are shown in Fig. 6(b), which was demonstrated by the authors to have the capability to describe changes of the evolving material state due to damage [72]. This metric was used to form a baseline based on which generalized distances were computed [71]. It should be noted that although features from multiple NDE methods could be used to form such metric values, for the results in this article only AE data were used based on previous work by the authors [71,72,74]. However, the method could be expanded to include additional NDE data sets.

Interestingly, by plotting the normalized version of the CMSD for two different experiments some differences are observed in Fig. 6(b). In this context, normalizing this metric has proved to be a good methodology for extracting a degradation criterion in the range between 0 and 1, which provides quantitative damage evolution trends which are consistent with similar finite element analysis modeling damage factors [72]. The normalized cumulative Mahalanobis squared distance clearly shows that although there might be differences in terms of damage initiation among experiments, distinct data trends can be obtained related to the progressive failure process in these composites which could be leveraged in the computational model updating, as described in Sec. 4.2.

Given the experimental investigation and NDE data sets related to progressive failure, the objective of the computational approach is to create a methodology that leverages such information into FEM models with progressive failure analysis capabilities. Specifically, Figs. 7(a) and 7(b) show that Model A (with no cohesive layers) fails significantly in being capable of predicting both damage initiation (compared to the experiment in Fig. 4) and the final failure point, even if it uses ply properties that are capable to match the elastic portion of the response of the material.

To explain the discrepancy in progressive and final failure predictions in Fig. 7, it should be noted that the experimentally determined mechanical behavior of CFRPs is affected by the voids, microcracks, and other manufacturing defects shown in Figs. 1 and 2. Such material imperfections do affect the overall composite properties in ways that cannot be tracked systematically by computational approaches unless a detailed testing and characterization approach is followed in which actual microstructure and material information is produced (e.g., using detailed X-ray micro-computer tomography data sets) to form inputs to modeling in digital twin type approaches. In addition, this discrepancy can also be attributed to the fact that free-edge delaminations are not considered although the results in Fig. 4 clearly show that they are dominant localized failure initiation points. Such free-edge delaminations occur due to interlaminar stresses generated because of the Poisson ratio mismatches and the degree of mutual influence of adjacent layers during progressive failure conditions.

To assess the effect of delamination, Model B is used and corresponding results are presented in Fig. 8. For Model B, the simulation was stopped (i.e., final failure was assumed to occur) when two adjacent cohesive layers were computed to be 85% damaged. The results shown in Fig. 8 demonstrate that even if the computational results for Model B are even closer to the corresponding experimental results of Fig. 4 based on the comparison of the elastic region, the activation of the different damage mechanisms shown in Fig. 8(b) creates discrepancies between the computational and experiments defined progressive failure.

To quantify this deviation, Fig. 9 provides a measure of the total damage computed by both Model A and Model B which is additionally compared with the normalized CMSD curve shown in Fig. 6 for the experimental results of Fig. 4. Based on the results shown in Fig. 9, it can be clearly seen that although both models provide progressive failure analysis capabilities, they differ from the provided experimental data set in terms of overall damage initiation and evolution.

This discrepancy between computational modeling and experimental results in terms of damage initiation and evolution is expected for the reasons mentioned earlier, and therefore, the research question at hand relates to the possibility that exists to leverage experimentally produced NDE data sets related to damage in CFRP to train models and therefore achieve asynchronous updating, which can then be used to obtain further insights, such as progressive failure analysis through thickness, as well as to assess the relative importance of one particular damage mechanism versus another, etc. More specifically to the type of inter-and intralayer failure modeling in this investigation, there are at least five independent parameters that must be provided by the user, namely, the longitudinal tensile strength, X_T, transverse tensile strength, Y_T, the fracture energy release rate during tension in the longitudinal direction, $GXTc$⁠, the fracture energy release rate during tension in the transverse direction, $GYTc$⁠, and finally, fracture toughness in different modes for cohesive material, $Gnc$⁠, $Gsc$⁠, and $Gtc$⁠, the change of which affects strain value for which damage initiates. It should be noted here the longitudinal shear strength, S_L, changes similar to Y_T since from a microscopic point of view, they are both proportional to matrix strength [75]. It is further assumed that $Gnc$⁠, $Gsc$⁠, and $Gtc$ will change by the same amount when adjustments are made because of model updating. Although specific experiments may be made to measure individual changes of these fracture toughness properties, such effort requires elaborate test data which are not available in the present investigation.

A more careful analysis of the results in Fig. 8 shows that the activation of delamination appears to initiate both earlier and contribute more compared to the longitudinal and transverse damage modes to the overall damage in Model B, which combines the Hashin criteria and the cohesive modeling approach for delaminations. Therefore, a possible updating of the cohesive parameters is expected to have a direct influence on the strain value at which damage initiation occurs in Model B. Practically, by adjusting the fracture toughness values in an iterative way a value can be found either by trial and error or by using some version of a search algorithm to adjust this parameter and achieve the needed agreement with a given experimentally obtained damage curve In this context, Fig. 10 shows results obtained by Model B with the use of original ply properties and 85% of the fracture toughness values, $Gnc,Gsc,andGtc$⁠, compared with the results shown in Fig. 8. Based on this form of updating, it can be seen that although the computed overall stress–stain curve in Fig. 10(a) is not significantly affected, the damage initiation point in the model approaches the corresponding one in the experiment as shown in Fig. 10(b). Interestingly, though, the computed damage curve past the initiation point continues to show deviations compared to the benchmark experiment chosen, demonstrating that this type of parameter sensitivity analysis requires several adjustments to match any given experimental data set.

To match the entire damage behavior as this is visualized by the cumulative damage versus strain curve in this investigation, the transverse (Y_T) strength (and consequently also the longitudinal shear (S_L) strength for the reasons mentioned previously) is reduced by 75%, in addition to the previous reduction of the fracture toughness values. It can be seen in Fig. 11 that this change makes the computed damage–strain curve in Fig. 11(b) approach even more the corresponding experimental curve. A further reduction of the fracture toughness in transverse direction $(GYTc)$ by 50% was found to cause additional changes, as shown in Fig. 12, where it can be seen that the damage curve by the model is very close to the experimental one.

The model updating approach presented in this paper is based on certain assumptions that could be further improved in subsequent steps which, however, are beyond the scope of the investigation presented herein. Specifically, as mentioned in the introduction, digital twin approaches are currently attempted which focus on representing the actual microstructure of CFRP with more details compared to the models used in this paper. Furthermore, several computational damage approaches have been proposed to address complex issues related to the damage behavior of CFRPs including damage initiation through thickness, interactions between damage modes that affect their relative activation and contribution to the overall progressive failure, and complex damage paths both at the ply level and though the different plies and more. Furthermore, the case selected in this investigation includes only one type of loading (tension) for one type of specimen (straight edge) and using a small number of NDE methods (acoustic emission, digital image correlation, and passive infrared thermography) while also using a specific method to produce a visualization of the activation and evolution of damage through the normalized cumulative Mahalanobis square distance leveraging only features of one of the available three NDE methods (in this case, only acoustic emission data sets were used). Moreover, the modeling approach used is based on certain parameters (strength as well as fracture toughness and energy release), which were updated in this investigation in a trial and error way, although more sophisticated methods can be used, as also mentioned previously.

The authors are aware of the issues mentioned earlier; however, to the best of their understanding, the approach presented in this paper presents for the first time a method to leverage NDE data in a statistically meaningful way (outlier analysis) to produce a visualization of overall damage evolution (i.e., more than defect size visualization typical with several NDE data processing methods) in relation to progressive failure analysis that could be leveraged by modeling for CFRP. The authors had previously developed a related NDE data-driven computational approach in the case of modeling damage in elastoplastic fracture of metal alloys [72]; however, the case of CFRP presented in this investigation poses unique challenges because of the direct influence of the layup and number of damage mechanisms which contribute to the overall damage at the specimen level.

In this context, it is further argued that the attempted model updating approach leveraging NDE data sets provides input to computational models that after being updated, asynchronously in this case, they could offer further insights into the progressive failure of CFRP that are not possible to be drawn neither by the NDE data sets nor by models available in this investigation separately. To elaborate on this point, Figs. 13 and 14 provide additional information related to the progressive development of the transverse damage and delamination modes, respectively, as a function of distance from the surface (i.e., through thickness) for a given strain increment (1.13%) which is near ultimate failure (results correspond to the model used in Fig. 12). It should be noted here that based on Fig. 1(b), loading occurs along the 90-deg fiber orientation.

The results in Fig. 13 demonstrate that all plies except the 90-deg ones (ply numbers 4, 8, 9, and 13 counting from left to right; the exact fiber orientation is indicated by the blue lines in the schematic that can be found at the bottom of each FEM contour plot and in accordance with the corresponding fiber designations in Fig. 3) show the activation of the transverse damage mode (similar to what is shown in Fig. 8(b)). This particular damage mode activation agrees practically with the fact that microcracks are expected to grow parallel to the fiber direction during loading in this specific direction. Hence, since the longitudinal direction (90 deg) is stronger than the transverse direction in this case, damage does not appear in 90-deg plies at the initial stages of loading and in fact up to the 1.13% strain. Hence both qualitative and quantitatively, the use of the NDE-informed model can provide additional information on the relative activation of given damage modes which cannot be practically computed by the normalized cumulative Mahalanobis square distance curve of Fig. 6 as the NDE measurements cannot be specific to neither given damage mechanisms nor a specific ply. Consequently, the combined use of NDE data and physics-based models provides a new paradigm on how incomplete (in terms of progressive failure analysis) experimental data can be used so that information related to the role of, for example, the stacking sequence can be extracted.

Moreover, Fig. 14 shows results related to the activation of delamination in relation to the 16-ply model with the additional 15 cohesive layers defined in Fig. 3(b) (the notation “COHIJ” is used here to declare the “I” and “J” plies on top and bottom, respectively, of each cohesive layers in this model. It should be noted here that based on Fig. 8(b) that shows the relative activation of delamination with respect to the Hashin criteria based damage modes, delamination appears to initiate first and cause the ultimate failure by the complete separation between the cohesive layer and the corresponding plies. Based on this information, the computed value for damage based on delamination is equal to “1” in a select number of cohesive layers as shown in Fig. 14.

The interlaminar forces and moment can be calculated at any point along the thickness of the laminate, z₀. That means that any point through thickness of the laminate is prone to delamination and not just the points between two adjacent plies. In this study, however, we just considered delamination between two adjacent plies. A matlab code was developed to compute the in-plane stresses (σ_y and σ_xy) based on classic laminate theory [66]. Then, the interlaminar forces and moment were determined.

This variable can also be compared with the quadratic nominal stress damage initiation criterion (QUADSCRT variable in ABAQUS). The comparison is valid for low strains, i.e., when cohesive damage evolution is not significant yet. Figure 15 provides the results of this new variable; it can be seen when the variable F_del reaches its maximum (i.e., COH12, COH23, COH56, COH67, COH1011, COH1112, COH1415, and COH1516), the contour plots show larger values which agree with the cohesive damage distribution through the thickness shown in Fig. 14 and were also found to agree with the analytical predictions using the classic laminate theory. Consequently, the use of the NDE-data assists in calibrating parameters of the damage mode activation criteria as explained in relation to the results shown in Fig. 12, which provide the FEM model that can be now used to further elaborate on the location and extend of damage through the thickness of the composite material.

Moreover, and to demonstrate damage evolution as a function of applied loading consistent with the need for progressive failure analysis, Fig. 16 shows the evolution of transverse damage mode for the second ply, and Fig. 17 the evolution of delamination in COH23, further demonstrating the type of PFA information that can be assessed by the presented approach.

The research approach in this investigation demonstrates a way to perform physics-based model updating of progressive failure in fiber-reinforced composite materials leveraging nondestructive evaluation data. The approach is based on a synthesis of two typically distinct elements, namely, sensing and computing that are both capable of providing descriptions of evolving material states caused by the complex failure processes in this type of materials. The key finding in this paper is that computational models of any degree of complexity in their mathematical and numerical descriptions of a given physical process (in this case damage) might fail to provide predictions comparable to experiments because actual sensing data were not used when calibrating them. Although such calibration process is case-specific and it does provide its own set of challenges, it is practically a necessary step if one considers an unknown state of any material or structure at the time that any sensing method is applied. In this context, the results described in this paper clearly show that when nondestructive evaluation data are used in conjunction with mechanical and other testing data, then the capability of the model to be updated based on performance characteristics being monitored, increase. Another key point in the investigation presented in this article is related to the concept of initiation versus progression/evolution of the failure process. Specifically, efforts were made in this paper to capture specific states of the evolving material behavior that are associated with the point of damage initiation. This choice has particular importance as in a number of physical processes, such as failure in this case, it is the transition between damage-free and damage states that is important. Hence, the approach presented shows that the model updating approach leveraging NDE data provides advantages in modeling the evolution of the failure process in complex material systems. In future work, the authors would attempt to both perform such model updating synchronously by leveraging internet of things (IoT) hardware/software integration as well as including a statistical/probabilistic approach to perform remaining useful life estimations, especially for fatigue test cases.

The authors acknowledge the financial support received by the SBIR contract W911QX-15-C-0045 funded by the U.S. Army as well as AlphaSTAR Corporation (Long Beach, CA) for discussions related to progressive failure analysis.

There are no conflicts of interest.

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.