Received date: January 18, 2017; Accepted date: February 04, 2017; Published date: February 08, 2017
Citation: Martínez Concepción ER, De Farias MM, Evangelista F (2017) XFEM Potential Cracks Investigation using Two Classical Tests. J Appl Mech Eng 6:251. doi: 10.4172/2168-9873.1000251
Copyright: © 2017 Martínez Concepción ER, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Applied Mechanical Engineering
Extended Finite Element Method (XFEM) is used in this work, first to perform the simulation of crack initiation and propagation mechanisms in plane models and then to determine the stress distribution singularities in the closest surroundings of a front fracture inserted in three-dimensional models. The essentials of XFEM is the well-known Finite Element Method (FEM) adding to degrees of freedom and enrichment functions, which serve to describe local discontinuities in the model. In XFEM, the fracture geometry is developed independent of the mesh, allowing it to move freely through the domain, without the need to adapt the mesh to discontinuity. In other words, the XFEM reproduces the discontinuity of the displacement field along the fracture, without discretizing this feature directly in the mesh. XFEM carry out the spatial discretization of two classic models in Fracture Mechanics: the single-edge-notch bending test (SEN (B)); and the disck-shaped compact tension test (CDT). The propagation criterion is based on the proportion of energy released and the stress intensity factors (SIF). The solutions provided by the XFEM numerical model indicated an excellent agreement with the results obtained from the experimental data.
Fracture mechanics; Fracture; Extended finite element method; XFEM; Enrichment; Stress intensity factor
The demand for a fracture analysis method has conducted important contributions since the 1960s. Research was initially distinguished by empirical, analytical, and semi-analytical fundamentals Gross et al.  Rice . These methods can be used in simple geometry problems and under specific boundary and loading conditions. For other more complex problems it is common to appeal to numerical methods due to the need to make several simplifications in the analytical models. Thus, several works were available using the traditional Finite Element Method (FEM) to analyze fracture toughness, but this resulted in a complex mesh which required to be adapted to the fracture surface and also be updated at each time-step to refine the elements size positioned in the surrounding area of the fracture tip. Recently, the Extended Finite Element Method (XFEM) has gained in popularity in its use by the scientific community, it allows developing strategies to analyze fractures without the need for a refinement of the mesh. The XFEM is considered an extension of the conventional FEM and is based on the unit partition concept (i.e., sum of the shape functions must be equal to unity), and were studied by Belytschko and Black  and Moës et al.  the developers of their initial working algorithm.
Since the publication of Belytschko and Black  it is usual to find fracture researches on two-dimensional field using XFEM and highlighted for the explicit description . In addition, the solution of fractures in three-dimensional models is constantly updated with news features for improve the approximations, from the purely explicit, implicit forms or coupling both as studied by Sukumar et al. ; Stolarska et al. ; Fries and Belytschko ; Baydoun and Fries .
In this paper, a numerical analysis is developed using the Extended Finite Element Method (XFEM) to study the initiation and propagation process in fractures. The validation of XFEM approximation by numerical simulations in two dimensions of fracture problems with known solutions and the numerically extraction of stress intensity factors ( K) in three-dimension solid is performed.
The XFEM was introduced first by Belytschko and Black  and Moës et al.  and incorporating enrichment functions and degrees of freedom in addition to the conventional approximation of finite elements in the region where the fracture is placed, to simulate discontinuities and singularities. The type of enrichment functions are described as asymptotic (capture the singularity at the crack tip) and discontinuous (representing the gap between the crack surfaces). The enriched area surrounding the crack tip and over of the fracture are exemplified in Figure 1 (a). The mathematical formulation to approximate the displacement field through an implicit-explicit description was introduce for Baydoun and Fries  as follows:
The first part of equation represents the classical FEM approximation, defined by a continuous shape function Ni( x ) and unknowns at nodal points . The enrichments are taken for the shape functions Ni*. I. is the set of all nodes in the domain. The discontinuity in the field of fracture displacements and the enrichment that captures the special behavior at the crack tip are considered, respectively, by the second and third terms. Two types of enrichment functions are implemented by the XFEM formulation, the Heaviside function, H(x), and the fracture tip asymptotic function B(x), as shown in Figures 1b and 1c.
Through process of investigation in two dimensions fracture problems, Moës et al.  explained exactly how discontinuity functions are added in the finite element approximation.
Denoted to Figure 1a the enhanced nodes belonging elements crossed by the fracture interface (set of nodes Iout), as well the nodes elements located in the crack tip (set of nodes I branch). Along the discontinuities, the nodes are enriched function degree entitled Heaviside; in this case the finite element boundaries are cut by the fracture, while the nodes in elements around the crack tip are enriched with the fracture tip functions named Branch. In Figure 1a encircled nodes are enriched with the Heaviside function and nodes with square symbols are enriched with crack tip functions.
The nodes set within a region around the discontinuity tips present a geometric enrichment. Since the origin of XFEM a geometric criterion was defined for the enrichment zone in order to determine the nodes that will be enriched by singularity functions. In particular, the strategy to include enrichment functions is useful for an efficient approximation of crack singularities and discontinuities as well as interface changes. The addition of discontinuous and asymptotic functions to the elements surrounding the crack tip allows to correctly capturing the singularity in this region Moës et al. . In conditions that the crack tip does not end in boundary elements, functions also describe the discontinuity on the fracture surfaces. Adopting a polar coordinates system, with origin on the crack tip and tangential coordinates on the trajectory propagation, as is shown in Figure 1a. Where is the shortest vector length extended from the crack tip and is the angle measured from rectangular to polar coordinates; the fracture tip enrichment functions for an elastic and isotropic material, were presented by Sukumar et al.  as follows:
The use of the enrichment function involves the addition of four degrees of freedom to all nodes in the improved region. For each freedom degree, is associated a term of the function. The pair (r, θ) represents local polar coordinates at the crack tip for the interval of -π ≤ θ≤ π. The expression is fundamental in the formulation of the function (x)for the reason that it describes the discontinuity on the fracture surfaces. Furthermore, this function is responsible for the representation of phenomenon along crack length, lies the fracture approximation occurs. The other functions of B(x) are used to improve the approximation of the solution in the zone proximately to the fracture tip.
For the discontinuity representation, the XFEM can make use a set of levels Sukumar et al. , endorsed for describe the crack interfaces geometry and without required to coincide with the element boundary, or that the surface of the interface matches the element faces. Before the element failure, the enrichment function degenerates into the conventional finite element; as soon as the element presents the damage, it activates the Level Set Method (LSM), which is based on enrichment functions that will assume the modeling of the discontinuities. Numerical simulation by XFEM and LSM incorporated allows modeling the movement of curves and surfaces in a fixed mesh. The mathematical difficulty to represents the fracture problem is assumed by the LSM. For the direction of fracture propagation, the XFEM arranged functions of set levels to track the fracture surface at each time step by Stolarska et al. .
Among other characteristics, in XFEM is implemented the approach known as phantom nodes by Song et al. ; Dassault Systèmes . This mathematical artifice is based on the internal duplication of each enriched element with the addition of phantom nodes.
The numerical simulation of two validation cases is target to test the XFEM accuracy through the comparison of the load-opening curve, computing the stress intensity factors and reproducing the direction of the fracture in the propagation development. The fracture cases analyzed are: disk-shape compact tension, CDT, (ASTM D7313, 2013), and the three-point notched bending tests, SEN (B), (ASTM E1820 . The ASTM standards describe the CDT and SEN (B) tests as important ways to measure fracture toughness. The experimental application determine the magnitude and singularity of the stress field originated in the volume surrounding the fracture front through the parameters, identified as stress intensity factors as exposed by Irwin . The crack problems were taken from experiments data on: diskshape compact tension available by Wagoner et al. and in concrete beams tested in flexion described in the research of Evangelista et al. . This researches describes the geometry and materials of the specimens, the laboratory test configuration and results (load-opening curves) obtained, and the analogous fracture energy values, statistically ensuring reliable results.
Disk-shape compact disk test
The first validation model is represented in Figure 2. The cylindrical sample extracted from a pavement and subjected to uniaxial tensile following the ASTM D7313 (Wagoner et al. in 2005; ASTM D7313, 2013). The specimen has the following dimensions: = 150, = 25, = 25, = 35 and = 110, and the length of the notch is + = 62.5, while the thickness is = 50. In Table 1, are summarized properties of the material used in the disk.
|Fracture Energy||GF= 328 N/m|
|Tensile Strength||ft' = 3.56|
|Young’s Module||E = 14.2|
|Poisson’s ratio||v = 0.35|
Table 1: Material properties for the compact disk (Wagoner et al., 2005).
Single-edge notched bending test
To obtain a load versus crack opening displacement that is used for determine fracture properties of the material, Evangelista et al.  performed this classic test, represented in Figure 3. The material parameters considered in the beams are exposed in Table 2.
|Fracture Energy||GF= 99 N/ m|
|Tensile Strength||ft' = 5.04 MPa|
|Young’s Module||E = 27 GPa|
|Poisson’s ratio||v= 0.19|
|Concrete Density||pconc = 2400 kg/m3|
Table 2: Material properties for the beam (Evangelista et al., 2013).
Discretization density analysis
The mesh density is analyzed thought in a function of stress intensity factor convergence for a simulation of CDT in the three-dimensional space, while in the SEN(B) model the mesh density convergence is investigated for the two-dimensional propagation of arbitrary crack tests, with base in the approximation of the peak load.
CDT model: For both 2D and 3D models, three mesh configurations were analyzed: coarse, intermediate and fine, as shown in Figure 4 for the plane model. The numerical results for CDT model in 2D and 3D, presented in the research, correspond to the fine discretization.
A convergence study as a function of mesh density, was performed in the 3D fracture models as shown in Figure 5 for the three types of discretization.
mesh analysis was started with approximately 5 ∙ 105 elements and is reached to 75 ∙ 105 elements. Accordingly increasing the mesh density in the simulations, that is, decreased size of the elements, it was observed that the stress intensity factor in the second mode (K11) was converging to a minimum, as shown in Figure 6. In the first simulation a K11 of 0.10MPa was obtained and in the sequential simulations the value decreased to approximately 0.040 MPa , that is, the increase in the number of elements suggests a better approximation in the magnitudes of the stress intensity factors. This statement is a consequence for how the test is conducted, where dominate a fracture propagation in a pure opening mode (Kl). After approximately 17∙105 finite elements the value of the stress intensity factor in Mode II of propagation begins to decrease in smaller increments, according to Figure 6, since this amount of elements was considered appropriate to carry out all the following simulations for the three-dimensional CDT fracture problem.
SEN (B) model: In the single-edge notched bending numerical simulation was considered a structured discretization with three different mesh densities, as can be observed in Figure 7. The type element use is the plane strain and the mesh was refined in the fracture region. The maximum loads obtained, for the same applied displacement, is shown in Table 3. In order to evaluate the mesh density performance, the four cases enumerated in Table 3, are offered the discretization configuration, the evolution of the relative error and the processing time for the different simulations.
|Model||Number of elements||Number of nodes||CPU time (min)||PmáxNodal(kN)||Approximation error to P = 3.56 kN (%)|
Table 3: Discretization analysis and simulation results for the SEN (B) model.
The error associated to the approximation of the peak load of the test was calculated according to the equation:
The Figure 8 shows the magnitude of the XFEM approximation error based on the maximum test load. In the abscissa axis is represented the number of nodes used in the model, variable used to characterize the mesh density, while the axis of ordinates represent the processing time consumed in the simulation, as well as the completed relative error through the XFEM approximation. As can be observed, the increase in the number of nodes implies in extended processing times and a significant reduction of the relative error, which influences are the practically insignificant value of 0.30%. It is reasonable, that XFEM errors in the numerical model are small, since the method implements enrichment functions around the fracture to obtain better approximations in the singularities of reaction forces field to a betterquality mesh model.
Opening displacements and stress intensity factors
The fracture openness results in the bi-dimensional model and the stress intensity factors in the three-dimensional models corresponds to the improved discretization of both tests. Plane models did not have to devote considerable time to shape the geometry and boundary conditions. In the three-dimensional model, the geometry and the variables to assign are in fact more complex. Firstly, contours are established to perform an evaluation of the corresponding contour integral and determine the stress intensity factor, In the simulations with inserted fractures it was possible to estimate for each off evaluation point, that is the stress intensity factor in modes I, II and III. Note that only the evaluation points along the fracture front are computed. The approximation is performed from the stress intensity factor value at each contour evaluation point by the following expression:
where is the number of rings excluded and is the total of contours requested in the contour integral. The variable specifies the number of these elements used in the contour integrals computation. The stress intensity factor is expected as the average value, generally the first contours are diverging values so omitted.
The deviation in the fracture path at the CDT test advises the presence of the pure opening mode. The Figure 9 shows the enriched elements that define the fracture faces.
However, Mode I propagation predominate in the bi-dimensional simulations, was found results which deviated fracture trajectory from the horizontal axis of the model. The effect shown in numerical modeling is related to the method capability to detect changes in the fracture energy. It was experimental that the Von Mises equivalent stresses at the Gauss mesh points take up maximum values around the notch vertex. The behavior response represented by the load-opening curve was compared with the experimental results in Figure 10, was observed a good correspondence between XFEM model and the research laboratory results. The curve slope for the loading period in the numerical simulation matches with the experimental results, at the maximum load obtained the XFEM estimates a slight higher than the results of Wagoner et al. and after the peak, the numerical results are intermediate in relation to the experimental curves.
The 3D model is used to investigate the opening process in a crack inserted to the numerical simulations and performed using the same material parameters of the two-dimensional model. The Figure 11 shows the distribution of stress intensity factors along the fracture length front. For the number of plotted points correspond a number of evaluation points used in the XFEM model, there is a predominance of Mode I fracture propagation.
SEN (B) model
Based on the experiment procedure, it is expected that the fracture simulation will exhibit Mode I dominated propagation, as can be perceived in Figure 12. The total fracture length was verified at the beam mid-span, achieve approximately the half distance between the notch tip and the model top. The symmetrical stresses are due to the geometry and contour conditions in the model, as well as the load position. The stress registered on the projection of the fracture tip at the beam top.
In a trajectory, close to the fracture line the development of the stress state was investigated, the stresses distribution in the beam midspan from the notch to the loaded surface in the interval of h (50 ≤ h ≤ 150) was considered. The final stress at the integration points for each crack interface element was plotted in Figure 13 and compared with the experimental results determined by Evangelista et al.  for the 0.5 pmax pós-peak loading stage. The difference in the 75-85 mm beam height sector is due to cohesive elements considerations in the numerical model used in the reference research. Despite this difference, was perceived a relationship in stress at the load point of the model.
The load-opening curves obtained by the XFEM model are compared with the experimental data in Figure 14. A maximum load achieved by XFEM model represents the reaction force acting at the top center node for a controlled displacement. During the loading path the numerical model present a premature change in the curve slope. This is a consequence numerical computations sensitivity, but in overall, the results are satisfactory. The distribution of the stress intensity factors is determined on a beam 3D model with an inserted vertical fracture, positioned over the notch in the model mid-span and represented by a plane. The 3D model was discretized with an accurate density of elements used in simulation approximately 1∙107 nodes. The scheming of the stress intensity factor computed by XFEM is shown in Figure 15. The theoretical virtual fracture length is the width of the beam and equal to 80 mm. As expected, was a propagation predominance in the pure opening mode, due to the experimental configuration, which agrees with the results obtained in the three-dimensional simulations [14,15].
Fracture problems were successfully solved using XFEM; eliminating the requirement to re-discretize the fracture surfaces, thus controlling computational costs and mesh projection errors associated with the conventional finite element method. Through the proper use of XFEM and refining the mesh in the fracture region were obtained excellent results in numerical approximated solutions to the experimental data. The two-dimensional model had an important role as an initial step to understand the problem simulation, and was used to learning the crack propagation due to the low processing time. In both cases of validation it was demonstrated that the maximum force reached in the two-dimensional numerical models is close to the experimental values. Convergence analysis over the mesh refinement returned results as accurately as possible. However, the fundamental discretization parameters must be carefully selected, to enhance the calculation speed and reduce the cost of monitoring. In threedimensional models with inserted fractures the stress intensity factors revealed the predominance of the pure opening Mode I.