*2.1. Method*

Fracture mechanics were widely used in crack analysis for their high accuracy, and the development of computer technology offers the convenience of calculating. Thus fracture mechanics are commonly used in pavement engineering. Using the finite element method for digital modeling can calculate the behavior of cracks. Responses of cracks could be summarized by modifying mechanical parameters of finite element model such as working conditions, structure combination, material composition.

Fracture mechanics divide cracks into three species (open crack—type I, sliding crack— type II, tear-open crack—type III) and stress intensity factors (*KI*, *KII*, *KIII*) were used to describe the tendency of cracking [12]. Linear superposition of three kinds of cracks can be used to describe the stress field at the tip of cracks. Normal stress and shear stress of the tip of cracks can be calculated using the following equations.

$$\sigma\_x = \frac{K\_I}{\sqrt{2\pi r}} \cos\frac{\theta}{2} \left(1 - \sin\frac{\theta}{2} \sin\frac{3\theta}{2}\right) + \frac{K\_{II}}{\sqrt{2\pi r}} \sin\frac{\theta}{2} \left(-2 - \cos\frac{\theta}{2} \cos\frac{3\theta}{2}\right) \tag{1}$$

$$\sigma\_{\theta} = \frac{K\_I}{\sqrt{2\pi r}} \cos\frac{\theta}{2} \left( 1 + \sin\frac{\theta}{2} \sin\frac{3\theta}{2} \right) + \frac{K\_{II}}{\sqrt{2\pi r}} \sin\frac{\theta}{2} \cos\frac{\theta}{2} \cos\frac{3\theta}{2} \tag{2}$$

$$\pi\_{xy} = \frac{K\_I}{\sqrt{2\pi r}} \sin\frac{\theta}{2} \cos\frac{\theta}{2} \cos\frac{3\theta}{2} + \frac{K\_{II}}{\sqrt{2\pi r}} \cos\frac{\theta}{2} \left(1 - \sin\frac{\theta}{2} \sin\frac{3\theta}{2}\right) \tag{3}$$

$$\pi\_{\Sigma\Gamma} = -\frac{K\_{III}}{\sqrt{2\pi r}} \sin\frac{\theta}{2} \tag{4}$$

$$\pi\_{yz} = \frac{K\_{III}}{\sqrt{2\pi r}} \cos\frac{\theta}{2} \tag{5}$$

where *r* represents the distance from the calculation point to the crack tip. *θ* represents the angle between the calculation point and the crack tip in the polar coordinate system and the polar axis. Analysis of the above formula shows that the key to determining the stress field at the crack tip is the stress intensity factor. These factors are mainly determined by the properties of external load, the shape of cracks, and the geometry of elastomer, which can be expressed by the equation below.

$$\mathcal{K}\_{\mathfrak{m}} = \varprojlim\_{|r| \to 0} \sqrt{2\pi r} Z\_{\mathfrak{m}}(r) \tag{6}$$

where *Zm* (*r*) is an analytic function related to the boundary conditions, which can be calculated using finite element method software ABAQUS. The stress intensity factor at the ultimate fracture of a material is called fracture toughness, usually expressed by *KIC*, which can be determined by experiments.

Compared with the stress intensity factor, J-integral is more suitable for analyzing fracture problems under elastoplastic conditions [13]. The J-integral is used to solve the situation that a certain range of plastic zone appears at the tip of the crack under the elastoplastic condition, which makes the problem very complicated. The J-integral can describe the stress and strain field strength in the crack tip area, and it can be easily determined by experiment. Use stress intensity factor and J-integral can reveal the crack behavior in pavement engineering. A plane crack problem was assumed when calculating J-integral. The integral is performed around the crack tip, starting from the lower crack surface, and stopping counterclockwise to the crack's upper surface, thus forming an integral loop. This integral value obviously has nothing to do with the integral path, so J-integral can reflect the intensity of the stress field near the crack tip. Figure 1 is a schematic diagram of the calculation of J-integral.

**Figure 1.** Calculation of J-integral.

J-integral is defined as:

$$J = \int\_{\Gamma} \left( \omega dy - \stackrel{\rightarrow}{T} \frac{\partial \stackrel{\rightarrow}{u}}{\partial x} ds \right) \tag{7}$$

→

where ω represents the strain energy density on the loop. *T* is the stress vector at any point on the loop. →*u* is the displacement vector at any point on the loop. In the case of linear elasticity (plane strain), the strain energy density *ω* is

$$
\omega = \frac{1}{2} \sigma\_{\overline{i}\overline{j}} \varepsilon\_{\overline{i}\overline{j}} = \frac{1+\mu}{2E} \left[ (1-\mu) \left( \sigma\_{11}^2 + \sigma\_{22}^2 \right) - 2\mu \sigma\_{11} \sigma\_{22} + 2\sigma\_{12}^2 \right] \tag{8}
$$

where *<sup>σ</sup>ij* and *εij* respectively represent the stress tensor and strain tensor at the crack. Substitute (8) into (7), we have

$$
\omega = \frac{K\_1^2}{2\pi r} \frac{1+\mu}{E} \left[ \cos^2 \frac{\theta}{2} \left( 1 - 2\mu + \sin^2 \frac{\theta}{2} \right) \right] \tag{9}
$$

If we take a loop Γ with the crack tip as the center and a radius of *r*,

→

*T*

$$
\int\_{\Gamma} \omega dy = \int\_{-\pi}^{\pi} \omega r \cos \theta d\theta = \frac{K\_1^2 (1 + \mu)(1 - 2\mu)}{4E} \tag{10}
$$

As for stress vector → *T*

$$=\sigma\_{\vec{\imath}\vec{\jmath}}\eta\_{\vec{\jmath}}\tag{11}$$

So

$$\int\_{\Gamma} T\_i \frac{\partial u\_i}{\partial \mathbf{x}\_i} ds = \int\_{-\pi}^{\pi} \left( T\_1 \frac{\partial u\_1}{\partial \mathbf{x}\_1} + T\_2 \frac{\partial u\_2}{\partial \mathbf{x}\_2} \right) r d\theta = -\frac{-K\_1^2 (1 + \mu)(3 - 2\mu)}{4E} \tag{12}$$

Substituting (10) and (12) into (7),

$$J = \frac{1-\mu^2}{E} K\_1^2 = \frac{K\_1^2}{E'} = G\_1 \tag{13}$$

That is, the J-integral is the crack propagation energy release rate *G* in the linear elastic state. This is a constant and can be obtained through experiments or theoretical calculations. Therefore, the J-integral is used to establish the basis for the determination of cracks and is consistent with the establishment of the *K* value.

When using the traditional finite element method to simulate static non-propagating cracks, it is necessary to consider the singularity of the crack tip stress field according to the principles of fracture mechanics. Singular meshes were used to solve the problem, such as setting the meshes of the crack tip into M font or refining meshes in the crack tip area. Difficulties in computing were generated due to a large number of meshes, especially in three-dimensional models. It is hard to compute in crack tip areas due to the low quantity of meshes. Finite element method software ABAQUS offered Extended Finite Element Method (XFEM) to solve the problem, which uses shape functions to simulate crack areas in the pavement structure. Thus, it does not need to remesh and improves the efficiency of calculating. Based on the thought of unit decomposition, XFEM added functions that can reflect the discontinuous characteristics of the crack area [14,15]. The thought of unit decomposition thinks that any functions *ψ*(*x*) can be expressed using the equations below.

$$\psi(\mathbf{x}) = \sum\_{I} N\_{I}(\mathbf{x}) q\_{I} \Phi(\mathbf{x}), \ \sum\_{I} N\_{I}(\mathbf{x}) = 1 \tag{14}$$

where *qI* is a parameter to be adjusted to make the expression reach the best approximation, *NI*(*x*) is a function that satisfies the element decomposition, and *Φ*(*x*) is an extended function. As for crack simulation, XFEM uses the equations below to simulate.

$$u^k = \sum\_I N\_I(\mathbf{x}) u\_I + \sum\_I N\_I(\mathbf{x}) q\_I \Phi(\mathbf{x}) \tag{15}$$

where *qJ* is the newly added degree of freedom and it has no meaning in physics. The only usage of *qJ* is to adjust the function *Φ*(*x*) to achieve the best approximation. Compared with the traditional finite element method, XFEM added more freedom degrees, thus improving its accuracy.

XFEM and fracture mechanics were used in finite element method software ABAQUS, which can simulate the crack behavior during its growth with considerable accuracy. Crack behavior can be calculated by modifying properties of models and materials, which can make predictions on pavement life and different working conditions, thus saves much time and cost compared with local experiments.
