Modelling of Fatigue Delamination Growth and Prediction of Residual Tensile Strength of Thermoplastic Coupons

1. Introduction

Composite materials and especially Carbon Fiber Reinforced Plastics (CFRPs) have become the major structural material in aerospace applications [1]. Although in most applications thermoset CFRPs have been used, thermoplastics find an increasing use due to specific characteristics such as the increased impact performance, the recyclability, and the weldability. However, the benefits of thermoplastics in the weight reduction and the lowering of recurring costs in aircraft production can be achieved only by the integration of several disciplines and the production of large integral parts. This means that focusing on small parts alone cannot achieve the full benefit of thermoplastic technology. In 2019, the Clean Sky 2 project STUNNING (Next Generation Multifunctional Fuselage Demonstrator) (https://www.clean-aviation.eu/media/results-stories/the-next-generation-multifunctional-fuselage-demonstrator-leveraging-thermoplastics-for-cleaner), aiming to build an 8-meter long narrow-body airliner fuselage from thermoplastic composites, was launched. After four years of development, a full thermoplastic fuselage part has been built using innovative manufacturing methods and incorporating novel joining techniques such as welding and co-consolidation. Characterization of the mechanical performance of the fuselage part is scheduled for the next phase of the program. In this process, the availability of efficient numerical models capable for simulating fatigue damage and predicting the residual strength of thermoplastics and their joints is of crucial importance as it could lead to a significant reduction of characterization time and cost.

Delamination, together with matrix cracking, are the primary fatigue failure modes in CFRPs, especially for high cycle fatigue. It is therefore very important to be able to simulate fatigue delamination propagation in CFRPs. Contrary to thermoset CFRPs, there have been reported only a few experimental and numerical works on delamination of thermoplastics; most of them have considered quasi-static loads. In the first effort, O’Brien [2] has conducted delamination tension fatigue tests on AS4/PEEK composite laminates. Uematsu et al. [3] have studied experimentally the behavior of delamination crack propagation of a carbon-fiber-reinforced thermoplastic polymer at high temperatures. Lachaud et al. [4] have conducted experimental and numerical studies of delamination caused by local buckling on thermoset and thermoplastic carbon-fibre composites. Zhou et al. [5] have described the mode I delamination behaviour of a unidirectional carbon-fibre/polyphenyleneetherketone (PEK-C) composite experimentally. Dube et al. [6] have performed an experimental investigation characterising the fatigue failure mechanisms of resistance-welded thermoplastic composites skin/stringer joints. Delamination has been found to be a major failure mechanism of the joints. Ruzek et al. [7] have conducted an experimental investigation of the effects of the loading frequency on the temperature change, fatigue behaviour, and failure mechanisms of carbon-fibre-fabric-reinforced polyphenylenesulphide (PPS) laminates, the thicknesses of which were varied by ply drops.

The numerical works are much less than the experimental ones and they have mainly considered quasi-static and impact loads. Keiichi et al. [8] have investigated the stress-strain response and the damage initiation/propagation mechanisms of the T700G/LM-PAEK material, in an open-hole configuration, experimentally and numerically. Liu et al. [9] have developed 3D FE models using the LS-DYNA program to study the effects of impact energy, ply angle and interfacial strength on the low-velocity impact performance of thermoplastic laminates. Sun et al. [10] have made comparisons between experimental and numerical studies of low-velocity impact damage for thermoplastic (IM7/PEEK) composites. Three LS-DYNA FE models (standard, continuum damage mechanics (CDM) and discrete) were developed, all using cohesive interface elements for delamination. Recently, Sioutis and Tserpes [11] have developed a fatigue interfacial crack growth model based on the cohesive zone modeling method which has been proven capable of efficiently simulating fatigue interfacial fracture of co-consolidated thermoplastic joints. Provided the co-consolidated interfaces have the same mechanical performance with the ply interfaces, the model of [11] could be used for simulating fatigue delamination of thermoplastic coupons.

The objective of the present paper is twofold: to examine whether the model of [9] could be used for simulating fatigue delamination propagation and to combine the fatigue model with quasi-static progressive damage models to develop a tool capable of predicting the residual tensile strength of thermoplastic coupons which when upscaled it may be used in the characterization process and the damage tolerant design of thermoplastic structures.

2. Experimental methods

The composite coupons were made from thermoplastic TC 1225 LM PAEK prepreg plies with a fiber volume fraction of 66%. The coupons were approximately 250 mm long, 25 mm wide and 2.24 thick. Tabs from the same material with dimensions of 50 mm x 25 mm x 2.24 mm were bonded to the coupons. The 0.14 mm-thick plies were stacked in a quasi-isotropic stacking sequence following a lay-up of [-45/0/45/90]2s [12].

Tension tests were conducted according to the ASTM D3039 standard [13] and fatigue tests according to the ASTM D3479 standard [14]. All tests were conducted using an MTS Universal Testing machine with a load cell of 100 kN and a data acquisition program. Tension tests were conducted in displacement control at a 0.5 mm/min loading rate. Fatigue tests were conducted in force control using a frequency of 5 Hz, a stress ratio of 0.1 and a maximum stress equal to 0.6 of the static strength.

The following test series were conducted:

• Five coupons were tested in tension to characterize the reference tensile behavior of the thermoplastic material.
• Five coupons were tested in fatigue to characterize the reference tension-tension fatigue behavior of the thermoplastic material.
• Six specimens were subjected to fatigue up to cycles equal to 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 of the average fatigue life of the coupons.
• The six fatigued coupons were tested in tension to measure the residual tensile stiffness and strength of the coupons.

3. Modelling approach

The adopted modeling approach is based on the work of Tserpes et al. [15]. It combines a quasi-static progressive damage model (PDM) [16,17], the cohesive zone modeling (CZM) method and a fatigue interfacial crack growth model [11]. The PDM is used to simulate intralaminar damage and the CZM method together with the fatigue crack growth model to simulate fatigue delamination. The execution of the models follows the experimental procedure described in Section 2. Initially, the PDM is used to simulate the tensile behavior of the thermoplastic material. Then, the fatigue crack growth model is used to simulate fatigue delamination. Finally, the PDM is used to predict the residual tensile strength of the thermoplastic coupons at various fatigue delamination states. At all stages, the numerical results are compared with the experimental ones for the models’ validation.

4. Numerical models

4.1. Quasi-static progressive damage model

To simulate the tensile behavior of the thermoplastic coupons, a progressive damage model (PDM) was used. The PDM was realized through the material model MAT-162 [18] of the LS-Dyna, which comprises a set of Hashin-type strain-based failure criteria [19] along with a damage mechanics-based material property degradation module proposed by Matzenmiller et al. [20]. The Hashin-type strain-based failure criteria considered are listed below.

Tension-Shear fiber mode 1

$$f_1-r_1^2={\left(\frac{E_a⟨{\epsilon }_a⟩}{S_{aT}}\right)}^2+\left(\frac{G_{ab}^2{\epsilon }_{ab}^2+G_{ca}^2{\epsilon }_{ca}^2}{S_{FS}^2}\right)-r_1^2=0,$$

(1)

Compression fiber mode 2

$$f_2-r_2^2={\left(\frac{E_a⟨{\epsilon }_a^{\prime }⟩}{S_{aC}}\right)}^2-r_2^2=0,{\epsilon }_a^{\prime }=-{\epsilon }_a-\frac{⟨-E_c{\epsilon }_c-E_b{\epsilon }_b⟩}{2E_a},$$

(2)

Crush mode 3

$$f_3-r_3^2={\left(\frac{E_c⟨-{\epsilon }_c⟩}{S_{FC}}\right)}^2-r_3^2=0,$$

(3)

Transverse compressive matrix mode 4

$$f_4-r_4^2={\left(\frac{E_b⟨-{\epsilon }_b⟩}{S_{bC}}\right)}^2-r_4^2=0,$$

(4)

Perpendicular matrix mode 5

$$f_5-r_5^2={\left(\frac{E_b⟨{\epsilon }_b⟩}{S_{bT}}\right)}^2+{\left(\frac{G_{bc}{\epsilon }_{bc}}{S_{bc0}+S_{SRB}}\right)}^2+{\left(\frac{G_{ab}{\epsilon }_{ab}}{S_{ab0}+S_{SRB}}\right)}^2-r_5^2=0,$$

(5)

Parallel matrix (delamination) mode 6

$$f_6-r_6^2=S_{delm}^2\left\{{\left(\frac{E_c⟨{\epsilon }_c⟩}{S_{cT}}\right)}^2+{\left(\frac{G_{bc}{\epsilon }_{bc}}{S_{bc0}+S_{SRC}}\right)}^2+{\left(\frac{G_{ca}{\epsilon }_{ca}}{S_{ca0}+S_{SRC}}\right)}^2\right\}-r_6^2=0$$

(6)

where ⟨ ⟩ are the Macaulay brackets, S_aT and S_aC are the tensile and compressive strengths in the fiber direction, and S_FS and S_FC are the layer strengths associated with the fiber shear and crush failure, respectively, S_bT and S_cT are the transverse tensile strengths of the corresponding tensile modes (ε_b > 0 or ε_c > 0), and S_ab0, S_bc0, and S_ca0 are the quasi-static shear strength values. Note that a scale factor Sdelm is introduced to provide a better correlation of the delamination area with the experiments. The scale factor S_delm can be determined by fitting the analytical prediction to the experimental data for the delamination area.

Under compressive transverse strain, ε_b < 0 or ε_c < 0, the damaged surface is considered to be “closed”, and the shear strengths are assumed to depend on the compressive normal strains based on the Mohr-Coulomb theory

$$S_{SRB}=E_b\tan \left(\phi \right)⟨-{\epsilon }_b⟩$$

(7)

$$S_{SRC}=E_c\tan \left(\phi \right)⟨-{\epsilon }_c⟩$$

(8)

where φ is a material constant as tan(φ) is similar to the coefficient of friction. The damage thresholds r_j, j = 1,…,6, have the initial values equal to 1 before the damage is initiated and are updated due to damage accumulation of the associated damage modes.

A set of damage variables ${\overline{\omega }}_i$ with _i = 1, …, 6; are introduced to correlate the onset and growth of damage to stiffness losses in the material. The compliance matrix $\left[S\right]$ is related to the damage variables as

$$\left[S\right]=\left[\begin{array}{cc}\begin{array}{ccc}\frac{1}{\left(1-{\overline{\omega }}_1\right)E_a}& -\frac{v_{ba}}{E_b}& -\frac{v_{ca}}{E_c}\\ -\frac{v_{ab}}{E_a}& \frac{1}{\left(1-{\overline{\omega }}_2\right)E_b}& -\frac{v_{cb}}{E_c}\\ -\frac{v_{ac}}{E_a}& -\frac{v_{bc}}{E_b}& \frac{1}{\left(1-{\overline{\omega }}_3\right)E_c}\end{array}& 0\\ 0& \begin{array}{ccc}\frac{1}{\left(1-{\overline{\omega }}_4\right)G_{ab}}& 0& 0\\ 0& \frac{1}{\left(1-{\overline{\omega }}_5\right)G_{bc}}& 0\\ 0& 0& \frac{1}{\left(1-{\overline{\omega }}_6\right)G_{ca}}\end{array}\end{array}\right]$$

(9)

The stiffness matrix [C] is obtained by inverting the compliance matrix, i.e., $\left[C\right]={\left[S\right]}^{-1}$. As suggested in [20], the growth rate of damage variables, ${\overline{\omega }}_i$, is governed by the damage rule of the form:

$${\dot{\overline{\omega }}}_i=max\left\{{\dot{\phi }}_j{q}_{ij}\right\},$$

(10)

where the scalar damage functions ${\dot{\phi }}_j$ controls the amount of growth and the vector-valued matrix $q_{ij}$ (i = 1,…,6, j = 1, …, 6) provide the coupling between the individual damage variables (i) and the various damage modes (j). The damage criteria of Eqs. (1) to (6) provide the damage surfaces in strain space for the unidirectional and fabric lamina models, respectively. Damage growth, ${\dot{\phi }}_j>0$, will occur when the strain path crosses the updated damage surface $f_j-r_j^2=0$ and the strain increment has a non-zero component in the direction of the normal to the damage surface, i.e., $\sum _i\frac{\partial f_j}{\partial {\epsilon }_i}{\dot{\epsilon }}_i>0$. Combined with damage growth functions ${\gamma }_j\left({\epsilon }_i,{\overline{\omega }}_i\right)$, ${\dot{\phi }}_j$ is assumed to have the form:

𝜑𝑗=𝑖𝛾𝑗𝜕𝑓𝑗𝜕𝜀𝑖𝜀𝑖

(11)

Choosing

$${\gamma }_j=\frac{1}{2}\left(1-{\phi }_j\right)f_j^{\frac{m_j}{2}-1}$$

(12)

And noting that

$$\sum _i\frac{\partial f_j}{\partial {\epsilon }_i}{\dot{\epsilon }}_i={\dot{f}}_j,$$

(13)

for the quadratic functions given by Equations (2.53) to (2.58) lead to:

$${\dot{\phi }}_j=\frac{1}{2}\left(1-{\phi }_j\right)f_j^{\frac{m_j}{2}-1}{\dot{f}}_j$$

(14)

Eq. (14) is the damage coupling matrix, and Eq. (15) and Table 1 explain how it is associated with the modulus reduction for the unidirectional lamina model.

$$q_{ij}^U=\left[\begin{array}{cc}\begin{array}{ccc}1& 1& 1\\ 0& 0& 1\\ 0& 0& 1\end{array}& \begin{array}{ccc}0& 0& 0\\ 1& 1& 0\\ 0& 0& 1\end{array}\\ \begin{array}{ccc}1& 1& 1\\ 0& 0& 1\\ 1& 1& 1\end{array}& \begin{array}{ccc}1& 1& 0\\ 1& 1& 1\\ 0& 0& 1\end{array}\end{array}\right]i=1,\dots ,6;j=1,\dots ,6$$

(15)

4.2. Fatigue crack growth model

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

To simulate fatigue delamination growth, a recently developed fatigue crack growth model for co-consolidated thermoplastics [11] was used. The model relies on a modified Paris’ law, which is used for the calculation of fatigue crack growth rate $\mathrm{d}a/\mathrm{d}N$ in form of

$$\frac{da}{dN}=c{\left(G_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)}^m$$

(16)

The intermediate mixed-mode c and m parameters are derived using the prediction model of Russel & Street [21], which requires as input only pure mode I and mode II fatigue data.

The numerical fatigue degradation has been implemented via the method of cumulative static and fatigue element damage parameter given by

$$d_{\mathrm{t}\mathrm{o}\mathrm{t}}=d_s+d_f,$$

(17)

where $d_{tot}$ is the total damage variable (0: undamaged, 1: failed), while $d_s$ and $d_f$ are the static and fatigue contributing damage variables, respectively. The stress state σ for each element in the fatigue activated zone has been degraded by

$$\sigma =\left(1-d_{\mathrm{t}\mathrm{o}\mathrm{t}}\right){\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$$

(18)

where ${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum traction of the affected element.

The fatigue loading is based on the loading envelope technique [22], where the actual sinusoidal spectrum is approached by a constant force which corresponds to the maximum fatigue loading value as can be seen in Figure 1. The iterative numerical procedure for the fatigue crack growth model is graphically presented in Figure 2.

The pure mode I and mode II modified Paris’ laws that have been used in the present application are:

$${\frac{da}{dN}}_I=0.0027G_{I,max}^{4.23}$$

(19)

$${\frac{da}{dN}}_{II}=0.00040G_{II,max}^{5.31}$$

(20)

where $\frac{da}{d{N}_I}$, $\frac{da}{d{N}_{II}}$ and $G_{I,\mathrm{m}\mathrm{a}\mathrm{x}}$, $G_{II,\mathrm{m}\mathrm{a}\mathrm{x}}$ are the fatigue crack growth rates and the maximum energy release rates per cycle for mode I and mode II, respectively. Figure 3 depicts the experimental $\frac{da}{d{N}_I}-G_{max}$ from which the modified Paris laws were derived. More details on the derivation of equations (19) and (20) and the respective graphs in Figure 3 can be found in [11].

4.3. FE models

The specimens modeled had the same dimensions as the coupons tested (see Section 2). 8-noded solid elements were at all FE models. The global damping for the models was specified at 800 for the entire analysis, while a Rayleigh damping coefficient of 0.1 was applied to achieve a higher structural rigidity. Furthermore, to eliminate non-physical hourglass instabilities, a Flanagan-Belytschko integration stiffness formulation was selected with a coefficient of 0.1. The mechanical properties of the thermoplastic layer assigned to the elements are listed in Table 2. The ply interfaces of the fatigue and residual strength FE models were modeled using intermediate cohesive layers. The cohesive elements were modeled using the UMAT 43c material model [18], designed for three-dimensional cohesive elements, and implemented as zero-thickness elements. The quasi-static tensile load was modeled by fixing the nodes at one end and applying an incremental displacement at the other end. The fatigue load was modeled by increasing the applied displacement until the maximum load is reached.

Figure 4. Mesh density of quasi-static tension model (a) front view and (b) side view.

Figure 4. Mesh density of quasi-static tension model (a) front view and (b) side view.

To predict the residual tensile strength of the thermoplastic coupons, the PDM and the fatigue model are combined by means of the numerical procedure described in the flowchart of Figure 5. For each fatigue state of interest, a dynain.ASCII keyfile was extracted, encompassing the geometry and the deleted cohesive elements, which delineate the delamination area of the state. Subsequently, the isolated deleted cohesive elements were imported into the initial fatigue model. Following this, the restoration of deleted elements to their original geometry was selected, resulting in the creation of double nodes at locations where delamination had occurred. The cohesive elements were then merged and removed, and tensile loading using MAT_162 (with properties identical to those in the model validation) was applied to extract the residual strength.

The failure mode of matrix cracking presents a challenge in precise modeling, as it does not conform neatly to established fatigue model or pre-existing damage state within residual strength prediction model. The information gleaned from experimental photos and C-Scan data, while valuable, remains insufficient to offer a clear and objective understanding of the evolution of matrix cracking. Consequently, MAT_162 has been chosen as the material model for both quasi-static models. This selection is based on its capability to delineate the evolution of matrix cracking and delamination at each layer, thereby providing a more comprehensive insight into the intricate mechanisms underlying these phenomena.

5. Numerical results

5.1. Tensile behavior

Figure 6 plots the experimental and numerical tensile force-displacement curves of the thermoplastic coupons. The observed scatter of the experimental curves is acceptable. The average tensile longitudinal strength is 917 MPa and the average longitudinal modulus is 59.75 Mpa. Regarding the numerical results, these represent very well the initial stiffness of the coupons, the maximum applied displacement, and the tensile strength. However, it does not capture accurately the decrease in the slope of the curves that takes place after the applied displacement of 1 mm, while the curvature of the experimental curve is due to visco-elastic behavior of the thermoplastic matrix and cannot be modelled precise, while knowledge about this behavior is not available from the manufacturer.

The main damage modes that drove the tensile failure of the coupons are matrix cracking and delamination. Matrix cracking initiated at 25 kN at the supported end of the coupon as can be seen in Figure 7. Specifically, it started at the middle 90˚ layers and propagated to the other 90˚ and ±45˚ layers. On the other hand, delamination initiated at 46 kN to a large extent and caused a sudden load-drop in the numerical curve, which does not appear in the experimental curves. Delamination started in the middle area of the coupon (Figure 8) at the middle 90˚/90˚ interface and propagated at the 90˚/±45˚ interfaces towards the supported end of the coupon.

5.2. Fatigue life

Figure 10 summarizes the results from the non-interrupted fatigue tests. The purpose of this test series was to measure the average fatigue life based on which the settings (number of intervals and number of cycles between intervals) of the interrupted fatigue tests will be determined. The average fatigue life is 106612 cycles, and the standard deviation is 25514 cycles. The relatively high standard deviation is attributed to the complexity of the fatigue failure mechanisms and the manufacturing defects which were present in the coupons. Figure 11 presents the fatigue lives of the coupons that have undergone the interrupted fatigue testing. In this case, the average fatigue life is 140037 cycles, and the standard deviation is 46640 cycles. The increase in the standard deviation compared to the non-interrupted tests is probably due to the processing errors that were introduced during the several stops and restarts of the testing process.

5.3. Fatigue delamination

The use of cohesive elements at each interface enables the simulation of delamination propagation as a function of the number of cycles in the entire coupon. For validation purposes, the predicted delamination propagation at the -45˚/90˚ interface is compared with C-scan images taken during interruptions of the fatigue test in Figure 12. As revealed, the model captures accurately both the initiation of delamination at the coupon’s edges at 10000 cycles and the propagation rate of delamination towards the middle of the coupon.

Figure 13 depicts the delamination propagation at all interfaces of the coupon for every 50000 cycles. The quasi-isotropic lay-up of the coupons has caused an unsymmetric and unequal delamination propagation between the interfaces. Delamination is mainly propagating to the -45˚/90˚ and the 90˚/45˚ interfaces at one side of the coupon and at the interfaces of the 90o layers, while it is almost absent from interfaces of the 0˚ layers. The different delamination propagation between the symmetric interfaces is confirmed by the experiments as can be seen in Figure 14 in the side-view photos taken from two coupons during the interruption of the tests. In the photos, the one-side extensive delamination at the neighboring 45˚/90˚ and 90˚ /45˚ interfaces is clearly visible. Also, at the later stages of fatigue loading, matrix cracking is also visible as for instance in Figure 14(a) at 100000 cycles.

5.4. Residual tensile strength

Combining the fatigue delamination model and the PDM, the residual tensile strength of the coupons was predicted. The numerical predictions are compared with the experimental measurements in Figure 15. As expected, both methods show a degradation of residual tensile strength with increasing fatigue cycles. The numerical values deviate from the experimental ones for the percentages of fatigue life of 0.2, 0.4 and 0.6. From 0.6 to 1.0 a very good agreement is observed between the two methods. The main difference between the model and the tests is that the model considers delamination, thus neglecting all other fatigue damage modes such as matrix cracking and fiber/matrix interface failure. However, this assumption should have led to an overestimation of residual tensile strength, which is not the case. The fact that the model underestimates the residual tensile strength for the percentages of fatigue life of 0.2, 0.4 and 0.6 reveals that the neglected failure modes have not accumulated significantly at those stages and that the modeling approach for delamination (nodes release) is too severe. Moreover, the good agreement at the high fatigue stages is because delamination in the tested coupons has become severe.

Regarding the predicted damage propagation in the coupons subjected to tension after fatigue, the pre-existing delamination is causing an extensive matrix cracking, which should have been already present if it has been considered by the fatigue model, and a fast propagation of delamination, especially close to the loaded ends. This has led to the fracturing of the coupons in two different areas located close to the two ends, as shown in Figure 16(a). This finding is verified by the tested coupons (Figure 16(b) and (c)).

6. Conclusions

In the present work, we have proposed a numerical modeling approach for simulating fatigue delamination growth and predicting the residual tensile strength of quasi-isotropic thermoplastic coupons. The development and implementation of the models was supported by quasi-static tension tests and fatigue tests. The model is capable of simulating fatigue delamination growth at all interfaces, thus giving a detailed insight of delamination, which is a crucial failure mechanism in fatigue of CFRPs. The respective results have been validated successfully against C-scan images taken on fatigued coupons at interruptions of fatigue tests. The fatigue model does not account for the other types of damage except for delamination. The combination of the fatigue model with the quasi-static PDM has led to accurate predictions of residual tensile strength of the coupons. This finding is an indication that the consideration of delamination alone is sufficient for capturing the degradation of residual strength in thermoplastic coupons. The proposed modeling approach if upscaled to be applicable to structural parts can be proved useful in the damage tolerant design and the development of structural health monitoring systems of composite structures.