Predicting the damage development in epoxy resins using an anisotropic damage model

Funding information Deutsche Forschungsgemeinschaft, Grant/ Award Number: 404502442 Abstract For fiber-reinforced plastics, the strain-rate dependent response is governed by the matrix behavior. In this work, the Goldberg model is considered for the epoxy matrix constitutive material model. Moreover, the strain-rate dependency is achieved by direct influence on the elastic modulus, the inelastic strain, and the material strain to failure. In addition, an anisotropic damage response is implemented and extended through a strain-rate dependent definition. Since the constitutive model relies on nonphysical parameters, a parameter study is further performed. Additional numerical investigations using a micro-mechanical model are performed. Tension and shear loading conditions are evaluated and the influence of different strain rates is explored. Furthermore, the implemented anisotropic damage model is compared and discussed against an isotropic damage model.


| INTRODUCTION
In the event of high strain-rate deformation or the socalled crash condition, fiber-reinforced plastic's (FRP) behavior deviates strongly in comparison to the quasistatic material response. This variation is primarily explained due to the damage spread and damage mechanism within the material. For this reason, the development of a suitable strain-rate dependent damage material model is required. It is well-known that within the FRP, the matrix is responsible for the overall strainrate dependency behavior. In particular, for epoxy materials, the Goldberg model has been extensively used and extended by details such as temperature dependence. [1][2][3] On the other hand, continuum mechanic approaches have been often used for damage modelling. For instance, Salas et al developed a gradual material stiffness degradation model. [4] Other approaches consider as well the damage direction as a feature in the so-called anisotropic damage model. This type of model finds mostly a direct application in concrete technology. To the authors best knowledge, only Pineda et al [5] have proposed an anisotropic damage model for FRP with epoxy matrix.
Within the scope of this work, an already existing Goldberg model for epoxy resins is extended by mean of an anisotropic damage feature. Subsequently, the implemented model is compared against the isotropic damage model proposed by Salas et al [4] under tensile and shear loading scenarios.

| MATERIAL MODEL
The key aspect of the Goldberg model is the description of the material's isotropic inelastic strain behavior as a function of the applied strain rate. [6] In its foundation, the Goldberg model shares the viscoplastic deformation of heated metals described by Bodner. [7] Since the material properties of plastics and metals differ in many respects, it seems contradictory to adopt this consideration. However, exist several properties in which plastics and heated metals resemble each other. For instance, stress saturation at high elongations is similar in both materials. [8] Furthermore, the Bodner model is supplemented so that effects observed in plastics, but not in metals, are mapped by the model. These effects include the influence of hydrostatic pressure on plastic elongation, and so on. [1] The Goldberg model has been extensively proved as a strain-rate dependent model for epoxy resin. [3,4,9,10] Within the model, the total strain rate is divided into an elastic part _ ε e ij and an inelastic part _ ε i ij . For the elastic part, the strain rate is calculated from the elastic modulus by mean of Hooke's law, whereas the inelastic part is obtained according to the following equation: Here, D 0 corresponds to the maximum inelastic strain rate, n is the degree of strain-rate dependence of the material, s ij is the deviatoric stress components, J 2 is the deviatory stress tensor, δ ij is the Kronecker delta, Z is a material constant that describes the internal molecular resistance in the material, and α is a material constant that determines the influence of the principal stresses. [1] Furthermore, σ e corresponds to the effective stress in the material, which is based on the Drucker-Prager criterion. Thus, σ e is calculated using Equation (2) where σ kk stands for the sum of the principal stresses.
Further information about the Goldberg model and the estimation of the parameters, such as Z and α, can be found in the literature. [6,11] Additionally, several modification for the Goldberg model has been proposed. For instance, Zheng [3] suggests a strain-rate dependent Young's modulus described by the following equation: Here, C corresponds to a material constant, _ ε 0 is the effective reference strain rate, _ ε is the current effective strain rate, and E 0 describes the elastic modulus measured at the effective reference strain rate. Additionally, Zheng [3] recommends a value of _ ε 0 = 1 1 s , since the elastic modulus for the evaluated epoxy resin can be assumed to be strain-rate independent. Thus, the effective strain rate is calculated as follows: with mean strain rate equal to: Gerlach et al [9] extended the linear approach of Zheng, [3] hereinafter referred to as E lin(t) , by a nonlinear portion E add(t) : Here, E diff(t) is the maximum modelable increase of the modulus of elasticity, while the material constants C 2 and C 3 are experimentally determined. Then, the elastic modulus is calculated as: Moreover, Table 1 summarizes the material parameters implemented in this study.

| Damage model
Damage mechanics describes the initiation of the damage within the material. Thus, for a load condition, the damage is initiated for a predefined failure criterion. For instance, a stress threshold. Then, if the damage is initiated, fracture mechanics models can be used to check whether the damage is spreading. [12] In the implemented model, the damage initiates after exceeding the strain at damage initiation ε f defined as: Here, ε in is the failure strain in the quasi-static loading case and ε diff is the difference of the strain in failure between the quasi-static and the high strain rate condition. _ ε is the acting strain rate. The strain rate is not further specified by Gerlach et al.
Salas et al [4] described a model in which the stress in an element is reduced by an increasing degradation factor Ω when the elongation of strength is exceeded. Ω is calculated from ε f , the current strain ε and the material constant m, which indicates the material brittleness. Complete degradation can lead to numerical instabilities, and therefore a maximum damage factor Ω max is introduced. [10] An experimental parameter determination is not described by Salas et al [4] and therefore m has to be adapted to results.
In the context of an isotropic damage model, Ω leads to a degradation of the stiffness matrix C in all directions. Therefore, all entries of C are multiplied by (1 − Ω) and the following relationship between stresses and strains is obtained. Notice that Ω corresponds to a scalar value: In the case of an anisotropic damage model, the stiffness matrix is not degraded in all entries. Anisotropic damage models were originally developed for the simulation of concrete and are still used successfully. [13][14][15] In the field of epoxy resins, the first successful attempts have already been made to implement an anisotropic damage model. [5,16] The difference compared to other models relies on the damage differentiation. Thus, when the damage criterion is satisfied, the material stiffness weakens is initiated. This degradation effect takes place exclusively perpendicular to the largest principal stress. The affected element is therefore rotated from the global x − y coordinate system into a local n − s coordinate system, where n points in the direction of the largest principal stress. The stiffness matrix in the local coordinate system C ns is then weakened in entries C 11 and C 33 to simulate the material failure. The simplest approach (cf. in ref. [11]) is to set the corresponding entries to zero and thus no stress transfer in the main stress direction is allowed. Also, a Poisson's ratio ν = 0 is assigned. [7] C ns = 0 0 0 Compared to experimental data, the complete failure does not correspond to the real material behavior. Besides, the sudden stress drop leads to numerical problems during simulation. [13] The shear retention factor β and the normal reduction factor μ are used to solve the problems. The shear retention factor represents the roughness of the fracture surfaces over which shear forces can still be transmitted despite the fracture. In most cases it is chosen as a constant, for example, μ = 0.2, but it can also be formulated as strain-dependent parameter. [13] On the other hand, the normal reduction factor μ simulates the increasing damage of the material and the gradually decreasing tensile strength, also called "tension-softening." μ is usually formulated as a function of elongation. In some publications, [13] the influence of ν is introduced again. Thus, the resulting stiffness matrix is shown in Equation (12). [7] In the scope of this work, the entries C 11 , C 55 , and C 66 are modified in the local coordinate system n − s. μ and β are replaced by (1 − Ω): Ongoing from the stiffness matrix, which is modified by the degradation factor Ω in the entries C 11 , C 55 , and C 66 the stresses for the next time step t i + 1 are then calculated using Equation (14): Additionally, Figure 1 shows the flowchart for the anisotropic model implementation in Abaqus. The userdefined material written as a VUMAT receives from Abaqus the current stress σ t i , current strain ε t i and the current strain increment Δε t i at the current time step t i . If ε t i of an element exceeds ε f , the value of Ω t i is calculated. To prevent healing of an element the evolution of Ω is restricted to be monotonically rising. With Ω t i being determined the maximum principal stress is calculated by using an Abaqus function. Based on the strain eigenvector e i, j , the transformation matrix T is formulated. Thus, Δε t i is transformed into the local coordinate system, while the product of C n,s : T − 1 : Δε t i ½ corresponded to the inverse transformation (eg, return into the global coordinate system). Then, the stress Δσ t i is computed. In the next step, the stress field for the new time increment is calculated using Equation (14).

| COMPARISON TO EXPERIMENTAL RESULTS
To simulate the tensile and compression tests, a cube model with a side length of 20 μm is implemented. For the calculation, the commercial finite element software Abaqus is used. The model consists of a single eight-node linear brick element with reduced integration (C3D8R). The top and bottom surface is shifted linearly by a displacement of 2 μm each. The degree of freedom of the lower surface is restricted to translation normal to the surface. The surfaces are pushed apart for the tensile tests and together for the compression tests. For the investigation of different strain rates, the time in which the displacement takes place is varied.
Tensile tests at strain rates of 0.001, 10 and 3800 seconds −1 are executed. Figure 2 shows the stress-strain curves of the tensile tests using the anisotropic damage model and the experimentally determined data from Gerlach et al. [9] Since the material behavior under tension and shear is exclusively investigated, only the strain rates 0.001, 10, and 3800 in tension are shown. For the simulation of the compression case, the material failed under the pressure of −122, −153, and −204 MPa for the strain rates 0.01, 10, and 3900 seconds −1 , respectively. In the experimental data, no failure is observed at the considered strain range (0 to −0.25). In contrast to the experiment, higher values for the stress are observed in the simulation. It is observed that this difference is proportional to the increase in the strain rate. Furthermore, the stress at failure in compression in the simulation at the strain rates of 0.01, 10, or 3900 seconds −1 is between 3.7 and 23% higher F I G U R E 1 Flowchart for the anisotropic damage model For the tensile tests, a difference between the lower (0.001 and 10 seconds −1 ) and higher (3800 seconds −1 ) strain rate is observed. At the lower strain rates, the stresses of the simulation up to the flattening of the curves are between 15 and 24% higher than those of the experiments.
Overall, the simulated data using in the anisotropic damage model are largely similar to the compression and tensile experimental data. Particularly in the tensile range up to strain rates of 110 seconds −1 , largely identical observations are made. In the compression range and at tensile strain rates from 1700, deviations between the simulation and the experiments are encountered. When comparing the results of the anisotropic damage model with the already validated model of Gerlach et al large similarities are observed. The strain rate-dependent tensile behavior can thus be regarded as verified. Due to the simple stress state, no differences between isotropic and anisotropic damage model is observed.
. The findings of Naik et al [17] showed that the shear properties of the evaluated resin depend on the applied shear strain rate. Therefore, further investigation for the understanding of the strain rate dependency in shear loading conditions should be addressed.

| INVESTIGATION ON THE NONPHYSICAL PARAMETERS OF THE MODEL
A parameter study is performed on the implemented anisotropic damage model from section 3. In total, the model contains 17 parameters (cf. Table 1). For better clarity, a defined selection of these parameters is examined in more detail in this section. The selection is based on how much information is already known about the parameter and its influence on material behavior. Table 2 shows five preselected parameters for closer examination. C 1 , C 2 , and C 3 are study since they are described by Gerlach et al [9] as nonphysical parameters. Therefore, their exact influence on the model is further evaluated. The parameter m can be described as brittleness of the degradation behavior while Salas et al [4] do not explain any significance. Due to this unspecific definition, the parameter m is selected for further investigation. Moreover, ε diff is described by Gerlach et al as the difference in elongation at break between quasi-static and crash loading. [9] However, the strain rate corresponding to the crash load is not specified in detail, which is why the parameter is examined in more detail.
To determine the maximum and minimum values of the parameter study, reference to literature values is made, if the case that information is available. For  [4] Salas et al do not publish the parameters used. A minimum or maximum value of 0.3 or 100 is defined for the parameter study. Only Gerlach et al use ε diff in their models, whereby both use ε diff = − 0.05. It is observed that in the study by Gerlach et al [9] this corresponds to the difference between the elongation at break of the highest strain rate and the elongation at break of the quasi-static strain rate. If this method is applied to the data of other publications, different values result for ε diff . Zheng [3] uses values between [−0.06, −0.035], Goldberg et al [1,11,18] values [−0.06, −0.04].
The parameter study is performed using the previous described cube-shaped model, which has already been T A B L E 2

Parameters
Unit Interval F I G U R E 2 Comparison of the simulation results using the anisotropic damage model with experimental data from Gerlach et al [9] used to verify the tensile behavior. A strain rate of 30 s −1 is evaluated. A total of 10 simulations are carried out, whereby the minimum or maximum value is used for each parameter, while the mean value from the respective minimum and the maximum value is used for the remaining four parameters. Figures 3 and 4 show the stress-strain curves for the parameter variations. It is observed that with an increase from C 1 from 0.01 to 0.3, the elastic modulus increment by 260% and thus the maximum stress of 115 MPa is achieved for a 45% lower elongation. The damage, recognizable by the stress drop, increases from C 1 from 0.01 to 0.03 later. By increasing C 2 from −1 Á 10 −8 to −1 Á 10 −12 , the modulus of elasticity is reduced by 2.5%. Otherwise, the stress-strain curve does not change. When determining the maximum and minimum values, it was not possible to fall back on literature values and the determination was arbitrary. With an increase from C 3 from 0.1 to 10, the modulus of elasticity increases by 48%, whereby the maximum stress of 115 MPa is achieved with a 20% lower elongation. The damage occurs at 46% lower elongation. The increase of m from 0.3 to 100 shows no change on the stress-strain diagram until the damage occurs. In the event of damage and for m = 100, the stress decreases 1600% faster. When increasing ε diff from −0.06 to −0.03, no change in the tensile and damage behavior is observed. Overall, it can be said that a variation within the given interval of the nonphysical parameter does not lead to significant changes in the stress-strain curve except for the parameters m and C 3 . Especially the variation of m shows that with adapting the parameter, the damage behavior of the material model can be adapted. Hence, m leads to the same degradation in the isotropic model and the anisotropic model the influence is not further investigated on the micro model.

| NUMERICAL INVESTIGATION ON THE MICRO MODEL
To understand the differences between the two damage propagation methods, simulations on a micro model are carried out. The model (cf. Figure 5) consists out of five randomly placed fibers. The model has a size of 19 × 19 μm 2 . The fiber has a diameter of 6.9 μm and the fiber volume content is 42%. The fiber is modeled as a F I G U R E 3 Stress-strain curve for the parameters C 1 , C 2 , and In two different load cases, the micro model is loaded in either tension or shear. The used model has periodic boundary conditions. A constant strain rate _ ε t is applied to both sides of the model into y-direction for the tension loading. For the shear loading, all four reference points are used to achieve shear and are subjected to _ ε s . As a result, the normalized displacement and reaction force of the reference point is going to be discussed. Furthermore, the influence of the strain rate is discussed.
Since no regularization method is implemented, no constant energy dissipation in the model variates is ensured. Thus, different mesh sizes are unpractical for comparison and therefore the influence of mesh size is not further investigated in this article. Nevertheless, further research in the field of regularization methods has to be done. One approach considers the fracture toughness to ensure that the energy dissipated by the material model is independent of the mesh size. [19] Hollmann and Hahn [20] found that the fracture toughness of the evaluated epoxy resin depends on the applied strain rate. Hence, an implementation with a strain-rate dependent fracture toughness can be further considered. Figure 6 shows the results of the simulations carried out for the load case tension. The force is normalized by the maximal force occurred in the graph. The displacement is normalized by the maximal displacement the model is subjected to. The force and displacement are measured at the reference point. Concerning the increasing strain rate, an increment in the slope of the forcestrain curve is observed. This effect can be explained due to the strain-rate dependency of the stiffness matrix. Furthermore, the maximum force increases as a function of the strain rate. This effect is grounded in the increasing strength of the matrix but more importantly, it is grounded in the increasing displacement to failure u d (cf. The stiffening of the matrix establishes a stress redistribution in the micro model and hence the resulting strains in the matrix are lower compared to a model subjected to a lower strain rate at the same displacement. Resulting from this effect, the failure strain ε f of the matrix is reached at a higher total displacement of the micro model.
Besides the strain rate effect, the damage model chosen has an impact on damage development of the micro model. Figure 6 shows the difference between the two models at different strain rates. In general, the damage develops faster in an anisotropic damage model. The anisotropic damage model is degraded only in the direction of the principal stress. This leads to a pronounced damage onset at a different location. Figure 7 shows the damage development of the two material models. It can be seen that the damage develops faster and the second damage front establishes at a lower strain. Figure 8 shows the results from the simulation with the micro model subjected to shear. The force is normalized by the maximal force occurred in the graph. The displacement is normalized by the maximal displacement the model is subjected to. In comparison to Figure 6, the strain rate dependence is similar. This means that the slope of the curve increases with an increasing strain rate, also the maximum force increments concerning the strain rate. Thus, noticeable differences in the damage model chosen are shown.
In comparison to tension loading, it seems that for higher strain rates, the differences in the two models can be neglected. Under shear loading at lower strain rates, the anisotropic damage model leads to a more pronounced damage development compared to the isotropic damage model. This effect attenuates for an increasing strain rate, so that at a strain rate of _ ε = 10, no significant difference is observed. Also, it can be depicted that the displacement where the two models differ increases with an increasing strain rate. Figure 9 shows the damage development at three different normalized displacements. It is found that for u 1 and u 2 , no significant difference in the damage development is present. The damage state at u 3 depicts a different appearance. Here, the two models differ from each other. This evolution is visible in the force-displacement curve (cf. Figure 7).

| CONCLUSIONS
A numerical model for the strain rate dependence of an epoxy resin with isotropic damage was implemented and extended with an anisotropic damage model. The influence of five nonphysical parameters (C 1 , C 2 , C 3 , m, ε diff ) on the stress curve for pure epoxy resin was investigated. It was found that not all parameters lead to significant changes in the stress-strain curve.
To show the influence of the strain rate and damage model, further investigation on a five fibers micro model was performed. The micro model was subject to either tension or shear loading. It was found that with an increment of the strain rate, the maximum force and the displacement to failure of the model increases. The increase in the displacement to failure is in contrast to the material behavior of the epoxy resin where a higher strain rate leads to a lower strain-failure. Besides, it was found that the selected damage model has a direct impact on the micro model's material behavior.
In general, the anisotropic damage model leads to pronounced damage propagation in the resin and hence the maximal forces and the displacement to failure are generally lower. Furthermore, a reduction of the differences between the two evaluated damage models under the shear loading condition was observed.