Numerical Investigation on Crack Propagation for a Central Cracked Brazilian Disk Concerning Friction

: This paper concerns the effect of friction on crack propagation for the centrally cracked Brazilian disk under diametric forces by using a modiﬁed ﬁnite element method. It shows that the mode II stress intensity factor decreases obviously with the increase of friction after the crack is closed, while friction has no inﬂuence on the stress intensity factor of mode I and T-stress. Meanwhile, there are some signiﬁcant inﬂuences on the crack propagation due to the change of the friction after the crack is closed with the appropriate loading angle and relative length of the crack. When T-stress is positive, the effect of friction becomes obvious and the crack propagation angle increases with a lager friction coefﬁcient. With increasing the friction, the deviation for the crack propagation trajectory increases and the curvature of path decreases, which may lead to the change of crack type. Additionally, the larger relative crack length can amplify the effect of friction, which is similar to the loading angle.


Introduction
Crack defects are often distributed in brittle media, including concretes and rocks. Under the action of complex external loads, these defects cause local cracking and form macroscopic cracks, resulting in the decrease of strength and stiffness of the medium, which may reproduce and extend to complete penetration. In order to keep the integrity and safety of structures, the study on crack initiation and growth is of guiding significance. Meanwhile, the stress intensity factors (SIFs) and T-stresses, as important fracture parameters of the linear elastic fracture mechanics (LEFM) used to evaluate the fracture performance of brittle materials, also have important practical significance in engineering fracture problems. Additionally, the central crack Brazilian disc has been widely used due to its convenience in changing the fracture mode of crack [1][2][3].
In any case, rough brittle materials with many closed cracks are subjected to the effect of the friction in actual underground engineering such as tunnel operation and gas exploration. Under these conditions, the fracture parameters, crack initiation and propagation of brittle materials under compressive-shear loading would be greatly affected. Therefore, in order to better understand the failure mechanisms and propagation behavior concerning the effect of friction of brittle materials such as rocks, some experimental tests were performed in this area. Erarslan et al. [4] carried out composite fracture loading by using the cracked chevron notched Brazilian disk (CCNBD) specimens made of three kinds of rock materials.
In the test, a displacement gauge was used to measure the displacement of radial compression in the crack tip. Moreover, this experimental result showed that in the case of compression-shear fracture, both the initiation angle and the initiation point of the crack would be affected. Nemat-Nasser and Horii [5] described compression-shear experiments of uniaxial compression concerning frictional contact. In order to control friction, the brass sheets or thin Teflon sheets were inserted into the crack. Their results showed that the specimen failure was caused by the axial splitting rather than shear failure and the crack propagated in the direction of the applied load, which was also demonstrated by the maximum tangential stress (MTS) criterion. Meanwhile, Cotterell [6] and Hoek et al. [7] found a similar phenomenon through a range of experiments. Furthermore, compressionshear experiments were reported as well on the law of crack propagation [8][9][10]. Although the experimental results could provide guidance in other directions, their implementation was usually limited by external conditions such as time-consuming and complex processes. Additionally, the cracks were not easy to contact due to the crack width [4] and some physical parameters such as friction coefficient could not be controlled easily. Therefore, theoretical method has attracted growing interest because of its convenience and reliability.
In the past few years, on the basis of previous reports by Dong [1] and Hua [11], Li [12] and Tang [13] analyzed the stress field inside the centrally cracked Brazilian disk (CCBD) specimen under different loading conditions when the crack surface contacted through the theoretical study. Furthermore, an analytical formula was obtained for calculating the Mode II SIFs in CCBD specimens under different loading conditions by using the weight function method (WFM). Meanwhile, prompted by some theoretical researches on T-stress [14,15], Tang [16] proved that the friction had no effect on T-stress in CCBD specimens subjected to diametric forces and confining pressure. Spagnoli et al. [17] discussed the effect of roughness between crack surfaces and friction coefficient on the dimensionless Mode II stress intensity factor and crack initiation angle using Distributed Dislocation Technique (DDT) and MTS criterion.
Currently, due to the complex boundary conditions, numerical approaches are increasingly used. For instance, Hua et al. [11,16] and Zhou et al. [18] accurately calculated the stress intensity factors and T-stresses for CCBD under concentrated force and confining pressure through conventional finite element method (FEM). Huang et al. [19] and Hou et al. [20] calculated those fracture parameters under different confining pressure coefficients based on the improved finite element method. Furthermore, various methods were reported for the fracture propagation under complex conditions. Lim et al. [21] developed a technique of automatic local meshing based on finite element method to simulate mixedmode fracture growth automatically. Erarslan et al. [4] and Haeri et al. [22] combined finite element method and MTS criterion to simulate the crack propagation from pure I to compression-shear fracture of Brazilian disk specimens. These results had a high matching degree with the experimental outcome. Meanwhile, numerical methods such as generalized particle dynamics, unconventional state-based peridynamics and phase field method, combining MTS and Mohr-Coulomb criterion, would simulate crack growth of different modes fracture [23][24][25]. Belytschko et al. [26] also proposed a convenient numerical method introducing an enrichment function to deal with complex practical problems, which was called as extended finite element method (XFEM). Based on XFEM, Dolbow et al. [27] combined traditional MTS criterion and LATIN method to explore the propagation of crack concerning frictional contact. The results showed that the simulated extension path was consistent with the experimental observations of Nemat-Nasser and Horii [5].
However, in the simulation of two-dimensional crack propagation problem, XFEM usually adopts first-order elements with weak analytical performance [28], which might result in inaccurate outcome when simulating compression-shear fracture. Moreover, there are few reports on fracture parameters and propagation of crack, considering the effect of friction for CCBD specimens through the general numerical methods. Therefore, prompted by an improved numerical method proposed by Huang et al. [19], a modified finite element method (MFEM) is developed. The MFEM involving the behavior of the frictional contact is based on FEM using user defined subroutine. Additionally, the purpose of this study is to investigate the influence of friction on fracture parameters and propagation of crack for CCBD specimens under radial compression forces when the crack surface is in contact, which is usually difficult to be studied through experimental tests.
The framework of this investigation is briefly summarized as follows. Section 2 mainly elaborates the principle of the interaction integral for MFEM and the generalized MTS (GMTS) criterion. Section 3 mainly illustrates when the crack surface is in contact, the analytical formula of SIFs and T-stresses concerning the friction by using the WFM for CCBD. Section 4 introduces some calculation examples for CCBD specimens about the effects of friction on fracture parameters (mainly SIFs and T-stresses) and crack propagation. Section 5 summarizes the conclusion of the whole work.

Two-Dimensional Interaction Integral with MFEM
The energy integral I, called the interaction integral, includes the contour Γ around the crack tip, which is a highly accurate method that has emerged in recent years to solve the two-dimensional fracture problem [29,30]. It is based on the quasi-static derivation process, assuming that the material has the characteristics of small deformation, linear elasticity and continuity.
where n j is the unit external normal vector of the contour Γ. In addition, the form of the strain energy density W for the interaction integral is as follows: In finite element theory, it is not easy to solve the infinitesimal contour integral near the crack tip. In order to facilitate the calculation in practical applications, the interaction integral formula is usually converted into an equivalent surface integral as shown in Figure 1. There is a new contour C set outside the contour Γ, which makes the region S closed. Then, through the divergence theorem, the line integral of Equation (1) can be further written as below [19,27,31] where q is defined as the variable weight function, which takes the value of 1 in region S and 0 in other places. In addition, there is a unit normal vector m j , which is perpendicular to profile C, as shown in Figure 1. Therefore, the stress intensity factors and T-stress in the real fields can be obtained by selecting the appropriate auxiliary fields. Because the selection of the stress intensity factors and T-stresses in the auxiliary fields is different for the interaction integral, the selection of auxiliary fields in this paper is consistent with that of the Refs. [19,[32][33][34].
When the integral path Γ approaches the tip of crack, the relationship between the stress intensity factors and the interaction integral I can be obtained as below factors and T-stresses in the auxiliary fields is different for the interaction integral, the selection of auxiliary fields in this paper is consistent with that of the Refs. [19,[32][33][34].
When the integral path Γ approaches the tip of crack, the relationship between the stress intensity factors and the interaction integral I can be obtained as below where E * is defined as follows: If the SIFs of the auxiliary fields satisfy special conditions [19], the SIFs of the actual fields can be derived as Similarly, T-stress can also be derived by interaction integral and its expression is as follows: where f is defined as a line load, which is applied on the crack propagation plane and along the line of the crack, as shown in Figure 2. Therefore, the stress intensity factors and T-stress in the real fields can be obtained by selecting the appropriate auxiliary fields. Because the selection of the stress intensity factors and T-stresses in the auxiliary fields is different for the interaction integral, the selection of auxiliary fields in this paper is consistent with that of the Refs. [19,[32][33][34].
When the integral path Γ approaches the tip of crack, the relationship between the stress intensity factors and the interaction integral I can be obtained as below where E * is defined as follows: If the SIFs of the auxiliary fields satisfy special conditions [19], the SIFs of the actual fields can be derived as Similarly, T-stress can also be derived by interaction integral and its expression is as follows: where f is defined as a line load, which is applied on the crack propagation plane and along the line of the crack, as shown in Figure 2.  In order to facilitate subsequent observation and comparison, the dimensionless forms of the SIFs and T-stress are adopted as [11,14]

Criteria for Crack Propagation
The crack initiation and propagation of brittle materials such as rocks and ceramics in engineering are often predicted by appropriate mixed mode fracture criterion. Currently, the common mixed mode fracture criteria have always been the focus of research, e.g., the maximum tangential stress criterion [35], the criterion of local symmetry [5,36], the minimum strain energy density criterion [37] and the maximum energy release rate criterion [38]. Additionally, the related literatures [5,22,27] show that the MTS criterion or the local symmetry criterion could be used to accurately simulate the compression-shear mixed mode propagation process consistent with the experimental outcome under appropriate conditions, which confirms the accuracy of the MTS criterion and extends its applicable range of MTS fracture criterion. However, prompted by related researches [15,16,25,[39][40][41], T-stress has a significant effect on propagation angle of the crack, while the traditional MTS criterion only contains the singular term, which would inevitably make the generalized MTS criterion considering the influence of T-stress more accurate in predicting the propagation path. Therefore, the GMTS criterion would be adopted to predict the initiation and propagation trajectory of the crack concerning the frictional contact in this work.
According to Hua et al. and Tang [16,41], the tangential stress component of the linear elastic stress field around the crack tip can be written as: where r and θ are the polar coordinates of the crack tip.
According to the theory of the criterion, the crack still propagates in the direction of the maximum tangential stress around the crack tip. In addition, the angle of the fracture initiation ( Figure 3) can be obtained from the following forms where r c is a material constant for brittle materials and represents a critical distance from the crack tip. In addition, r c can be expressed as where K IC is the fracture toughness of the pure mode I; σ t is the tensile strength of the brittle material.

The WFM Concerning Frictional Contact
It is assumed that the studied Brazilian disc obeys the assumptions of small deformation, isotropic and linear elasticity. According to the research of Dong [1] and Hua [11,14], the expressions of the stress intensity factors and T-stress for CCBD under diametrical forces without involving the friction can be obtained.  Meanwhile, the equivalent stress intensity factor K eq can be adopted and the crack begins to expand when it is greater than or equal to K IC [42]. The equivalent stress intensity factor can be expressed as below:

The WFM Concerning Frictional Contact
It is assumed that the studied Brazilian disc obeys the assumptions of small deformation, isotropic and linear elasticity. According to the research of Dong [1] and Hua [11,14], the expressions of the stress intensity factors and T-stress for CCBD under diametrical forces without involving the friction can be obtained.
where h 1 (x, a), h 2 (x, a), h 3 (x, a), obtained by the boundary configuration method, are the weight function independent of the loading condition. σ r , σ θ and τ rθ are the stress components for the uncracked disk. In addition, the detailed formulas mentioned above can be found in the Refs. [43,44]. When the crack changes from pure mode I to mixed mode, the stress intensity factors and T-stress of CCBD can be further written as follows: where σ is the normal stress, where σ = P/(πBR). The relative length of crack α equals a/R. f ji and A ji (j = 1, 2 and i = 1, 2, . . . , n) are the coefficients [1]. However, the situation is much more complicated when the crack changes from pure mode II to compression-shear state so that the crack surface is closed. At this time, with the increase of the loading angle, the boundary of the crack surface may overlap or even cross, which would cause the normal pressure and friction on the crack surface [2]. In addition, when the crack surface is closed, CCBD is equivalent to a crack-free disk. The central crack of the CCBD is assumed to be a linear crack that is in contact under radial forces. The shear stress of crack surface changes due to the friction generated by the contact [45]. Assuming that the friction coefficient is µ, through the Coulomb friction law, the changed shear stress τ rθ can be obtained as follows [46]: where shear stress τ rθ should be greater than the friction, otherwise the crack would be in the state of cementation. Meanwhile, the formula for the mode II stress intensity factor considering friction can be further derived from Equation (16) as bellow According to the above analysis, the friction generated by contact will not affect the stress components (σ r and σ θ ) except the shear stress (τ rθ ). Therefore, the formulas for mode I stress intensity factor (K I ) and T-stress (T), involving the behavior of the friction, are consistent with those without considering the friction. The compression-shear crack surfaces cannot be intercalated due to the consideration of contact. Therefore, the K I in compression-shear fracture is taken as 0 [17]. Therefore, the mode I SIF concerning the frictional contact can be sorted into the following form Then, the dimensionless intensity factors and T-stress can be respectively written as bellow: The derived Equations (24)- (26) are the basis of GMTS criterion, which could establish the mathematical model for crack propagation.

Results and Discussion
The fracture test under radial forces was numerically performed using the MFEM. The advantage of the MFEM, implemented through the ABAQUS and PYTHON, is that it is based on the FEM using user defined subroutine, combining the GMTS criterion with the interaction integral method. During the simulation (Figure 4), the propagation of the studied crack is obtained by the elongation of many small cracks and its successful implementation also depends on the overall remeshing technique [47,48] and the method of incremental crack propagation [21,49]. In addition, considering the influence of mesh size h and increment of the crack propagation length i on the calculation accuracy, both h and i adopted the recommended values (generally, half of the internal parameters is chosen as the h and h/2 ≤ i < h) [19,25,49]. The combination of these technologies can make the MFEM have high accuracy and versatility for two-dimensional linear elastic problems under various loads. Meanwhile, the two-dimensional and isotropic CCBD numerical model ( Figure 5), concerning the friction contact on the crack surfaces, was adopted, which prevented the elements of crack surface from infiltrating each other. In addition, around the crack tip of the model, the second order plane stress grids of six-node element (CPS6) and eight-node element (CPS8), which are set with singularity, are used to ensure the accuracy of the calculation. In addition, due to the well crack closed state and crack growth paths obtained by inserting the metal sheets, the model parameters (Table 1), which come from the mixture of fine sand, moderate pozzolanic cement and water performed by Haeri et al. [22], are used for simulation. For the sake of simplicity, only the friction coefficient, loading angle and relative crack length were changed in the following work, while other model parameters remained fixed. The detailed calculation results and discussions of several case studies about the effects of friction contact on fracture parameters and crack propagation in this section are as follows.
obtained by inserting the metal sheets, the model parameters (Table 1), which come from the mixture of fine sand, moderate pozzolanic cement and water performed by Haeri et al. [22], are used for simulation. For the sake of simplicity, only the friction coefficient, loading angle and relative crack length were changed in the following work, while other model parameters remained fixed. The detailed calculation results and discussions of several case studies about the effects of friction contact on fracture parameters and crack propagation in this section are as follows.   obtained by inserting the metal sheets, the model parameters (Table 1), which come from the mixture of fine sand, moderate pozzolanic cement and water performed by Haeri et al. [22], are used for simulation. For the sake of simplicity, only the friction coefficient, loading angle and relative crack length were changed in the following work, while other model parameters remained fixed. The detailed calculation results and discussions of several case studies about the effects of friction contact on fracture parameters and crack propagation in this section are as follows.

SIFs and T-Stress Concerning Frictional Contact
The purpose of this section is to study the effect of friction on stress intensity factors and T-stresses at the onset of fracture. Figure 6 clearly shows the dimensionless K and T, calculated by the WFM and MFEM, versus the loading angle with different friction coefficients (µ = 0, 0.2, 0.4, 0.6 and 0.8) and relative lengths of the crack (α = 0.3, 0.5 and 0.7). By comparison, the results between MFEM and WFM have strong consistency with a maximum error within 3%, which shows that the proposed MFEM is reliable.  First, some variation rules of the stress intensity factors and T-stress with different loading angles and relative crack lengths can be observed, which is similar to the research conclusions of some scholars [14,16,20,25]. Meanwhile, when the crack is closed, the setting the frictional contact makes the normalized mode I SIF K I * = 0, while the K I * and T* remain the same trends no matter what the friction coefficient is. These results indicate that the SIF of mode I and T-stress are indeed independent of friction, which is consistent with the conclusions predicted by WFM.
However, it can be seen from the Figure 6, the friction has a significant influence on the mode II stress intensity factor when the crack is closed. In addition, as increasing the friction coefficient µ, the K II * decreases obviously when relative length of crack α and loading angle β have the certain values. For example, with the friction coefficient increases from 0 to 0.8, for the loading angle of 60 • , the K II * decreases 120.7% when α = 0.5. In addition, with the larger α and β, the effect of friction is more significant as the amount of increase in friction coefficient is constant. For example, as the friction coefficient increases from 0 to 0.6, the K II * decreases 14.6% when α = 0.3 and β = 40 • , while the K II * decreases 58.5% for α = 0.5 and β = 50 • . It is well known that the crack is in a viscous state when the K II * is less than 0 [50]. Figure 6 also shows that when µ is equal to 0 (the crack surface is smooth enough), the crack would not be in a sticky state from beginning to end with the loading angle increases from 0 to 90 • . However, with a larger friction coefficient, the loading angle corresponding to K II * = 0 decreases gradually, which indicates that the crack is in a viscous state earlier under the condition. At the same time, the larger loading angle and relative crack length will make this phenomenon more obvious.
In general, the foregoing results show that the friction has a significant effect on the mode II SIF after the crack is closed, which indicates that it is necessary to consider influence of the friction in the real contact state, but has no influence on the stress intensity factor of mode I and T-stress.

Crack Initiation
The effects of friction force, relative crack length and loading angle on crack initiation and propagation of CCBD specimens are studied in this section. For some of the reasons illustrated in Section 2.2, the generalized MTS criterion, considering the influence of Tstress, would be used to predict the initiation and propagation of the crack. It can be seen from the conclusion in Section 4.1 that a large friction coefficient will cause the crack to be in the adhesive zone earlier [50], which makes the normal propagation of the crack gradually difficult [16]. Meanwhile, the phenomenon is more obvious under larger relative length of crack α and loading angle β. However, many literatures have also shown some experimental results of the instability of compression-shear cracks, e.g., (1) the shifting phenomenon of crack initiation point [4,22], which indicates that the crack initiation point would shift to the center of the crack at a larger loading angle; (2) the multi-branching phenomenon of the crack [51], which generally occurs in compression-shear fracture under a certain loading angle. Therefore, in order to explore the law of crack initiation and propagation stably, the friction coefficient µ is equal to 0, 0.3 and 0.6 respectively and the relative crack length α = 0.3 and 0.5 with the loading angle β ≤ 60 • .
The crack propagation angle θ 0 , which is defined as positive in the counterclockwise direction, versus the loading angle β is compiled in Figure 7. It can be seen that there is a certain amount of the deviation between the crack propagation angle predicted by GMTS criterion and that predicted by MTS criterion in most regions. That is because K I , considering friction contact, is always 0 after the crack surface is closed and T-stress has a significant effect on propagation angle of the crack [40]. Meanwhile, there are some regions of loading angle with gradually decreasing the above deviation, because the T-stress at this stage changes from negative to positive, as shown in Figure 6. On the whole, the crack propagation angle predicted by GMTS increases with the increase of β. When Tstress is positive, the effect of friction becomes obvious and the crack propagation angle increases with the increase of friction. Further, a larger relative crack length magnifies the above effect. criterion and that predicted by MTS criterion in most regions. That is because KI, considering friction contact, is always 0 after the crack surface is closed and T-stress has a significant effect on propagation angle of the crack [40]. Meanwhile, there are some regions of loading angle with gradually decreasing the above deviation, because the T-stress at this stage changes from negative to positive, as shown in Figure 6. On the whole, the crack propagation angle predicted by GMTS increases with the increase of β. When T-stress is positive, the effect of friction becomes obvious and the crack propagation angle increases with the increase of friction. Further, a larger relative crack length magnifies the above effect.

Crack Propagation Trajectory
In this section, the dynamic effect would not be taken into account. First, this section compares the fracture trajectories obtained by MFEM and the experimental results [22],

Crack Propagation Trajectory
In this section, the dynamic effect would not be taken into account. First, this section compares the fracture trajectories obtained by MFEM and the experimental results [22], as shown in Figure 8. It can be found that, in spite of the setting of the contact, when the friction coefficient is 0, the crack would eventually deviate from the original crack line and extend towards the loading direction. Under this condition, the increase of relative crack length and loading angle will increase the deviation of the crack, which could not change the type of crack. According to the previous research [8], it can be found that this kind of crack is a tensile wing crack. However, as the friction increases, the crack propagation is significantly affected. Appl. Sci. 2021, 11, x FOR PEER REVIEW 13 of 16          Considering the symmetry of the CCBD, the propagation path extending from the right tip of the central crack is taken as the research object. In addition, the effects of different friction coefficients on growth trajectories of the crack are shown in the Figures 9 and 10. The crack propagation path is greatly affected, which is manifested that the crack propagation path deviate upward. At the same time, as increasing the friction coefficient, the amount of the upward deviation increases and the curvature of the path decreases. In addition, there are some influences on crack type due to the change of friction coefficient. It can be found that when the loading angle β is equal to 30 • , no matter what the friction coefficient is, the crack propagation path of the CCBD sample does not change significantly, which makes only the tensile wing crack occur. However, as β = 60 • , tensile crack or mixed tensile-shear crack would begin to occur. For example, when β = 45 • and µ = 0.3, the crack propagation path of the crack is like a straight line, without bending and extends in the direction of loading, which is called the tensile crack. Similarly, under the condition of this loading angle, mixed tensile-shear crack can be observed when µ = 0.6 and α = 0.5. According to the previous literature [50], when the values of α and β for the CCBD are small, the effect of friction is not obvious due to the incomplete contact of the crack surface. In addition, these conclusions can be embodied in the Figures 9 and 10. Therefore, based on the above conclusions, the appropriate loading angle and relative crack length should be required, which makes the friction have a non-negligible effect on the crack surface.

Conclusions
The present work fits into the stream of study and analysis of the friction impact on fracture parameters (SIFs and T-stresses) and crack propagation of CCBD specimens under diametric forces. Considering the frictional contact, the MFEM is developed to investigate the above effects. In addition, the following conclusions can be drawn from the simulation results of several examples: (1) The friction has a significant effect on the mode II SIF after the crack is closed, but has no influence on the stress intensity factor of mode I and T-stress. When the crack is closed, the mode II SIF decreases obviously as increasing the friction. (2) There are significant effects on the crack propagation angle and crack propagation trajectory due to the change of the friction after the crack is closed with the appropriate relative crack length and loading angle. (3) When T-stress is positive, the influence of friction becomes obvious and the crack propagation angle increases with a lager friction coefficient. (4) After the crack is closed, as increasing the friction, the amount of the deviation increases and the curvature of the path decreases. Furthermore, the crack type is easier to change with the increase of friction. Funding: This paper received no external funding.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.