Analysis of crack growth problems using the object-oriented program bemcracker2D

This paper presents an application of the boundary element method to the analysis of crack growth problems in linear elastic fracture mechanics and the correlation of results with experimental data. The methodology consists of computing stress intensity factors (SIFs), the crack growth path and the estimation of fatigue life, via an incremental analysis of the crack extension, considering two independent boundary integral equations, the displacement and traction integral equations. Moreover, a special purpose educational program for simulating two-dimensional crack growth based on the dual boundary element method (DBEM), named BemCracker2D, written in C++ with a MATLAB graphic user interface, has been developed and used to verify the adopted methodology. The numerical results are compared with those of the finite element method (FEM) and correlated with experimental data of fatigue crack-growth tests for twodimensional structural components under simple loading, aiming to demonstrate the accuracy and efficiency of the methodology adopted, as well as to evaluate the robustness of the BemCracker2D code.


INTRODUCTION
he boundary element method (BEM) is a powerful analysis tool for elasticity problems, particularly in the field of fracture mechanics.The incremental crack growth numerical techniques proposed by Blandford et al. [1] and Portela et al [2] are the most popular in the field.The first authors adopted the sub-regions approach where an artificial boundary connects to the crack and divides the structure into two new sub-regions, which is disadvantageous in T THEORY he numerical technique adopted for the simulations is the dual boundary element method (DBEM) [11], which considers two independent boundary integral equations, for the displacements and tractions.Both integral equations use the same integration path for each pair of matching source points but are applied to distinct boundaries along the crack, giving rise to independent algebraic equations.These governing equations are described in more detail below.The displacement boundary integral equation, in the absence of body forces and assuming continuity of displacements at a boundary point x', is given by, where i and j are Cartesian components; ij(x',x) and Uij(x',x) represent Kelvin fundamental solutions for traction and displacement, respectively; CPV denotes Cauchy principal value, and c ij (x') are geometric coefficients given by  ij /2 for a smooth boundary at point x', where ij is the Kronecker delta.
The integrals in eq. ( 1) are regular in case the distance r between the source point x' and the field point x is non-zero.
When r tends to zero, the fundamental solutions present strong singularities of order 1/r for ij , and weak singularities of order ln(1/r) for U ij .
In the absence of body forces and assuming continuity of displacements and tractions at a boundary point x' on a smooth boundary, the stress components  ij are given by,

T
In the above equation, HPV denotes Hadamard principal value.The tensors Sijk(x',x) and Dijk(x',x) contain derivatives of Tij(x',x) and Uij(x',x), respectively.Similarly to eq. ( 1), when r tends to zero, Sijk presents a hipersingularity of order 1/r 2 while D ijk has a strong singularity of order 1/r.Thus, we can write the traction components t j on smooth parts of the boundary in the form where ni represents the ith component of the unit outward normal to the boundary at point x'.Eqns.( 1) and ( 3) form the basis of the DBEM and are applied here by using the standard BEM approach.After discretising the boundary into elements, assuming a quadratic variation of the displacements and tractions within each element, and applying the discrete version of Eqns.( 1) and ( 3) at a number of collocation points equal to the number of boundary nodes, a linear system of algebraic equations is obtained which allows the calculation of the initially unknown displacements and tractions at the boundary nodes, using the system Hu=Gt ( 4) where H and G contain the integrals of T ij and U ij , respectively, given by Eqn.(1), or the integrals of S ijk and D ijk given by Eqn.
(3), respectively.The vectors t and u contain the traction and displacement components at the boundary nodes, respectively.
The system described by Eqn. ( 4) can be assembled as, Ax=By=f where x is a vector containing the unknown values of ti and ui, while y is a vector containing the boundary conditions  and  .The matrices A and B are obtained by rearranging the matrices H and G in the conventional BEM manner.

Determination of Stress Intensity Factors
The elastic stress field at the crack tip can be determined by a correlation between the size, geometry and a Stress Intensity Factor (SIF), which defines their magnitude and also is used in the propagation angle prediction and growth increments.Therefore, a material can resist for crack growth without brittle fracture while the SIF falls below a critical value KIc, which is the material fracture toughness.
The stress intensity factors are generally obtained using a displacement or stress extrapolation technique, which requires sufficient mesh refinement at the crack tip in conventional FEM approaches [12].An alternative is to use methods relying on an energy approximation, therefore avoiding singularities close to the crack tip.The evaluation of SIFs using a pathindependent integral is particularly efficient with the BEM as displacements, displacement derivatives and stresses at internal points are directly obtained from the boundary integral equations.
In this paper, the J-integral technique for obtaining SIFs is employed.According to [2], this is a very efficient postprocessing technique for solving crack problems in general and is based on integrals along a path-independent boundary [13].Here, the decomposition technique for uncoupling the SIFs with the J-integral was implemented considering a circular contour path around each crack tip, as illustrated in Fig. 1, whose reference system has its origin at the tip of a traction-free crack.
According to [13], the relationship between the J-Integral and the SIFs is given by, ' where K I and K II are the SIFs at the crack tip for modes I and II, respectively, E' is the Young's modulus, with E'=E for plane stress and E'=E/(1-v 2 ) for plane strain, where v is the Poisson ratio.
Considering the SIFs uncoupling, the J-integral can be given by the sum of the two integrals, The decomposition of Eqn.(7) leads to, 2 2 , ' ' The integration process takes place along a circular path around the crack tip, as illustrated in Fig. 1.A simple trapezoidal rule can be adopted along the circular path.
Figure 1: Circular contour path for calculation of the J-integral [12].

Prediction of Fatigue Crack Growth
According to [14], there are several criteria to compute the crack-extension direction in the linear elastic regime, as the Maximum Principal Stress (MPS) criterion.This criterion postulates that the crack growth will occur in a perpendicular direction to the maximum principal stress, that is, the crack propagates in the direction that maximizes the circumferential stress in a closed area at the crack tip.Here, this criterion will be used since it describes the local direction for the crack propagation considering mixed mode fracture mechanics, with the SIFs calculated by the J-integral technique.
According to the MPS criterion, the local direction of crack-extension t is determined by the expression, sin (3cos 1) 0 where  t is an angular coordinate with centre at the crack tip and measured from the crack axis ahead of the tip.
In an incremental analysis, the procedure used to set the direction of the nth crack extension increment requires a correction angle β in the tangential direction θ t(n) of the crack path due to the continuity criterion of Eqn.(9), as shown in Fig. 2. Accordingly, where the length increment of the crack propagation, a, tends to zero, the angle t (n + 1) tends to zero as well, and therefore the angle  is fixed meaning that, in the limit, the crack-extension direction tends to the tangential direction of the continuous crack path.This correction angle  is given by  = t(n+1)/2, where t(n+1) corresponds to the direction of the next crack-extension increment, which is calculated by the MPS criterion.
The fatigue crack propagation analysis presented here stems from the results extracted from the incremental analysis of the crack growth, considering the stress intensity factors calculated for each increment.Thus, based on the fundamental postulate of LEFM, a fracture criterion involving the SIFs is given by, where K Ieq is the stress intensity factor (mode I), defined in a mixed-mode analysis, and K Ic is the stiffness fracture plane strain, or a critical SIF value defined as a material property.Thus, we consider a load cycle with constant amplitude, for simplicity, described by a static load level and stress amplitude ratio R, given by, where  min and  max correspond to the minimum and maximum applied stresses, and K Ieqmin and K Ieqmax are the stress intensity factors.
Formulating an analysis at a maximum stress level, we can define an intensity factor range, KIeq, which is given by,   1 I e q I e qm a x I e qm i n I e qm a x Applying the Paris Law [15], we can calculate the crack growth rate as where a is the crack length, N is the number of load cycles, C and m are empirical constants of the material.Eqn. ( 13) expresses the fatigue life and stands for the variation of the number of load cycles required to propagate the crack, or the length advanced at each increment.

MATERIAL AND METHODS
ere, the automation of the modelling process and calculations by the BemCracker2D code, together with BEMLAB2D GUI, as well as the experimental methodology and the FEM models due Miranda [8] compose the materials and methods applied in this work.

Automatic Crack Growth
Automating a crack propagation analysis means describing the whole geometry of the two-dimensional model of the cracked structure and assigning a boundary element discretization for it, physical attributes such as loads, support conditions and material properties.This description consists of a set of geometric continuous and discontinuous quadratic boundary elements that defines the contour of the regions and the cracks, respectively, as well as boundary conditions associated to the elements and domain parameters associated to the regions.

H
The analysis initially set the pre-processor and then continues with the incremental analysis process of the crack growth, whose propagation path is also discretized by short straight segments or elements and, for each increment, the BEM is applied to perform the stress analysis.At each increment, the J-integral is used to compute the SIFs through Eqn.(8), which are then evaluated by Eqn. ( 9) to adjust and define the crack path growth direction and compute the fatigue life as a function of the crack length using Eqn.(13).Finally, the post-processing analysis allows the visualization of the results.The automation process described above has been done through the implementation of two codes: the BEMLAB2D GUI for pre-and post-processing, written in MATLAB, and the BemCracker2D program for boundary element analysis, written in C++ and based on Object Oriented Programming (OOP) [16].The following subsections present a brief review of both codes, whose automation scheme is represented by the flowchart in Fig. 3.

BemLab2D GUI
The BEMLAB2D program is a graphical user interface (GUI) for pre-and post-processing, written in MATLAB and focused on two-dimensional mesh generation and visualization, as well as for visualizing the results of elastostatic analyses produced using the BemCracker2D program.BEMLAB2D GUI is based on user-defined actions through the tool buttons, mouse and dialogs, which modules and features are described below and illustrated in Fig.

BemCracker2D Program
The BemCracker2D program is based on our previous work [17,18], with respect to the standard BEM modelling (displacement equation and continuous quadratic elements) and the incremental analysis strategy due to Aliabadi [11], with respect to the traction equation for crack modelling.BemCracker2D takes advantage of OOP and interprets the attributes, or data generated in the BEMLAB2D code as objects, and properly assigns to each method the form specific classes sorting out the complex data structure that exists in the modelling process.It also assigns characteristics and abstract behaviors.This separation between data and methods grouping is known as encapsulating, a key concept of OOP [19].Tab. 1 shows the main concepts of OOP.In OOP, a class is a description of the set of attributes and methods that define the objects group.In C++ language, classes can protect the access structure for the other data entities through three access attributes: private, public and protected.Each attribute provides an occultation level for class members and is used to ensure that only identified codes groups have access to parts of the pre-selected class.
In this context, the BemCracker2D program is based on the class diagram shown in Fig. 6, where BemCrk_BEMSYS with its attributes and methods is the driving class program and main link with the BEMLAB2D interface.Tab. 2 summarizes the purpose of each class outlined in the class diagram in Fig. 6, and their respective instances/objects, aiming at understanding the incremental analysis which was implemented according to the sequence diagram in Fig. 7.

Experimental Data
Based on the experimental results obtained in Miranda [8], where crack propagation fatigue tests were performed in 2D geometries under simple loadings, this section presents the experimental methodology adopted and the crack propagation models tested.Moreover, a brief description is given of the FEM codes Quebra2D [8,9] and ViDa [10] used to corroborate the numerical analysis with the results obtained by BemCracker2D (BC2D).
The methodology adopted by Miranda [8] was aimed at testing real specimens pre-shaped by finite elements, by means of the Quebra2D program for obtaining the SIF values along the propagation path, and then adjusting them by a suitable analytical function f(a/w), which is an input data for the ViDa program, in order to validate their numerical results.For this, all tests were based on the preparation of fatigue crack propagation tests, as follows:  Crack test preparation: modelling of specimens (CP) by the FEM, making the CP and remodeling of the real CP;  Loading application: simple sinusoidal loads of constant amplitude;  System test: servo-hydraulic machine of high speed;  Measure of the crack length: optical microscopy, image analysis system and the flexibility variation of the method;  Propagation test preparation: standard CP, obtaining the expression of K I (a) to specify the test loading program.
For curved cracks, standard CPs has been modified by the introduction of holes positioned in a way as to bend the crack according to the desired path.For this, the CP was initially modelled with the FEM and a hole was inserted to obtain the crack path.This process was done using the Quebra2D program and once designed the CP was made numerically with the hole in the predetermined position and then remodeled with FEM for any manufacturing defect correction.
With the adjusted numerical model, the incremental crack propagation was simulated, making it grow in small steps under simple loads and according to the criteria mentioned in section 2. The K I values were obtained along the propagation path and adjusted by an analytical function f(a/w), were

Experimental Specimens
Five specimens were tested under quasi constant SIF: 1 (one) SENB (Single Edge Notched Bend) specimen under four points bending load including a hole (Fig. 8) and ( 4   As described previously, the tests were performed with quasi constant SIF, because a constant amplitude cycled force (20 Hz) is applied on specimen, and as soon as the crack propagates, the force is reduced to maintain the same level of SIF.A hydraulic system used to carry out the tests is composed of a servo-controlled machine with load capacity of 250 kN and stiffness of 0.43 GN/m.Values of force, displacement and deformation were read constantly to control the force applied to the specimen.

RESULTS
he following results of crack propagation fatigue analysis under simple loading, obtained with the BC2D program, are compared with those obtained in Miranda [8,9]:  Crack path mesh and KI and KII values;  Graphics of the normalized crack propagation path: f(a/w) x a/w;  Fatigue life diagrams.For all considered models, the BEMLAB2D code was used in the crack model generation and boundary element mesh, and the following modeling strategy being adopted:  The crack boundary is modelled with quadratic and discontinuous boundary elements, with ratio over the crack 0.4-0.3-0.2-0.1 and the element size containing the initial crack tip as equal to 0.25mm;  Quadratic semi-discontinuous boundary elements are used at the intersection between the crack and the edge;  The remainder of the boundary is modeled with continuous quadratic elements.T Initially, a FEM analysis with ABAQUS [20] for K I evaluation in a beam with a centered hole as illustrated in Fig. 10, was carried out aiming to calibrate the BC2D program.Three simulations, varying the crack size and the mesh dimension, were carried out to compare the SIF values when the crack tip is located at a distance of the main boundary hole equal to 2.5, 6.4 and around 12 mm, respectively.It is important to note that this procedure have been used in the other examples, where the final meshes presented were obtained after convergence for a given SIF value.The results obtained with the BC2D code using different numbers of elements are compared with those obtained by ABAQUS and are presented in Tab. 3. Bold values correspond to the final convergent values for the respective meshes.For the crack nearest to the hole (≈12 mm), the maximum size allowed by ABAQUS (minimum radius of the J-Integral) to perform the analysis was 12.1 mm, while with BC2D was 12.2 mm.The results in Tab. 3 show that both FEM and BEM produced similar KI values when the crack tip is distant from the inclusion (2.5 and 6.4 mm), although the discrepancy between the numbers of elements containing the crack tip with FEM and BEM should be noted.On the other hand, as the crack approaches the hole (12.1 and 12.2 mm), the values obtained with the FEM are smaller, which justifies the need for a better FEM mesh refinement in the hole region and at the crack tip.

Results for modified SEN
Fig. 11 shows the initial mesh generated by BC2D with an enlarged detail of the initial crack, where the BEM discretization has 124 nodes and 62 quadratic elements, 24 continuous on the hole boundary and 4 discontinuous on each crack surface.Fig. 14 shows the results obtained with BC2D for the beam with a hole.The crack propagation size is 1 mm with 13 increments and 6-12-24-36 quadratic boundary elements in the hole.In this, it is observed that the insertion of the hole has a significant influence on the f(a/w) value for both methods, with a difference which extends along the path to closer to the hole (last increment for BC2D).The BEM mesh becomes stabilized with between 24 and 36 elements in the hole with f(a/w) = 2.29, while for the FEM this value is 1.91 using the same data a/w = 0.5167.This difference is due to the number of load increments and the initial crack length used in the FEM, whose f(a/w) value in the last increment was 2.13 and, therefore, very close to that found with BC2D.Fig. 16 shows the initial mesh for the CTS01 test generated by the BC2D program, with an enlarged detail of the initial crack.The other CTS models have the same mesh size and will not be shown here for simplification.This mesh has 168 nodes and 84 quadratic elements, 16 continuous in the main hole and 4 discontinuous on each surface crack.Fig. 19 show the results for f(a/w) for the modified CTS (with a hole) obtained by FEM (a = 6.5 mm) and BC2D (a = 9 mm) compared to results from the literature [22] for a CTS with no hole.The results for f(a/w) for the CTS01 and CTS03 models are similar until values of a/w between 0.5 and 0.55, and thereafter exhibited lower values than the CTS without a hole, but with the BC2D the values are larger than FEM and far closer to the standard curve, which was expected, since the propagation paths deviate from the main hole, as seen in Fig. 17.For the CTS02 and CTS04 models, the f(a/w) curves are also similar up to a/w ≈ 0.4 and, from there, with higher values than the curve values corresponding to the case without the hole.Note also that the CTS04 values are higher than CTS02 due to the greater proximity to the main hole.Importantly, due to the initial crack size being different between BC2D and FEM, the CTS02 and CTS04 models have also a different final a/w value, as well as two increments over the FEM.

DISCUSSION AND CONCLUSION
n this work, the dual boundary element method (DBEM) has been applied to two-dimensional crack propagation modelling.The prediction of the crack path, involving an incremental analysis with the BEM incorporating both the displacement and traction equations, has been performing with the aim of assessing the DBEM formulation and validate their numerical results by the BC2D program with those obtained experimentally and using a conventional FE formulation by the Vida program.The automation of the modelling process and the incremental load analysis were performed with the BEMLAB2D and BemCracker2D programs.The second code is invoked by BEMLAB2D GUI to perform a stress analysis with the DBEM and, within each load increment, compute the stress intensity factors using the J-Integral technique, the crack growth direction with an adjustment of the maximum principal stress criterion and lastly, the fatigue life diagrams using Paris law.

I
Results of SIF, crack path and lives have been compared to DBEM.About SIF solutions, there were small difference between FEM and DBEM due to the different numerical methods.Crack path gave similar direction with almost none visual difference.Lives gave also small differences, but they are similar to the same differences with FE solutions.These life differences are due to the difference of SIF solutions obtained by methods.In summary, numerical proposed results have been found to be in good agreement with experimental ones.They also have shown that the DBEM provides high solution accuracy and greatly simplifies the modeling, especially when included in the OOP paradigm and associated with C++, because the objects interact with each other, searching, retrieving and exchanging necessary information with BEMLAB2D GUI, that provides flexibility to the BemCracker2D program.

Figure 3 :
Figure 3: Flow chart of the automatic crack growth.
4a.  GEOMETRY (Module I): This module is independent and has two purposes, the 2D model construction using drawing tools (POINTS, LINES, ARCS and ZONES), and model identification without / with cracks;  MESH (Module II): This module allows the mesh generation, for the model created in module I, for three different numerical methods (BEM, FEM and Meshless); however, for the last two methods, the interface is limited by the generation, visualization and storage of the mesh geometry;  BOUNDARY CONDITIONS (Module III): This module is specific for analysis with the BemCracker2D program.It has three boundary condition types: Displacements, Tractions and Unknown;  ELASTOSTATIC ANALYSIS (Module IV): This module is specific for analysis with BemCracker2D program and it is responsible for the three analysis types: STANDARD BEM (no cracked structures), WITH NO CRACK GROWTH and WITH CRACK GROWTH (cracked structures).For the latter analysis model, the dialog in Fig. 4b is opened requesting information such as number of increments, increment size (integer multiplier element that contains the crack-tip size) and Paris parameters;  GRAPHICAL RESULTS (Module V): This module is specific for analysis with the BemCracker2D program and it is responsible for previewing the following graphical results: DEFORMED MESH, STRESS MESH, STRESS INTENSITY FACTORS, CRACK GROWTH PATH, FATIGUE LIFE and CRACK GROWTH REALTIME.Fig. 5 illustrates the interface functionality BEMLAB2D hierarchy, which includes all steps involved in the analysis procedure.

Figure 6 :
Figure 6: Class Diagram of the BemCracker2D program.

Figure 7 :
Figure 7: Sequence diagram of the BemCracker2D program.
) four CT (Compact Test) specimens, including holes in different positions as shown in Fig. 9.The holes included in specimens force the crack paths move out of a straight direction.The modified SEN model is shown in Fig. 8 where a rectangular beam of dimensions 125x30x10 mm with load application points s=50 mm and r=25 mm from the beam center, and a hole radius R=5.2 mm positioned 9.3 mm from the left of the initiator notch.Fig. 9 present the modified CTS model for different position of the main hole.The diameter of the main hole is 7mm, positioned at horizontal and vertical distances A and B, respectively, from the notch root.The CPs were made from SAE 1020 carbon steel with Young modulus E = 205 GPa, yield strength Sy = 285 MPa, tensile strength SU = 491 MPa and reduction in area RA = 53.7%.The Paris equation da/dN = 8.59x10 -14 .K 4.26 was adjusted with average load R = 0.1.

Figure 8 :
Figure 8: Geometry of the modified CP SEN.

Figure 9 :
Figure 9: Geometric details of the modified CTS.

Figure 11 :
Figure 11: Initial mesh for CP modified SEN with Bemcracker2D program.

Fig. 12
Fig.12shows the crack path predicted by programs BemCracker2D with BEM and Quebra2D with FEM.The initial mesh of the FEM model comprises 1995 elements and 4185 nodes and, in the final propagation step, this has increased to 2585 elements and 5467 nodes, while the BEM increased from 62 to 114 elements and 124 to 228 nodes.

Fig. 13
Fig. 13 shows the results obtained with BC2D compared to analytical ones from the literature [21] for the SEN with no hole.

Figure 13 :
Figure 13: Results for SEN standard (no hole): Stress Intensity Factor K I and f(a/w) expression.

Fig. 15
Fig. 15 shows the BEM fatigue life results compared with the ViDa program and Experimental data.

Figure 15 :
Figure 15: Fatigue life for CP SEN modified.

Fig. 17
Fig.17shows the crack path mesh obtained by the BC2D and Quebra2D programs, comparing the BEM and FEM results for the four modified CTS.The initial FEM mesh consists of about 1300 elements and 2300 nodes and, in the final crack propagation step, it has increased to about 2200 elements and 5500 nodes, while the BEM mesh increased from 84 to 142 elements and 168 to 284 nodes.

Fig. 18
Fig.18shows the results obtained with BC2D compared to the literature [22] for the CTS with no hole.

Figure 18 :
Figure 18: Results for CP CTS standard (no main hole): Stress Intensity Factor KI and f(a/w) expression.
Fig. 20 shows fatigue life diagrams for the CTS models, obtained with the ViDa program and Experimental data, both compared with BC2D program.

Figure 20 :
Figure 20: Fatigue life and residual strength for CTS01 and CTS02.

Table 1 :
objects to encapsulate data, functionality and behavior Member A variable or method (operation, service) within an object EncapsulationThe binding of data and methods into a class of objects Inheritance Means that derived classes inherit the variables and methods from their base classes, while possibly redefining or adding new variables and methods; this creates a hierarchy of ancestor classesPolymorphismThe ability of classes in a hierarchy to share names for methods that behave appropriately to the particular class for which they were designedConstructor/destructor Methods to create/eliminate objects Key Concepts of OOP.

Table 2 :
Main classes, instances and purpose.

Table 3 :
KI evaluation with FE and BE methods.