# Fatigue Life Prediction of the Zirconia Fixture Based on Boundary Element Method

Authors: Y. W. Wang; J. J. Ye; W. P. He; G. G. Cai; B. Q. Shi
Journal: Shock and Vibration
 (2020)

---

## Abstract

Zirconia grinding fixtures have been widely used in semiconductor industry to improve the quality and precision of the products. For maximizing the service life and minimizing the risks of accidental damage, it is critical to have a better understanding of the fatigue life of zirconia grinding fixtures. To this end, a boundary element method is developed in this paper to investigate their crack growth and fatigue life. To validate the proposed method, the stress intensity factor of a typical plate structure with initial cracks is considered. On this basis, Paris Law is employed in the boundary element model to further study the crack growth and stress distributions in the zirconia fixture under cyclic loads. Numerical results show that stress concentration occurs at the pillar of the fixture, and crack growth is perpendicular to the loading direction.


---

## Body

## 1. Introduction

In semiconductor industry, both grinding and polishing technologies play important roles in the procedures of product manufacturing. In order to improve the quality and precision of products, it is critical to choose an ideal auxiliary fixture to control their surface profiles [1, 2]. Due to the properties of high strength, low density, and high wear and corrosion resistances, zirconia has been one of the most chosen materials to manufacture the fixture [3, 4]. However, a zirconia fixture is very susceptible to microscopic cracks and damage due to its high brittleness and low fracture toughness, which will further affect the quality of grinding or polishing surfaces [5]. To insure the product quality and improve service life of zirconia fixtures, investigations of failure mechanisms and fatigue life prediction have drawn great attention in scientific community.Finite element method (FEM) is considered to be an effective method in investigating the fracture properties [6–9]. The basic idea of a FEM approach is to discretize a continuous solution domain into many subdomains that are connected by nodes. Moreover, unknown variables in the solution domain are expressed by approximate functions. Many researchers have employed FEM to solve fracture problems [10], including fatigue crack growth [11–13] in engineering structures. Chen et al. studied the impact of incision angle and depth on the strength [14] and crack propagation [15] of zirconia ceramic. The results show that the effect of notch angles on crack propagation increased with notch depth. El-sayed et al. [16] introduced a three-dimensional explicit FEM combined with the critical plane approach and multiaxial low cyclic fatigue mode to predict crack orientation and fatigue crack extension. Xue [17] employed FEM to analyze the crack growth mechanisms in the process of zirconia ceramics cutting, which was used for optimization design of process parameters. However, when a FEM is used, crack singularity is always a challenging issue and causes inefficiency in numerical simulations. This is due to the fact that an extreme refined mesh in crack regions is required to acquire a high calculation accuracy, which is extremely computationally expensive and may be impracticable.An alternative approach to study fracture of materials is to use boundary element method (BEM). Compared with the traditional FEM, the BEM provides a higher computational efficiency and accuracy in solving stress concentration problems. BEM has been extensively applied for the solution of stress intensity factor (SIF) [18, 19] and crack singularity problems [20, 21]. Ming and Yong-yuan [22] developed a time-domain boundary element method in studying dynamic SIF of a 3D crack. Excellent numerical results were acquired by employing isoparametric elements and singular elements. Peixoto et al. [23] presented a new boundary element formulation for quasi-brittle failure analysis. The boundary element method was extended to explore the quasi-brittle material fracture. Wang and Yao [24] proposed a fast BEM for the thermal analysis of fiber-reinforced composites. Numerical results clearly demonstrated the validity of the proposed model and its potential for large-scale analysis. Hilgendorff et al. [25] investigated irreversible damage accumulation in high-cycle fatigue regime by using BEM. Results show that the simulation of slip bands was in good agreement with experimental observations.In this paper, an extension of the BEM is presented to investigate fatigue crack propagations of the zirconia grinding fixture used in semiconductor industry. For each example, the open crack is considered. Stress intensity factor, which is related to the crack size and structural geometric characteristics, is accurately calculated. Moreover, a highlight of the crack evolutions under the fatigue load is studied. The outline of this paper is as follows: Section2 introduces the advantages of BEM in studying the crack singularity problem. Section 3 presents the theory of boundary element method and the calculation procedures of the SIF. In Section 4, two typical plate structures with cracks are investigated, and numerical results are compared with finite element method and theoretical solutions to verify the proposed method. In Section 5, a 3D model of zirconia fixture is established to study its fatigue life. The conclusions are summarized in Section 6.


## 2. Crack Singularity Problem

Some classic mechanical methods can be employed to solve crack problem [26–29]. For the grinding fixtures, many transnational corporations and research institutes have invested much capital in studying their design, materials, and manufacturing technique to improve their quality. There is no doubt that the experimental-based method is only considered to be an auxiliary method to study the mechanical properties of grinding fixtures due to their variety in kind and high experimental fees.Singularity problems, such as structural crack, can be easily found in grinding fixtures under multiple cyclic loads. The crack singularity leads to a deterioration of the performance in the procedure of finite element solution due to the singularity property of1/r, where r is the vector in column-coordinate of crack tip [30]. To acquire the numerical solutions with a high accuracy at crack tip, researchers need to employ the denser meshes or higher-order element near the crack region. To overcome this problem, wavelets have been applied to finite element analysis to identify crack location and depth of pipe [31], as well as rotor system [32]. Meanwhile their studies mentioned above are limited to detect surface cracks. In addition, the meshes will be redivided in studying the crack growth, which leads to the reduction of calculation accuracy and efficiency in some extent.Compared with other solution methods based on the partial differential equations, the BEM reduces the dimension of the crack problem. In addition, it is easier for the researchers to discretize the boundary rather than the solution domain. Relatively simple elements are needed to simulate the boundary shape, and the linear algebraic equations with lower order will be acquired. Therefore, the BEM is considered to be a higher effective method to solve the crack singularity problem compared to the FEM. Herein, the BEM is employed to effectively solve the SIF, and the BEM is extended to investigate the mechanical properties and predict their service life.


## 3. Modeling Procedures for Crack Problems Based on the BEM

### 3.1. Scheme of the BEM

The boundary integral equation in a 3D elastic domain without consideration of body forces can be written as follows [33]:(1)CijPUjP+∫SUjQTij∗P,QdSQ=∫STjQUij∗P,QdSQ,where Uj and Tj represent the displacement and traction components, respectively. Uij∗ and Tij∗ are their corresponding fundamental solutions. P and Q are any two points in the solution domain. S and Cij are the boundary surface and geometrical coefficients, respectively. It should be noted that the coefficients Cij are closely dependent on the location of the point P. When P is located in the solution domain, Cij=δij. When P is on a smooth boundary, Cij=1/2, and when P is outside the solution domain, Cij=0.In each element, the local coordinateXi can be represented by the node global coordinate Xil in the element by employing the interpolation function Nlξ,η, that is,(2)Xi=∑l=1MNlξ,ηXil,where M indicates the number of nodes in each element. ξ and η are the local coordinates as shown in Figure 1. Xil indicates the global coordinates of the node l relative to the element i.Figure 1
The relationship between global coordinate system and local coordinate system.In each element, the traction forcesTj and displacements Uj in local coordinate system can be expressed as(3)Uj=∑l=1MNlξ,ηUjl,Tj=∑l=1MNlξ,ηTjl,where Ujl and Tjl indicate the displacement and traction components along the coordinate axis j in global coordinate system.After discretizing into surface elements, the boundary integral equation can be expressed as follows:(4)CijPUjP=∑e=1Ne∑l=1M∫−11∫−11NKξ,ηUij∗P,QJdξdηTjl−∑e=1Ne∑l=1M∫−11∫−11NKξ,ηTij∗P,QJdξdηUjl,where J is Jacobian matrix. The boundary S is discretized into Ne elements with Nb boundary nodes. According to the Gauss numerical integration formula, it can be written as(5)∫−11∫−11NKξ,ηUij∗P,QJdξdη=G^ijelP,∫−11∫−11NKξ,ηTij∗P,QJdξdη=H^ijelP.Substituting equation (5) into equation (4), it can be written as(6)CijPUjP=∑e=1Ne∑l=1MG^ijelPTjl−∑e=1Ne∑l=1MH^ijelPUjl.By employing the transformation matrixZjel, node values of Ujl and Tjl in global coordinate system are written as follows:(7)Ujl=ZjelU,Tjl=ZjelT,where U and T are displacement and traction column vectors of the boundary nodes. The transformation matrix Zjel is the diagonal matrix. Substituting equation (7) into equation (6), the equations can be easily expressed in matrix format, that is,(8)CU=GT−H′U,where(9)GijP=∑e=1Ne∑l=1MG^ijelPZjel,Hij'P=∑e=1Ne∑l=1MH^ijelPZjel.


### 3.2. Stress Intensity Factor

In fracture mechanics, the open mode or type I crack (as shown in Figure2) is one of the most common cracks. Herein, the SIF of the crack is studied by BEM to lay a foundation for predicting the crack evolutions under cycling loadings.Figure 2
Open crack mode.The SIF, which relates to many factors, such as crack size and structural geometric characteristics, reflects the strength of elastic stress field at crack tip. In the modeling procedures by the BEM, longitudinal stresses (σx and σy) and shear stress τxy at the crack tip in global coordinate system as shown in Figure 3 can be expressed as follows:where the coordinate system r−θ indicates the polar coordinate of point A near the crack tip as shown in Figure 3(b).(10)σx=KI2πrcosθ22−sinθ2sin3θ2,σy=KI2πrcosθ21−sinθ2sin3θ2,τxy=KI2πrsinθ2cosθ2cos3θ2,Figure 3
Stress status near the crack tip. (a) Diagram of a three-dimensional crack. (b) Crack coordinate systems.
(a)(b)Displacement components at the crack tip can be written as(11)uv=21+υKIEr2πcosθ21−2υ+sin2θ2sinθ22−2υ+cos2θ2,where the parameters u and v denote the displacement in horizontal and vertical directions, respectively. E is the elastic modulus; υ is Poisson’s ratio.If the angleθ is equal to 180o, the displacement components at crack region can be simplified as(12)v=41+υ1−υEr2πKI.By employing the displacement components derived from the BEM, the SIF can be expressed as(13)KI=E1+v1−v2πrvBEM.


## 3.1. Scheme of the BEM

The boundary integral equation in a 3D elastic domain without consideration of body forces can be written as follows [33]:(1)CijPUjP+∫SUjQTij∗P,QdSQ=∫STjQUij∗P,QdSQ,where Uj and Tj represent the displacement and traction components, respectively. Uij∗ and Tij∗ are their corresponding fundamental solutions. P and Q are any two points in the solution domain. S and Cij are the boundary surface and geometrical coefficients, respectively. It should be noted that the coefficients Cij are closely dependent on the location of the point P. When P is located in the solution domain, Cij=δij. When P is on a smooth boundary, Cij=1/2, and when P is outside the solution domain, Cij=0.In each element, the local coordinateXi can be represented by the node global coordinate Xil in the element by employing the interpolation function Nlξ,η, that is,(2)Xi=∑l=1MNlξ,ηXil,where M indicates the number of nodes in each element. ξ and η are the local coordinates as shown in Figure 1. Xil indicates the global coordinates of the node l relative to the element i.Figure 1
The relationship between global coordinate system and local coordinate system.In each element, the traction forcesTj and displacements Uj in local coordinate system can be expressed as(3)Uj=∑l=1MNlξ,ηUjl,Tj=∑l=1MNlξ,ηTjl,where Ujl and Tjl indicate the displacement and traction components along the coordinate axis j in global coordinate system.After discretizing into surface elements, the boundary integral equation can be expressed as follows:(4)CijPUjP=∑e=1Ne∑l=1M∫−11∫−11NKξ,ηUij∗P,QJdξdηTjl−∑e=1Ne∑l=1M∫−11∫−11NKξ,ηTij∗P,QJdξdηUjl,where J is Jacobian matrix. The boundary S is discretized into Ne elements with Nb boundary nodes. According to the Gauss numerical integration formula, it can be written as(5)∫−11∫−11NKξ,ηUij∗P,QJdξdη=G^ijelP,∫−11∫−11NKξ,ηTij∗P,QJdξdη=H^ijelP.Substituting equation (5) into equation (4), it can be written as(6)CijPUjP=∑e=1Ne∑l=1MG^ijelPTjl−∑e=1Ne∑l=1MH^ijelPUjl.By employing the transformation matrixZjel, node values of Ujl and Tjl in global coordinate system are written as follows:(7)Ujl=ZjelU,Tjl=ZjelT,where U and T are displacement and traction column vectors of the boundary nodes. The transformation matrix Zjel is the diagonal matrix. Substituting equation (7) into equation (6), the equations can be easily expressed in matrix format, that is,(8)CU=GT−H′U,where(9)GijP=∑e=1Ne∑l=1MG^ijelPZjel,Hij'P=∑e=1Ne∑l=1MH^ijelPZjel.


## 3.2. Stress Intensity Factor

In fracture mechanics, the open mode or type I crack (as shown in Figure2) is one of the most common cracks. Herein, the SIF of the crack is studied by BEM to lay a foundation for predicting the crack evolutions under cycling loadings.Figure 2
Open crack mode.The SIF, which relates to many factors, such as crack size and structural geometric characteristics, reflects the strength of elastic stress field at crack tip. In the modeling procedures by the BEM, longitudinal stresses (σx and σy) and shear stress τxy at the crack tip in global coordinate system as shown in Figure 3 can be expressed as follows:where the coordinate system r−θ indicates the polar coordinate of point A near the crack tip as shown in Figure 3(b).(10)σx=KI2πrcosθ22−sinθ2sin3θ2,σy=KI2πrcosθ21−sinθ2sin3θ2,τxy=KI2πrsinθ2cosθ2cos3θ2,Figure 3
Stress status near the crack tip. (a) Diagram of a three-dimensional crack. (b) Crack coordinate systems.
(a)(b)Displacement components at the crack tip can be written as(11)uv=21+υKIEr2πcosθ21−2υ+sin2θ2sinθ22−2υ+cos2θ2,where the parameters u and v denote the displacement in horizontal and vertical directions, respectively. E is the elastic modulus; υ is Poisson’s ratio.If the angleθ is equal to 180o, the displacement components at crack region can be simplified as(12)v=41+υ1−υEr2πKI.By employing the displacement components derived from the BEM, the SIF can be expressed as(13)KI=E1+v1−v2πrvBEM.


## 4. Crack Propagation Analysis of the Typical Structures

To validate the proposed BEM method,KI of the plate structures made by zirconia materials are compared with the FEM and theoretical results. On this basis, the crack propagation is further discussed according to the maximum circumferential stress criterion [34]. The zirconia material parameters are as follows: elastic modulus is 220 GPa, and Poisson’s ratio is 0.25.


### 4.1. Rectangular Plate with Bilateral Cracks

A circular hole ofΦ2mm located at the geometric center of the plate, whose dimensions are 50 mm length, 10 mm width, and 1 mm thickness, is considered as shown in Figure 4. In order to easily and clearly observe crack evolutions, the bilateral precracks with length of 0.1 mm are implanted at the maximum stress regions. The constraints and external loadings are shown in Figure 4.Figure 4
Plate structure with bilateral cracks.For the open crack mode as shown in Figure2, the theoretical solutions of KIt can be written as [31](14)KIt=Fσπa,where the parameters σ and a indicate the stress at crack tip and the distance between crack tip and the center of circular hole, respectively. F indicates an effective coefficient and its value 0.84 is employed according to [35]. It can be easily found that the SIF is closely dependent on crack length. The value of KIt is equal to 15.6153MPa⋅m1/2 when the initial crack length is equal to 0.1 mm.To discuss the relations between SIF and crack length, the boundary element model of the plate as shown in Figure4 is established. When the multiregion algorithm is used, the model is divided into two regions by the crack surface along crack direction. In the modeling procedures, the quadrilateral elements are employed to discrete the plate structure. In the stress concentration region near the hole, more density triangle elements are required to solve the stresses with a high accuracy. For a comparison, numerical results acquired by the FEM are also investigated under an identical condition. Numerical results KI of the plate with four different crack lengths (0.1 mm, 0.2 mm, 0.3 mm, and 0.4 mm) are shown in Figure 5(a). Here, the horizontal axis and vertical axis in this figure indicate the number of discrete points and theoretical results of the KI acquired by the BEM, respectively. It can be easily found that the crack growth will sharply increase the SIF at crack tip. Meanwhile it is interesting to mention that KI is slightly dependent on quantities of discrete points. In other words, it takes relatively less discrete points to solve the SIF with a high accuracy by using the BEM. To further validate the method, numerical results by the FEM and theoretical solutions by experience equations are both considered. The error analysis can be seen in Figure 5(b). It can be concluded that the BEM agrees well with the FEM when the crack growth is from 0.2 mm to 0.4 mm. If a small crack of 0.1 mm is considered, the BEM provides a higher accuracy than the FEM, and the relative error is equal to 8.22%. A higher calculation accuracy of the SIF will be acquired by the BEM when the crack growth from 0.3 mm to 0.4 mm is considered in the case.Figure 5
(a)KI value along the crack front. (b) Relative error comparisons between the FEM and the BEM.
(a)(b)Figure6 indicates stress distributions and crack evolutions of the plate by employing the BEM. It can be easily found that the maximum stress is located at the crack tip, and crack growth is perpendicular to the loading direction. In addition, the stress distribution variations with respect to crack evolutions can be easily discerned. The maximum stresses at the crack tip are sharply increased with the crack evolution. In detail, the maximum stress is approximated to 35.05 MPa once the crack length of 0.5 mm is considered as shown in Figure 6(a). Once when the crack extends to 4.0 mm as shown in Figure 6(h), the maximum stress at crack region increases rapidly to 541.0 MPa. In addition, the identical crack propagations will be found in both sides of the circle hole due to the structural symmetry of the studied plate.Figure 6
Mises stress distributions of the plate with different crack lengths: (a) 0.5 mm, (b) 1.0 mm, (c) 1.5 mm, (d) 2.0 mm, (e) 2.5 mm, (f) 3.0 mm, (g) 3.5 mm, and (h) 4.0 mm.
(a)(b)(c)(d)(e)(f)(g)(h)


### 4.2. Rectangular Plate with a Unilateral Crack

A rectangular plate made of zirconia materials with a unilateral crack is also considered here to further validate the established model. The dimension parameters are as follows: length, 10 mm; width, 3 mm; and height, 30 mm. A rectangular precrack with a length of 1 mm is located at the stress concentration region to investigate the crack growth. The constraints and external loadings are applied on the structure as shown in Figure7.Figure 7
Plate structure with a unilateral crack.For the fracture mode of this plate with a unilateral crack as shown in Figure7, theoretical solutions of the SIF KIα can be acquired according to the manual of the SIF [31]; that is,(15)KIα=Fσπa,where the parameters F and a indicate the related coefficient and crack length, respectively. When the initial crack length a=1mm is considered, the parameter F is equal to 1.2. Taking these parameters into equation (15), the theoretical result of KIα is 20.3042MPa⋅mm1/2.For investigating theKI variations with respect to crack growth, the BEM model is established. Quadrilateral elements are used to discretize the structure. With a full consideration of stress concentration in the crack region and improving the calculation efficiency, one-quarter discretization ratio of the elements away from the crack region is employed in establishing the model. Figure 8(a) presents the numerical results of KI at crack tip, and five different lengths (1 mm, 2 mm, 3 mm, 4 mm, and 5 mm) are considered in the example. Similar to the results of a plate with bilateral cracks, it can be easily found that the unilateral crack growth will also sharply increase KI. In addition, it is difficult for researchers to discern the effect of the number of discrete points on calculation accuracy when KI is considered. To further validate the established model, results by the FEM and theoretical method are also considered, and relative errors are shown in Figure 8(b). In general, numerical results by the BEM agree well with the theoretical results. In addition, it can be concluded that the BEM provides a higher accuracy than the FEM. In detail, the errors acquired by the BEM and FEM are 1.75% and 7.65% when crack length of 1 mm is considered. The maximum relative error acquired by the BEM is 0.64% when crack extends to 3 mm.Figure 8
KI values at the crack tip. (a) Relations between discrete points and KI. (b) Relative error comparisons between the FEM and the BEM.
(a)(b)In order to execute the crack propagation analysis, a precrack with length of 1 mm is also implanted into the boundary element model as shown in Figure9(a). Figures 9(b)–9(i) indicate crack evolutions of the plate with the increasing of external loading. It can be seen that the maximum stress is located at the crack tip, and the crack growth is perpendicular to the loading direction. With the crack growth, the maximum stress will be further increased. In detail, the maximum stress is 18.01 MPa when an initial crack of 1 mm is considered. With an increase of crack length, the maximum stress of the plate at crack tip is 256.1MPa when the crack length of 7 mm is considered in the example.Figure 9
Mises stress distributions of plate with different crack lengths: (a) 1 mm, (b) 2 mm, (c) 3 mm, (d) 4 mm, (e) 5 mm, (f) 6 mm, (g) 7 mm, (h) 8 mm, and (i) 9 mm.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 4.1. Rectangular Plate with Bilateral Cracks

A circular hole ofΦ2mm located at the geometric center of the plate, whose dimensions are 50 mm length, 10 mm width, and 1 mm thickness, is considered as shown in Figure 4. In order to easily and clearly observe crack evolutions, the bilateral precracks with length of 0.1 mm are implanted at the maximum stress regions. The constraints and external loadings are shown in Figure 4.Figure 4
Plate structure with bilateral cracks.For the open crack mode as shown in Figure2, the theoretical solutions of KIt can be written as [31](14)KIt=Fσπa,where the parameters σ and a indicate the stress at crack tip and the distance between crack tip and the center of circular hole, respectively. F indicates an effective coefficient and its value 0.84 is employed according to [35]. It can be easily found that the SIF is closely dependent on crack length. The value of KIt is equal to 15.6153MPa⋅m1/2 when the initial crack length is equal to 0.1 mm.To discuss the relations between SIF and crack length, the boundary element model of the plate as shown in Figure4 is established. When the multiregion algorithm is used, the model is divided into two regions by the crack surface along crack direction. In the modeling procedures, the quadrilateral elements are employed to discrete the plate structure. In the stress concentration region near the hole, more density triangle elements are required to solve the stresses with a high accuracy. For a comparison, numerical results acquired by the FEM are also investigated under an identical condition. Numerical results KI of the plate with four different crack lengths (0.1 mm, 0.2 mm, 0.3 mm, and 0.4 mm) are shown in Figure 5(a). Here, the horizontal axis and vertical axis in this figure indicate the number of discrete points and theoretical results of the KI acquired by the BEM, respectively. It can be easily found that the crack growth will sharply increase the SIF at crack tip. Meanwhile it is interesting to mention that KI is slightly dependent on quantities of discrete points. In other words, it takes relatively less discrete points to solve the SIF with a high accuracy by using the BEM. To further validate the method, numerical results by the FEM and theoretical solutions by experience equations are both considered. The error analysis can be seen in Figure 5(b). It can be concluded that the BEM agrees well with the FEM when the crack growth is from 0.2 mm to 0.4 mm. If a small crack of 0.1 mm is considered, the BEM provides a higher accuracy than the FEM, and the relative error is equal to 8.22%. A higher calculation accuracy of the SIF will be acquired by the BEM when the crack growth from 0.3 mm to 0.4 mm is considered in the case.Figure 5
(a)KI value along the crack front. (b) Relative error comparisons between the FEM and the BEM.
(a)(b)Figure6 indicates stress distributions and crack evolutions of the plate by employing the BEM. It can be easily found that the maximum stress is located at the crack tip, and crack growth is perpendicular to the loading direction. In addition, the stress distribution variations with respect to crack evolutions can be easily discerned. The maximum stresses at the crack tip are sharply increased with the crack evolution. In detail, the maximum stress is approximated to 35.05 MPa once the crack length of 0.5 mm is considered as shown in Figure 6(a). Once when the crack extends to 4.0 mm as shown in Figure 6(h), the maximum stress at crack region increases rapidly to 541.0 MPa. In addition, the identical crack propagations will be found in both sides of the circle hole due to the structural symmetry of the studied plate.Figure 6
Mises stress distributions of the plate with different crack lengths: (a) 0.5 mm, (b) 1.0 mm, (c) 1.5 mm, (d) 2.0 mm, (e) 2.5 mm, (f) 3.0 mm, (g) 3.5 mm, and (h) 4.0 mm.
(a)(b)(c)(d)(e)(f)(g)(h)


## 4.2. Rectangular Plate with a Unilateral Crack

A rectangular plate made of zirconia materials with a unilateral crack is also considered here to further validate the established model. The dimension parameters are as follows: length, 10 mm; width, 3 mm; and height, 30 mm. A rectangular precrack with a length of 1 mm is located at the stress concentration region to investigate the crack growth. The constraints and external loadings are applied on the structure as shown in Figure7.Figure 7
Plate structure with a unilateral crack.For the fracture mode of this plate with a unilateral crack as shown in Figure7, theoretical solutions of the SIF KIα can be acquired according to the manual of the SIF [31]; that is,(15)KIα=Fσπa,where the parameters F and a indicate the related coefficient and crack length, respectively. When the initial crack length a=1mm is considered, the parameter F is equal to 1.2. Taking these parameters into equation (15), the theoretical result of KIα is 20.3042MPa⋅mm1/2.For investigating theKI variations with respect to crack growth, the BEM model is established. Quadrilateral elements are used to discretize the structure. With a full consideration of stress concentration in the crack region and improving the calculation efficiency, one-quarter discretization ratio of the elements away from the crack region is employed in establishing the model. Figure 8(a) presents the numerical results of KI at crack tip, and five different lengths (1 mm, 2 mm, 3 mm, 4 mm, and 5 mm) are considered in the example. Similar to the results of a plate with bilateral cracks, it can be easily found that the unilateral crack growth will also sharply increase KI. In addition, it is difficult for researchers to discern the effect of the number of discrete points on calculation accuracy when KI is considered. To further validate the established model, results by the FEM and theoretical method are also considered, and relative errors are shown in Figure 8(b). In general, numerical results by the BEM agree well with the theoretical results. In addition, it can be concluded that the BEM provides a higher accuracy than the FEM. In detail, the errors acquired by the BEM and FEM are 1.75% and 7.65% when crack length of 1 mm is considered. The maximum relative error acquired by the BEM is 0.64% when crack extends to 3 mm.Figure 8
KI values at the crack tip. (a) Relations between discrete points and KI. (b) Relative error comparisons between the FEM and the BEM.
(a)(b)In order to execute the crack propagation analysis, a precrack with length of 1 mm is also implanted into the boundary element model as shown in Figure9(a). Figures 9(b)–9(i) indicate crack evolutions of the plate with the increasing of external loading. It can be seen that the maximum stress is located at the crack tip, and the crack growth is perpendicular to the loading direction. With the crack growth, the maximum stress will be further increased. In detail, the maximum stress is 18.01 MPa when an initial crack of 1 mm is considered. With an increase of crack length, the maximum stress of the plate at crack tip is 256.1MPa when the crack length of 7 mm is considered in the example.Figure 9
Mises stress distributions of plate with different crack lengths: (a) 1 mm, (b) 2 mm, (c) 3 mm, (d) 4 mm, (e) 5 mm, (f) 6 mm, (g) 7 mm, (h) 8 mm, and (i) 9 mm.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 5. Fatigue Crack Propagation Analysis of the Zirconia Grinding Fixture

Figure10 is a typical 3D model of grinding fixtures made by zirconia materials, which is widely used in semiconductor processing. For this fixture, since it is used to replace the old fixture of steel material, it had been special toughening treatment to improve property. The components C and C' indicate the fixed pin holes, which are used to restrain the grinding fixture. The symbols A1–G1 indicate the force adjusting heads of the fixture. In the manufacturing process, the strip-shaped workpiece is constrained on the surface of grinding fixture. The structural dimensions of the fixture are as follows: length of 220 mm, width of 50 mm, and thickness of 5 mm.Figure 10
The model of grinding fixture.Due to the high hardness and brittleness of zirconia materials, machined surfaces of the fixture are prone to produce microcracks, and the cracks will immediately accumulate. In service, periodical cycling loads will be applied on the holes of the adjusting heads to obtain a desired bottom surface profile. It is critical for researchers to grasp crack status and fatigue life to maximize their service life. To this end, the 3D model of the fixture based on the BEM is established. In the fatigue life analysis, the surfaces of the fixed pin holes C and C′ are fully restrained. A cyclic loading is supposed to apply on the hole of force adjusting heads A1–G1. An external load of 0.1 N is applied on the surfaceK to grasp stress distributions of the fixture.The flow chart of fatigue crack predictions is shown in Figure11; it can be summarized as follows:(1)
The Object Solid Modeler (OSM) is employed to establish the 3D BEM of the grinding fixture.(2)
Meshing the surface and inputting the material properties, boundary conditions. The discretized model based on the BEM is established by employing the Fracture Analysis Code (FRANC3D).(3)
Before the fatigue analysis, a small external load is applied on the fixture, and the stress and deformation cloud diagrams will be acquired. According to the stress distributions, an initial crack is embedded into the maximum stress area of the fixture.(4)
Applying cyclic loading on the fixture, the crack will accumulate and extend. At each crack increment, the SIF will be recalculated.(5)
According to the history data of the SIF and the corresponding crack length, the fatigue life of fixture will be determined.The process of crack propagation can be summarized as follows:(1)
Based on the boundary element method (BEM), the displacement field of crack surface is solved, and the distribution of SIF along the current crack front is solved by crack opening displacement (COD) fitting method(2)
Based on classical Paris Erdogan Law, the propagation size and position of each point on the crack front edge are calculated(3)
The new crack front is obtained by fitting the new position after propagation by least square method(4)
Repeat the above steps for the specified simulation accuracyFigure 11
The flow diagram of FRANC3D.


### 5.1. Stress Distributions and Crack Propagations

In the mechanical analysis of the fixture, it will be found that the maximum stress is located at the root of the pillarP as shown in Figure 12.Figure 12
Mises stress distributions of the fixture. (a) Fixture. (b) Partial enlarged figure.When calculating the stress intensity factors of three-dimensional cracks by FRANC3D, the cracks are usually simplified to round, ellipse, and rectangle according to the actual situation. According to the experimental observation, the fatigue crack in the pillar of the fixture is basically rectangular, which is continuously expanding forward. Herein, we simplify the crack model into a rectangle in order to further investigate the fatigue crack propagation, and a unilateral crack of 0.1 mm is implanted into the fixture as shown in Figure13(a). Figures 13(b)–13(f) indicate the relations between crack lengths (0.5 mm, 1.0 mm, 1.5 mm, 2.0 mm, and 2.5 mm) and stress distributions. Similar to the conclusions acquired by the plate structure, the maximum stress in each case is located at the crack tip, and the crack growth is perpendicular to the loading direction. In detail, when the crack length is 0.1 mm as shown in Figure 13(a), the maximum stress and KI at crack tip are 162.0 MPa and 61.26MPa⋅mm1/2, respectively. With the crack growth, the maximum stress at crack region will be sharply increased as shown in Figures 13(b)–13(f). For instance, the maximum stress and KI are 776.2 MPa and KI=394.4MPa⋅mm1/2 when the crack extends to 1.5 mm, respectively.Figure 13
Mises stress distributions of the fixture with a unilateral crack. (a) 0.1 mm. (b) 0.5 mm. (c) 1.0 mm. (d) 1.5 mm. (e) 2.0 mm. (f) 2.5 mm.
(a)(b)(c)(d)(e)(f)


### 5.2. Fatigue Crack Propagation Analysis

In order to determine the fatigue life of the fixture under cyclic loads, herein the classical Paris Erdogan Law is employed. The relationship between crack growth rate and SIF range can be written as [36, 37](16)dadN=CΔKeffm,where N and a indicate the number of applied loadings and crack length, respectively. da/dN describes the crack growth per cycle; it has a power law relationship with the stress intensity factor amplitude. ΔKeff is a parameter related to the SIF. C and m are material constants. Take C=2.5×10−14 and m=4.5.According to the crack propagation law studies mentioned above, an initial crack is supposed and located at the bottom of the pillar. The crack evolution trajectory in the fixture can be confirmed as shown in Figure14. The SIF at the crack tip along the crack trajectory E (Figure 14(b)) can be first calculated. When dealing with the three-dimensional crack propagation problem, the maximum circumferential stress criterion is used to judge the propagation direction of the crack front edge under the specified length.Figure 14
Crack in the fixture. (a) Cross section of the crack. (b) Partial enlarged figure of the crack trajectory.Stress intensity factor (SIF) obtained from each crack growth step (crack length) can be determined by maximum circumferential stress. As shown in Figure15, it can be easily seen that the relations between KI at crack tip and crack length show an evidently nonlinear property. When the cracks extend to 1.2 mm and 1.3 mm, maximum stress intensity factors KI of the fixture at crack tip are 295.899MPa⋅mm1/2 and 382.171MPa⋅mm1/2, respectively. According to the fracture toughness of 300MPa⋅mm1/2 of the fixture acquired by Wum et al. [38], the crack length of 1.2 mm is considered to be finial failure length of this fixture.Figure 15
KI values of the fixture at crack tip with respect to crack length.It can be seen from this that the correspondingda and ΔKeff in classical Paris Erdogan Law can obtain the corresponding dN after determining the coefficient C and m related to the material. Figure 16 indicates the load cycle number versus crack length of the fixture obtained by the boundary element method. It can be found that the crack growth is at a relatively slow rate if the crack length is less than 1.2 mm. Accordingly, the fatigue life of zirconia grinding fixture is approximated to 1.6×104 when crack limit length is confirmed to be 1.2 mm. According to the usage experience, the cyclic number of the fixtures is about 100 times per day. To maximize production safety, the number of cycles is supposed to be 1.5×104. In other words, the fatigue life can be considered as 15000/100 = 150 days. Therefore, it is recommended that a monthly inspection of the appearance of the fixture should be executed during the production manufacturing process. In addition, once the crack length in the pillar reaches 1.2 mm, the fixture should be replaced to avoid a sudden fracture of the fixture.Figure 16
Relationships between cycle number and crack length.


## 5.1. Stress Distributions and Crack Propagations

In the mechanical analysis of the fixture, it will be found that the maximum stress is located at the root of the pillarP as shown in Figure 12.Figure 12
Mises stress distributions of the fixture. (a) Fixture. (b) Partial enlarged figure.When calculating the stress intensity factors of three-dimensional cracks by FRANC3D, the cracks are usually simplified to round, ellipse, and rectangle according to the actual situation. According to the experimental observation, the fatigue crack in the pillar of the fixture is basically rectangular, which is continuously expanding forward. Herein, we simplify the crack model into a rectangle in order to further investigate the fatigue crack propagation, and a unilateral crack of 0.1 mm is implanted into the fixture as shown in Figure13(a). Figures 13(b)–13(f) indicate the relations between crack lengths (0.5 mm, 1.0 mm, 1.5 mm, 2.0 mm, and 2.5 mm) and stress distributions. Similar to the conclusions acquired by the plate structure, the maximum stress in each case is located at the crack tip, and the crack growth is perpendicular to the loading direction. In detail, when the crack length is 0.1 mm as shown in Figure 13(a), the maximum stress and KI at crack tip are 162.0 MPa and 61.26MPa⋅mm1/2, respectively. With the crack growth, the maximum stress at crack region will be sharply increased as shown in Figures 13(b)–13(f). For instance, the maximum stress and KI are 776.2 MPa and KI=394.4MPa⋅mm1/2 when the crack extends to 1.5 mm, respectively.Figure 13
Mises stress distributions of the fixture with a unilateral crack. (a) 0.1 mm. (b) 0.5 mm. (c) 1.0 mm. (d) 1.5 mm. (e) 2.0 mm. (f) 2.5 mm.
(a)(b)(c)(d)(e)(f)


## 5.2. Fatigue Crack Propagation Analysis

In order to determine the fatigue life of the fixture under cyclic loads, herein the classical Paris Erdogan Law is employed. The relationship between crack growth rate and SIF range can be written as [36, 37](16)dadN=CΔKeffm,where N and a indicate the number of applied loadings and crack length, respectively. da/dN describes the crack growth per cycle; it has a power law relationship with the stress intensity factor amplitude. ΔKeff is a parameter related to the SIF. C and m are material constants. Take C=2.5×10−14 and m=4.5.According to the crack propagation law studies mentioned above, an initial crack is supposed and located at the bottom of the pillar. The crack evolution trajectory in the fixture can be confirmed as shown in Figure14. The SIF at the crack tip along the crack trajectory E (Figure 14(b)) can be first calculated. When dealing with the three-dimensional crack propagation problem, the maximum circumferential stress criterion is used to judge the propagation direction of the crack front edge under the specified length.Figure 14
Crack in the fixture. (a) Cross section of the crack. (b) Partial enlarged figure of the crack trajectory.Stress intensity factor (SIF) obtained from each crack growth step (crack length) can be determined by maximum circumferential stress. As shown in Figure15, it can be easily seen that the relations between KI at crack tip and crack length show an evidently nonlinear property. When the cracks extend to 1.2 mm and 1.3 mm, maximum stress intensity factors KI of the fixture at crack tip are 295.899MPa⋅mm1/2 and 382.171MPa⋅mm1/2, respectively. According to the fracture toughness of 300MPa⋅mm1/2 of the fixture acquired by Wum et al. [38], the crack length of 1.2 mm is considered to be finial failure length of this fixture.Figure 15
KI values of the fixture at crack tip with respect to crack length.It can be seen from this that the correspondingda and ΔKeff in classical Paris Erdogan Law can obtain the corresponding dN after determining the coefficient C and m related to the material. Figure 16 indicates the load cycle number versus crack length of the fixture obtained by the boundary element method. It can be found that the crack growth is at a relatively slow rate if the crack length is less than 1.2 mm. Accordingly, the fatigue life of zirconia grinding fixture is approximated to 1.6×104 when crack limit length is confirmed to be 1.2 mm. According to the usage experience, the cyclic number of the fixtures is about 100 times per day. To maximize production safety, the number of cycles is supposed to be 1.5×104. In other words, the fatigue life can be considered as 15000/100 = 150 days. Therefore, it is recommended that a monthly inspection of the appearance of the fixture should be executed during the production manufacturing process. In addition, once the crack length in the pillar reaches 1.2 mm, the fixture should be replaced to avoid a sudden fracture of the fixture.Figure 16
Relationships between cycle number and crack length.


## 6. Conclusions

In this study, the BEM was developed to predict fatigue life of the zirconia fixture. To validate the proposed method, the stress intensity factors of plate structures with a unilateral precrack and bilateral precracks are investigated. Compared with the theoretical solutions and the FEM, numerical results by the BEM provide a high accuracy in predicting the SIF. The results show that the crack propagation will sharply increase the SIF at crack tip. On this basis, the crack propagation and service life of the zirconia grinding fixture under cyclic load are further studied. The results indicate that the maximum stress is located at the root of the pillar, and the crack length of 1.2 mm is suggested to be the finial failure length of the fixture to avoid a sudden fracture for the fixture.



---

*Source: 8863379-2020-08-12.xml*