An Improved Cohesive Zone Model for Interface Mixed-Mode Fractures of Railway Slab Tracks

Abstract

The interface crack of a slab track is a fracture of mixed-mode that experiences a complex loading–unloading–reloading process. A reasonable simulation of the interaction between the layers of slab tracks is the key to studying the interface crack. However, the existing models of interface disease of slab track have problems, such as the stress oscillation of the crack tip and self-repairing, which do not simulate the mixed mode of interface cracks accurately. Aiming at these shortcomings, we propose an improved cohesive zone model combined with an unloading/reloading relationship based on the original Park–Paulino–Roesler (PPR) model in this paper. It is shown that the improved model guaranteed the consistency of the cohesive constitutive model and described the mixed-mode fracture better. This conclusion is based on the assessment of work-of-separation and the simulation of the mixed-mode bending test. Through the test of loading, unloading, and reloading, we observed that the improved unloading/reloading relationship effectively eliminated the issue of self-repairing and preserved all essential features. The proposed model provides a tool for the study of interface cracking mechanism of ballastless tracks and theoretical guidance for the monitoring, maintenance, and repair of layer defects, such as interfacial cracks and slab arches.

1. Introduction

Chinese railway track systems (CRTS) have successfully served for more than 10 years in China’s high-speed railway (CRH) and have performed well during the period. However, with increasing operation time and the influence of complex temperature and environmental conditions, hundreds of interfacial cracks (as shown in Figure 1) between track slab and cement asphalt mortar (CA mortar) have appeared on the high-speed railway tracks [1,2]. Under extremely high temperatures in summer, defects of slab arching also occur.

Typical interlayer defects, such as slab arch [3], are closely related to the interfacial cracks between the track slab and the under-layer. During operation, the track directly undertakes the effects of the cyclic loads from the high-speed train and environmental temperature, which increases the possibility of interface cracking. As a vertical multilayer and longitudinally heterogeneous structure, the slab ballastless track has weak parts between the new and old concrete interface and composite connection surface. Therefore, a reasonable simulation of interlayer interactions is the key to studying the defects of track structures.

Cohesive zone model (CZM), an effective and favored crack model in interface fracture mechanics, has been widely used to simulate crack initiation and propagation in various materials, such as metals [4,5,6], polymers [7], ceramics [8], concrete [9,10,11], and fiber-reinforced composites [12]. In CZM, material failure is characterized by a traction–separation law, which relates the traction across the crack to the corresponding separation [13]. The approach ensures that CZM maintains the continuity conditions mathematically and removes the singularity present in Linear Elastic Fracture Mechanics (LEFM) [14,15,16]. With the development of CZM over the past six decades since it was proposed by Barenblatt [17] and Dugdale [18], a large variety of traction–separation laws have been established. The most popular are bilinear [19,20,21,22], trapezoidal [23,24,25], exponential [26,27,28], polynomial [29,30,31], and so on.

Various cohesive zone models may have different applicable conditions due to different initial assumptions. For instance, the trapezoidal cohesion zone model proposed by Tvergaard [23] could not consider the situation where the mode I fracture energy is not equal to the mode II fracture energy. The exponential cohesive zone model proposed by Xu and Needleman [27] could consider the difference values of normal and tangential fracture energies, but when the two fracture energies are different, there is a “self-repairing” problem at the crack tip under mixed-mode loading and unloading.

The Park–Paulino–Roesler (PPR) model is a kind of polynomial traction–separation law for mixed-mode fractures that was proposed by Park et al. [32] in 2009. This model is versatile because it can consider different fracture energies with respect to fracture modes and can be applied to represent various material softening responses, i.e., ductile, brittle, and quasi-brittle, due to the controllable softening given by the shape parameters [13,32]. More significantly, the model guarantees the consistency of the cohesive constitutive relationship under mixed-mode conditions [30,33,34].

Due to the above advantages and the convenient implementation in commercial software ABAQUS as a user subroutine [34,35,36], the PPR model has been utilized to investigate a wide range of failure phenomena and cited in many papers. The model was found to still have limitations that need to be improved. Nguyen et al. [37] indicated that due to the different cohesive interaction regions between the normal and tangential tractions when fracture energies are different, one traction component might become zero while the other traction component had not yet vanished. This situation does not conform to reality in which normal and tangential tractions typically fail simultaneously when a fracture happens.

In addition, Spring et al. [38] noted that the unloading/reloading relationship, which was commonly utilized in conjunction with the PPR model, produced self-healing behavior when the crack underwent unloading/reloading. To address this issue, a new coupled unloading/reloading relationship, which maintained the thermodynamic consistency of the PPR cohesive model, was developed [38]. More recently, the research by Gilormini et al. [39] showed that the new unloading/reloading relationship prevented the questionable features that might appear when the original model [34,35] was used, but also bred a new issue regarding damage initiated from the very beginning of the loading process. This model ignores the initial elastic region.

In this paper, an alternative simplified PPR traction–separation law and an improved unloading/reloading relationship are developed and validated using multiple cases that could effectively eliminate the above issues and preserve all essential features of the original one. The modeling method of connections between the layers of the slab track as proposed in this paper can contribute to the mechanism of high-speed railway (HSR) interlayer defects, on-site monitoring, inspection, and maintenance.

This paper is organized as follows. The review of the original PPR model (traction–separation law) and unloading/reloading relationship are presented in Section 2. Section 3 shows the modification of the original PPR model and the comparison of the modified model and original through example cases. Section 4 introduces the improvement of the unloading/reloading relationship and demonstrates that the improved one is effective with the example used in [39]. Then, Section 5 presents the application of the proposed model to analyze interface damage of railway slab track. Finally, the paper is summarized in Section 6.

2. Original Models

The PPR model was designed for pure loading conditions and did not contain a built-in unloading/reloading relationship [38]. To simulate the fracture submitted to the general loading conditions, such as loading, unloading, and reloading, the PPR model was combined with an unloading/reloading relationship [34]. The original PPR model and unloading/reloading relationship are introduced shortly in the following subsections.

2.1. Original PPR Model

The fundamental issue in cohesive zone modeling is the definition of traction–separation law, which gives the constitutive behavior of the fracture. The original PPR model defines the traction–separation law by taking the derivative of the cohesive fracture potential. The potential consists of polynomials formulated in terms of a normal separation (${\mathsf{\Delta }}_n$) and a tangential separation (${\mathsf{\Delta }}_t$), and it is expressed as [32]:

$$\begin{gathered} \Psi \left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=\min \left({\varphi }_n,{\varphi }_t\right)+\left[{\Gamma }_n{\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha }{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^m+\langle {\varphi }_n-{\varphi }_t\rangle \right] \\ \left[{\Gamma }_t{\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta }{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^n+\langle {\varphi }_t-{\varphi }_n\rangle \right] \end{gathered}$$

(1)

Therefore, the traction–separation law is calculated

$$\begin{gathered} T_n\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=\frac{\partial \Psi }{\partial {\mathsf{\Delta }}_n}=\frac{{\Gamma }_n}{{\delta }_n}\left[m{\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha }{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{m-1}-\alpha {\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha -1}{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^m\right] \\ \left[{\Gamma }_t{\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta }{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^n+\langle {\varphi }_t-{\varphi }_n\rangle \right] \end{gathered}$$

(2)

$$\begin{gathered} T_t\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=\frac{\partial \Psi }{\partial {\mathsf{\Delta }}_t}=\frac{{\Gamma }_t}{{\delta }_t}\left[n{\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta }{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{n-1}-\beta {\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta -1}{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^n\right] \\ \left[{\Gamma }_n{\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha }{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^m+\langle {\varphi }_n-{\varphi }_t\rangle \right]\frac{{\mathsf{\Delta }}_t}{\left|{\mathsf{\Delta }}_t\right|} \end{gathered}$$

(3)

where $\langle ·\rangle$ is the Macaulay bracket, i.e., if $x\le 0$, then $\langle x\rangle =0,$ otherwise $\langle x\rangle =x$.

There are eight basic parameters (${\varphi }_n$, ${\varphi }_t$, ${\sigma }_{max}$, ${\tau }_{max}$, $\alpha$, $\beta$, ${\lambda }_n$, and ${\lambda }_t$) involved in the PPR model [32]. The PPR model considers different normal and tangential fracture energies (${\varphi }_n$ and ${\varphi }_t$), different cohesive strengths (${\sigma }_{max}$ and ${\tau }_{max}$), and controls the shape of the traction–separation law using the parameters $\alpha$, and $\beta$ and the initial slope indicators ${\lambda }_n$, and ${\lambda }_t$. The influence of $\alpha$, $\beta$, ${\lambda }_n$, and ${\lambda }_t$ on the material softening response were introduced in detail in [32].

These eight parameters could be obtained by fitting the interface stress—displacement relation measured in the splitting and shearing model test of concrete and mortar bonded composite specimens [40]. From these eight parameters, the following quantities can be deduced, which are used in (1), (2), and (3):

$$m=\frac{\alpha \left(\alpha -1\right){\lambda }_n^2}{\left(1-\alpha {\lambda }_n^2\right)},n=\frac{\beta \left(\beta -1\right){\lambda }_t^2}{\left(1-\beta {\lambda }_t^2\right)}$$

(4)

$${\Gamma }_n={\left(-{\varphi }_n\right)}^{{\varphi }_n-{\varphi }_t/\left({\varphi }_n-{\varphi }_t\right)}{\left(\frac{\alpha }{m}\right)}^m,{\Gamma }_t={\left(-{\varphi }_t\right)}^{{\varphi }_t-{\varphi }_n/\left({\varphi }_t-{\varphi }_n\right)}{\left(\frac{\beta }{n}\right)}^n\mathrm{for}\left({\varphi }_n\ne {\varphi }_t\right)$$

(5)

$${\Gamma }_n=-{\varphi }_n{\left(\frac{\alpha }{m}\right)}^m,{\Gamma }_t={\left(\frac{\beta }{n}\right)}^n\mathrm{for}\left({\varphi }_n={\varphi }_t\right)$$

(6)

$${\delta }_n=\frac{{\varphi }_n}{{\sigma }_{max}}\alpha {\lambda }_n{\left(1-{\lambda }_n\right)}^{\alpha -1}\left(\frac{\alpha }{m}+1\right){\left(\frac{\alpha }{m}{\lambda }_n+1\right)}^{m-1}$$

(7)

$${\delta }_t=\frac{{\varphi }_t}{{\tau }_{max}}\beta {\lambda }_t{\left(1-{\lambda }_t\right)}^{\beta -1}\left(\frac{\beta }{n}+1\right){\left(\frac{\beta }{n}{\lambda }_t+1\right)}^{n-1}$$

(8)

where ${\delta }_n$ and ${\delta }_t$ are the normal final crack opening width and the tangential final crack opening width, respectively. If ${\mathsf{\Delta }}_n\ge {\delta }_n$ or ${\mathsf{\Delta }}_t\ge {\delta }_t$, the tractions $T_n$ and $T_t$ are set to zero. Therefore, the traction–separation law is only valid in a region. To keep things simple, the separations (${\mathsf{\Delta }}_n$, ${\mathsf{\Delta }}_t$) are assumed to be positive here. Then, the region can be expressed as $\left[\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)|0\le {\mathsf{\Delta }}_n\le {\delta }_n,0\le {\mathsf{\Delta }}_t\le {\delta }_t\right]$. Considering the region, the normal and tangential cohesive tractions of the PPR model are plotted in Figure 2 with different fracture energies (e.g., ${\varphi }_n=100 \mathrm{N}/\mathrm{m}$, ${\varphi }_t=200 \mathrm{N}/\mathrm{m}$, and other cases), cohesive strengths (e.g., ${\sigma }_{max}=40 \mathrm{MPa}$, ${\tau }_{max}=30 \mathrm{MPa}$), shapes (e.g., $\alpha =5$, $\beta =1.3$), and initial slope indicators (e.g., ${\lambda }_n=0.1$, ${\lambda }_t=0.2$).

The normal cohesive traction (on the left in Figure 2) illustrates the fracture behavior of a typical quasi-brittle material, while the tangential cohesive traction (on the right in Figure 2) describes a plateau-type behavior. If ${\varphi }_n<{\varphi }_t$ (Figure 2a,b), the tangential cohesive traction was properly defined in the rectangular region corresponding to the final crack opening widths (${\delta }_n$,${\delta }_t$) as mentioned above, while in the same region, the normal cohesive traction $T_n\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)$ existed as negative (Figure 2a), which is contradictory to the nature of cohesive tractions. Similarly, if ${\varphi }_n>{\varphi }_t$, the normal cohesive traction was properly defined in the rectangular region, while the tangential cohesive traction was negative in some areas, as illustrated in Figure 2c,d. If ${\varphi }_n={\varphi }_t$ (Figure 2e,f), the normal and tangential tractions were non-negative in the same region.

To prevent the unphysical response, Park et al. [32] redefined the region by narrowing it to make the cohesive traction non-negative in new region, and the traction was set to zero if it was out of the new region. Taking ${\varphi }_n<{\varphi }_t$ as an example, the change of region for the normal traction is demonstrated in Figure 3 (separations are assumed positive here). The parameter ${\overline{\mathsf{\delta }}}_{\mathrm{t}}$ in Figure 3 is the tangential conjugate final crack opening width, and it is the single root of ${\Gamma }_t{\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta }{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^n+\langle {\varphi }_t-{\varphi }_n\rangle =0$ between 0 and ${\delta }_t$ [32].

For the new cohesive interaction region (on the right in Figure 3), one border of the new region is the normal final crack opening width ${\delta }_n$. The other border is the tangential conjugate final crack opening width ${\overline{\mathsf{\delta }}}_{\mathrm{t}}$. Due to ${\overline{\mathsf{\delta }}}_{\mathrm{t}}<{\delta }_t$, the new region was smaller than the original one $\left[\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)|0\le {\mathsf{\Delta }}_n\le {\delta }_n,0\le {\mathsf{\Delta }}_t\le {\delta }_t\right]$ (on the left in Figure 3), whereas the region of the tangential traction is the original one as shown in Figure 2b when ${\varphi }_n<{\varphi }_t$. This means the cohesive interaction regions of the normal and tangential tractions are different, and the tangential traction may still be large, while the normal traction has vanished in some regions. In other words, when a fracture happens, the normal and tangential tractions will not fail simultaneously. This is unrealistic for most interfaces encountered in engineering practice.

2.2. Unloading/Reloading Relationship

The original unloading/reloading relationship, which was commonly used with the PPR model, was linear to the origin [35], and expressed as follows

$$T_n^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_n\left({\mathsf{\Delta }}_n^{max},{\mathsf{\Delta }}_t\right)\frac{{\mathsf{\Delta }}_n}{{\mathsf{\Delta }}_n^{max}}$$

(9)

$$T_t^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_t\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t^{max}\right)\frac{{\mathsf{\Delta }}_t}{{\mathsf{\Delta }}_t^{max}}$$

(10)

where ${\mathsf{\Delta }}_n^{max}$ and ${\mathsf{\Delta }}_t^{max}$ are the largest values of ${\mathsf{\Delta }}_n$ and ${\mathsf{\Delta }}_t$ reached so far. If ${\mathsf{\Delta }}_n<{\delta }_n^{peak}$ (respective ${\mathsf{\Delta }}_t<{\delta }_t^{peak}$), with ${\delta }_n^{peak}={\lambda }_n{\delta }_n$ (resp. ${\delta }_t^{peak}={\lambda }_t{\delta }_t$), then ${\mathsf{\Delta }}_n^{max}=0$ (resp. ${\mathsf{\Delta }}_t^{max}=0$), and $T_n^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_n\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)$ (resp. $T_t^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_t\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)$). If ${\mathsf{\Delta }}_n\ge {\delta }_n^{peak}$ (resp. ${\mathsf{\Delta }}_t\ge {\delta }_t^{peak}$), then ${\mathsf{\Delta }}_n^{max}={\mathsf{\Delta }}_n$ (resp. ${\mathsf{\Delta }}_t^{max}={\mathsf{\Delta }}_t$). That is to say, the original unloading/reloading relationship is activated when the normal or tangential separation is past the peak cohesive strength.

Spring et al. [38] found that the original unloading/reloading relationship was not thermodynamically consistent and produced self-healing behavior. To address this issue, a new coupled unloading/reloading relationship was proposed.

$$T_n^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_n\left({\mathsf{\Delta }}_n^{max},{\mathsf{\Delta }}_t^{max}\right)\frac{{\mathsf{\Delta }}_n}{{\mathsf{\Delta }}_n^{max}}$$

(11)

$$T_n^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_n\left({\mathsf{\Delta }}_n^{max},{\mathsf{\Delta }}_t^{max}\right)\frac{{\mathsf{\Delta }}_n}{{\mathsf{\Delta }}_n^{max}}$$

(12)

where ${\mathsf{\Delta }}_n^{max}$ and ${\mathsf{\Delta }}_t^{max}$ are updated as soon as ${\mathsf{\Delta }}_n>0$ and ${\mathsf{\Delta }}_t>0$. This means the linear unloading/reloading response applies even before any peak has been passed.

Gilormini et al. [39] compared the two unloading/reloading relationships. They demonstrated that the new unloading/reloading relationship performed better than the original one and did not have the above questionable features. However, they also indicated that the new one did not include an initial elastic region, since the energy was dissipated by increasing the damage from the very beginning of the loading process. To address this issue, our paper improves the unloading/reloading relationship (see Section 4).

3. Simplified PPR Traction–Separation Law

The traction–separation law of the PPR model is adjusted here to avoid the issues mentioned in Section 2.1. The modifications of the traction–separation law are interpreted below. Then, based on previous studies [32], the path dependence of work-of-separation is investigated with respect to proportional and non-proportional paths to demonstrate the consistency of the simplified PPR traction–separation law. Finally, the simplified model was verified by simulating a mixed-mode bending test and comparing with the original model.

3.1. Modification

From Figure 2, we concluded that the cohesive interaction regions for the normal and tangential tractions were the same only if ${\varphi }_n={\varphi }_t$. Substituting ${\varphi }_n={\varphi }_t$ into Equations (2) and (3), we obtain the traction–separation law as follows.

$$\begin{gathered} T_n\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=\frac{{\varphi }_n}{{\delta }_n}{\left(\frac{\alpha }{m}\right)}^m{\left(\frac{\beta }{n}\right)}^n{\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta }{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^n \\ \left[\alpha {\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha -1}{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^m-m{\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha }{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{m-1}\right] \end{gathered}$$

(13)

$$\begin{gathered} T_t\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=\frac{{\varphi }_n}{{\delta }_t}{\left(\frac{\alpha }{m}\right)}^m{\left(\frac{\beta }{n}\right)}^n{\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha }{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^m \\ \left[\beta {\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta -1}{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^n-n{\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta }{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{n-1}\right]\frac{{\mathsf{\Delta }}_t}{\left|{\mathsf{\Delta }}_t\right|}. \end{gathered}$$

(14)

The traction–separation law only depends on mode I fracture energy ${\varphi }_n$. To account for different values of ${\varphi }_n$ and ${\varphi }_t$, the mode II fracture energy ${\varphi }_t$ is substituted for ${\varphi }_n$ in the equation for the tangential traction. Therefore, the final form of simplified PPR traction–separation law is given by

$$\begin{gathered} T_n\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=\frac{{\varphi }_n}{{\delta }_n}{\left(\frac{\alpha }{m}\right)}^m{\left(\frac{\beta }{n}\right)}^n{\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta }{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^n \\ \left[\alpha {\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha -1}{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^m-m{\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha }{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{m-1}\right] \end{gathered}$$

(15)

$$\begin{gathered} T_t\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=\frac{{\varphi }_t}{{\delta }_t}{\left(\frac{\alpha }{m}\right)}^m{\left(\frac{\beta }{n}\right)}^n{\left(1-\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^{\alpha }{\left(\frac{m}{\alpha }+\frac{{\mathsf{\Delta }}_n}{{\delta }_n}\right)}^m \\ \left[\beta {\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta -1}{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^n-n{\left(1-\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{\beta }{\left(\frac{n}{\beta }+\frac{\left|{\mathsf{\Delta }}_t\right|}{{\delta }_t}\right)}^{n-1}\right]\frac{{\mathsf{\Delta }}_t}{\left|{\mathsf{\Delta }}_t\right|} \end{gathered}$$

(16)

The simplified PPR traction–separation law is similar to the original PPR model and can also consider different fracture energies, cohesive strengths, and various material softening behaviors. The noteworthy merits of the simplified model are that the energy constants ${\Gamma }_n$ and ${\Gamma }_t$ are omitted (other parameters are the same as the original model), and the formulas are unified regardless of what the fracture energies are. Taking ${\varphi }_n=100 \mathrm{N}/\mathrm{m}$ and ${\varphi }_t=200 \mathrm{N}/\mathrm{m}$ as an example, the normal and tangential cohesive tractions of the simplified model are plotted in Figure 4. Figure 4 shows that the normal and tangential tractions are both properly defined in the same regions as expected. In the following section, the applicability of the simplified model is demonstrated using multiple cases.

3.2. Path Dependence of Work-of-Separation

The analysis of work-of-separation is a way to study the behavior of a coupled cohesive zone model [13,32,41]. In this paper, we compare the work-of-separation of the simplified PPR traction–separation law (SPPR) with the original PPR model for proportional separation paths and non-proportional paths. The fracture parameters in [32] were utilized in this investigation: ${\varphi }_n=100 \mathrm{N}/\mathrm{m}$, ${\varphi }_t=200 \mathrm{N}/\mathrm{m}$, ${\sigma }_{max}=3 \mathrm{MPa}$, ${\tau }_{max}=12 \mathrm{MPa}$, $\alpha =3$, $\beta =3$, ${\lambda }_n=0.01$, and ${\lambda }_t=0.01$.

3.2.1. Proportional Separation

The proportional separation path is shown in Figure 5. The variable $\theta$ in Figure 5 is the separation angle between the path direction and tangent, and ${\mathsf{\Delta }}_r$ is the separation for the proportional path. With the increase in ${\mathsf{\Delta }}_r$, the interface gradually debonds. The work-of-separation is calculated with the following expression [32].

$$W_{sep}={\int }_0^{{\delta }_r}{T}_n\left({\mathsf{\Delta }}_{r}sin\theta ,{\mathsf{\Delta }}_{r}cos\theta \right)sin\theta d{\mathsf{\Delta }}_r+{\int }_0^{{\delta }_r}{T}_t\left({\mathsf{\Delta }}_{r}sin\theta ,{\mathsf{\Delta }}_{r}cos\theta \right)cos\theta d{\mathsf{\Delta }}_r$$

(17)

where ${\delta }_r=\sqrt{{\delta }_n^2+{\delta }_t^2}$. The first term in the work-of-separation expression is the work conducted by the normal traction ($W_n$), and the second term in the expression is the work conducted by the tangential traction ($W_t$). $W_{sep}=W_n={\varphi }_n$ when the separation angle $\theta$ is 90°. When $\theta =0°$, the work-of-separation $W_{sep}$ and $W_t$ are the same as the mode II fracture energy ${\varphi }_t$.

Figure 6 illustrates the variation of $W_{sep}$, $W_n,$ and $W_t$ with respect to the separation angles. The results for the PPR model are on the left and for the SPPR model are on the right. The changing laws of $W_{sep}$, $W_n,$ and $W_t,$ with respect to the separation angles for different models, are the same. Especially, when the separation angle is 0° or 90°, the curves for the SPPR model are exactly the same as the PPR model.

If $\theta$ is equal to 0°, $W_{sep}$ and $W_t$ increase from 0 to the mode II fracture energy (200 N/m) with the increase in ${\mathsf{\Delta }}_r$, while $W_n$ remains zero. When $\theta$ is equal to 90°, $W_{sep}$ and $W_n$ reach the mode I fracture energy (100 N/m), and $W_t$ stays at zero. For the intermediate angles ($0°<\theta <90°$), the $W_{sep}$, $W_n,$ and $W_t$ of both models change monotonically with respect to the increase in the separation angle $\theta$. These verify that the PPR model and SPPR models both guarantee the consistency of the cohesive constitutive model.

There is a difference between the PPR model and the SPPR model. When $0°<\theta <90°$, the work conducted by the normal traction $W_n$ for the PPR model only has a small change with increases in the separation angle. In contrast, the SPPR model has a more obvious and uniform change within the whole separation angles. This is due to the fact that the cohesive interaction region for normal traction of the PPR model is smaller than the SPPR model here (${\varphi }_n<{\varphi }_t$), leading to a smaller $W_n$ for the PPR model under mixed-mode fracture conditions.

3.2.2. Non-Proportional Separation

The non-proportional separation path is shown in Figure 7. Path 1 is that the interface is loaded in the normal direction until ${\mathsf{\Delta }}_n={\mathsf{\Delta }}_{n,max}$; then, complete tangential separation occurs. Accordingly, path 2 is that the interface is first loaded in shear up to ${\mathsf{\Delta }}_{t,max}$, and then completely broken in the normal direction [41]. The expressions of the work-of-separation for the two paths were given by [32]:

$$W_{sep}={\int }_0^{{\mathsf{\Delta }}_{n,max}}{T}_n\left({\mathsf{\Delta }}_n,0\right)d{\mathsf{\Delta }}_n+{\int }_0^{{\delta }_t}{T}_t\left({\mathsf{\Delta }}_{n,max},{\mathsf{\Delta }}_t\right)d{\mathsf{\Delta }}_t$$

(18)

$$W_{sep}={\int }_0^{{\mathsf{\Delta }}_{t,max}}{T}_t\left(0,{\mathsf{\Delta }}_t\right)d{\mathsf{\Delta }}_t+{\int }_0^{{\delta }_n}{T}_n\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_{t,max}\right)d{\mathsf{\Delta }}_n$$

(19)

For the first path (Figure 7a), ${\mathsf{\Delta }}_{n,max}=0$ represents the pure mode II fracture, while ${\mathsf{\Delta }}_{n,max}={\delta }_n$ describes the pure mode I fracture. Similarly, for the second path (Figure 7b), when ${\mathsf{\Delta }}_{t,max}$ is zero, the separation path illustrates the pure mode I failure, while ${\mathsf{\Delta }}_{t,max}={\delta }_t$ represents the pure mode II fracture. The change of ${\mathsf{\Delta }}_{t,max}$ from 0 to ${\delta }_t$ (resp. ${\mathsf{\Delta }}_{n,max}$ from 0 to ${\delta }_n$) demonstrates the gradual change of the mode mixity from the mode I fracture to the mode II fracture (resp. from the mode II fracture to the mode I fracture). Based on Equations (18) and (19), the work-of-separation may change with the increasing of ${\mathsf{\Delta }}_{n,max}$ or ${\mathsf{\Delta }}_{t,max}$. If the work-of-separation has a monotonic variation from one fracture mode to the other fracture mode, this demonstrates the consistency of the cohesive constitutive model [32,41].

Figure 8 shows the variation of $W_{sep}$, $W_n,$ and $W_t$ with respect to the two paths, under the condition of ${\varphi }_n<{\varphi }_t$. The results for PPR model are on the left and for the SPPR model are on the right. $W_{sep}$, $W_n$, and $W_t$ all change monotonically for both models. For path 1 (Figure 8a,b), the curves of $W_{sep}$, $W_{n,}$ and $W_t$ for the SPPR model are exactly the same as the PPR model.

Figure 8a,b show that the work conducted by the tangential traction $W_t$ gradually decreases from ${\varphi }_t$ to 0, while the work conducted by the normal traction $W_n$ increases from 0 to ${\varphi }_n$. The work-of-separation $W_{sep}$ is the sum of $W_n$ and $W_t$, and this monotonically varies from the value of ${\varphi }_t$ to the value of ${\varphi }_n$ by increasing ${\mathsf{\Delta }}_{n,max}$ from 0 to ${\delta }_n$. For path 2 (Figure 8c,d), the change rules of $W_{sep}$, $W_n,$ and $W_t$ are the exact opposite of those in path 1. There is a kink point on the curves of $W_n$ and $W_{sep},$ as shown in Figure 8c, but not in Figure 8d.

The separation at the kink point corresponds to the border ${\mathsf{\Delta }}_t={\overline{\delta }}_t$ of the original PPR model, where ${\overline{\delta }}_t$ is the tangential conjugate final crack opening width as previously described in Section 2.1. When ${\mathsf{\Delta }}_t$ is smaller than ${\overline{\delta }}_t$, the normal cohesive interaction is obtained based on Equation (2). When ${\mathsf{\Delta }}_t$ is greater than ${\overline{\delta }}_t$, the normal traction is set to zero. The normal cohesive interaction is then not smooth but piece-wise continuous at ${\mathsf{\Delta }}_t={\overline{\delta }}_t$. As a result, the $W_n$ and $W_{sep}$ also have the kink point at the same location. In contrast, the curves of $W_{sep}$, $W_n,$ and $W_t$ for the SPPR model, as shown in Figure 8d, are continuous and smooth. This is because both the normal and tangential cohesive interactions for the SPPR model are continuous and smooth in the region $\left[\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)|0\le {\mathsf{\Delta }}_n\le {\delta }_n,0\le {\mathsf{\Delta }}_t\le {\delta }_t\right]$. This indicates that the SPPR model describes the mixed-mode fracture better.

Additionally, the same conclusion can be reached when the mode I fracture energy is greater than the mode II fracture energy as shown in Figure 9. For the PPR model, the kink point occurs in path 1 because the tangential cohesive interaction is piece-wise continuous while being continuous and smooth for the SPPR model.

3.3. Mixed-Mode Bending (MMB) Test Verification

The simplified PPR traction–separation law is verified here and compared to the original PPR model by simulating the mixed-mode bending (MMB) test. The MMB test has been widely used to validate the applicability of CZM for mixed-mode fracture [37]. The configuration of the test is shown in Figure 10. Following the geometry parameters of the MMB test, specimens were considered: L = 51 mm, h = 1.56 mm, a₀ = 33.7 mm, c = 60 mm, and B = 25.4 mm.

Numerical simulations of the mixed-mode fracture were implemented using the commercial software ABAQUS with a user-defined element (UEL) subroutine, and such a subroutine for the PPR model was given in the work of [34]. In this paper, user element subroutines (UEL) were also utilized to implement the simplified PPR traction–separation law. Since the mesh, FE element type, boundary conditions, as well as the solving method are all the same as in [34], those items are not covered again here.

In this study, two cases were tested, one with the same fracture energy (${\varphi }_n={\varphi }_t=1 \mathrm{N}/\mathrm{m}$) and another with different fracture energies (${\varphi }_n=1 \mathrm{N}/\mathrm{m} and {\varphi }_t=2 \mathrm{N}/\mathrm{m}$). The cohesive strength ${\sigma }_{max}={\tau }_{max}=200 \mathrm{MPa}$, shape parameter $\alpha =\beta =3$, and the initial slope indicator ${\lambda }_n={\lambda }_t=0.02$ were the same for both cases. The numerical results were compared to the analytical solution given in [32].

For the same fracture energy, the computational results for different models are illustrated in Figure 11a. The results for the SPPR model and PPR model were the same, and coincided with the analytical solutions. For the case of different fracture energies (Figure 11b), the computational results for the SPPR model were in better agreement with analytical solution compared with the PPR model under the same conditions. The results for the PPR model were relatively small, as shown in Figure 11b. The reason is that the effective region for the PPR model is smaller than for the SPPR model when ${\varphi }_n\ne {\varphi }_t$, leading to the smaller tractions and energies under mixed-mode fractures as mentioned before.

4. Improved Unloading/Reloading Relationship

Previous studies [38,39] demonstrated that the original unloading/reloading relationship was not thermodynamically consistent and produced self-healing behavior. In addition, the new unloading/reloading relationship proposed by Spring et al. [38] did not include the initial elastic region. To prevent these issues, an improved unloading/reloading relationship was developed. The modifications of the unloading/reloading relationship are interpreted below. Then, the comparison of the three models is presented in Section 4.2. For convenience in the presentation of the results, the original unloading/reloading relationship is referred to as model (i) here. The new unloading/reloading relationship developed in [38] is referred to as model (ii), while the improved one proposed in this paper is referred to as model (iii).

4.1. Modification

The reason why model (ii) has a lack of an initial elastic region is that the variables ${\mathsf{\Delta }}_n^{max}$ and ${\mathsf{\Delta }}_t^{max}$ in Equations (11) and (12) are updated at the very beginning. Referring to the definition of model (i), ${\mathsf{\Delta }}_n^{max}$ and ${\mathsf{\Delta }}_t^{max}$ should not be updated unless certain conditions are met, for example, the peak cohesive strength should be passed. Therefore, how to determine the peak becomes the key. For model (i), ${\delta }_n^{peak}$ and ${\delta }_t^{peak}$ are used. However, ${\delta }_n^{peak}$ and ${\delta }_t^{peak}$ are the separations corresponding to the peak cohesive strength under mode I and mode II fractures, respectively. Under the conditions of mixed-mode fracture, the separations are not ${\delta }_n^{peak}$ and ${\delta }_t^{peak}$, as illustrated in Figure 12. Figure 12 also shows that the peaks change with the variation of mode mixing. Thus, the separations corresponding to the peak under mixed-mode fractures are not convenient to obtain. For this, an alternative method is presented here to estimate the peak, which is based on the gradients of the tractions.

The improved unloading/reloading relationship (model (iii)) is expressed as

$$T_n^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_n\left({\mathsf{\Delta }}_n^{\chi },{\mathsf{\Delta }}_t^{\gamma }\right)\frac{{\mathsf{\Delta }}_n}{{\mathsf{\Delta }}_n^{\chi }}$$

(20)

$$T_t^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_t\left({\mathsf{\Delta }}_n^{\chi },{\mathsf{\Delta }}_t^{\gamma }\right)\frac{{\mathsf{\Delta }}_t}{{\mathsf{\Delta }}_t^{\gamma }}$$

(21)

where ${\mathsf{\Delta }}_n^{\chi }$ and ${\mathsf{\Delta }}_t^{\gamma }$ are state variables, and ${\mathsf{\Delta }}_n^{\chi }={\mathsf{\Delta }}_n$ and ${\mathsf{\Delta }}_t^{\gamma }={\mathsf{\Delta }}_t$ by default. This indicates that $T_n^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_n\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)$ and $T_t^{\upsilon }\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)=T_t\left({\mathsf{\Delta }}_n,{\mathsf{\Delta }}_t\right)$, until the following conditions are met: the gradients of tractions $\frac{T_n\left({\mathsf{\Delta }}_n^i,{\mathsf{\Delta }}_t^i\right)-T_n\left({\mathsf{\Delta }}_n^{i-1},{\mathsf{\Delta }}_t^{i-1}\right)}{{\mathsf{\Delta }}_n^i-{\mathsf{\Delta }}_n^{i-1}}\le 0$ or $\frac{T_t^{\upsilon }\left({\mathsf{\Delta }}_n^i,{\mathsf{\Delta }}_t^i\right)-T_t^{\upsilon }\left({\mathsf{\Delta }}_n^{i-1},{\mathsf{\Delta }}_t^{i-1}\right)}{{\mathsf{\Delta }}_t^i-{\mathsf{\Delta }}_t^{i-1}}\le 0$. Then, ${\mathsf{\Delta }}_n^{\chi }$ and ${\mathsf{\Delta }}_t^{\gamma }$ become the largest values of ${\mathsf{\Delta }}_n$ and ${\mathsf{\Delta }}_t$ reached so far. Since ${\mathsf{\Delta }}_n^i>{\mathsf{\Delta }}_n^{i-1}>0$ and ${\mathsf{\Delta }}_t^i>{\mathsf{\Delta }}_t^{i-1}>0$, ${\mathsf{\Delta }}_n^{\chi }={\mathsf{\Delta }}_n^i$ and ${\mathsf{\Delta }}_t^{\gamma }={\mathsf{\Delta }}_t^i$. Once one of above conditions is satisfied, both ${\mathsf{\Delta }}_n^{\chi }$ and ${\mathsf{\Delta }}_t^{\gamma }$ are updated. That is to say, the improved unloading/reloading relationship (model (iii)) is activated when one of above conditions is met.

4.2. Comparison

In this section, comparisons of the three models are drawn using the example in [39], where the following set of parameters is used: ${\varphi }_n=100 \mathrm{N}/\mathrm{m}$, ${\varphi }_t=300 \mathrm{N}/\mathrm{m}$, ${\sigma }_{max}=2 \mathrm{MPa}$, ${\tau }_{max}=4 \mathrm{MPa}$, $\alpha =3$, $\beta =5$, ${\lambda }_n=0.20$, and ${\lambda }_t=0.25$. The loading process consists of three steps. First, a proportional mixed-mode loading where ${\mathsf{\Delta }}_n={\mathsf{\Delta }}_t$ is applied up to a predefined value $\mathsf{\Delta }$. Then, a proportional mixed-mode unloading where ${\mathsf{\Delta }}_n={\mathsf{\Delta }}_t$ is carried out down to 0. Finally, a mode I reloading (keeping ${\mathsf{\Delta }}_t=0)$ is conducted.

Unlike the work in [39], the simplified PPR traction–separation law proposed in this paper is used here instead of the PPR model. Therefore, ${\overline{\mathsf{\delta }}}_{\mathrm{t}}$ is not used, and the fracture occurs for either ${\mathsf{\Delta }}_n={\delta }_n$ or for ${\mathsf{\Delta }}_t={\delta }_t$. Based on the parameters above, ${\delta }_n=0.099 \mathrm{mm}$, ${\delta }_t=0.171 \mathrm{mm}$, ${\delta }_n^{peak}=0.020 \mathrm{mm},$ and ${\delta }_t^{peak}=0.043 \mathrm{mm}$.

Figure 13 shows how the dissipated energy for the three models change with the increase in the proportional loading amplitude $\mathsf{\Delta }$. The three models give the same energy values at the beginning and end. The variation of the energy value for model (i) is quite different from model (ii) and model (iii), while model (ii) and model (iii) are almost the same. The dissipated energy given by model (iii) is identical to model (i) when $\mathsf{\Delta }<0.017\mathrm{mm}$, which is a constant equal to the mode I fracture energy of ${\varphi }_n=100N/m$.

When $\mathsf{\Delta }\ge 0.017\mathrm{mm}$, the change law of the energy value for model (iii) is exactly the same as model (ii). There are three discontinuities for model (i), at $\mathsf{\Delta }={\delta }_n^{peak}=0.020 \mathrm{mm}$, $\mathsf{\Delta }={\delta }_t^{peak}=0.043 mm$, and $\mathsf{\Delta }={\delta }_n=0.099 \mathrm{mm}$, whereas model (ii) gives a smooth continuous curve, and model (iii) only has one discontinuity at $\mathsf{\Delta }=0.017 \mathrm{mm}$. These differences are explained below by analyzing the changes of the traction components during the loading process.

First, consider the proportional loading amplitudes $\mathsf{\Delta }$ around 0.017 mm. Figure 14 presents the variations of the traction components during the loading process for a proportional loading amplitude of $\mathsf{\Delta }=0.016 \mathrm{mm}$. The computation results for model (i) and model (iii) are identical; thus, both are presented using Figure 14a. As can be observed in Figure 14a, during unloading, both traction values back up along the same curves that they followed during loading, when the peak of tractions was not reached.

Consequently, the final energy is only dissipated in pure mode I reloading, which is equal to ${\varphi }_n=100 \mathrm{N}/\mathrm{m}$. In contrast, both unloading curves follow straight lines with the model (ii), as shown in Figure 14b. This is due to no peak value being needed to start using the linear response. That is say, damage is assumed to occur at the beginning, and the initial elastic region is ignored. Due to the assumed damage, the energy dissipated in pure mode I reloading for model (ii) is smaller than for model (i) and model (iii). As a result, the total dissipated energy (97.0 J/m²) for model (ii) is lower than ${\varphi }_n$.

When $\mathsf{\Delta }=0.017 \mathrm{mm}$, the peak of normal traction is reached under the mixed-mode loading (Figure 15a,b), and therefore, the linear unloading response of model (iii) is activated. As a consequence, the variations of the traction components during the whole loading process for model (iii) become the same as for model (ii), both presented in Figure 15b. They still apply for larger proportional loading amplitudes $\mathsf{\Delta }>0.017 \mathrm{mm}$.

Thus, in the following analysis, the results for model (ii) and model (iii) are all displayed with the same diagrams. Due to the variation of response, the dissipated energy for model (iii) changes from 100 J/m² (${\varphi }_n$) at $\mathsf{\Delta }=0.016 \mathrm{mm}$ to 97.0 J/m² at $\mathsf{\Delta }=0.017 \mathrm{mm}$. Such an energy discontinuity is inherent to any CZM model that obeys a curved line in the reversible range and an unloading straight line when irreversibility has appeared, which was discussed in detail by Gilormini et al. [39]. Therefore, the small energy jump is accepted.

When $\mathsf{\Delta }=0.019 \mathrm{mm}$, the peak of normal traction is exceeded, as shown in Figure 16a,b. However, because of $\mathsf{\Delta }<{\delta }_n^{peak}=0.020\mathrm{mm}$, the normal traction for model (i) still returns along the loading path during the unloading process, leading to a questionable response that the traction increases with the decrease in separations. When $\mathsf{\Delta }=0.020 \mathrm{mm}$, the ${\delta }_n^{peak}$ value is reached, and thus the model (i) is activated. Similar to that for model (iii) at $\mathsf{\Delta }=0.017 \mathrm{mm}$ mentioned above, there is a small energy jump for model (i) due to the change from an elastic region to a softening region. In contrast, for model (ii) and model (iii), the dissipated energy varies continuously.

Consider now the proportional loading amplitudes $\mathsf{\Delta }$ around ${\delta }_t^{peak}=0.043 \mathrm{mm}$. Figure 17a, for $\mathsf{\Delta }=0.042\mathrm{mm}$, shows that tangential traction component $T_t$ still returns along the loading path and increases significantly during the unloading process. When $\mathsf{\Delta }=0.043 \mathrm{mm}$ (Figure 17c), the ${\delta }_t^{peak}$ value is reached, and therefore, the tangential unloading component in model (i) is activated as well. On account of the added energy that is dissipated by $T_t$ during the proportional loading/unloading process, the total energy given by model (i) increases sharply, which induces a jump at $\mathsf{\Delta }=0.043 \mathrm{mm}$ in Figure 13. In contrast, model (ii) and model (iii) have a smooth evolution of dissipated energy, as can be observed in Figure 17b,d.

Finally, consider the proportional loading amplitudes $\mathsf{\Delta }$ around ${\delta }_n=0.099 \mathrm{mm}$. When $\mathsf{\Delta }=0.098\mathrm{mm}$ (Figure 18a,b), which is slightly below the critical value ${\delta }_n$, both tractions are near 0 at the end of loading. For model (i), there is an increase in $T_t$ during unloading, and, due to that, ${\mathsf{\Delta }}_n$ in the tangential term of model (i) varies during proportional unloading. When the proportional loading amplitude $\mathsf{\Delta }=0.099\mathrm{mm}$, the critical value ${\delta }_n$ is reached, and hence fracture is complete (Figure 18a,b). As a result, the unloading and reloading phases no longer exist, and the dissipated energy for the three models becomes the same, equal to 155.2 J/m².

From the above analysis, the original unloading/reloading relationship (referred to as model (i) here) may induce some questionable responses, such as increasing traction during unloading. The new unloading/reloading relationship proposed by Spring et al. [38] (referred to as model (ii) here) may bring damage at the beginning and ignore the initial elastic region. The improved unloading/reloading relationship (referred to as model (iii) here) proposed in this paper combines the merits of the above two models, which prevents the issues mentioned above, and defines an elastic region before a softening regime.

5. Application

Interface damage, which even happens in the construction phase, has become a major problem for China Railway Track System (CRTS-II) slab track. To reveal the behavior of the slab track under the difference of temperatures, the effect of daily changing temperature on the curling behavior and interface stress of slab track in the construction stage was researched by the authors [2]. As a follow-up study, the interface damage of slab track under daily changing temperature is analyzed by implementing the improved cohesive zone model in this section.

The CRTS-II slab track consists of precast slab, CA mortar, and concrete base, as shown in Figure 19. All of these components are modeled according to actual size. The dimensions, material properties, mesh, FE element type, and boundary conditions of each component are all the same as in [2]; those items are not covered again here.

Interface cracks usually occur between the track slab and CA mortar, as shown in Figure 1. The interlaminar cracking is modelled based on the constitutive model proposed in this paper, by using the commercial software ABAQUS with a user-defined interaction (UINTER) subroutine. The validated interface parameters [42] are ${\varphi }_n=2.6\mathrm{N}/\mathrm{m}$, ${\varphi }_t=4\mathrm{N}/\mathrm{m}$, ${\sigma }_{max}=0.015\mathrm{MPa}$, ${\tau }_{max}=0.015\mathrm{MPa}$, $\alpha =2$, $\beta =2$, ${\lambda }_n=0.1$ and ${\lambda }_t=0.1$. Due to symmetry of the geometry and loading conditions, only a quarter of the slab track is established. The 3-D finite element model of CRTS-II slab track is built as presented in Figure 20.

Based on the proposed model, the interface damage is simulated under gravity load and measured temperature. The measured temperature was input into the model as temperature load using user-defined subroutines named UTEMP [2]. In the analysis, the summer temperature (Figure 9 in [2]) is taken as an example. As the initial stress field has an influence on the stress history and stress level, the time of 14:30 with the maximum temperature difference is selected as the starting time.

Figure 21 shows the interface crack opening (COPEN) distribution of the slab track system as a result of temperature change. It is found that the damage at four corners is the most obvious. Such damage mode is exactly the same as that observed in high-speed railway lines.

The normal and the two shear stresses of the interface at the slab corner are shown in Figure 22. It can be observed that the interface damage is mainly caused by the presence of normal and lateral shear stresses. It is worth noting that the stresses change smoothly for the model proposed in the paper and PPR model, while piece-wise continuously for the cohesive zone model in ABAQUS. Moreover, the problem of self-repair for the PPR model is found in Figure 22b,c. The cause was mentioned before.

Figure 23 shows the interface normal stresses (CPRESS) between the slab and CA mortar layer when interface cracking happens. It is found that the stress distribution for the model proposed in the paper and PPR model is almost the same, and continuously changes with location. However, that for the cohesive zone model in ABAQUS is rugged and unreasonable. For example, the tensile and compressive stresses occur simultaneously around slab corner. This may be due to the stress oscillation [33].

6. Conclusions

A simplified cohesive zone model combined with an improved unloading/reloading relationship was proposed in this paper to overcome certain shortcomings of the original model, and was validated using multiple cases.

First, the traction–separation laws of the PPR model under different conditions of fracture energies were compared. We concluded that the cohesive interaction regions for the normal and tangential traction components were different, when the mode I fracture energy was not equal to the mode II fracture energy. This may lead to an undesired response where one traction component is still very large while the other traction component has vanished, which is unrealistic for most interfaces encountered in civil engineering practice. To address this issue, the simplified PPR model was developed based on the original model. We found that the simplified model had unified formulas and cohesive interaction regions regardless of the fracture energies. The investigations of the path dependence of work-of-separation and the simulation of the mixed-mode bending test both demonstrated that the simplified model guaranteed the consistency of the cohesive constitutive model and had better performance in modeling the mixed-mode fracture.

When a loading/unloading/reloading process was applied, we observed that the original unloading/reloading relationship, which was commonly utilized with the PPR model, induced questionable responses, such as increasing the traction during unloading. The new unloading/reloading relationship proposed by Spring et al. [38] ignored the initial elastic region. By conducting an analysis of the above issues and the causes, the unloading/reloading relationship was improved based on the gradient of traction. We verified that the improved unloading/reloading relationship prevented the above issues and defined an elastic region before a softening regime.

The proposed model provides a tool for the research of the interface cracking mechanism of ballastless tracks. After the above analysis and verification, the proposed model solves the problem of ”self-repair” in the existing models and can correctly simulate the interface damages and cracking process under reciprocating loads. By using the UINTER platform of ABAQUS/Standard user interface constitutive subroutine, the module of interlaminar cracking analysis based on the constitutive model proposed in this paper could be constructed.

After coupling the module with the main structure model of ballastless track, a nonlinear finite element model of multilayer slab ballastless track system that could accurately simulate the interlayer compound mode cracking was constructed. Based on the model, the mechanism of interface cracking can be analyzed in detail [42]. The results of the research on the defect mechanism of the ballastless track can provide a scientific basis for the maintenance of the defects of ballastless tracks and guide the research of the monitoring of track service status, such as monitoring point placement and data analysis.

The proposed model could model the initiation and propagation of interface cracks under a coupled thermo-mechanical operating condition; however, it does not take into account the time/temperature dependency of the interfacial fracture parameters, which is regarded as our future work.