Numerical Simulation of Fracture Behavior of Brittle Solids under I/III Loading

+is study used numerical analysis to carry out a large number of numerical model calculations based on a new semicircular bend (SCB) model. Instead of existing numerical computation methods (J-integral), the M-integral method was used to calculate the mixed mode fracture parameters (KI, KII, and KIII), investigate the influence of the geometry and material parameters on the fracture behavior under mixed mode I/III loading, and predict the fracture path. +e results revealed that the limitations in the scope of the mixity parameter M in the previous studies can be overcome to a certain extent. +e range of M was established under different Poisson ratios and can be used as a reference for actual material testing. +e simulation path is in good agreement with the experimentally obtained fracture path, and the proposed method can be used to simulate the fracture path under mixed mode I/III loading.


Introduction
Fracture mechanics focuses on the mechanical properties of materials and structures with cracks. e cracks can be inherent in the material or generated during manufacturing, and the existence and growth of these cracks decrease the bearing capacity of the structure or even cause its failure. e semicircular bend (SCB) is a classic fracture mechanics test that is widely used to investigate solid materials owing to its inherent favorable properties, such as its simplicity, minimal machining requirement, testing convenience, and ease of reaching tensile failure as reported by Mahinda and Kuruppu [1]. e SCB specimen was proposed by Chong [2][3][4] and has been gradually improved. Additionally, it has been widely used for determining the fracture parameter K IC , and its potential application in the mixed mode fracture test has been investigated [5][6][7][8][9]. e investigation of the SCB specimen mainly includes fracture toughness testing and determining the fracture parameters [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28], discussing and revising the limitations of the extension criteria [11,13,16,25,29,30], and numerically predicting the fracture path [14,20,23,27,[31][32][33]. Owing to the limitation of the specimen shape, the abovementioned studies focused on the fracture parameters of mode I, mode II, and mixed mode I/II. e out-of-plane shearing mode (mode III) has not been as extensively investigated, although it exists in actual situations. Additionally, existing testing methods are highly dependent on the machine and test fixture and mainly focus on metal materials, ceramics, and other artificial materials [34][35][36][37][38].
Because engineering accidents, such as rock fracture and pavement rupture, occur frequently, the investigation of mode III has attracted a great amount of interest. Moreover, observations have shown that flat cracks tend to reorient themselves to oblique planes during propagation and can grow under mixed-mode I/III conditions [39,40]. erefore, test methods and specimens under the abovementioned conditions have been proposed for rocks, concrete, asphalt, and other materials. Berto et al. applied the strain energy density (SED) criterion to pure modes and mixed mode I/III loading conditions [41][42][43][44][45]. Aliha et al. [46][47][48][49] proposed a new test configuration called the ENDB specimen for investigating the fracture behavior under the mixed mode I/ III. Additionally, this configuration can achieve an arbitrary composition of mode I and mode III. Notably, however, more variables need to be controlled in the experiment. Linul et al. [50,52] conducted experiments with foam materials and obtained fracture toughness data for the mixed modes I/II and I/III. Subsequently, Pirmohammad [53] proposed a new SCB specimen type, whose mixed mode I/III fracture behavior is easier to investigate using the proposed specimens with a tilted crack and further investigated the fracture behavior of asphalt concrete at low temperature.
In previous studies [49,53], the change of the geometric parameters and material parameters had a significant effect on the mixed mode I/II loading of SCB specimens. However, only few studies have investigated whether the abovementioned conditions affect the fracture behavior of mode I/ III. In this study, a relatively accurate numerical analysis method was used to investigate the feasibility of realizing mixed mode I/III loading and the effects of the specimen's geometric parameters and Poisson's ratio on these loading conditions. Finally, the fracture path was simulated and analyzed based on the fracture parameters.

3D Numerical Analyses
In the previous numerical studies [41,53,54], the fracture parameters were directly calculated using the J-integral method in the Abaqus software. However, according to Li, the J-integral is more suitable for pure mode I fracture problems [55]. For the mixed mode fracture, the equivalent domain equation of the interaction integral (M-integral) [56,57] is the most accurate method for calculating the fracture parameters. e SCB macrostructural model was established in Abaqus; then, the mesh was remeshed using adaptive grid technology in franc3D, and the bottom inclined crack was inserted. Finally, different K I , K II , and K III values were obtained by the M-integral. e schematic diagram of the semicircular geometry is shown in Figure 1; the radius of the SCB specimen (R), specimen thickness (t), and half of the loading span (S) in the SCB test were determined by modeling in Abaqus. Figure 2(a) shows the Abaqus grid model with C3D10 cells. Figure 2(b) shows the mesh model after defining the crack length (a) and the angle of crack (α). Figure 3 shows the element distribution and crack tip types. A singular element was inserted at the crack tip, and conservation integral evaluation was carried out around the two element rings at the crack front. e inner ring of a 15-node singular wedge-shaped element and the outer ring of a 20-node hexahedral element with symmetric meshes were used to reduce the local discrete errors. Notably, the fracture parameters (K I , K II , and K III ) calculated at the surface contact zone are not reliable, because the local plane strain condition is no longer maintained. e element at the crack tip may be severely distorted, as shown in Figure 2(c), and the grid symmetry may be lost. Moreover, Figure 2(c) shows that the transition region at the crack tip has pyramidal elements, while the global region is a tetrahedral element type.
In the created models, radius of SCB specimen (R), specimen thickness (t), and applied load (P) are invariants, respectively, 75 mm, 32 mm, and 40.3437 kN. Geometrical and material parameters such as the half of loading point span ratio (S/R), crack length ratio (a/R), crack angle (α), and Poisson's ratio (v) were varied in the following ranges: 5,10,15,20,25,30,35,40,45,48,50,52,55, 60}, and v � {0.1, 0.2, 0.3, 0.4}. When α is 0°, it is pure mode I fracture. With the increase of α, mixed mode I/III fracture can be obtained. e fracture parameters (K I , K II , and K III ) can be directly extracted using the following expressions. For linear analysis, we can add two valid solutions, and the result is a valid solution, as follows: (1) Let us consider the corner mark (1) solutions as the Abaqus results and the corner mark (2) solutions as the solutions that we can select. ese can be substituted into the expression for the J-integral as follows: According to Betti's reciprocal theorem, the following relationship holds: By collecting terms, the following relationship can be obtained: with e crack-tip energy release rates can be determined from Irwin's crack closure integral for small scale yielding assuming plane strain conditions, as follows: 2 Advances in Civil Engineering By equating the two definitions for the M-Integral, we can obtain the following relationship: where We used the Abaqus results for solution (1) and selected the three simple auxiliary solutions (2a), (2b), and (2c). Table 1 presents the values of the asymptotic solutions for the pure mode cracks.
From the analytical expressions for the crack-front fields, we can obtain the following: By substituting the value of K (2) , reported in Table 1, in equation (9), we can obtain three equations for the unknown K (1) 's, as follows: where v is Poisson's ratio, Γ is the integral curve away from the tip of the crack, δ 1j is the Kronecker delta, q is a function that is equal to one at the crack tip and zero on the boundary of the integration domain and can be interpreted as a virtual crack extension, q t is the value of the q function along the crack front, and L is the length of the cylindrical domain along the crack front.
ij is the solution provided by the Abaqus results.
K I , K II , and K III can be expressed by the dimensionless parameters Y I , Y II , and Y III , as follows: To verify that the numerical simulation method can accurately realize mixed I/III mode fracturing, the K I , K II , and K III of the crack tip along the specimen thickness were adopted after normalization processing as shown in Figure 4. Here, K n (n � I, II, III) is the fracture parameter and K Im is the maximum value of mode I. e crack length ratio (a/R) of the model is 0.26, one half of the loading point span ratio (S/R) is 0.8, and the crack angle (α) is 0°, 20°, 40°, and 60°, respectively. Figure 4(a) shows that, when α was 0°, fracture occurred under the pure mode I condition; K I along the direction of the specimen thickness remained approximately unchanged within a certain range (0.1 < z/t < 0.9) and declined in the free surface, which is consistent with the numerical simulation results reported by Aliha and Saghafi [54]. Moreover, z/t was considered as 0.5 for the K I value of the specimen. Figures 4(b)-4(d) shows that, as α increased, the effects of modes I and III dominated within a certain range (0.2 < z/ t < 0.8) along the thickness direction of the specimen. e proportion of mode II was small and could be ignored when z/t was 0.5, which is consistent with previously reported numerical analysis results [46][47][48][49]. Considering that the surface contact zone calculation of K I , K II , and K III is not reliable, and according to the traditional definition of the singular stress field close to the free surface, it is widely accepted that the crack tip stress field is different, as discussed by Bažant and Estenssoro [58]. Hence, K I , K II , and K III close to the free surface should not be considered, and this numerical method is valid and feasible for testing the mixed mode I/III fracture. Notably, in mixed mode I/III loading, the crack tip along the thickness of the specimen completely corresponds to the point where z/t is 0.5. erefore, subsequent calculations considered this point as the value of K I , K II , and K III of the specimen.  e mixity parameter M e can be used to describe the relative contributions of mode I and mode III, as expressed by equation (15). Under the condition of pure mode I, the value of M e is equal to one, and decreases as the contribution of mode III increases: To quantify the role of mode I and mode III in the numerical simulation, we propose to substitute equation (16) for equation (15), and the parameter |M e | can vary from 0 to 1, corresponding to the pure mode III and mode I loading conditions, respectively:

Numerical Calculation Results and
Analysis. e crack length ratio (a/R) remained unchanged at 0.267 and a � 20 mm. e influences of the crack angle (α), half of the loading point span ratio (S/R), and Poisson's ratio (v) for Y I , Y III , and |M e | were investigated. Figures 5 and 6 show the variations of Y I , Y III , and |M e |, with a crack inclination angle (α) for different values of S/R and v in the SCB specimen.  Advances in Civil Engineering e following conclusions were drawn from image analysis: (a) Y III remained at zero, only when α was 0°. For the same α, the increase of S/R and v promoted the Y I . But for the Y III , the increase of S/R and v led to its decrease and increase, respectively.  Advances in Civil Engineering (α) for different values of a/R and v in the SCB specimen. Figure 9 shows the |M e | min obtained with the change of a/R and v. e following conclusions were drawn from image analysis: (d) e value of α was 60°, corresponding to the |M e | min of all SCB specimens. As v and a/R became smaller, a smaller |M e | min value was obtained and the minimum value was 0.39. In this case, the contribution of model III was obviously higher compared with that of model I. Additionally, for smaller a/R, v only had a minor influence on |M e | min ; however, when a/R was larger, |M e | min was more sensitive to the variations of v.
Notably, the rule discussed in the previous section is not described.

Fracture Path Simulation.
Using the K I , K II , and K III values of the crack tip and the maximum tensile stress criterion (MTS), three model groups (S/R � 0.8; a/R � 0.267; α � 0°, 30°, 45°, and 60°) were obtained for the crack failure and growth path, as shown in Figures 10 and 11. To facilitate the comparison, the failure mode after the image was flipped horizontally is in good agreement with the previous results [53].
Under mixed mode I/III loading conditions, the crack surface growth morphology was different from the inplane fracture. Figure 12(a) shows that, when an in-plane fracture (mode I, mode II, and mixed mode I/II) occurred, the crack surface remained in the same plane despite a certain deflection. However, Figure 12(b) shows that, when the mixed mode I/III fracture occurred, the crack surface was twisted and deviated from the direction of the crack surface by forming an angle called the kink angle.
Considering the direction of the crack kink angle (ϕ), the normal stress and shear stress components on the new plane can be obtained by carrying out coordinate transformation as follows: where x, y, and z are the original axes and x′, y′, and z′ are the axes of the crack plane after deflection. e transformed stress intensity factor can be expressed as follows: e relationship between the release rate of the mechanical energy and the angle is expressed as follows: By superposing the contributions of mode III, the corresponding transformed stress intensity factor can be expressed as follows.
By substituting equation (18) into equation (19), the expression of the release rate of the mechanical energy changing with the angle can be obtained for the mixed mode I/III. Figure 13 shows the normalization of the mechanical energy release rate along with the change of the crack kink angle (ϕ).
As can be seen in the graph, in the case of pure mode I fracturing (K III /K I � 0), G(ϕ)/G(0) became maximum at the kink angle (ϕ) of 0°. erefore, the pure mode I fracture extended in the direction of the original when the kink angle was 0°. When the superposition of the mode III contribution was added (K III /K I � 0.25), as shown in Figure 13, the crack kink angle deviated from the original crack plane by an angle that can be considered as a correction by the contribution of mode III. Finally, the crack moved up in the direction perpendicular to the maximum principal tension and changed from mixed mode I/III to pure mode I. e numerical simulation results presented in Figure 10 reveal that the extension surface gradually returned to a plane from a curved surface close to the loading point and finally fractured close to the loading point.
is phenomenon is consistent with the fact that the fracture surface is divided into three regions [53]. Table 2 presents the kink angle between 0.2 and 0.8 at the crack front with different α, as extracted from Franc3D. e SCB specimen parameters were S/R � 0.8, a/R � 0.267, and v � 0.2. Additionally, it was observed that the kink angle of the crack is related to α, and a larger α resulted in a greater kink angle.

Conclusions
e following conclusions were drawn from this study: (a) Numerical modeling combined with Abaqus and Franc3D analysis was used to investigate the influence of half of the loading point span ratio (S/R), crack length ratio (a/R), crack angle (α), and Poisson's ratio (v) on the behavior of mixed mode I/ III loading. (b) e increase of S/R, a/R, and v promoted the dimensionless parameters Y I . But for the Y III , only the increase of a/R and S/R led to its decrease. As α increased, the sensitivity of Y I decreased with the changes of v, while the sensitivity of Y III became maximum at |Y III | max . (c) e obvious approach toward reducing the mixity parameter |M e | is to decrease α. Additionally, it is possible to reduce |M e | by changing S/R or a/R. Moreover, the SCB specimen materials with small v exhibited the best reduction effect. (d) e |M e | min values of the SCB specimens were obtained under different a/R and v, which can be useful as a references for SCB specimens of different materials subjected to mixed mode I/III loading tests. (e) e numerical parameters obtained by the M-integral were used to simulate the crack growth and failure path, which are in good agreement with previous experimental results. e relationship between α and the kink angle (ϕ) was quantitatively analyzed, and their positive correlation was confirmed. e proposed method can be used to simulate the fracture path under mixed mode I/III loading.

Data Availability
All data used to support the findings of this study are included within the article.

Conflicts of Interest
ere is no conflicts of interest regarding the publication of this paper.