Microsoft Word - numero_45_art_6 G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 67 Analysis of crack growth problems using the object-oriented program bemcracker2D Gilberto Gomes University of Brasilia, Brasil ggomes@unb.br, http://orcid.org/0000-0002-8385-9042 Antonio C. O. Miranda University of Brasilia, Brasil acmiranda@unb.br, http://orcid.org/0000-0002-5225-7428 ABSTRACT. 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 two- dimensional 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. KEYWORDS. Dual Boundary Element Method; Crack Growth; Stress Intensity Factors; Fatigue Life; C++ code. Citation: Gomes, G., Miranda, C. O. A., Analysis of crack growth problems using the object-oriented program bemcracker2D, Frattura ed Integrità Strutturale, 45 (2018) 67- 85. Received: 07.03.2018 Accepted: 18.04.2018 Published: 01.07.2018 Copyright: © 2018 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 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 http://www.gruppofrattura.it/VA/45/8.mp4 G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 68 an incremental analysis due to the repeated introduction of an artificial boundary leading to unwanted approximations along the interface between the sub-regions and the formation of larger systems of algebraic equations. The second authors proposed a technique based on the application of two boundary integral equations, for displacements and tractions, each applied to one crack surface. The latter technique, named DBEM, has the advantage of not requiring remeshing in the neighbourhood of the crack-tip, as required by the FEM and the sub-regions BEM. Other pioneering works in the field are those of Ingraffea et al. [3], Becker [4], Myiazaki et al. [5], Chen and Farris [6] and Bush [7], all of which are characterised by the application of the displacement boundary integral equation and the sub-regions approach. This paper deals with two-dimensional crack growth with the DBEM, using an in-house developed computer code called BemCracker2D for incremental analysis as well as for the calculation of stress intensity factors (SIFs) using the J-integral, the crack growth path and fatigue life. Two structural components under simple loading have been analysed, derived from experimental tests of fatigue crack propagation due to Miranda [8, 9], including their numerical results with the FEM codes Quebra2D [8, 9] and ViDa [10], for simulation of two-dimensional crack growth and the prediction of fatigue life, respectively, to support the BEM analysis. The functionality of these programs, as well as their advantages and disadvantages, will not be discussed here. Moreover, the whole treatment of the crack modelling and the visualisation of the crack path, the SIF diagrams, the f(a/w) curves and the fatigue life diagrams were conducted with the code BEMLAB2D GUI, a pre- and post-processing program that interacts with BemCracker2D. One main contribution of this work concerns the development of computer software for educational use that allows a complete numerical analysis of elastostatic engineering problems, ranging from CAD modelling (BEMLAB2D GUI) to the visualisation of crack propagation results including the calculation of fatigue life through the BemCracker2D code, validated by experimental results and their FEM correlations This paper is organised as follows: initially, a brief review of the DBEM formulation for crack problems is presented, as well as the techniques used to obtain the SIFs, the crack-extension direction and fatigue life. This is followed by a discussion on the automation of the modelling process and calculations by the BemCracker2D code, together with visualisations by using the code BEMLAB2D GUI. The paper then presents the experimental methodology and the FEM models of Miranda [8], followed by numerical results with BemCracker2D and correlations between the DBEM and both experimental and FEM results. 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,    ' ' ( ', ) (x)d (x) ( ', ) (x)d (x)ij j ij j ij jc x u x CPV T x x u U x x t        (1) 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 cij(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 Uij. 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,              1 ' ', ', 2 ij ijk k ijk kx HPV S x x u x d x CPV D x x t x d x        (2) T G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 69 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/r2 while Dijk has a strong singularity of order 1/r. Thus, we can write the traction components tj on smooth parts of the boundary in the form                  1 ' ' ', ' ', 2 j i ijk k i ijk kt x n x HPV S x x u x d x n x CPV D x x t x d x        (3) 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 Tij and Uij, respectively, given by Eqn. (1), or the integrals of Sijk and Dijk 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 (5) 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 path- independent 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 post- processing 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, 2 2 ' I IIK K J E   (6) where KI and KII 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-v2) 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, I IIJ J J  (7) G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 70 The decomposition of Eqn. (7) leads to, 2 2 , ' ' I III IIK K J J E E   (8) 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) 0I t II tK K    (9) 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, Ieq IcK K (10) where KIeq is the stress intensity factor (mode I), defined in a mixed-mode analysis, and KIc 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, Ieq minmin max Ieq max K R K     (11) G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 71 where min and max correspond to the minimum and maximum applied stresses, and KIeqmin and KIeqmax are the stress intensity factors. Figure 2: Schematic of incremental crack-extension direction [12]. Formulating an analysis at a maximum stress level, we can define an intensity factor range, KIeq, which is given by,  1Ieq Ieq max Ieq min Ieq maxK K K K R     (12) Applying the Paris Law [15], we can calculate the crack growth rate as  m Ieq da C K dN   (13) 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 G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 72 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. Figure 3: Flow chart of the automatic crack growth. 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. 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. G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 73 Figure 4: (a) BEMLAB2D GUI lay-out; (b) Dialog for crack growth analysis. Figure 5: Functionality hierarchy interface BEMLAB2D. 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. G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 74 Terms Definition Object An instance of a class Class A template for objects to encapsulate data, functionality and behavior Member A variable or method (operation, service) within an object Encapsulation The 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 classes Polymorphism The ability of classes in a hierarchy to share names for methods that behave appropriately to the particular class for which they were designed Constructor/destructor Methods to create/eliminate objects Table 1: Key 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. Figure 6: Class Diagram of the BemCracker2D program. G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 75 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. Class Instances/object Purpose BemCrk_BEMSYS bem Reading, Analysis, Solution and Printing the crack model BemCrk_GENRL grl Assembling the system Ax = By = f BemCrk_ELEMENT elem Assembling the continuous (C) and discontinuous (D) quadratic elements BemCrk_MESH mesh Element mesh assembly (C/D) BemCrk_QUADJAC quad Compute Gauss points, Jacobian and shape functions (C/D) BemCrk_CRACKS crk Assembly of crack growth increment BemCrk_JINTGR jint Compute Stress Intensity Factors BemCrk_SYS sys Arrays allocation and LU Solver Table 2: Main classes, instances and purpose. Figure 7: Sequence diagram of the BemCracker2D program. G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 76 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 KI (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 KI values were obtained along the propagation path and adjusted by an analytical function f(a/w), were  I a K a a f W         (14) 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) 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.K4.26 was adjusted with average load R = 0.1. Figure 8: Geometry of the modified CP SEN. G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 77 Figure 9: Geometric details of the modified CTS. 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. Figure 10: Beam test for SIF KI evaluation. T G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 78 Initially, a FEM analysis with ABAQUS [20] for KI 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. Crack-tip distance (mm) ABAQUS BemCracker2D Elements KI Elements KI Crack-tip Total Crack Total 2.5 66 90 156 981 1075 1218 1865 1885 1964 2.12 1.93 1.83 1.80 1.61 1.47 4/face 54 62 74 1.46 6.4 156 1049 1389 1489 2.76 2.73 2.62 4/face 54 62 74 2.63 12.1/12.2 156 932 1058 1091 10.16 9.94 9.628 4/face 62 74 86 98 9.65 9.71 9.75 9.77 Table 3: KI evaluation with FE and BE methods. 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. Figure 11: Initial mesh for CP modified SEN with Bemcracker2D program. G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 79 Fig. 12 shows 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. Figure 12: Mesh crack growth with BemCracker2D (BEM) and Quebra2D (FEM). Fig. 13 shows the results obtained with BC2D compared to analytical ones from the literature [21] for the SEN with no hole. Figure 13: Results for SEN standard (no hole): Stress Intensity Factor KI and f(a/w) expression. G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 80 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. Figure 14: Results for CP modified SEN (holed) and standard (no hole). Fig. 15 shows the BEM fatigue life results compared with the ViDa program and Experimental data. Figure 15: Fatigue life for CP SEN modified. G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 81 Results for modified CTS 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. Figure 16: Initial mesh for modified CTS (CTS01) by Bemcracker2D program. Fig. 17 shows 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. Figure 17: Crack growth: BemCracker2D (BEM) and Quebra2D (FEM). G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 82 Fig. 18 shows the results obtained with BC2D compared to the literature [22] for the CTS with no hole. Figure 18: Results for CP CTS standard (no main hole): Stress Intensity Factor KI and f(a/w) expression. Figure 19: Results for modified (holed) and standard (no hole) CP CTS. 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 G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 83 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. 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: Fatigue life and residual strength for CTS01 and CTS02. 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 G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 84 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. ACKNOWLEDGMENT he authors are grateful to the Brazilian National Research Council (CNPq) and to the Brazilian Coordination for the Improvement of Higher Education (CAPES) for the supporting funds for this research. The authors also thank the Graduate Programme in Structural Engineering and Civil Construction in the Department of Civil and Environmental Engineering at the University of Brasilia. REFERENCES [1] Blandford, G. E., Ingraffea, A. R. and Liggett, J. A. (1981). Two-dimensional stress intensity factor computations using the boundary element method, Int. J. Numer. Methods Eng., 17, pp. 387-404. [2] Portela, A., Aliabadi, M.H., and Rooke, D.P. (1992). The dual boundary element method: Effective implementation for crack problems, Int. J. Numer. Methods Eng., 33, pp. 1269-1287. [3] Ingraffea, A. R., Blandford, G. E. and Ligget, J. A. (1983). Automatic Modelling of Mixed-Mode Fatigue and Quasi- Static Crack Propagation Using the Boundary Element Method, Proc. of Fracture Mechanics: Fourteenth Symposium, ASTM STP 791, J. C. Lewis and G. Sines (Editors), ASTM, I, pp. 407-426. [4] Becker, A. A. (1986). The Boundary Integral Equation Method in Axisymmetric Stress Analysis Problems, Springer- Verlag, Berlin. [5] Myiazaki, N., Ikeda, T., and Munakata, T., (1989). Analysis of stress intensity factor using the energy method combined with the boundary element method, Comput. Struct., 33, pp. 867-871. [6] Chen, S.Y and Farris, T.N. (1994). Boundary element crack closure calculation of axisymmetric stress intensity factors, Comput. Struct., 50, pp. 491-497. [7] Bush, M.B. (1999). Simulation of contact-induced fracture, Eng. Anal. Boundary Elem., 23, pp. 59-66. [8] Miranda, A.C.O., Meggiolaro, M.A., Castro, J.T.P., Martha, L.F. and Bittencourt, T.N. (2002). Applications of Automation Technology in Fatigue and Fracture Testing and Analysis, in: M.P. Braun AA, Lohr RD (Eds.), Fatigue crack propagation under complex loading in arbitrary 2D geometries, ASTM STP, 1411, pp. 120-146. [9] Miranda, A.C.O., Meggiolaro, M.A., Castro, J.T.P., Martha, L.F. and Bittencourt, T.N. (2003). Fatigue life and crack path predictions in generic 2D structural components, Eng. Fract. Mech., 70, pp. 1259-1279. [10] Meggiolaro, M.A., Castro, J.T.P., (1998). ViDa 98 – a visual damagemeter to automate fatigue design under complex loading, Revista Brasileira de Ciencias Mecânicas, 20, pp. 666-685. [11] Aliabadi, M.H. (2002). The Boundary Element Method – Applications in Solids and Structures, Wiley, Chichester, 2. [12] Chan, S.K., Tuba, I.S., Wilson, W.K. (1970). On the finite element method in linear fracture mechanics, Eng. Fract. Mech., 2, pp. 1-17. [13] Rice, J.R. and Tracey, D. M. (1973). Computational Fracture Mechanics, in Numerical and Computer Methods in Structural Mechanics, Academic Press, New York. [14] Castro, J.T.P. and Meggiolaro, M.A. (2016). Fatigue Design Techniques: Vol. III - Crack Propagation, CreateSpace Independent Publishing Platform. [15] Paris, P.C. (1962). The Growth of Fatigue Cracks Due to Variations in Load. PhD Thesis Lehigh University, USA. [16] Booch, G. (1994). Object-Oriented Analysis and Design with Applications, The Benjamin/Cumming Publishing Company. [17] Gomes, G., and Noronha, M.A.M. (2012). A B-Rep data structure and object GUI programming to implement 2D boundary elements. International Journal of Applied Mathematics and Computation, 4, pp. 369-381. [18] Gomes, G., Delgado Neto, A.M. and Wrobel, L.C., (2016). Modeling and View of 2D Cracks Using Dual Boundary T G. Gomes et alii, Frattura ed Integrità Strutturale, 45 (2018) 67-85; DOI: 10.3221/IGF-ESIS.45.06 85 Integral Equation. Proceedings of the XXXVII Iberian Latin-American Congress on Computational Methods in Engineering, ABMEC/RIPE, Brasília, DF, Brazil, 2(6), pp. 120-133. [19] Stroustrup, B. and Ellis, M.A., (1990). The Annotated C++, Reference Manual, AT&T Bell Laboratories, Murray Hill, New Jersey, USA. [20] ABAQUS Software Documentation, Version 6.9 by SIMULIA, (2009). [21] Anderson, T.L. (1995). Fracture Mechanics: Fundamentals and Applications. 2nd ed, CRC Press, Boca Raton, FL, USA. NOMENCLATURE i, j Cartesian components ij(x’,x) Kelvin fundamental solutions for traction Uij(x’,x) Kelvin fundamental solutions for displacement CPV Cauchy Principal Value HPV Hadamard Principal Value ij the Kronecker delta Sijk(x’,x), Dijk(x’,x) Stress and Displacements Tensors SIF Stress Intensity Factor KI, KII SIFs at the crack tip for modes I and II E Young’s modulus v Poisson ratio t Angular coordinate with centre at the crack tip  Correction angle IeqK Equivalent SIF (mode I) IcK Stiffness fracture IeqK    Equivalent intensity factor range da/dN Crack growth rate << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles true /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Error /CompatibilityLevel 1.4 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.0000 /ColorConversionStrategy /CMYK /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams false /MaxSubsetPct 100 /Optimize true /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments true /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile () /AlwaysEmbed [ true ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 300 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.50000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 300 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.50000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.50000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile () /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName () /PDFXTrapped /False /CreateJDFFile false /Description << /ARA /BGR /CHS /CHT /CZE /DAN /DEU /ESP /ETI /FRA /GRE /HEB /HRV (Za stvaranje Adobe PDF dokumenata najpogodnijih za visokokvalitetni ispis prije tiskanja koristite ove postavke. Stvoreni PDF dokumenti mogu se otvoriti Acrobat i Adobe Reader 5.0 i kasnijim verzijama.) /HUN /ITA /JPN /KOR /LTH /LVI /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken die zijn geoptimaliseerd voor prepress-afdrukken van hoge kwaliteit. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.) /NOR /POL /PTB /RUM /RUS /SKY /SLV /SUO /SVE /TUR /UKR /ENU (Use these settings to create Adobe PDF documents best suited for high-quality prepress printing. Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.) >> /Namespace [ (Adobe) (Common) (1.0) ] /OtherNamespaces [ << /AsReaderSpreads false /CropImagesToFrames true /ErrorControl /WarnAndContinue /FlattenerIgnoreSpreadOverrides false /IncludeGuidesGrids false /IncludeNonPrinting false /IncludeSlug false /Namespace [ (Adobe) (InDesign) (4.0) ] /OmitPlacedBitmaps false /OmitPlacedEPS false /OmitPlacedPDF false /SimulateOverprint /Legacy >> << /AddBleedMarks false /AddColorBars false /AddCropMarks false /AddPageInfo false /AddRegMarks false /ConvertColors /ConvertToCMYK /DestinationProfileName () /DestinationProfileSelector /DocumentCMYK /Downsample16BitImages true /FlattenerPreset << /PresetSelector /MediumResolution >> /FormElements false /GenerateStructure false /IncludeBookmarks false /IncludeHyperlinks false /IncludeInteractive false /IncludeLayers false /IncludeProfiles false /MultimediaHandling /UseObjectSettings /Namespace [ (Adobe) (CreativeSuite) (2.0) ] /PDFXOutputIntentProfileSelector /DocumentCMYK /PreserveEditing true /UntaggedCMYKHandling /LeaveUntagged /UntaggedRGBHandling /UseDocumentProfile /UseDocumentBleed false >> ] >> setdistillerparams << /HWResolution [2400 2400] /PageSize [612.000 792.000] >> setpagedevice