Simulation of the Gri th’s Crack Using Own Method of Predicting the Crack Propagation

The paper presents the results of numerical simulation of brittle material fracture initiated by Griffi th’s crack situated at any angle with respect to the load direction. The simulations were performed with the Simulia Abaqus FEA system using the X-FEM method. The direction of the crack propagation was determined using the user’s subroutine written in Fortran, implemented in the Abaqus system. The authors programmed several conditions for the direction of the fracture propagation: maximum principal stress criterion, Ottosen-Podgórski criterion, three criteria based on displacements around the crack tip, and the MTS criterion, which is based on stress intensity factors. For the purposes of this study, two shapes of the crack tip were analyzed sharp and blunted. The relationship between the direction of the initial fracture and the direction of crack propagation was found. The simulations were carried out for the angles of the initial crack from 0° to 90°, every 10°. A theoretical analysis of the implemented criteria was also carried out and compared to numerical results. The infl uence of the shape of the crack tip on the simulation result was analyzed. It was shown that the simulation results closest to the theoretical results were obtained for the own implementation of the stress-based criteria. It has also been proven that the shape of the crack tip has little eff ect on the result if the fi nite element mesh is properly densifi ed. However, the sharp crack has a disadvantage for small initial crack angles. The Author’s method of predicting the crack propagation has been proven correct and eff ective. This work is also the basis for the further development of the described method.


INTRODUCTION
This paper presents the results of fracture simulation using the user's own implementation of the crack propagation criteria with user subroutines in the Abaqus FEA system. This method was described in the previous author's work [1]. The analyzed task is the propagation of the Griffi th crack located in an infi nite 2D area in a quasi-brittle material, subjected to a pair of uniformly distributed tensile loads. Six criteria of the fracture propagation angle have been programmed. The fracture was simulated using the extended fi nite element method (X-FEM) developed by Moes et al. [2], which is one of the popular methods of simulating a discrete crack in the Finite Element Method. This method consists of modifying the element shape function by adding an enrichment function (Heavyside function), which is responsible for dividing a fi nite element at any point without the need to modify the element mesh.
Unfortunately, this method has some drawbacks. In the Finite Element Method, the stress values at the crack tip are not available, but only at the surrounding integration points, so the direction of the crack propagation calculated by the program is often incorrect. However, the Abaqus system allows to introduce own criteria for determining the direction of the crack propagation using the socalled user subroutines. The user subroutine to defi ne the damage initiation criterion (UDMGINI) is the subroutine responsible for determining the direction of the crack propagation. This subroutine, written in Fortran, reads the appropriate input data Simulation of the Gri th's Crack Using Own Method of Predicting the Crack Propagation (most often these are user-specified values and values read during the simulation, such as stresses, displacements, coordinates) and as a result, sends the crack propagation vector and the criterion value to the Abaqus solver. The criterion value, in the simplest case, is the ratio of tensile stress to tensile strength. If the value of the criterion exceeds 1, the program continues the crack using the X-FEM method to the next element. This algorithm is repeated in each iteration step of calculations.
There are few publications in the literature describing the crack simulation with the use of the UDMGINI subroutine in the Abaqus system. Most of the works, however, describe the cracking of steel [3,4], composites [5,6], and polymers [7]. However, no studies describe their own criteria in quasi-brittle materials such as concrete or rocks, which will be analyzed in this paper.
Griffith's crack simulations are also an interesting research topic. In the literature, there are publications by Yin et al. [8] or Dewapriya et al. [9] where the simulation of symmetrical fracture in graphene was performed with the use of Molecular Dynamic Simulations, or, for example, the work by Wang et al. [10], in which the kinked Griffith's crack was analyzed using the Discrete Element Method. The topic of kinked crack has appeared in the literature for a long time [11,12]. In this paper, the angle of propagation of the Griffith's kinked crack will be verified by computer simulations using own implementations of the crack propagation criteria.

General description of the problem
The analyzed task is Griffith's crack shown in Figure 1. It is a crack located inside an infinite twodimensional area stretched by a pair of uniformly distributed loads. In the figure, the angle β is the angle between the initial crack and the horizontal axis (perpendicular to the load), and the angle α is the predicted angle between the initial crack line and the actual crack formed under load. In this section, the relationship between the initial crack angle and the fracture angle will be described. There are many criteria that determine the direction of the crack propagation. The following criteria were used in this study: the maximum principal stresses criterion, the Ottosen-Podgórski criterion, three criteria based on displacements around the crack tip, the maximum circumferential tensile stress (MTS) criterion, and the maximum energy release rate criterion (MERR) which are based on stress intensity factors. The listed criteria are described in the following subsections.

Stress-based criteria
The criterion of the maximum principal stresses assumes that the fracture will be guided in the direction in which the maximum principal stress (tension) at a constant radius around the crack tip locally reach the minimum value, which corresponds to the maximum value of the stress gradient when viewed from the material side. The values of the stresses around the crack tip in mode I of fracture are described by the following formulas in a polar coordinate system [13]: where: r is the distance between the crack tip and the analyzed point, θ is the angle between the crack line and the analyzed point, K I is the stress intensity factor in mode I. This situation is illustrated in Figure 2.
Similarly, the stresses around the crack tip for mode II are shown in the following formulas: Stress intensity factors in mode I and II for theoretical kinked Griffith's crack can be described as: where: K I0 is the stress intensity factor in mode I for β = 0°. This allows to apply the formula for the maximum principal stress as the maximum of the two formulas below: Maximum principal stresses around the Griffith's crack for radius r = 1 are presented in Figure 3. It can be seen that for mode I (initial crack at angle β = 0°) the local minimum is obtained for angle 0°, for mode II (initial crack at angle β = 90°) is for angle 70.529°, and for mixed-mode, these should be values in between. Based on these results, it was found that the fracture angle α is the value of this minimum. The relationship between the initial angle β and the fracture angle α for all criteria will be presented at the end of this chapter.
The Ottosen-Podgórski (O-P) criterion by Podgórski [14,15] is expressed as a relationship of three alternative stress tensor invariants:  where: P(J) is a function describing the cross-section of the boundary surface by a deviator plane σ 0 = const., and is given by the equation P(J) = cos(⅓arccosξJ -φ), J here stands for alternative stress tensor invariant and ξ, φ, C 0 , C 1 , C 2 are parameters that can be calculated from material constants [15,16]. In this case, instead of the maximum principal stresses around the crack tip, the values of the material effort function ω is calculated as: where: σ 0 and τ 0 are the normal and tangential octahedral stresses determined at the tested point around the crack tip in the polar system, respectively; and σ f and τ f are the values of critical stresses proportional to σ 0 and τ 0 , respectively.
By inserting stresses σ x , σ y and τ xy from Eq. (4) as data around the crack tip, the dependence of the angle around the crack tip θ and the material effort ω can be found. It is very similar to the graph of the maximum principal stresses in Figure 1, so ω also takes a local minimum in the direction of the crack propagation.

Displacements-based criteria
Similarly to the stresses, the displacement formulas around the tip of the Griffith's crack for modes I and II can also be found [13]: where: = � 3 − 4 for plane strain (3 − ) 1 + for plane stress (9) where: μ is the shear modulus, ν is the Poisson's ratio (for all calculations assumed as 0.2), u x and u y are horizontal and vertical displacements at any point in the polar system around the crack tip, respectively.
The diagram of these displacements is shown in Figure 4. The displacements can be also rotated to the local coordinate system: where: u r will be called as displacements along the radius and u θ as displacements perpendicular to the radius. u are resultant displacements.
These displacements are shown in Figure 5. According to the diagram in mode I, three different conditions can be deduced. The crack will occur in the direction 0°, in which: • the displacements u r reach a local minimum, • the first derivative of the displacements u θ reaches the local minimum, • the resultant displacements u reaches the local minimum.
In mode II, these conditions are met for the following angles: • 45.275° for minimum displacements u r , • 52.776° for a minimum of the first derivative of the displacements u θ , • 0° for minimum magnitude displacements u.
As it can be seen, these values are significantly different from these obtained by the stress-based criteria. The relationship between the initial angle β and the fracture angle α will be presented later.

Criterions based on stress intensity factors
The next analyzed criterion is the maximum circumferential tensile stress (MTS) criterion by Erdogan and Sih [11]. This is a criterion based only on the stress intensity factors, where, the fracture occurs in the direction for which the circumferential stress σ θθ described by the formula below reaches the local maximum: Analyzing the graph in Figure 6 it turns out that the maximum of this function is for θ = 0° in mode I of fracture and 70.529° in mode II.
The last criterion is the maximum energy release rate criterion (MERR) developed by Hussain et al. and Wu [17,18]. It defines the crack angle where the energy release rate takes the maximum: (12) Where: where: E -is the young modulus assumed 1 for simplicity.
The maximum energy release rate is obtained for θ = 0° in mode I and 70.529° in mode II. Figure  7 shows a diagram of the relationship between the initial crack angle and the fracture direction. As it can be seen, all criteria based on stresses and stress intensity factors give practically the same results, while criteria based on displacements give very different results. In the first of these criteria, for an initial angle less than 60°, the fracture will be at an angle greater than 60° (right-down direction) and for an initial angle greater than 60°, the fracture will be at an angle less than 60° (right-up direction). The crack is horizontal for an initial angle of 0° and 60°.

NUMERICAL METHOD OF DETERMINING THE DIRECTION OF THE CRACK The built-in method in Abaqus
The main goal of the authors of this study was to verify their own procedure for determining the crack path in the Abaqus FEA system. In the beginning, it was also decided to test the effectiveness of the built-in method of the Abaqus system. The default criterion that works with the X-FEM method is the maximum principal stress criterion. The method of its operation is presented in this subsection based on [19].
Abaqus in the X-FEM method in a plane stress state allows the use of only four-node CPS4 elements, with four integration points near each node. However, it is possible to select elements with reduced integration CPS4R (with one integration point in the middle of the element). CPS4 elements have been selected for all simulations. When a fracture reaches a finite element, its shape function is enriched. Abaqus reads the stresses at the integration points (Gauss points) and then calculates the principal stresses from them. These stresses are averaged from four points and fracture is conducted when these stresses exceed the tensile strength specified by the user. The direction of the crack is at the same angle as the averaged angle α of rotation of the stress tensor components to the principal stresses. This situation is shown in Figure 8.
This is a very simple method but has many disadvantages. Only the stresses at single integration points are always considered. If an element with four integration points is selected, the values from these four points are averaged. Thus, the stresses in the finite elements adjacent to the element containing the crack tip are still omitted. Another disadvantage is that the stresses are taken from the integration points. A more correct approach would be to take the stress values from a point close to the crack tip.
This algorithm is repeated for every iteration of the calculation to find the solution. If the solution is found the crack is lead to the next element and the next load increment is calculated.

Author's method
A detailed description of the own method of predicting the direction of the crack propagation can be found in the author's previous publication [1]. Here only its assumptions will be presented.
In each load increment for the maximum principal stress criterion and the Ottosen-Podgórski criterion the algorithm is as follows: • The coordinates of the crack tip are searched, • Stresses at several dozen closest integration points around the crack tip are read (50-100 points were adopted to ensure that the result was sufficiently accurate and so the calculations did not last unnecessarily long), • From these stresses, the maximum principal stresses or the material effort ω for the Ottosen-Podgórski criterion are determined (Fig. 9a), • These values are reduced to the unit radius (Fig. 9b), • The obtained values are analyzed for their dependence on the angle θ around the crack tip, • The 6th order polynomial is fitted to these points by the least square method, • The local minimum of the polynomial closest to the angle from the previous load increment is found by the bisection method (Fig. 9c).
The reduction of the value to the unit radius means that the values are correspondingly decreased if the distance of their point from the crack tip is less than 1 or increased if the distance is greater than 1. This assumption is derived from the formula (1) in which it can be seen that the stresses are decreasing at a rate of 1/ √r. This algorithm was programmed in Fortran as Abaqus User Subroutine and the code is compiled during each simulation.
The main advantage of this method is that stress from multiple integration points is used in the calculation, not just stress at a single finite element. Moreover, in the built-in method, the tangential stresses, which are often disturbed near the crack tip, have a big influence on the results. In the own method, the stresses at a greater distance from the crack tip, which are already more stabilized, are taken into account. In addition, the tangential stresses are of less importance. This makes the crack path much more correct and smooth, as shown in the previous publication [1].
Criteria based on displacements have been implemented very similarly. Displacements are read at element nodes, not at integration points. So, there are fewer points, which is associated with lower accuracy of the results. Displacements in the horizontal and vertical directions have to be properly corrected. In the case of stresses, theoretically, they grow to infinity at the crack tip, while the displacements according to the formula (8) tend to 0. In fact, in various simulations, depending on the model, the displacements at the crack tip will not be equal to 0, hence the displacements at all points need to be decreased by values obtained at the crack tip. Then the displacements along the radius u r , the displacements perpendicular to the radius u θ, and the resultant displacements u are calculated as in equation (10). In the same way, as in the previous criteria, they are placed on the graph, reduced to a unit radius, according to the assumption in the formula (8) that the displacements increase at the rate √r. Again, the polynomial is matched to these values and depending on the selected criterion, the minimum displacements u r , the minimum derivative of displacements u θ, or the minimum displacements u is searched.
In the MTS criterion, the direction of the crack propagation is determined by the minimum of the function (11). In this case, it was necessary to find the stress intensity factors at the crack tip during simulation. For this purpose, the direct method was used: (15) where: u x and u y are the differences in displacements on both sides of the crack along the crack and perpendicular to the crack, respectively (u x = u 1r -u 2r , u y = u 1θ -u 2θ as in Fig. 10).
To find only the direction of the crack propagation, only the ratio between the stress intensity factors is enough, which, based on formula (15) is: (16) This means that the calculation algorithm in this criterion is very similar to the displacementbased criteria. In this case, instead of searching the angle for which the minimum displacement values were obtained, we find the displacements u r and u θ for the angles -180° and 180° and these are the values u 1r , u 2r , u 1θ , and u 2θ .  Figure 11a. It was ensured that the size of the model was large enough in relation to the size of the initial crack. The model mesh is shown in Figures 11b-c. The entire model had dimensions of 100×100, and the length of the initial crack was 8, while the mesh size ranged from 0.02 near the crack tip to 0.5 around the initial crack, to 5 for the rest of the model. The load was applied as a distributed force on the upper and lower edges of the model. The value of the load does not matter, because only the direction of the crack propagation is needed. 10 models were made, with the initial crack at an angle β from 0° to 90°, every 10°.
One problem with this task is that it is point symmetric. In fact, the fracture should appear symmetrically on both sides of the initial crack. Unfortunately, the X-FEM method in the Abaqus system, and especially the described own method, is not suited to simulate more than one crack. For this reason, it was decided to dense the mesh in only one crack tip to force the fracture there. The main purpose of the following calculations is to find the initial crack propagation angle.
After a crack appears, the angle of propagation will probably change with the distance from the initial crack. Therefore, a big problem is properly interpreting the calculated initial angle. A crack through the first finite element is always incorrect because there are no crack tip coordinates defined by the X-FEM algorithm. Before fracture, the coordinates of the integration point with the greatest principal stresses are taken as the coordinates of the hypothetical crack tip. Unfortunately, this raises a certain problem presented in Figure 12a. The program finds the integration point of the highest principal stresses and conducts the fracture at an angle similar to the one obtained in the theoretical analysis. The figure shows that the fracture cannot start from the model edge, which also misleads subsequent calculations, because then there are two crack tips.
To eliminate this problem, a blunted crack can be used, but this one also has its drawbacks -in close proximity, the crack path shape is not consistent with the theoretical assumptions of this task. For this reason, two pairs of all simulations were performed -with a sharp crack and blunted crack.
Because the fracture in the first element is always incorrect the crack needs to stabilize. It was decided to measure the crack angle as an average for the first 2-3 cracked finite elements.

RESULTS
Simulations for sharp and for blunted crack were performed for all criteria. However, simulations for the criterion of minimum magnitude displacements u were omitted, as the results of the theoretical analysis completely differ from the results for the other criteria (Fig. 7). Also, the MERR criterion is not implemented at the moment. Figure 12b shows the map of maximum principal stresses obtained in the simulations around the tip of the initial crack at the angle of 0°. Lower stresses in the horizontal direction are clearly visible here, which are responsible for the local minimum as in the theoretical approach to this task. According to the assumptions of the own method of predicting the crack propagation direction, the crack should be led towards this local minimum. Figure 13 shows the diagrams of the relationship between the angle of the initial crack and the predicted direction of the fracture. Here, the results for sharp and for blunted crack were compared with the directions obtained by the theoretical analysis as well as the results of the procedure built into the Abaqus system. The comparison of the implemented procedures with the results obtained in the laboratory and "in situ" tests is given in the previous authors paper [1].
It can be seen from these figures that the propagation angles obtained in the simulations are closest to the theoretical solutions for the stress-based criteria (Fig. 13a-b). In these simulations, a change in the crack propagation angle after exceeding the angle of 60° of the initial crack is also noticed, as it was for theoretical calculations. It turns out that the criterion of maximum principal stresses built into the Abaqus system (Fig. 13f) is less precise than the above criteria implemented by the authors.
The displacement-based criteria give results similar to the stress-based criteria (Fig. 13c-d), however, they are far from theoretical calculations. In fact, in this case, the crack runs close to the horizontal direction. This may be related to a slightly different distribution of displacements in the simulation results compared to the theoretical displacements from the formulas (8). It may also be related to the lack of point symmetry after the initiation of the crack in the first element. Another reason may be the incorrect assumptions of these criteria themselves.
The MTS criterion also gives incorrect simulation results (Fig. 13e), especially for large angles of the initial crack. In this criterion, the stress intensity factors were calculated based on displacements which, as mentioned earlier, may not be correct. In addition, the displacements important for determining the stress intensity factors are read for the angles -180° and 180° relative to the angle from the previous increment, which means that these are displacements near the initial crack (see Fig. 14) which are subject to disturbances.
The above problems with the MTS criterion make the crack path predicted by this criterion very inaccurate (see Fig. 15b). For comparison, Figure 15a shows a smooth fracture line obtained from the criterion of maximum main stresses. The MTS criterion for large initial fracture angles gives completely incorrect results as shown in Figure 15c (see also Fig. 13e).
Finally, Figure 16 shows the comparison of the crack path obtained in the simulations and the theoretical analysis for the maximum principal stresses criterion, with the initial crack every 10°. As it can be seen, despite some differences visible in Figure 13a, the direction of the reality, however, the stress-based criteria implemented using the own method of predicting the crack propagation direction give the best results.

CONCLUSIONS
The paper presents the results of simulation of the Griffith crack propagation inclined at any angle with respect to the direction of the tensile load. The direction of the crack propagation appearing at the tip of the initial crack was analyzed. The X-FEM method in the Simulia Abaqus FEA system was used to simulate the crack propagation. The following criteria were analyzed: the maximum principal stresses criterion, the Ottosen-Podgórski criterion, the criterion of minimum displacements along the radius, the criterion of the minimum derivative of displacements perpendicular to the radius, and the maximum circumferential tensile stress criterion (MTS), all implemented by the Abaqus User Subroutines. The obtained results were confronted with the criterion built into the Abaqus system. The results were compared with the theoretical calculations. It turned out that the results obtained for the stress-based criteria are the closest to reality. Unfortunately, the own implementation of criteria based on displacement and stress intensity factors has some drawbacks, leading to incorrect results. Perhaps for the MTS criterion, it is necessary to find another way to determine the stress intensity factors.
It is also worth noting the slight influence of the shape of the crack tip. A sharp crack may make it difficult to find the correct fracture of the first element, while a blunted crack does not have the defined exact tip coordinates. To overcome this problem, the best solution would be to use triangular elements that have integration points at the vertices of the triangle, and therefore precisely at the model nodes. Unfortunately, Abaqus does not allow the use of triangular elements, as they currently do not work with the X-FEM method.
In the future, it is planned to thoroughly investigate and improve criteria based on displacement and stress intensity factors. The long plans include creating own FEM code, in which it would be possible to use triangular elements and in which the authors would have full control over the algorithm. Currently, only an incomplete modification of the algorithm is possible, which is allowed by the Abaqus system with the user subroutines.