A computational damage model, which is driven by material, mechanical behavior, and nondestructive evaluation (NDE) data, is presented in this study. To collect material and mechanical behavior damage data, an aerospace grade precipitate-hardened aluminum alloy was mechanically loaded under monotonic conditions inside a scanning electron microscope, while acoustic and optical methods were used to track the damage accumulation process. In addition, to obtain experimental information about damage accumulation at the laboratory scale, a set of cyclic loading experiments was completed using three-point bending specimens made out of the same aluminum alloy and by employing the same nondestructive methods. The ensemble of recorded data for both cases was then used in a postprocessing scheme based on outlier analysis to form damage progression curves, which were subsequently used as custom damage laws in finite element (FE) simulations. Specifically, a plasticity model coupled with stiffness degradation triggered by the experimentally defined damage curves was used in custom subroutines. The results highlight the effect of the data-driven damage model on the simulated mechanical response of the geometries considered and provide an information workflow that is capable of coupling experiments with simulations that can be used for remaining useful life (RUL) estimations.

## Introduction

Damage is a multiscale, spatially distributed, stochastic process bridging nucleation and growth, that varies significantly among different material types. The damage process spans four distinct stages including incubation, microstructurally small initiation, physically small growth, and extensive growth processes [1,2]. In fact, it has been reported that the incubation, nucleation, and microstructurally small growth regions consume most of the material's life particularly in high cycle fatigue loading schemes [1]. During this significant portion of material life, local microstructural features, which, for example, in metals, include grain size and orientations, inclusions, voids, etc., play a dominant role which might not be detrimental; however, it may lead to conditions that favor the subsequent development of known catastrophic failure modes such as cracks. From an engineering perspective, this fact has been recently leveraged by investigating innovative sensing methods [3–5], as well as new materials design approaches [6] in order to make a transition from damage monitoring to damage precursor identification and understanding [7–10].

Damage models discussed in the literature generally focus on defining material degradation due to damage using phenomenological approaches that target the analytical description of the particular role of individual mechanisms. In several cases, a damage variable, *D*, is defined as the ratio between the damaged and the total area of a reference volume. If damage is isotropic a scalar damage variable is defined; anisotropic damage leads to the adoption of a second-order damage tensor. The damage parameter has a value of zero in the undamaged state and a value of one when final material failure has occurred [11]. In this context, Kuna and Wippler [12] used a unified Chaboche model derived from a free energy potential which was enhanced to account for void growth. A thermodynamic approach was used in Richard et al. [13] to model damage in quasi-brittle materials. Other efforts have brought a more physics-based approach to the definition of damage. For example, Liu and Zheng [14] provided a review of damage modeling techniques in composite laminates including the use of phenomenological damage tensors in continuum damage models and multiscale finite element methods that derive their damage parameters via stochastic analysis of the underlying micromechanisms. Moreover, Rinaldi et al. [15] introduced a microstructure-aware damage parameter in two-dimensional lattice models that is related to the coherence length between microcracks. For an overview of microstructure-sensitive computational models with an application to fatigue crack growth, the authors refer to McDowell and Dunne [1] and references therein.

Several characterization and experimental methods have been employed to quantitatively monitor damage from microscale samples [4,16–22] to full-scale structures [23,24]. In the case of precipitate-hardened alloys, which is the material used in this article, damage has been reported to initiate at intermetallic inclusions resulting from the alloying process which lead to early crack nucleation followed by growth [25,26]. In addition to image-based characterization, nondestructive evaluation (NDE) techniques have been proven useful in monitoring the damage processes at the material level. For example, digital image correlation (DIC) methods are currently routinely used in experimental mechanics to track full field deformations; localizations observed in such datasets have been successfully linked to ductile fracture [24,27], twinning [19,28–30], and particle fracture [4,31–33]. Additionally, acoustic emission (AE) has proven useful in detecting microscale damage such as dislocation motion and particle fracture [4,34,35], as well as macroscale damage such as crack growth [27,36–40]. Additional NDE methods reported to monitor material damage including infrared thermography [41–43], electromechanical impedance [44–46], and ultrasonic testing [47–49].

By combining experimental and simulation methods, a different class of data-driven damage models has been reported. For example Barbero et al. [50] proposed a second-order damage tensor that is formed and evolved with model parameters defined through experimental data. In addition, Kompalka et al. [51] used successive model updating to refine an initial polynomial approximation of a damage parameter used to define the damage evolution at a number of integration points in finite element method (FEM) simulations. Moreover, Zárate et al. [52] predicted fatigue crack growth from AE data in a probabilistic model that provides predictions of the stress intensity factor. More recently, Loutas et al. [53] used a machine learning protocol to incorporate AE data features into a damage model for composites, that is, agnostic to the particular damage mechanisms for the specific material used in experiments.

Given this state of the art, this paper presents a novel damage model that leverages microstructural level material information along with mechanical behavior and NDE data to define damage evolution curves. A computational workflow is then developed to use such data-defined damage information into simulations of both monotonic and cyclic behavior.

## Experimental and Computational Approach

### Material Characterization Data.

An aerospace grade precipitate-hardened aluminum alloy (Alcoa, Pittsburgh, PA), Al 7075-T6351, rolled to 2.5 mm thickness was used in the as-received condition. The primary alloying element is zinc (5.1–6.1%) with traces of silicon (0.4%), iron (0.5%), copper (1.2–2.0%), manganese (0.3%), magnesium (2.1–2.9%), and chromium (0.18–0.28%) according to literature [54] and verified by energy dispersive spectroscopy. The addition of these elements increases the yield and ultimate strength of the material; however, these same strengthening elements have been shown to be nucleation points of damage [4,55], particularly at the larger particles that form from the alloying process [2,56]. A scanning electron microscope (SEM) (FEI XL30 ESEM) equipped with an EDAX electron backscatter diffraction camera was used to obtain the microstructural information shown in Fig. 1. The data plotted in Figs. 1(a) and 1(b) reveal a nearly equiaxed grain structure with an average grain size of 60 *μ*m. The presence of a rolling texture was also observed in the pole figures in Fig. 1(c), where the [100] axis is aligned with the rolling direction in a 1.5 × 1.5 mm^{2} area.

### Experimental Setup and Mechanical Data.

Experiments were conducted using a microtensile stage (GATAN) capable of applying up to 2100 N inside an SEM microscope as shown in Figs. 2(a) and 2(b). Additionally, NDE data including full field strain and acoustic emission activity were recorded similar to previous work by the authors [4,5]. Specifically, AE monitoring was achieved inside the SEM microscope by use of a BNC (Bayonet Neill–Concelman) pass-through on the custom SEM load stage door shown in Fig. 2(a). In this investigation, two broadband PICO sensors with a frequency range from approximately 150 to 750 kHz, shown in Fig. 2(c), were used. The sensor frequency response is shown in Fig. 2(d). Signals were recorded from the sensors with a sampling rate of 10^{7} samples per second using a peak definition, hit definition, and hit lock out time of 100, 500, and 500 *μ*s, respectively. A 25 dB threshold was used for the test based on eliminating ambient signals collected while no load was applied. Figure 3 shows the specimen geometry used, which has a reduced area gage section to ensure the specimen fractures inside the monitored region, while providing sufficient space to attach the AE sensors.

The 1.3 × 2.5 × 3 mm gage section ensures that there are approximately 10,000 grains in the gage and 200 grains in the cross section which is sufficient to obtain the bulk material behavior. Prior to testing, specimens were ground and polished to a 0.05 *μ*m surface finish to minimize the effects of surface defects in the damage evolution of the material. In addition, ink toner, with speckle size of approximately 8 *μ*m, was applied to the surface to provide sufficient contrast for DIC imaging in a field of view (FOV) of 1.3 × 2.5 mm, which is indicated by the gray field in the center of the specimen in Fig. 3. The specimen was then loaded at a rate of 0.1 mm/min until failure while images were captured via live video monitoring using a four-frame average at 30 frames per second. Figure 4(a) shows the load versus measured axial strain curve along with full field deformation maps at the six locations marked in Fig. 4(a).

In order to obtain experimental information in a fatigue setting, a set of cyclic loading experiments was also completed. For these experiments, a prismatic beam of dimensions 200 × 50 × 19 mm was loaded under three-point bending tensile cyclic conditions ($R=0.1$) for 50 cycles. During the experiment, AE and DIC information was collected and a similar procedure to the monotonic case was followed.

### Acoustic Emission Data and Microstructural Evidence.

Real time recorded AE activity was passed through a frequency filter that matched the range of the sensors. Thirty-two AE features were extracted; selected AE features were correlated with the load in Fig. 5, revealing a spike of AE activity at the onset of yielding. Figure 5(a) shows that most of the recorded AE signals have amplitudes below 40 dB with the exception of the activity near the yielding point and the final failure of the specimen. Further, a clear distinction between signals with higher peak frequencies, above 400 kHz, and lower peak frequencies, below 200 kHz, is observed in Fig. 5(b). Moreover, the amount of energy within the higher frequency range (400–600 kHz), called partial power 3, is shown in Fig. 5(c) to be decaying after yielding. This observation agrees with the trends seen in Fig. 5(d), where the majority of count activity happens around the yield point where the cumulative hits and energy start to increase. Specifically, the energy appears to instantly jump at yielding, while the hits gradually increase.

Examining the specimen postmortem allows for the identification of damage mechanisms as shown previously by the authors [4,27,37] such as cracking and particle fracture. Specifically, Fig. 6 shows the presence of several damage mechanisms observed by ex situ microscopy and while the specimen was still in the load frame inside the SEM both near (i.e., within 1 mm) and further away (further than 1 mm) from the crack tip. These micromechanisms which include precipitate fracture, microcracks and slip lines are possibly related to the large amount of AE activity, particularly, the higher amplitude hits in Fig. 5(a) and the jump in AE energy in Fig. 5(d) that appear at the onset of yielding.

### Outlier Analysis.

A normalized subset of the extracted AE features was determined by using the correlation and variance within all 32 features available. A complete list of the available AE features that can be extracted in recorded AE activity and their definition is given in ASTM E1316-17a. The correlation between each feature and all others used the data from the experiment and a correlation matrix. The mean correlation value for each feature was then calculated and the features with the lowest values below a correlation value of 0.25 where selected. Choosing the least correlated features ensures that most of the information in the recorded AE signals is preserved while using the smallest number of features to parametrically represent them. In this case, five features (namely, the decay angle, reverberation frequency, initiation frequency, fast Fourier transform peak frequency, and peak frequency) were determined to be sufficient to represent the AE data.

where **X** is, in this case, a *n* × 5 matrix where each column represents each feature and each row represents a single AE signal, **S** is the covariance matrix of the feature values, and **μ** is defined as the baseline defined by noise AE signals (i.e., by using those signals recorded before loading was applied). The diagonal of the resulting matrix represents the MSD value and the normalized in a scale between 0 and 1 cumulative sum of all the points produces the damage curve shown in Fig. 7, which is plotted as a function of the total strain reached in the experiments. Interestingly, the MSD curve presents a sharp increase near the onset of the plateau region in the stress–strain curve of Fig. 7, while it continues to grow with a decreasing rate until final failure. It should be noted also that at peak stress the MSD shows a ∼60% remaining useful life (RUL) demonstrating that this parameter can provide a quantitative way to estimate it based on input data from the NDE setup. Furthermore, the fact that the MSD curve was defined as a function of measured strain allows its connection to model parameters as explained next.

The same procedure used to examine the AE and DIC data for the in situ monotonic tension case was used to process the data obtained during the three-point bending fatigue test. Figure 8 shows the raw AE data in the form of amplitude, peak frequency, cumulative hits and cumulative energy overlaid with the evolution of the maximum longitudinal strain value observed by DIC. Similar behavior occurs in this condition as for the in situ monotonic loading case where a concentration of hits with peak frequencies above 400 kHz occurring early in the specimen's life and decaying as more cycles are applied. A similar trend is observed in amplitude while the cumulative hits and absolute energy show a linear behavior with the exception of a spike in energy at the same location as the large cluster of “high”- amplitude and “high”-frequency hits occurs. Interestingly, there is no significant increase in any of the AE features plotted once the onset of strain localization is observed. Figure 8(d) shows the location of this localization and its size at the end of the test.

Similar to the results in Fig. 7, the MSD curve was obtained for the cyclic data of the three-point bending tests and the results are shown in Fig. 9. Here, the damage parameter is plotted with respect to the structure's life fraction where the value of 1 corresponds to the presence of a crack such as the one shown in Fig. 8(d).

### Computational Approach.

The experimental information obtained through the multiphysics experiments and the subsequent outlier analysis was implemented in a computational framework in order to study the damage accumulation trends in specimens of various geometries. Specifically, an isotropic hardening scheme was selected to describe the plasticity effects under monotonic and cyclic conditions in the tensile regime. This was coupled with a continuum local damage law based on stiffness degradation that characterizes void nucleation and growth, and implements the MSD information into the finite element (FE) domain. To this aim, a user subroutine (UMAT) was developed and the MSD curve previously derived was used in a user hardening (UHARD) subroutine similar to previous work by the authors [57].

The elastic strain and stress increments were also determined and all field quantities (stress components, effective plastic strain) were updated at the end of each increment.

*D*is an overall (scalar) damage variable. The stiffness and load reductions alter the local stress–strain response. The response of the scalar damage parameter

*D*with respect to the applied strain to the tested specimen is shown in Fig. 7. For applied strain up to 0.1, the behavior corresponds to the elastic regime of the stress–strain curve, and therefore, damage accumulation was not considered for this strain range. Stiffness degradation was initiated after yielding and then followed the response of Fig. 7. The damage variable evolution (Fig. 7) was then mathematically approximated as a function of accumulated strain given by the polynomial expression

Given this development, the computational procedure developed is described as a pseudocode in Fig. 10. This workflow is used to find values of stress, elastic and plastic strain, stiffness components, and damage parameter for every time increment at each integration point in the model.

*i*th increment, the hardening is evaluated as

where $\sigma yi$, $\epsilon ip$, and $\sigma yi\u22121$, $\epsilon i\u22121p$ are pairs of yield stress and plastic strain values at the current and previous increments, respectively. In the next step, the total strain at any integration point is compared to the damage initiation criterion. If the criterion is met, the current value of the damage variable *D* is calculated according to Eq. (7) and the element stiffness and load are updated using Eq. (6). If the damage criterion is not yet met, no changes are made to the stiffness or stress tensors. Finally, the material tangent stiffness matrix is calculated and returned to the finite element solver.

This numerical procedure was successfully implemented in abaqus [58] as a user subroutine, suitable for both two-dimensional and three-dimensional problems. It is noted that due to the nonlinearities introduced, the problem is either solved statically or dynamically with the aid of (artificial) damping. Two three-dimensional geometrical configurations were used to study the incorporation of this method, a dogbone specimen under simple tension and a beam under three-point bending loading. The boundary conditions for both geometries are shown in Fig. 11. The height of the dogbone specimen was 152.54 mm, its thickness equal to 3.18 mm and the reduced section width was 13.54 mm. The top and bottom ends were assumed rigid in accordance with the actual experimental conditions. Hence, the sample was fixed at the top rigid area while displacement was applied at the bottom (Fig. 11(a)). Accordingly, the beam had a total length of 200 mm, while its width and height were 50 mm and 19 mm, respectively. A half longitudinal symmetry model of the beam was used as shown in Fig. 11(b). Both models were discretized using structured meshes of linear eight-node brick elements (C3D8). The dogbone model consisted of ∼320,000 degrees-of-freedom, with an average mesh size of 0.5 mm and the beam comprised ∼240,000 degrees-of-freedom in a biased mesh of a minimum size equal to 1 mm (in the center) and 5 mm at the ends.

## Results and Discussion

### Coupon-Level Analysis Under Monotonic and Cyclic Loading.

The data-driven damage model was first tested on a dogbone specimen subjected to monotonic loading to explore the effect of damage on the bulk mechanical response. Two cases were considered; in the first, the stress–strain response was computed using only an isotropic hardening plasticity scheme, while in the second, the stiffness degradation step according to the process described in Fig. 10 was added. The results are shown in Fig. 12(a), where the “undamaged” and “damaged” designation was adopted to describe the cases without and with the damage model, respectively. Also, the damage parameter output is plotted for the “damaged” case. To obtain the data for these curves, stress and strain values were averaged in a region defined by two parallel lines set at a distance of 10 mm from the midpoint, where deformation localizations were expected to occur. It can be observed that the numerical response of the “undamaged” specimen matches exactly the experimental data, as given in Fig. 5, which was expected as the actual stress–strain curve obtained experimentally was used to define the plasticity in these models.

On the other hand, when the stiffness degradation scheme that incorporates the MSD curve of Fig. 7 was added, the global mechanical behavior was affected significantly. The response follows the undamaged material behavior up to a strain level of $\epsilon \u22484.6%$, after which the stress gradually reduces with further increase in strain. It is worth noting that at an element level, damage initiates at earlier strain stages (after yielding); however, this local effect is not manifested immediately in the macroscopic stress–strain response. The rate of change in the subsequent softening branch follows the behavior given by the MSD (see Fig. 7), i.e., the damage accumulates at a rapid rate up to $\epsilon =12%$ and then saturates until the ultimate specimen failure.

The same dogbone specimen was also subjected to tension–tension cyclic conditions (using a load ratio of $R=0.1$) for 15 cycles and the computed stress–strain responses for both undamaged and damaged cases are plotted in Fig. 12(b). The undamaged specimen follows the corresponding trend of the monotonic loading case and has generally constant drops in stress at every cycle. A softening effect is noticed after the sixth cycle at $\epsilon \u22487.5%$, which is attributed to significant necking that occurred in the field of view due to large deformations in this model. The response obtained from this simulation can be well compared with other attempts in the literature (e.g., see Ref. [59]).

As a comparison, the response of a specimen where the damage curve was added is given in the same graph. The results of this simulation follow the corresponding stress–strain curve of the undamaged specimen for the first three cycles. After an average strain of $\epsilon \u22484.8%$, a softening behavior is obtained since the stress level achieved with further increase in strain are much lower than that of the undamaged specimen. It is noted that damage initiation is observed at earlier strain stages at a finite element level, as shown in the damage parameter output. After the third loading cycle, a drop of almost constant slope is noticed. The limited capacity of the specimen to undertake load is also noted since after the fourth loading cycle the unloading branches do not coincide with the corresponding ones of the undamaged specimen, i.e., the strain range between cycles is less in this case. As a result, this simulation is concluded at a lower overall strain level. It should be also mentioned that the damage parameter shows a similar trend to the monotonic loading case. A plateau in the damage accumulation is noted at all unloading points.

Moreover, it is worth investigating the effect of the damage model application in the overall deformation mechanism of the specimen. In Fig. 13, full-field results (displacement *u _{y}*, von Mises stress, and strain

*ε*) for the undamaged and damaged specimens under monotonic loading are plotted at the same loading stage (the response is similar in the cyclic loading case). In both cases we observe localized deformation around the midpoint of the specimen's height; however, this localization is symmetric in the undamaged specimen case and antisymmetric in the damaged one. In the first case, the deformation mode corresponds to tensile necking, i.e., a reduction of the specimen's cross-sectional area in a small region. On the other hand, the use of a local damage model leads to shear deformation mode, as discussed by the authors [57]. This is mostly evident in the von Mises stress plot where a diagonal zone of high intensity across the specimen's width is formed. It should be also noted that the localization region is affected by the mesh characteristics (element type, mesh density, and discretization type).

_{yy}### Application to Three-Point Bending Modeling.

A beam under three-point bending loading conditions was investigated next to examine the differences caused by the damage model in this geometry that resembles a structural component. It is reminded that similar to the previous analysis, the experimental information in the form of stress–strain input and damage parameter evolution is applied everywhere in the structure at the element level. Results for monotonic loading conditions are presented in Fig. 14. Specifically, contour plots for the displacement *u _{y}*, von Mises stress, and strain

*ε*are given for the case of plasticity (Fig. 14(a)) and plasticity coupled with damage (Fig. 14(b)) at the same loading stage. The use of fixed colormaps emphasizes the shift in the response between the two solutions. The maximum beam deflection (at the center) and the normal strain on the

_{yy}*y*-direction are higher when the damage input is considered while the von Mises stress is lower as the load carrying capacity of the structure reduces due to damage. A quantitative comparison of the maximum values in the three plots yields a 27% increase in the vertical displacement, a 23% reduction in the von Mises stress and an 8% increase in the normal strain, when damage is added in the simulation.

Moreover, the full-field damage evolution at different loading stages is presented in Fig. 15. The damage is accumulated at the center of the beam where the load is applied in a symmetric way with respect to the beam's thickness. The maximum damage value for the four plots are 27%, 52%, 73%, and 100%. In snapshot (d), it is clear that several elements have reached ultimate damage levels. As a reference, it is worth evaluating the ratio of maximum displacement over the displacement at failure (case d) for each of the cases (a) to (c). The corresponding values are 23%, 42%, and 79%.

A beam of the same geometry was also studied under (tensile) cyclic three-point bending loading conditions both using the plasticity law only and employing the stiffness degradation approach. In this series of simulations, the stress amplitude per number of fatigue cycles information (*S*–*N* curve) was utilized as extracted from the literature [60] (Fig. 16). A static initial analysis under the maximum load of the given *S*–*N* curve was performed first to define the stress at each node of the structure. Using the stress distribution from this analysis, the fatigue life for each integration point was calculated using the *S*–*N* diagram. Then, the simulation was run at discrete number of cycles setting the number of loading cycles as a function of computational time. For any given number of cycles, the life fraction was calculated, and the associated value of the damage parameter was introduced in the analysis based on the MSD curve presented in Fig. 9 In this way, the stiffness was degraded at each integration point depending on the loading cycle the analysis was running on.

The results of this numerical investigation are shown in Fig. 17, which is organized in a similar manner to the monotonic case discussed earlier (see Fig. 14): on the top row full-field plots for the displacement *u _{y}*, von Mises stress, and strain

*ε*are given for the case of plasticity only while on the bottom row the same quantities are presented for the case of plasticity coupled with damage, at the same point of fatigue life. Under these loading conditions it is observed that all quantities are significantly increased when damage information is added in the analysis due to the reduction in stiffness. Specifically, the value of the maximum deflection is 1.75 times higher, the maximum von Mises stress 3.41 times higher, and the normal strain is increased by 4.96 times.

_{yy}The effect of this increase is further explained in Fig. 18, where the full-field damage evolution (*a*) and the maximum beam deflection for the two approaches (*b*) at different stages of fatigue life is shown. The accumulation of damage is already noticed at 1024 cycles while an ultimate failure value is reached in some finite elements at 16,384 cycles. In Fig. 18(b), it is shown that the maximum beam deflection remains constant when damage is not incorporated in the model. On the other hand, the maximum displacement increases significantly with respect to the number of loading cycles. At 4096 cycles, the values of the two simulations practically coincide, while for 16,384 the increase is 3.84%, for 65,536 cycles 16.4% and for 524,288 cycles 74.5%. Therefore, the degradation of stiffness is essential so as the outcome of the cyclic loading simulation becomes realistic.

## Concluding Remarks

A data-driven continuum level local damage model based on multiphysics experimental information was proposed in this work. The material was tested inside an SEM and its deformation was monitored using nondestructive acoustic (AE) and optical metrology (DIC) techniques. Data collection was followed by outlier analysis and a microstructurally sensitive damage law was generated. This experimental information was incorporated in numerical models based on custom subroutines implemented in a finite element method that combine a plasticity model with a stiffness degradation approach. Monotonic and cyclic numerical simulations of different geometries were completed in order to study the damage evolution process. The results showed the influence of a microstructure-based damage model in the macroscopic material response. The behavior deviated significantly when compared to analysis based on a macroscopic plasticity law. This effect was emphasized in high-cycle fatigue simulations where the use of a damage model is essential to predict a realistic response. In conclusion, this work contributes to the development of physics-based data-driven damage models that aim to enhance our understanding on material failure and identification of damage precursors.

## Acknowledgment

The results reported were obtained by using computational resources supported by Drexel's University Research Computing Facility and AlphaSTAR Corporation GENOA Software.

## Funding Data

Army SBIR Health Conscious Structures for Zero-Maintainence Rotorcraft Platforms Phase 1 (Contract No. W911QX-17-C-0015).