Global cracking elements: a novel tool for Galerkin-based approaches simulating quasi-brittle fracture

Following the so-called Cracking Elements Method (CEM), recently presented in \cite{Yiming:14,Yiming:16}, we propose a novel Galerkin-based numerical approach for simulating quasi-brittle fracture, named Global Cracking Elements Method (GCEM). For this purpose the formulation of the original CEM is reorganized. The new approach is embedded in the standard framework of the Galerkin-based Finite Element Method (FEM), which uses disconnected element-wise crack openings for capturing crack initiation and propagation. The similarity between the proposed Global Cracking Elements (GCE) and the standard 9-node quadrilateral element (Q9) suggests a special procedure: the degrees of freedom of the center node of the Q9, originally defining the displacements, are"borrowed"to describe the crack openings of the GCE. The proposed approach does not need remeshing, enrichment, or a crack-tracking strategy, and it avoids a precise description of the crack tip. Several benchmark tests provide evidence that the new approach inherits from the CEM most of the advantages. The numerical stability and robustness of the GCEM are better than the ones of the CEM. However, presently only quadrilateral elements with nonlinear interpolations of the displacement field can be used.

For evaluation of the ease of implementation, the reliability, robustness, and the efficiency of numerical methods for simulation of quasi-brittle fracture, the following criteria are listed in form of questions with "yes" or "no" answers: 1. Does the new method show great differences comparing to the other continuum-based methods? 2. Does it provide results that strongly depend on the discretization? 3. Does it require reprocessing such as remeshing? 4. Does it introduce extra degrees of freedom other than displacements? 5. Does it need to predict the crack path (position as well as orientation)? 6. Does it require more computing effort than a traditional finite-element method (FEM)-based method such as the smeared-crack approach?
A good method is one with predominantly negative answers to these questions. In this article, a novel numerical approach is presented, with "No" as the answer to all of the above questions. It is named the global cracking elements method (GCEM). Its characteristics are as follows: 1. The GCEM is developed on the basis of the original cracking elements method (CEM), presented in References 1, 2, 12, and 16. It is embedded in the conventional continuum-based framework. It may be considered as a conventional FEM approach with a special type of elements. 2. Unlike the smeared-crack approach, 17 the discretization dependency of the GCEM is negligible. To demonstrate this advantage, most of the numerical examples in this work are characterized by irregular meshes which were automatically built by Gmsh. [18][19][20] Besides, the mesh used in the GCEM can be relatively coarse. 3. For putting it more precisely, the GCEM is a type of strong-discontinuity-embedded approach (SDA), [21][22][23] which obviously does not need remeshing. 4. Unlike phase field methods, 24,25 extended finite-element methods (XFEMs), 26,27 or numerical manifold methods, [28][29][30] the GCEM does not require a precise description of the stress state at crack tips and does not need extra degrees of freedom. * 5. Unlike the traditional SDA and the XFEM, the GCEM does not need a crack-tracking strategy. 31,32 Disconnected element-wise cracking segments, passing through the elements, are used for representing crack paths. The orientation of the crack is determined locally, providing self-propagating cracks. 1 6. The global cracking elements used in the framework of the GCEM, are formally similar to the nine-node quadrilateral elements. They provide a symmetric and sparse stiffness matrix. 16,33 Thus the numerical efficiency is reasonably high.
The main features of the GCEM are as follows: • The designed special type of global cracking elements is formally similar to the Q9 element. The shape functions, the local stiffness matrix, and the residuals will be provided in this work.
• The domain is initially discretized by Q9 elements. Once a crack appears, the standard Q9 element will become a global cracking element by "borrowing" the displacement degrees of freedom to indicate local crack openings.
• Since the crack openings become global unknowns, the GCEM is numerically more stable and more efficient than the CEM. Fewer iteration steps than in the CEM are needed. At the same time, the self-propagating property of cracking elements is maintained.
The article is organized as follows: In Section 2, the adopted mixed-mode softening law, the kinematics, and the designed cracking element are described. Numerical examples are presented in Section 3, including an L-shaped panel test and disk tests with initial slots. In most of the examples, irregular meshes, with different sizes of the elements, are used. Section 4 contains concluding remarks. *Nonetheless, a small trick is used insofar as the degrees of freedom of the center node of the Q9 are "borrowed" for indicating crack openings, see the following sections for details.

Traction-separation law
A mixed mode traction-separation law 34,35 is used in this work, with n and t denoting the crack openings in the normal and the parallel direction, respectively. The equivalent crack opening is defined as The traction components in the normal and the parallel direction, that is, T n and T t , are obtained as where f t is the uniaxial tensile strength, G f is the fracture energy, G f ,0 is the threshold value of G f , with G f ,0 = 0.001G f assumed herein, 0 is the threshold opening, given as 0 = 2G f ,0 ∕f t , and mx denotes the maximum opening, which the crack has ever experienced. Its value is updated at the end of every loading step if mx > 0 . T mx = L 2 ( mx ) is the corresponding traction. Figure 1 shows T eq ( eq ) . From Equations (1) and (2), the following relations can be obtained These relations will be used for computing the stiffness matrix and the residual of the cracking element.

Kinematics
The deduction of the displacement field u(x) and of the corresponding strain field (x) of the cracking element can be found in References 1 and 16. The result for u(x) is given as where u(x) is the regular part, H s (x) is the Heaviside function across the discontinuity surface, with H s (x) = 1 on one side of the subdomain and H s (x) = 0 on the other side, and (x) is a differentiable function with (x) ∈ [0, 1]. The result for (x) is obtained as where, unit vectors corresponding to n and t are denoted as n and t, respectively. In the presented symmetric formulation, † ∇ ||n. Then, based on energy conservation, 16 ∇ = n∕l c is obtained, with l c = V∕A. In this formula, V is the volume of the cracking element and A is the area of an equivalent crack surface, that is, the surface of a parallel crack, passing through the center point of the element, as shown in Figure 2. Hereby, l c corresponds to the classical characteristic length. It only depends on the shape of the element and on the orientation of the crack, but not on its position. For the FEM approximation of the eight-node quadrilateral element (Q8), Equation (5) is obtained as where (⋅) S denotes the symmetric part of the tensor 36 and (⋅) (e) is the value of the respective quantity of the Q8 element e. For example, N (e) i is a shape function. Furthermore, in References 1, 2, and 16, an idea referred to as "center representation strategy", is presented. Once a crack appears in an element, the stress/strain state at the center Gauss point is used for representation of the stress/strain state of the whole element. Considering the equivalence of forces of discrete and embedded models, 16 the following relations exist not just at the center point, but in the whole cracking element: where C (e) is the elasticity tensor. This requires a comment. For elastic analysis, nine Gauss points are used for Q9. In one element, the stress/strain states at different Gauss points can be very different, especially when using a coarse mesh. Some Gauss points can experience cracking whereas some others may not. 37 Thus, center representation strategy may be inaccurate. However, presently † In this approach (x) does not have to be known, while in the XFEM and in some other types of SDA (x) = ∑ N i i is assumed. Therefore, the proposed approach is easier to implement, and it has fewer assumptions. no alternative is available. The authors of Reference 16 obtained ∇ = n∕l c and l c = V∕A based on energy conservation, considering the stress/stain state at the center point. In other words, determination of l c (see Figure 2) is unavailable when assuming multiple cracks appearing at different Gauss points in a single element. In appendix A of Reference 16, numerical studies also indicated local stress locking can happen if the center representation strategy is not used.

Global cracking element
In this subsection, the Voigt notation is used for representation of second-and fourth-order tensors with corresponding vector and matrix forms. 38 For the Q8 element, the displacement vector is given as The B matrix is given as where The unit vectors n (e) and t (e) are given as Equation (6), a special matrix B (e) is introduced as follows:  (7) can be rewritten as follows: Next, for convenience we assume that B (e),1 denotes B (e) for the center Gauss point (center representation strategy). Thus, Equation (6) gives The elastic energy E (e) is obtained as where ∫ (⋅)d(e) is the integral of e ⋅ (e) in the element e and C e denotes the matrix form of C (e) . The elastic energy is used for checking whether the equilibrium iteration by means of the Newton-Raphson (N-R) method converges, see the next section for the details.
For the iteration step l within load step i in the framework of the N-R iteration, 1,35 the following relation is obtained: where Δ (⋅) denotes an increment with respect to the corresponding value at the preceding load step, i − 1, and ΔΔ (⋅) stands for an increment with respect to the value at the last N-R iteration step, l − 1. The N-R iteration process was described in Reference 2.
Considering the balance equation ∇ = F for quasistatic loading, where F denotes the loading force, and making use of the equivalence relations according to Equation (10), after linearization the global balance equation is obtained as follows: where, (14) is not useful for numerical analysis, because neither K (e) s nor K (e) is symmetric. As a remedy of this problem, Equation (14) is slightly changed to which shall be used to obtain the element stiffness matrix and the residue of the global cracking elements. However, also K (e) s,new is not symmetric. Moreover, this matrix has many zero elements, rendering the iteration numerically unstable. Fortunately, K (e) s,new can be replaced by the following symmetric matrix: which is the final element stiffness matrix of a global cracking element in the framework of a standard Galerkin weak form. It is worth mentioning that the matrix K (e) new , see Equation (15), appearing as the final term on the right-hand side of Equation (15) remains the same. Otherwise serious stress locking would occur (see appendix A of Reference 16).
Finally, the nodal displacement vector U (e) pseudo-Q9 and the matrix B (e) pseudo-Q9 are defined as follows: They are similar to the displacement vector and the B matrix of the Q9 element. The original displacement of the center node of this element has been replaced by the element-wise crack openings (e) n and (e) t . By introducing ΔΔ (e) as a global unknown quantity, the relationship between ΔΔU (e) and ΔΔ (e) is established with the help of Equation (15). As compared to the CEM, where (e) i−1 + Δ (e) l−1 is obtained from U (e) i−1 + ΔU (e) l−1 by an additional iteration on the element level, the N-R iteration in the framework of the GCEM is more efficient, because fewer iteration steps are needed. This will be verified in the numerical investigation. In summary, the formulation of the cracking element, presented in Reference 1, is converted to matrix form, resulting in a pseudo-Q9 formulation. With this in mind, the presented formulation can be easily changed to an enriched form of the crack opening. 39 In other words, the extra degrees of freedom for indicating crack openings are only introduced if cracks appear.

Initiation and propagation of cracks
The local criterion, presented in Reference 1, is used for determining the orientation n (e) of the crack aŝ , which permits determination of n (e) based on Equation (18). A disadvantage of Equation (18) is that n (e) rotates during loading. This is even the case within one loading step if, that is, a new crack appears. However, since n (e) depends on̂( e) but not (e) , this rotation is generally small when in comparisons to some other classical rotating crack models such as the ones in References 36 and 40 where the cracks may rotate by 90 • during the loading process, as was reported in late 1980s. 41-43 § Similar to the procedure described in Reference 1, the uncracked elements of the domain are separated into two regions: the propagation region and the new root-search region. This is done with a simple criterion, checking whether the element shares at least one boundary with a cracking element, see Figure 3. The following strategy is used to find the next cracking element: where the propagation region is always at first searched for the existence of a propagation region ¶ . If a new cracking element appears in the propagation region, this region will consequently expand, and the search will be continued. If no new cracking element is found in the propagation region, the new root-search region will be checked. If that a new cracking element appears in the new root-search region, the propagation region will further expand and the search in the ‡ Using the Voigt notation,̂( e) = B (e),1 U (e) is obtained. However,̂( e) in Equation (18)  propagation region will again be continued. If still no new cracking element is found in the new root-search region, the next loading step will start. With this strategy, both the initiation and the propagation of cracks can be captured. The algorithm of the described procedure is illustrated in Figure 4. The procedure is simpler than the one used in the framework of the CEM. 1 The domain is firstly discretized by Q9 elements. Once a new crack appears, the standard Q9 element is converted to a pseudo-Q9 cracking element, as described in the previous section. Because of the center representation strategy, the criterion in Equation (21) is only considered for the center point of the elements. In other words, if (e) RK > 0 only appears at noncenter Gauss points, the element is not recognized as a cracking element. This causes oscillations of force-displacement curves when using a very coarse mesh, as reported in References 12 and 35.
As mentioned before, the degrees of freedom of the center node of the Q9 element are "borrowed" for indicating element-wise crack openings. Therefore, at least theoretically, these degrees of freedom can be returned, meaning that the displacements of the center point are returned to the displacement field. This occurs if, for example, a normal Q9 element switches to a cracking element. Then the calculated eq is very small, that is, almost equal to zero. Herein, the criterion is very strict insofar as only in case of (e) eq = 0 and if the cracking element is not neighbor of any other cracking elements, the mentioned degrees of freedom will be returned to the center node, indicating its displacement. In all of our numerical examples, no degrees of freedom were returned. In other words, in all of the numerical examples the transformation from the standard Q9 to the pseudo-Q9 cracking element was never followed by its reversal.

NUMERICAL EXAMPLES
Plane-stress assumptions hold for all the numerical examples presented in this section. Although relatively coarse as well as fine meshes are considered in the examples, generally the meshes are finer than typical meshes used for elastic analysis with high-order elements. This was done in order to avoid oscillations of force-displacement curves reported in References 12 and 16. Some adaptive approaches are still under developing.

L-shaped panel test with regular meshes
The L-shaped panel test, investigated in this article, was studied before in References 1 and 16. The set-up and the material parameters are shown in Figure 5. The displacement increment is chosen as Δd = 10 μm. The regular meshes are firstly considered, as shown in Figure 6. The force-displacement curves are shown in Figure 7, generally indicating good agreement with the results obtained by the XFEM. 34 It can be found mesh III produced slightly higher force responses which is caused by the different cracking patter. Figures 8 and 9 show the deformation and (element-wise) crack-opening plots and the maximum principal stress at Gauss points, respectively, for d = 450 μm. The stress plot shows that for cracking elements with large values of the crack opening, the stress is fully released. However, for very few noncracking (normal Q9) elements, the maximum principal stress of which at the center Gauss point, does not exceed the tensile strength, while the maximum principal stress at some other Gauss points may temporarily exceed the tensile strength. This results in the oscillation behavior of force-displacement response in case of using very coarse mesh, as pointed out in Reference 16. A local smeared crack approach 44 could be helpful in the prospective implementation, which is still ongoing. The deformation and crack-opening plots and the maximum principal stress at Gauss points for d = 600 μm are shown in Figures 10 and 11, indicating a situation similar to the one in Figures 8 and 9.

L-shaped panel test with irregular meshes
Then, irregular meshes automatically built with Gmsh are considered, as shown in Figure 12. Figure 13 show the force-displacement curves, again indicating good agreement with the results obtained by the XFEM 34 and by the CEM. In the calculations, the total elastic energy, given in Equation (12), was used for checking whether the iteration converges. Thus, if  Figure 14. It is seen that the GCEM needs fewer iteration steps than the CEM. However, generally the numbers of iterations are high. This is because (i) the tolerance in Equation (22) is relatively small and (ii) the flowchart in Figure 4 indicates that once a new global cracking element appears, the N-R iteration will be rerun.

Brazilian disk tests
Brazilian disk tests with single or double slots were studied experimentally in Reference 46. Numerical simulations were conducted by the phase field method, 47 the peridynamics theory, 48 and by the lattice element model. 49 Similar to the GCEM, these methods do not need crack tracking. However, it will be shown that in case of the GCEM, good results are obtained with coarser meshes. The models of the tests are shown in Figure 19. The slot width is chosen as 1 mm. The disk test is a typical brittle test, where a snap-back behavior can be expected once the first crack appears. However, arc-length methods that can account for this behavior are not considered in this work. Still, for the displacement increment of Δd = 0.5 μm, no convergence problem occurred. The meshes are shown in Figure 20. They are automatically generated by Gmsh which changes slightly with different analysis of inclination, , of the crack.

Results of intact disk
The model of intact disk is firstly simulated, considering two values of the loading width "w". For intact disks, the relationship between the ultimate loading per unit thickness, F max , and the tensile strength f t is given as which is used for obtaining normalized peak loads. The force-displacement curves and the crack-opening plots are shown in Figure 21. For the two chosen values of "w", the obtained peak loads agree very well with the analytical solution. However, after cracking, the force does not drop to zero. This is the consequence of the center representation strategy. As mentioned before, the element will be recognized as a cracking element only if the center point fulfills the cracking criteria, summarized in Equation (21). In disk tests, some elements close to the loading position can experience local stress concentrations. In that case, the stress states at some noncenter Gauss points fulfill the cracking criteria but the stress state at the center point remains as a relatively low value. Such elements are considered as normal Q9 elements for very long period. This is a drawback of the proposed model. It can be partially mitigated by using a small value of "w" and a fine mesh. In the following tests in this section, "w = 4 mm" is taken.

Results of disk tests with a single slot
For disk tests with a single slot, once the force-displacement curve begins to drop, the disk is considered to be completely damaged and the simulation is terminated. The force-displacement curves are shown for different meshes in Figure 22, indicating similar peak loads and load histories. The normalized peak loads are shown in Figure 23 for the GCEM, the phase field method, 48 and the peridynamics theory. 48 The results provide evidence of the reliability of GCEM.

Results of disk tests with double slots
For disk tests with double slots, the disks are continuously loaded until several major cracks appear. The force-displacement curves and the normalized peak loads are shown in Figures 27 and 28. Generally, they are similar to the results published in Reference 47. Unfortunately, in most cases the forces did not drop to zero at the final step. This is partly caused by the drawback brought by the center representation strategy mentioned in Section 3.3.1. The other reason is that geometric nonlinearity (large displacement effects) has not been considered in this work. The results for the crack openings are shown in Figures 29 to 32, indicating that the GCEM is capable of capturing complex crack growth.

CONCLUSIONS
In this article, for extending earlier work on the CEM, a novel standard Galerkin-based numerical approach, named GCEM was presented. The GCEs used in the GCEM represent a special type of elements, formally similar to Q9 elements. This feature makes it easy to implement the elements into the FEM. Moreover, unlike the conventional XFEM and SDA, the GCEM does not need nodal/element-wise enrichment or a crack-tracking strategy. Unlike the phase field method, a relatively coarse mesh can be used in the GCEM to obtain reliable results with negligible mesh dependency. Unlike methods based on equivalent-type theories, the GCEM does not introduce bonds/lattices/links and accompanied additional assumptions. It stays in the conventional continuum-based framework. The numerical tests have indicated the reliability, efficiency, and robustness of the GCEM. This method is capable of capturing both the initiation and the propagation of cracks and of obtaining reliable results of mixed-mode crack openings. A disadvantage of the GCEM is that presently only quadrilateral elements with nonlinear interpolation of the displacement field are available. Such elements are obviously more complex than elements with linear interpolation of the displacement field.