*Article* **Modeling Stiffness Degradation of Fiber-Reinforced Polymers Based on Crack Densities Observed in Off-Axis Plies**

**Matthias Drvoderic, Martin Pletz and Clara Schuecker \***

Chair of Designing Plastics and Composite Materials, Department of Polymer Engineering and Science, Montanuniversitaet Leoben, 8700 Leoben, Austria; matthias.drvoderic@unileoben.ac.at (M.D.); martin.pletz@unileoben.ac.at (M.P.)

**\*** Correspondence: clara.schuecker@unileoben.ac.at

**Abstract:** A model that predicts the stiffness degradation in multidirectional reinforced laminates due to off-axis matrix cracks is proposed and evaluated using data from fatigue experiments. Off-axis cracks are detected in images from the fatigue tests with automated crack detection to compute the crack density of the off-axis cracks which is used as the damage parameter for the degradation model. The purpose of this study is to test the effect of off-axis cracks on laminate stiffness for different laminate configurations. The hypothesis is that off-axis cracks have the same effect on the stiffness of a ply regardless of the acting stress components as long as the transverse stress is positive. This hypothesis proves to be wrong. The model is able to predict the stiffness degradation well for laminates with a ply orientation similar to the one used for calibration but deviates for plies with different in-plane shear stress. This behavior can be explained by the theory that off-axis cracks develop by two different micro damage modes depending on the level of in-plane shear stress. It is found that besides influencing the initiation and growth of off-axis cracks, the stiffness degradation is also mode dependent.

**Keywords:** crack detection; fiber-reinforced polymers; fatigue damage model; composite fatigue; off-axis cracks

### **1. Introduction**

Components made from multidirectional fiber-reinforced composite laminates experience several distinct damage mechanisms when exposed to fatigue loads. The macroscopic damage mechanisms are matrix cracks, delamination, and fiber failure. This sequence of damage mechanisms during fatigue loading can be categorized into characteristic damage states [1]. The first fatigue-damage state of multi-axial laminates is matrix cracking in off-axis plies where multiple matrix cracks develop and grow in number and length. These so-called off-axis cracks typically span the whole thickness of a ply and propagate along the fiber direction. One of the main effects of off-axis cracks is a significant stiffness reduction of the laminate but not immediate failure of the component [2–8].

Since multiple similar cracks form under fatigue loading, the crack density is used as a measure for the amount of damage in the material in many progressive damage models [9–17]. These models describe the evolution of off-axis cracks as well as their effect on the stiffness of a laminate. Therefore, off-axis crack densities are often used in the development and calibration of these models or to compare their predictions to experimental data. Transmitted or transilluminated white light imaging (TWLI) can be used for transparent composites like glass fiber-reinforced polymers (GFRPs). It is an efficient, reliable and relatively simple method to capture off-axis cracks [2,10,18–20]. TWLI uses a light source placed behind a transparent specimen and a camera on the other side. An undamaged specimen appears bright as it only absorbs a small portion of the light. Cracks, on the other hand, scatter the light and therefore appear as dark lines in the images. For non-transparent laminates, more sophisticated techniques like computed tomography

**Citation:** Drvoderic, M.; Pletz, M.; Schuecker, C. Modeling Stiffness Degradation of Fiber-Reinforced Polymers Based on Crack Densities Observed in Off-Axis Plies. *J. Compos. Sci.* **2021**, *6*, 10. https://doi.org/ 10.3390/jcs6010010

Academic Editor: Stelios K. Georgantzinos

Received: 7 December 2021 Accepted: 23 December 2021 Published: 29 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

may be used [21]. Up to now, a few methods have been developed to compute the crack density from TWLI images. The simplest is to count all cracks along a straight line normal to the fiber direction and divide the number of counted cracks by the length of the line. This method only takes the number of counted cracks and not their length into account and the results are influenced by the selection of the path. As shown in Refs. [17,22,23], including the crack length results in an improved description of the average damage state of the material. A better approach that includes some information of the crack length is the weighted crack density used by Quaresimin et al. [7]. It clusters the cracks into eight groups of crack lengths and then computes a weighted average. Still, manually counting and categorizing the cracks is labor intensive and prone to human errors. An automated algorithm that takes the images as input and detects all cracks in a given direction has been developed by Glud et al. [24]. Based on this algorithm, we have developed *CrackDect*, an open-source package for evaluating crack densities [25]. With this package, even large fatigue test series can be evaluated efficiently. Figure 1 qualitatively shows the evolution of off-axis cracks and the associated stiffness degradation as well as examples of images taken during fatigue tests. CrackDect takes these images and computes the crack density.

**Figure 1.** Evolution of off-axis cracks during fatigue tests. The stiffness stays constant as long as there are no cracks (1). After the onset of matrix cracking (2), the cracks grow in number and size and the stiffness starts to decrease due to the damage. (3) The last stage of matrix damage in composites is crack saturation and finally total failure usually due to other damage mechanisms like fiber failure and delamination.

Many fatigue models have been developed to describe the effect of fatigue damage in composite laminates. Degrieck and Van Paepegem [26] sorted them into three categories: fatigue life models, phenomenological models, and progressive damage models. To accurately describe the effect of distinct damage modes, progressive damage models are the most promising candidates since they take the actual cause of the degradation of a mechanical property—the damage—into account. Usually, this is done in a two step approach. In the first step, a damage model describes the evolution of a damage variable with respect to the undamaged material. The second step computes the effect of this damage on mechanical properties. Many models have been developed that establish a connection between off-axis cracks and laminate stiffness [9–17,27]. Often, the Finite Element Method (FEM) is used to compute the elastic response of a laminate. One widely-used approach is to model a representative volume element of the laminate with cracks in one or more plies [17,19]. Resolving the laminate into its plies and computing the effect of cracks on the ply-level has the advantage that the elastic response of any laminate can be computed. The drawback of FEM is that it is time-consuming when implemented as a sub-model for evaluation of large composite structures. Carraro et al. [9] further compared several modeling techniques and developed an analytical model based on shear lag analysis with the capability to compute the elastic response of any symmetric laminate. One drawback of this model is that the laminate must be symmetric and the elastic response is computed for cracks in both symmetric plies. This means that bending loads which yield unsymmetrical stress distributions in the laminate cannot be accounted for. Schuecker et al. [27] proposed a static damage model that computes the stiffness degradation based on the Mori-Tanaka method. The Mori-Tanaka method is a micro-mechanical approach to compute homogenized material properties of a material consisting of a matrix and an embedded inclusion [28]. This allows to compute the stiffness degradation for each ply separately and then combine all plies using classical laminate theory to compute the overall stiffness of the damaged laminate. In Schuecker's approach, the effect of crack-like void inclusions on the stiffness of a ply is computed. Even though this model has been developed for static load cases, it should also be applicable to fatigue since the stiffness degradation is only dependent on the amount of cracks. The main advantages of this model are:


Naturally, the model has limitations and prerequisites arising from classical laminate theory and the Mori-Tanaka method, but it has proven to provide an overall relatively simple yet effective approach to compute stiffness degradation of laminates due to cracks. Also, it does not consider delaminations or other damage mechanisms of composite laminates [27,30].

In this work, Schuecker's degradation model is combined with crack detection by replacing the evolution function for static loads by the crack evolution detected from experimental data. Based on the off-axis crack density from experiments, the stiffness degradation computed by the degradation model is compared to experimental stiffness data. Opposed to other works like [31], where crack densities of similar specimens have been averaged, we take the crack density of each individual test and compute the resulting stiffness degradation. This enables a comparison of experimental stiffness data and predictions based on experimental crack density data for individual specimens. Also, a mostly-automated procedure to fit experimental crack density data from the automated crack detection is presented.

### **2. Methods**

### *2.1. Experimental Fatigue Data*

Experimental data from fatigue tests of ±*θ<sup>s</sup>* GFRP laminates from [32] is used for comparison with computations of the damage model. The specimens with a gauge length of 100 mm, a width of 20 mm, and a thickness of 2 mm (12 layers) had been cut from GFRP laminate plates produced from a unidirectional glass fiber weave and epoxy resin. The plates had been produced by vacuum pressing manually-impregnated glass fiber layers. The stacking sequence of the laminates is [+*θ*3/ − *θ*3]*s*. For the fatigue tests, stresscontrolled sinusoidal load cycles with an R-ratio of 0.1 and a frequency of 5 Hz had been

periodically interrupted to perform displacement-controlled quasi-static tensile tests. The servo-hydraulic material testing system MTS 810 by MTS Systems Corporations had been used for all tests with an optical displacement measurement system (CV-X100 by Keyence) to provide a free field of few for the images taken during the tests. This procedure allows to track the change in stiffness as a function of the number of cycles and take images for the crack detection. For a more detailed description of the experimental setup, the reader is refered to [32]. In addition to the stiffness of the laminates, the crack density is evaluated from TWLI images with crack detection. For this work, the results of the fatigue tests for ±45°, ±60°, and ±75° laminates are used. Laminates with a ply orientation of less than ±45° show delamination as the main damage mode and only little off-axis cracks. The unidirectional 90° laminates show hardly any cracks before final failure.

For each laminate type, two fatigue tests had been conducted. Table 1 lists the static stiffness, transverse strength *R*2, and in-plane shear strength *R*<sup>12</sup> of the material. In the referenced data from [32], a miscalculation had happened in the evaluation of the in-plane shear strength and in-plane shear modulus *G*12, which is corrected in this work. Additionally, the static ply properties are corrected with respect to the fiber volume fraction of the individual specimens. The procedure is described in Appendix A. The aforementioned miscalculation does not effect the validity of the tests since only the evaluation had to be redone. The load level of the fatigue tests, which is the ratio of the maximum load in the fatigue test to the static strength of the laminate, is computed by the Puck failure criterion [33]. The exact computation is given in Appendix B. The load levels of the tests are 75% for ±45°, 78% for ±60°, and 74% for ±75° laminates. The slight differences between load levels arise from the corrections done in the evaluation.

**Table 1.** Elastic constants of the composite ply material from static tests.


### *2.2. Crack Detection*

The Python package *CrackDect* is used to to automatically evaluate the crack density from the TWLI images [25]. This package includes a sightly modified crack detection algorithm compared to [24]. Example pictures of a specimen at the beginning, during, and at the end of the fatigue test prior to image processing are shown in Figure 2. The processing pipeline of the images is as follows:


The exact procedure of this processing pipeline is explained in [25] where the opensource code of all functions can be obtained. The input parameters are listed in Table 2 for each test series. It is observed that the crack width of the first major visible cracks varies slightly with the fiber orientation. Cracks in the ±45° specimens appear to be thinner than in ±75° specimens. Therefore, the average width of cracks that should be detected by the crack detection is set individually for each test series to get comparable results between the test series. To avoid artifacts in the crack detection or false detection due to inherent noise in the images, cracks of less than 50 pixels (0.7 mm) in length are excluded from the evaluation. The crack density is defined as

$$\rho\_c = \frac{\sum\_{i=1}^{n} L\_i}{A} \quad , \tag{1}$$

with *Li* as the length of crack *i* and *A* as the area of the region of interest (evaluation area). Since the crack detection only computes the crack density based on pixels, the conversion from pixel to millimeter is also listed in Table 2.

**Table 2.** Input parameters for CrackDect. The coordinates of the region of interest *x*<sup>0</sup> to *y*<sup>1</sup> are given in pixel.


**Figure 2.** Example of TWLI images taken from a ±45° specimen with a load level of 75% before the test (**a**), after 254 cycles (**b**), and and after 4013 cycles (**c**). In (**a**), the evaluation area is shown for the crack density (region of interest) marked with the blue rectangle. In (**b**), a typical crack pattern for a ±45° is laminate shown and (**c**) shows the last image taken before failure. The shift of the specimen can easily be observed by the drift of the black line at the bottom from (**a**–**c**). In (**c**), delamination between the plies can be spotted as dark areas.

Instead of using the extracted densities directly, a crack density function is defined by fitting a cumulative Weibull distribution function to the experimental crack density data. The Weibull function is used since it has a form similar to the experimental crack density plotted over the number of cycles in logarithmic space. A direct fit of such a crack density function to the experimental data resulted in convergence problems, even with non-linear least squares algorithms (scipy.optimize.curve\_fit) [34]. Therefore, a two-step approach was used to achieve a satisfactory quality of the fit. This fitting process is qualitatively illustrated in Figure 3. The first step is a linear regression in the region where the crack density increases linearly. For experiments that reach crack saturation, the linear fit is done from 30% to 85% of the maximal crack density. If crack saturation is not reached because the specimen fails prior to that, the region for the linear regression extends to 100%. In the second step, the following three-parameter cumulative Weibull distribution function is fitted to the linear regression

$$\rho\_{\varepsilon}^{fit}(n) = \left[1 - \exp\left(-\left(\frac{n - n\_0}{\lambda}\right)^k\right)\right] \cdot \rho\_{\varepsilon}^{sat} \quad , \tag{2}$$

with *λ* and *k* as scale and shape parameters respectively, and *n*<sup>0</sup> to shift the fitted curve along the *x*-axis. Note that *n*<sup>0</sup> is a fitting parameter and does not correlate with the cycles to damage initiation *ninit* defined later in Section 3.1. Since the Weibull distribution function has a span of 0–1, the fitted crack density function is scaled with the saturation crack density *ρsat <sup>c</sup>* .

**Figure 3.** The experimentally obtained crack densities are fitted with a two-step approach. The first step is a linear regression from 30% to 85% of the saturation crack density. In the second step, a cumulative Weibull distribution function is fitted to the linear regression.

### *2.3. Damage Model*

To compute the stiffness degradation of a laminate for a given damage state, Schuecke's damage model is used [27]. Like most progressive damage models, it consists of a part that describes the evolution of a damage variable and a part that computes the effect of this damage on the stiffness of the material. The first part computes the damage variable as a function of the loads for each ply. Then, in the second part, the degradation of stiffness based on the damage variable is computed using the Mori-Tanaka method. Here, only the second part of the model is used since the evolution of the damage variable is obtained by calibration to the fitted crack density functions. The damage variable in the model represents the volume fraction of crack-like inclusions in the Mori-Tanaka formulation, which for void inclusions, is given by

$$E\_{MTM, \text{void}} = E^{(m)} \left[ I + \frac{V}{1 - V} (I - \mathcal{S})^{-1} \right]^{-1} \,, \tag{3}$$

where *V* is the inclusion volume fraction, *I* is the identity tensor, *E*(*m*) is the initial stiffness tensor of the ply, and *S* is the Eshelby tensor [35,36]. To compute the Eshelby tensor for transversely isotropic materials, the numerical computation scheme by Gavazzi et al. [37] is used.

The Eshelby tensor for an inclusion depends on the surrounding material and the shape of the inclusion. It is assumed that the inclusion geometry is the same for all cracks and independent of the orientation of the ply. Often, an extremely sharp or disk-like inclusion geometry is used when the effect of cracks in a material is computed by the Mori-Tanaka method. Experimental evidence shows that off-axis cracks are often not straight but have a crooked path since the crack has to find the way of least resistance between the fibers. Cracks sometimes even split and merge on the way through the ply [20,38]. Here, the introduction of an orientation tensor to give idealized penny shaped cracks an orientation distribution similar to the crooked paths of the real cracks is avoided. Instead, one oblate ellipsoidal pore that represents the homogenized effect of these cracks qualitatively is used. It has been reported in [39,40] that this approach gives satisfactorily results. In this work, an aspect ratio of 100,1,10 in the 1,2,3-direction of the ply is chosen. The reasoning is as follows: The cracks are substantially longer in fiber direction than in out of plane direction so the 1-direction is set to 10 times the 3-direction. Also, cracks are relatively flat compared to the thickness of the ply so the 3-direction is 10 times the 2-direction. A schematic representation of this idealized inclusion is shown in Figure 4a and its effect on the ply properties is shown in Figure 4b. The curves of *E*1, *E*2, *ν*<sup>12</sup> and *G*<sup>12</sup> as function of to the inclusion volume show that *E*<sup>2</sup> is reduced the most relative to its initial value. The stiffness in fiber direction *E*<sup>1</sup> is reduced only slightly up to an inclusion volume fraction of 0.1. This behavior qualitatively agrees well with the stiffness degradation of a ply due to off-axis cracks. The degradation

of *E*1, *E*2, *ν*<sup>12</sup> and *G*<sup>12</sup> could be fine-tuned by adjusting the aspect ratios if additional data from static tensile tests were available that allows to relate the engineering constants of a ply directly to off-axis cracks.

**Figure 4.** Representation of off-axis cracks in one ply (**a**) looking in fiber direction. The crack path is often not straight and can branch or merge. Therefore, an oblate ellipsoid is chosen as a single inclusion shape to compute the homogenized effect of these cracks. The effect of this inclusion is shown in (**b**) with the aspect ratios of 100,1,10 in 1,2,3-direction, respectively.

### *2.4. Calibration*

The Mori-Tanaka method uses inclusion volume fractions to compute the effect of an inclusion in a surrounding matrix material, but since the crack density is analyzed, the model needs to be calibrated to experimental data. A link between the crack density and the inclusion volume fraction for the analyzed stiffness degradation must be established. For this, the ±45° tests are used to calibrate the model. The following equation links the crack density (*ρc*) to the inclusion volume fraction (*V*) with a simple correlation factor (*μ*).

$$V = \mu \cdot \rho\_{\mathbb{C}} \tag{4}$$

The calibration process is illustrated in Figure 5. First, the experimental stiffness data is smoothed to reduce scatter using a lowess filter [41] with a window length of 40% the range of cycles. Then, the stiffness degradation relative to its initial value up to crack saturation (see Table 3) is calculated from this smoothed curve. Parallel to this, a model of the laminate is built with classical laminate theory. For each ply, the stiffness is reduced according to the damage model (see Equation (3)) as a function of the inclusion volume fraction. The inclusion volume fraction to reach the experimental stiffness degradation is optimized with the minimization algorithm from SciPy (*scipy.optimize.minimize*) [34]. From this, the correlation factor *μ* can easily be computed from Equation (4). This procedure is carried out for the two ±45° tests and the average is used as the global calibration constant for the material.

**Table 3.** Results of the crack detection. The cycle number to damage initiation *ninit*, saturation *nsat*, and crack density growth rate in semi-logarithmic space *dρc*/*d*(*log*(*n*)) are listed to compare the laminates.


**Figure 5.** A schematic representation of the correlation process. The stiffness and crack density is evaluated in fatigue tests. The needed inclusion volume fraction is optimized to yield the same stiffness degradation as observed in the experiment. With Equation (4), the correlation constant is then computed from the experimental saturation density.

### **3. Results and Discussion**

### *3.1. Crack Detection Results*

The results of the crack detection are shown in Figure 6. An example of the cracks detected in +60° direction for the ±60° T1 specimen after 759 load cycles is shown in Figure 6a. At this state, the crack density is 0.8 mm/mm2. The bigger cracks are detected well while cracks smaller than 50 pixel or 0.7 mm in length are filtered out. Since the laminate consists of 12 plies, cracks in the bottom layer will appear fainter and not as sharp as cracks in the top layer. Therefore, the detection is less reliable for cracks in the bottom layer. Cracks in negative fiber-direction are not included because the noise, blur and overlap with cracks in positive fiber direction from the top plies resulted in too many false detections. Since the plies in negative direction develop approximately the same amount of cracks as in positive direction (see Figure 2), only the positive direction is analyzed.

**Figure 6.** Results of the crack detection. An example of the detected cracks for ±60° T1 is shown at 759 cycles (**a**). Cracks are only detected in the chosen direction of +60°. (**b**) shows the crack density results for all tests and the crack density functions.

Figure 6b shows the evolution of the crack density over loading cycles and the corresponding fitted crack density functions for all specimens. In the ±45° laminates, off-axis cracks initiate earlier than in ±60° and ±75° laminates. The ±45° laminates also show the highest crack saturation density. The ±75° laminates do not reach the state of crack saturation because they fail before. The growth rate of the crack density seems to be approximately constant for all tests when plotted on a logarithmic scale.

The results of the crack detection are listed in Table 3. The point of damage initiation *ninit*, crack saturation *nsat* and the crack density growth rate of the linear region are used to compare the crack detection results. The saturation density *ρsat <sup>c</sup>* is taken manually at the point where the crack accumulation reaches a plateau. Since the ±75° laminates do not reach crack saturation, the maximum crack density of ±75° T1 is used for the fitted curves. In the ±45° tests, the crack density increases again after reaching a first plateau due to delaminations that cause problems with the crack detection. Therefore, the saturation density is taken at the first plateau. The point of damage initiation *ninit* is defined as the first point with a crack density above 0.1 mm−<sup>1</sup> and the point of crack saturation *nsat* are the cycles needed to reach the saturation crack density *ρsat <sup>c</sup>* . For ±75° tests, crack saturation is not reached. The crack density growth rate is the slope of the linear regression, the first step of the fitting process.

These results show that the automated crack detection is suited to efficiently obtain off-axis crack densities of multidirectional GRFP laminates. Damage initiation is easily detected and the crack density function can be modeled by a three-parameter cumulative Weibull distribution by fitting it to the linear regression of the region of constant crack density increase in semi-logarithmic space (see Figure 6b). The number of cycles until the saturation density is reached are also captured by the crack detection. It should be noted that the accuracy of the crack detection decreases when approaching saturation due to delamination and merging of cracks (see Figure 2c). From our experience, almost all cracks are detected up to approximately 80 percent of the saturation level. From there on, the merging of side-by-side cracks into one black line and delamination leads to misses and false detections. To improve the accuracy of the crack detection near saturation, image differencing techniques that allow to see only changes from one image to the next could be used (see [24]). The prerequisite to this is an extremely precise shift correction. This can be achieved with position markers on the images that allow a precise tracking of image shift and distortion. Since our images did not have these markers, a simpler version of the shift correction had to be used.

### *3.2. Stiffness Degradation Model*

The calibration of the damage model for the material from [32] yields a correlation constant of 0.011 mm for the chosen inclusion aspect ratios of 100,1,10. Note that the correlation constant depends on the chosen aspect ratios. The results of the computed stiffness degradation are shown in Figure 7 along with experimental data for all laminates. The comparison between model and experiments for the ±45° laminates and ±60° T2 shows good agreement. Since the calibration constant is computed from the ±45° tests, good agreement of the stiffness drop is expected for ±45° specimens. Nonetheless, the curves follow the same shape as the experiments which is determined by the inclusion geometry. This indicates that the chosen inclusion geometry describes the effect of off-axis cracks for this laminate and our assumptions on the inclusion geometry are reasonable. For ±60° T1, the detected crack density is higher than for ±60° T2. Therefore, the computed stiffness degradation is also higher compared to ±60° T2. The trend of experimental data and the shape of the curves from the degradation model for ±45° and ±60° laminates is similar. Contrary to the ±45° tests, the experimentally observed stiffness degradation does not stop at the saturation of the off-axis cracks for the ±60° laminates (see Figure 7). As listed in Table 3, crack saturation is reached at around 2000 cycles for the ±60° tests. The computed stiffness of the ±75° laminates drops earlier than the experimental curves. It seems that off-axis cracks in ±75° laminates do not have the same effect at the local ply coordinate system as for ±45° laminates.

The error of the model increases from the ±45° laminates, where the calibration is carried out to the ±75° laminates. Carraro et al. [42] have shown that the macroscopic damage initiation is driven by two damage modes that depend on the level of in-plane shear stress. For plies with mostly in-plane shear stresses *σ*<sup>12</sup> in the ply coordinate system shown in Figure 4a, the driving force for damage evolution at the micro scale is local maximum principal stress (LMPS). For plies with mostly positive in-plane transverse stress *σ*22, local hydrostatic stress (LHS) is the driving force. The shift between LMPS and LHS in GFRP occurs at a fiber direction of around 60° [42]. Fractographic images also show different crack patterns for off-axis cracks depending on the in-plane shear stress. It has also been found that shear stress significantly reduces the number of cycles for damage initiation [7,43]. With differences in the micro-structure of the fracture plane between the two damage types, the effect of cracks on the ply stiffness can also be expected to be different for the LMPS/LHS damage types. For our model, this would require a separate correlation for the ±45° and ±75° tests, since these tests correspond to LMPS and LHS type damage, respectively. The presence of these two separate damage types would also explain the large difference in the number of cycles to damage initiation from ±45° to ±75° laminates (see Figure 6). The ±45° and ±60° tests show damage initiation at less than 500 cycles with ±45° being a bit lower than ±60°. On the other hand, the ±75° tests show damage initiation at more than 50,000 cycles, although the fatigue load level for all tests has been set to about 75% of the static strength. These findings back the theory of two distinct microscopic damage types controlled by in-plane shear stresses.

**Figure 7.** Stiffness degradation for the individual tests and damage model results. Comparison for (**a**) ±45° tests, (**b**) ±60° tests, and (**c**) ±75° tests.

In Figure 8, the experimental stiffness is plotted over crack density up to saturation. At the beginning, a small drop in stiffness without an increase in the crack density is visible for some specimens. After this initial drop in stiffness, all curves except ±75° T1 show a linear correlation up to saturation. Note that this does not conclude that the degradation of stiffness is linear since the crack density follows a S-shaped curve. For ±75° T1, a distinctive kink compared to the linear regression is visible. Also, the scatter of the experimental data for ±75° T1 (see Figure 7) is higher at the beginning. At around 200,000 cycles, this scatter nearly vanishes. This could be an indication of a problem in the evaluation or experimental procedure for this specimen. The linear relation between crack density and stiffness shows that the crack density is a good choice as a damage parameter.

**Figure 8.** Relationship of the experimentally determined stiffness and crack density up to saturation. For the ±45° and ±60° specimen, the stiffness degradation correlates linearly with the crack density. Note that the fact that the stiffness for the ±75° specimen is higher compared to ±60° is a result of the given ply material properties and in accordance with classical laminate theory.

### **4. Conclusions**

In this work, the effect of off-axis cracks in GFRP laminates is studied by using test results of the crack density and computing the effect of these cracks on the stiffness of the laminate. It is shown that crack detection can be used to efficiently evaluate images of off-axis cracks from fatigue tests and the fatigue crack density function can be modeled by a three-parametric cumulative Weibull distribution function. A mostly-automated scheme is presented for the calibration of the crack density functions from experimental data. Furthermore, the stiffness degradation model for multidirectional fiber-reinforced polymer laminates from Schuecker, which uses Mori-Tanaka homogenization on the plylevel, is tested against experimental fatigue data. The results suggest that it is necessary to distinguish the effect of off-axis cracks on the stiffness of a ply depending on the microscopic crack type. This requires a separate calibration for cracks formed under LMPS and LHS conditions, respectively. This observation agrees well with the theory by Carraro and Quaresimin which also distinguishes the evolution of off-axis fatigue cracks based on the micro-damage mechanisms driven by LMPS and LHS. A distinct jump in cycles to damage initiation is also found where the micro-damage mechanisms change. A new test campaign is under way to further test the damage model for both microscopic damage types.

**Author Contributions:** Conceptualization, C.S., M.P. and M.D.; methodology, M.D., C.S. and M.P.; software, M.D., validation, M.D.; formal analysis, M.D.; investigation, M.D.; data curation, M.D. and C.S.; writing—original draft preparation, M.D.; writing—review and editing, C.S. and M.P.; visualization, M.D.; supervision, C.S. and M.P.; project administration, C.S.; funding acquisition, C.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** Part of this work has been performed within the COMET-project Experimental and numerical analysis of the damage tolerance behavior of manufactured induced defects and bonded repairs in structural aerospace composite parts (project-no.: VI-3.04) at the Polymer Competence Center Leoben GmbH (PCCL, Austria) within the framework of the COMET-program of the Federal Ministry for Transport, Innovation and Technology and the Federal Ministry for Digital and Economic Affairs with contributions by Montanuniversität Leoben (Chair of Designing Plastics and Composite Materials) and MAGNA Powertrain Engeneering Center Steyr GmbH CO KG. The PCCL is funded by the Austrian Government and the State Governments of Styria, Lower Austria and Upper Austria.

**Data Availability Statement:** The data presented in this study are available upon request from the corresponding author. The algorithms of the crack detection are available at GitHub (https: //github.com/mattdrvo/CrackDect, accessed on 28 December 2021).

**Acknowledgments:** Special thanks go to René Rieser and the Chair of Materials Science and Testing of Polymers, Montanuniversitaet Leoben for providing experimental data and images for the crack detection.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **Appendix A. Quasi-Static Material Parameters**

The basic characterization of the quasi-static material parameters is done with UD0° specimens for *E*<sup>1</sup> and *ν*12, UD90° specimens for *E*<sup>2</sup> and ±45° specimens for *G*12. In the referenced data set [32], a miscalculation has happened in the evaluation of the in-plane shear modulus *G*<sup>12</sup> which is corrected here according to DIN EN ISO 14129. Additional measurements of the fiber volume fraction revealed differences between the specimens. Therefore, the ply stiffness of each laminate is corrected to a fiber volume fraction of 45%. For this, the semi-empirical Chamis model is used to approximate the engineering constants as a function of the fiber volume fraction since it showed good agreement with experimental data for GRFP [44]. This modified rule of mixture (ROM) replaces the fiber volume fraction in the iso-stress model (Reuss) for *E<sup>c</sup>* <sup>2</sup> and *<sup>G</sup><sup>c</sup>* <sup>12</sup> with the root of the fiber volume fraction. The iso-strain model (Voigt) for *E<sup>c</sup>* <sup>1</sup> and *<sup>ν</sup><sup>c</sup>* <sup>12</sup> is not altered. For isotropic matrix and fibers, this leads to the following set of equations:

$$\begin{aligned} E\_1^\varepsilon &= (1 - V^f)E^m + V^f E\_1^f\\ E\_2^\varepsilon &= \left[ \frac{1 - \sqrt{V^f}}{E\_m} + \frac{\sqrt{V^f}}{E\_2^f} \right]^{-1} \\ \nu\_{12}^\varepsilon &= (1 - V^f)\nu^m + V^f \nu^f \\ G\_{12}^\varepsilon &= \left[ \frac{1 - \sqrt{V^f}}{G\_m} + \frac{\sqrt{V^f}}{G\_f} \right]^{-1} \end{aligned} \tag{A1}$$

For the correction, the elastic constants of the matrix *E<sup>m</sup>* and *ν<sup>m</sup>* must be known. Since the elastic constants for each laminate at a certain fiber volume fraction are tested, the elastic constants of the fibers can be estimated. This estimate is then reinserted in the same equations to obtain the elastic constants of the composite as a function of the fiber volume fraction. In Table A1, the fiber volume fractions for each laminate and the tested elastic constants are listed. With this data, the elastic constants at a fiber volume fraction of 45% in Table 1 are computed. In the degradation model, the stiffness is corrected for each laminate individually.


**Table A1.** Fiber volume fractions of the laminates and the quasi-static elastic constants in plycoordinates determined from tests.

### **Appendix B. Fatigue Load Level**

The fatigue load level of laminates is computed by Puck´s failure criterion [33]. A load level of 75% means that the laminate is loaded in the fatigue tests up to 75% of its static strength. In Figure A1, the stress space for Puck mode A with the stress vectors of the tests is shown. In the ±45° laminates used to test the in-plane shear strength according to DIN EN ISO 14129, significant transverse stresses *σ*<sup>22</sup> are present. The shear strength *R*<sup>12</sup> for the material is therefore corrected using Puck´s failure criterion based on the strength of the ±45° and UD90° laminates. With this corrections, the fatigue level of the tests are 75% for ±45°, 78% for ±60° and 74% for ±75° laminates. As a comparison, the stress vector for a typical carbon fiber-reinforced laminate is also shown. The higher ratio of fiber stiffness to transverse stiffness yields nearly no transverse stress so no correction is necessary for carbon fiber-reinforced composites.

**Figure A1.** Puck mode A failure area (red) with stress vectors of the tests. In the ±45° laminate, significant transverse stress *σ*<sup>22</sup> is present. Therefore, the in-plane shear strength *R*<sup>12</sup> is corrected with Pucks failure criterion (black arrow). As a comparison, the stress vector for a typical ±45° CFRP laminate is also shown (red arrow).

### **References**

