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

Zirconia grinding ﬁxtures 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 ﬁxtures. 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 ﬁxture under cyclic loads. Numerical results show that stress concentration occurs at the pillar of the ﬁxture, and crack growth is perpendicular to the loading direction.


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][7][8][9]. e 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][12][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. e results show that the effect of notch angles on crack propagation increased with notch depth. El-sayed et al. [16] introduced a threedimensional 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. is 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. e 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 largescale 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. e outline of this paper is as follows: Section 2 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. e conclusions are summarized in Section 6.

Crack Singularity Problem
Some classic mechanical methods can be employed to solve crack problem [26][27][28][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.
ere 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. e crack singularity leads to a deterioration of the performance in the procedure of finite element solution due to the singularity property of 1/ � r √ , where r is the vector in columncoordinate 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. erefore, 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.

Modeling Procedures for Crack Problems
Based on the BEM 3.1. Scheme of the BEM. e boundary integral equation in a 3D elastic domain without consideration of body forces can be written as follows [33]: where U j and T j represent the displacement and traction components, respectively. U * ij and T * ij are their corresponding fundamental solutions. P and Q are any two points in the solution domain. S and C ij are the boundary surface and geometrical coefficients, respectively. It should be noted that the coefficients C ij are closely dependent on the location of the point P. When P is located in the solution domain, C ij � δ ij . When P is on a smooth boundary, C ij � 1/2, and when P is outside the solution domain, C ij � 0.
In each element, the local coordinate X i can be represented by the node global coordinate X l i in the element by employing the interpolation function N l (ξ, η), that is, where M indicates the number of nodes in each element. ξ and η are the local coordinates as shown in Figure 1. X l i indicates the global coordinates of the node l relative to the element i.
In each element, the traction forces T j and displacements U j in local coordinate system can be expressed as where U (l) j and T (l) j 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: where J is Jacobian matrix. e boundary S is discretized into N e elements with N b boundary nodes. According to the Gauss numerical integration formula, it can be written as Substituting equation (5) into equation (4), it can be written as By employing the transformation matrix Z el j , node values of U (l) j and T (l) j in global coordinate system are written as follows: where U and T are displacement and traction column vectors of the boundary nodes. e transformation matrix Z el j is the diagonal matrix. Substituting equation (7) into equation (6), the equations can be easily expressed in matrix format, that is, where 3.2. Stress Intensity Factor. In fracture mechanics, the open mode or type I crack (as shown in Figure 2) 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. e 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 Displacement components at the crack tip can be written as 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 180 o , the displacement components at crack region can be simplified as By employing the displacement components derived from the BEM, the SIF can be expressed as

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]. e zirconia material parameters are as follows: elastic modulus is 220 GPa, and Poisson's ratio is 0.25.

Rectangular Plate with Bilateral Cracks.
A circular hole of Φ2 mm 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. e constraints and external loadings are shown in Figure 4.
For the open crack mode as shown in Figure 2, the theoretical solutions of KI t can be written as [31] 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. e value of KI t is equal to 15.6153 MPa · m 1/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 Figure 4 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. e 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 6 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.
e 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. Shock and Vibration

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. e 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. e constraints and external loadings are applied on the structure as shown in Figure 7.
For the fracture mode of this plate with a unilateral crack as shown in Figure 7, theoretical solutions of the SIF KI (α) can be acquired according to the manual of the SIF [31]; that is, 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.3042 MPa · mm 1/2 .
For investigating the KI 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   Shock and Vibration errors acquired by the BEM and FEM are 1.75% and 7.65% when crack length of 1 mm is considered. e maximum relative error acquired by the BEM is 0.64% when crack extends to 3 mm.
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 Figure 9(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.1 MPa when the crack length of 7 mm is considered in the example. Figure 10 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. e components C and C' indicate the fixed pin holes, which are used to restrain the grinding fixture. e 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. e structural dimensions of the fixture are as follows: length of 220 mm, width of 50 mm, and thickness of 5 mm.

Fatigue Crack Propagation Analysis of the Zirconia 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 surface K to grasp stress distributions of the fixture. e flow chart of fatigue crack predictions is shown in Figure 11; it can be summarized as follows: 35   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. e 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

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 pillar P as shown in Figure 12.
When calculating the stress intensity factors of threedimensional 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 Figure 13(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.26 MPa · mm 1/2 , respectively. With the crack growth, the maximum stress at crack region will be sharply increased as shown in Figures 13(b)    relationship between crack growth rate and SIF range can be written as [36,37] da 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. ΔK eff 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. e crack evolution trajectory in the fixture can be confirmed as shown in Figure 14. e 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.
Stress intensity factor (SIF) obtained from each crack growth step (crack length) can be determined by maximum circumferential stress. As shown in Figure 15, 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.899 MPa · mm 1/2 and 382.171 MPa · mm 1/2 , respectively. According to the fracture toughness of 300 MPa · mm 1/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.
It can be seen from this that the corresponding da and ΔK eff 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 × 10 4 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 × 10 4 . In other words, the fatigue life can be considered as 15000/100 � 150 days. erefore, 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.

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. e 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. e 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.
Data Availability e raw/processed data required to reproduce these findings cannot be shared at this time due to technical and time limitations.

Conflicts of Interest
e authors declare no conflicts of interest. e first author is grateful to the Engineering Department, Lancaster University, for the support he has received during the course of his visit.