Multiphysics Based Simulation of Damage Progression in Composites

The long-term properties of continuous fiber reinforced composite materials are increasingly important as applications in airplanes, cars, and other safety critical structures are growing rapidly. Although a clear understanding has been established for initiation, growth and accumulation of damage, it is still unclear when and how the interactions of these local events lead to the development of a “critical” fracture path resulting in a sudden change of global properties and possible rupture. In the present paper, we simulate damage development in a neat polymeric resin using X-FEM analysis, and conduct concomitant dielectric response analysis with a COMSOL simulation model to study the collective defect structure as it develops in a model system. Our studies reveal inflection points in the predicted global dielectric response vs. strain that are related to changes in local damage growth rates and modes that clearly indicate impending fracture and capture the progressive change in material state.


Introduction
Composites are designed to accommodate distributed damage development consisting of various types of defects and even multiple breaks in the same reinforcing fiber.However, local changes in the material state also have significant effects on the prognosis of the future behavior of a composite material system.
Reifsnider et al. [1] [2] have described how the crack coupling process that takes place creates a fracture plane through the thickness of a fibrous composite material leading to rupture of the specimen at the global level through the formation of a "critical path".Under various applied field conditions, heterogeneous materials degrade progressively.To evaluate such changes in material state Electrochemical Impedance Spectroscopy (EIS) and Broadband Dielectric Spectroscopy (BbDS) are robust tools to extract the material-level information, including the morphology changes caused by micro-defect generation and the orientation and connectivity of those defects.Various models are well established to simulate the dielectric response of heterogeneous materials such as effective medium theories, bounding methods, percolation theory, random walk, hopping models, Fourier expansion, finite difference time domain methods (FDTD), RC network methods, and so on [3].A framework developed by Reifsnider et al. [4] has postulated that these techniques capture the change in material state and are sensitive to damage at the ending stages of life with the capability to detect imminent fracture.
Vadlamudi et al. [5] [6] [7] [8] [9] considered the variation of dielectric properties with large deformations to study the influence of change in shape of the specimen on the dielectric study.Haider [10], Raihan et al. [11] [12] [13], Bekas et al. [14] have demonstrated the capability of this technique to characterize and detect damage in composite material systems.Most recently, Reifsnider et al. [15] have postulated the concept of a "critical" path generated by coalescence of various defects and mechanisms.

Dielectric Study of Material Systems
Heterogeneous material systems are strongly dielectric materials because of their ability to store charge at various interfaces between constituent elements, voids, defects, cracks, etc.The framework postulated in [4] uses this behavior to capture the change in material state as shown in Figure 1 damage progresses in the material system we see more interfaces being created and hence more charge accumulation and higher permittivity.At high frequencies in the range of MHz to GHz we can observe orientational or ionic polarization due to the presence of moisture or other polar molecules in the form of impurities.Also polarization density can be estimated using the relation ( ) where P is the polarization density (C/mm 2 ), E is the vector electric field intensity (V/mm) and 0 ε and r ε are the dielectric Permittivity of vacuum and the modeled material, respectively.This implies that the polarization relies heavily on the dielectric permittivity of the material and on the impurities that are accumulated in the crack openings.Hence the working environment of the system is crucial and its effect has been shown previously [4] [11].
To calculate the dielectric variables we solve the Poisson's equation for the electric potential from the continuity equation as shown in Equations ( 2)-( 5) [16].
( ) where J is the Current Density (A/mm 2 ), ( ) is the current per unit volume (A/mm 3 ), D is the Electric Displacement Field, (C/mm 2 ), σ is the dielectric Conductivity (S/mm), and φ is the Electric Potential (V).Once, the potential distribution is solved, we then calculate the current from which parameters such as impedance and capacitance are calculated using the following relations [16].
where Z is the dielectric Impedance of the system (Ω), Y is the dielectric Admittance of the system (S), C is the dielectric Capacitance of the system (F), ω is the angular frequency (rad), A is the area between the electrodes (mm 2 ), d is the distance between the electrodes (mm).

Experimental Results and Motivation
Fazzino et al. [4] performed end loaded bending tests on plain weave glass/epoxy composite coupons to introduce damage as shown in Figure 2 Using in-situ tensile testing, Raihan et al. [11] induced damage in an off-axis plain weave E-glass fabric combined with epoxy resin matrix and performed dielectric studies on the samples to get real time response with damage as shown in Figure 3(a).Also, a novel edge replication technique was used to observe the damage progression with increasing load as shown in Figure 3 So, these experimental evidences have given clear indication that the degradation created by an applied mechanical field is captured by the material response to an applied vector electric AC field, and that the measured response also captures different stages in degradation mechanisms that are not generally obtained from a stress-strain curve.This also suggests that mechanical and electric fields can be coupled to get additional information about the "beginning of the end of life"; for example, Reifsnider et al. have developed the concept of "Heterogeneous" fracture mechanics that showed the direct association of the dielectric response of material systems with the mechanical loading on the systems [15] [18].
In this paper we induce damage in the specimen and capture the change in material state using the concept of "Heterogeneous" fracture mechanics and predict the "beginning" of end of life.

Current Framework
The current research presents a combined mechanical-dielectric study of the growth of a defect through the thickness of a coupon sample as a control experiment to determine how the observable changes in dielectric response define the details of the actual growth and how they can be used to predict the onset of a "critical path" or other abrupt change in observable behavior.
In the present study, a 2D model of a neat polymeric (epoxy) resin coupon with a hole (to simulate the cross section of a coupon specimen for example) is considered and is modeled in ABAQUS TM .The coupon is deformed in Mode-I using displacement control until failure.Damage is modeled using X-FEM based on the Virtual Crack Closure Technique (VCCT) criteria in ABAQUS TM which propagates the crack when the strain energy release rate exceeds the critical strain energy release rate (G c ) [19].There are many other advanced schemes for crack propagation e.g., [20] [21], but the intention of current research is to obtain a crack path driven by consistent mechanics and to perform dielectric studies of that path development in a model material.
After mechanical loading, the deformed mesh was exported into COMSOL TM at various increments of damage development in the coupon using the a/w ratio to define each analysis event.A vector electric field was applied through the thickness of the coupon perpendicular to loading direction as shown in Figure 4 and global dielectric properties were plotted at various a/w ratios.The framework for the analysis is shown in Figure 4. with a Poisson's ratio (ν = 0.35) and stiffness (E = 3450 MPa).To simulate damage, the damage initiation strength was assumed to be about 60 MPa and the critical energy release rate of (G Ic = 0.177 N/mm) was used [20].A mechanical 1% strain was applied along the right edge of the coupon using an initial minimum displacement of 0.015 mm and a maximum displacement of 0.075 mm per increment.However, the applied displacement per increment is automatically controlled by ABAQUS TM and varies it as per the convergence of the analysis.
The displacements along the left edge were constrained in the x direction and at the midpoint on the left edge, displacement in the y direction was constrained like a pin support as shown in Figure 5.The analysis was performed using Plane Strain assumptions and the mesh was generated using 4 node Continuum Plane Strain (CPE4) quadrilateral elements.The hole in the coupon acted as the stress concentrator that drove the damage initiation mechanism.

Dielectric Study of the Open Hole Tension Model
The deformed mesh was exported into COMSOL TM at various increments of damage development.The Electric Currents Shell (ECS) interface of COMSOL TM was used.The dielectric properties of epoxy are widely available in the literature; for this study the coupon had a dielectric permittivity ( r ε = 3.15) and dielectric conductivity (σ = 1.64 E-17 S/mm) [22].The imported undeformed model with an electric potential φ = 1V is applied on the top edge, the bottom edge is grounded as shown in Figure 6.Remaining edges are insulated prescribed by the condition shown in Equation ( 10) [16].

Results and Discussion
In the deformation process at around 0.39% strain the coupon was completely  absorb the moisture from the environment or lose moisture to environment [23] and in this diffusion process the probability of moisture settling in the newly formed defect surfaces is high, thus in this analysis the authors assume that during that phase of damage initiation moisture will migrate to the crack boundaries.For moisture the values of Dielectric permittivity ( r ε = 20) and dielectric conductivity (σ = 5 E-7 S/mm) are well documented in the literature [24].Figure 8 shows an example (a/w ratio = 0.91) of the deformed mesh exported from ABAQUS TM into COMSOL TM (zoomed near the vicinity of the hole).
During the dielectric study, impedance of the material system as a function of frequency for various a/w ratios was plotted using Bode Plots in Figure 9(a).
From impedance plots, it is observed that before damage initiation the response looks like a parallel plate capacitor, which is expected, and as damage initiates and progresses the resistance to charge displacement decreases (due to the heterogeneity in crack regions) and finally when the critical fracture path forms the impedance of the system further decreases indicating that conductive behavior is occurring.To better understand the degradation of properties of a material system it is recommended to have a reference and to normalize the properties with respect to that reference, which then indicates the relative magnitude or the intensity of property change with damage progression.Hence in the current study, the dielectric properties of the material system are initially calculated at the undeformed stage of the coupon and these are used to normalize the response at different stages of life.A sample of such a procedure is shown below.
Consider, the impedance at the initial undeformed state as Z 0 , the normalized impedance will be calculated as For the initial undeformed state norm Z will be 1; the variation in normalized impedance with damage progression at various excitation frequencies can be seen in Figure 9(b).
From Figure 9(b), it is observed that at low frequencies the normalized impedance changes definitively with the property degradation.As we move towards higher frequencies, the normalized impedance curve starts shifting up as seen from 1000 Hz to 1 MHz and becomes insensitive to any material state change.This can be validated by the phase angle plots.The phase angle (δ) of a system is calculated by the arctangent of the ratio of Im(Z) to Re(Z).If phase angle is negative implying voltage lags current, the system acts as a capacitor, where as if the phase angle is positive this implies an inductive behavior.For the current material system, the phase angle is negative implying the capacitive nature.Figure 10(a) shows the variation of phase angle with damage progression at different frequency spectra.
This implies in the absence of applied electric fields above a particular frequency the dipoles are no longer influenced, i.e. they have relaxed.This characteristic time is known as the relaxation time which can be estimated using the From Figure 10(b), it is observed that up to a/w = 0.49 there is one relaxation frequency in the range of 100 Hz to 1KHz; after a/w = 0.49 there are two relaxation frequencies in the range of 100 Hz to 10K Hz.The former relaxation corresponds to interfacial polarization and the latter corresponds to orientational polarization due to dipolar moisture molecules.It is to be noted that the orientational relaxation starts after a certain damage growth implying the increasing moisture content corresponding to increasing damage growth.Although the Cole-Cole scheme is not generally an accurate way of representing the data for dielectric materials, it can be used to make a rough estimate of the minimum time required for dipole relaxation.The relaxation frequency can be estimated using the relations shown in Equations ( 12)-( 13).

( ) ( )
( ) ( ) where ε * is the complex permittivity of the system, ε ′ is the real part of com- plex permittivity of the system, ε ∞ is the permittivity of the system at high frequency, s ε is the permittivity of the system at low frequency, 0 τ is the debye relaxation time.
Using the Cole-Cole plot for non-dielectric data, one can observe the plot looks like a semicircle with center on the ε ′ axis with the peak at the center of the semicircle, average of ε ∞ and s ε , and where 0 ωτ is 1 from which 0 τ can be estimated [25].To evaluate the relaxation frequency of the system let us consider the permittivity data at a/w = 0.49 shown in Table 1.
For the data in Table 1, the relaxation frequency was estimated to be close to 350 Hz.This is a rough estimate for the given data since the Cole-Cole plot is not a good fit for the dielectric data recorded.Hence, at higher frequencies one can observe a shift in the normalized impedance curve and relaxation in the normalized permittivity curve.This implies that in order to predict material state change with this technique it is required to operate in the low frequency spectra.
Table 1.Permittivity of the system at a/w = 0.49.
Frequency, Hz ε', Re(ε*) ε'', Im(ε*)   Since, this is a homogeneous model system with no reinforcing phases, it can be observed that once damage initiates there will be continuous growth of the crack up to fracture as there were no crack arresting reinforcements or energy absorbing mechanisms.The same was depicted by the Force-Displacement data shown in Figure 12.Typically, this response provides no warning or indication at the onset of damage initiation or before failure, and once again it is to be noted that in this system the damage initiation is the crucial stage that one would like to predict. Figure 12 shows a plot of the Force-Displacement and Capacitance-Displacement response of the system just before fracture at 1 mHz frequency.From Figure 12 it is clear that the dielectric response provides an indication of the onset of damage initiation which is the beginning of the end of life for the current material system while the mechanical response provides no warning of the onset of damage.
To find indicators of unsteady crack growth it is important to study the rate at which the dielectric variables are changing with damage.The variation in dielectric permittivity from Figure 10(b) (page 9) depicts the rate at which the normalized permittivity varies with damage progression.A sudden increase at damage initiation; followed by gradual increase up to 35% life; followed by steady acceleration from 35% to 63%; and then unsteady acceleration till 90% and sudden drop at failure is observed.The same trend was observed in the crack growth rate with displacement as shown in Figure 13(a).
From Figure 13(a), it is observed that up to 63% of life the crack growth rate was steady, and from there on the growth became unsteady and reached its end at a/w = 1.For further understanding, when we plotted the growth rate of normalized crack length and rate of change of normalized capacitance with applied displacement; both of the responses have an identical behavior as shown in Figure 13(b).The data for the plot are given in Table 2.In Figure 13(b), we have classified the regions as three stages, based on the crack growth rate and rate of change of normalized capacitance with displacement.It is clear that an increase in the normalized dielectric capacitance under the influence of a vector Electric field corresponds to the damage growth induced by a mechanical displacement field, and that an inflection point was determined at 63% of life which is the precursor for the unstable crack growth.

Conclusions and Future Work
A homogenous model material system was considered to study a coupled Electrical-Mechanical response during damage development.Mode-I loading was applied on an epoxy resin and crack propagation was modeled using an XFEM technique with a VCCT criterion in ABAQUS TM .The deformed model was exported into COMSOL TM at various increments of damage in order to perform the dielectric study by applying a vector electric field normal to the loading di- Figure 1(b) [9].At low frequencies in the range of mHz one can capture the data related to the charge trapped around material interfaces because of interfacial polarization, As (a), and performed BbDS analysis to conduct a dielectric study of the samples from the be-V.Vadlamudi et al. ginning to end of life.Fazzino's data clearly showed the change in material state at end of life as shown in Figure 2(b).Majumdar et al. [17] performed 3D X-ray microscopy on those samples to study the development of fracture paths through the thickness, as shown in Figure 2(c) and Figure 2(d).
(c).The edge replicas were able to show different mechanisms of initiation, growth, and coupling as shown by the sequence in Figure 3(b).

Figure 3 .
Figure 3. (a) In situ Tensile testing setup; (b) Dielectric study response; (c) Edge Replication images of different stages of damage progression.

Figure 4 .
Figure 4. Framework of current study

V.Figure 5 .
Figure 5. Framework of current study.

Figure 6 .Figure 7 .
Figure 6.2D Model of coupon for dielectric response with boundary conditions.

Figure 10 .
Figure 10.(a) Variation of phase angle with damage progression at different frequencies; (b) Variation of normalized Permittivity with damage progression.

From Figure 9 &
3.535794991 ( ε ∞ ) 5.97E-04 V. Vadlamudi et al.When the material completely fractures, the normalized impedance plots clearly showed a conductive behavior as in Figure9(b) (page 8); this implies that there is a change in state of the material and hence a change in behavior of the material.This is validated by the phase angle plots which clearly show a shift in phase at fracture in the low frequency spectra.When the phase angle becomes zero, this implies that both the voltage and current are at the same phase and the material doesn't exhibit a capacitive nature, which was observed in normalized permittivity plots shown in Figure10(b) (page 9).The important outcome of this study is that the frequency spectra dictates what are predicted in the material behavior, and hence for the remaining dielectric plots the authors consider only low frequency spectrum behavior to capture the material state changes.Figure10(page 8 & page 9 respectively) it was observed that the change in material state (property degradation) was captured in this study.However, to convince ourselves and to better understand the physics, we reran the analysis with same properties as the neat resin in the open crack region and with air in the open crack regions.The dielectric properties of air are dielectric permittivity ( r ε = 1.0005) and dielectric conductivity (σ = 1 E-18 S/mm)[26].When the same properties as neat resin were used, it was observed that there was no significant variation in the predicted dielectric response, and whatever variations were observed were due to changes in shape (geometry) as shown in Figure11.With air in the cracked regions, the global permittivity initially increased and then started decreasing as shown in Figure11.This makes sense because the presence of a poor conductor in the cracked region increases the resistance to charge displacement leading to higher impedance and lower capacitance.Hence, the effect of working environment dictates the details of the change in global response as observed in[4] [11].

Figure 11 .
Figure 11.variation of global permittivity with various materials in crack regions.

Table 2
. Rate of norm.Crack growth (norm a) & rate of change of norm.Capacitance (norm C). a/w ratio Displacement (u), mm d(norm a)/du d(norm C)/du