Cohesive Discontinuities Growth Analysis using a Nonlinear Boundary Element Formulation

Fracture mechanics and crack propagation phenomenon have been widely studied by the scientific community in recent decades. Fracture mechanics contains a group of theories that allows the accurate modelling of the mechanical behaviour of bodies subjected to crack propagation. Therefore, this knowledge domain is robust and accurate to explain and to represent the failure process of structures. In recent years, this subject has been particularly studied with a focus on the development of numerical formulations and mathematical techniques capable of modelling the crack growth processes in complex systems and bodies [1,2].


Introduction
Fracture mechanics and crack propagation phenomenon have been widely studied by the scientific community in recent decades. Fracture mechanics contains a group of theories that allows the accurate modelling of the mechanical behaviour of bodies subjected to crack propagation. Therefore, this knowledge domain is robust and accurate to explain and to represent the failure process of structures. In recent years, this subject has been particularly studied with a focus on the development of numerical formulations and mathematical techniques capable of modelling the crack growth processes in complex systems and bodies [1,2].
In spite of linear elastic fracture mechanics be applied successfully in several types of structures, the results provided by this approach are not consistently when quasi-brittle materials are considered. For such materials, the fracture process zone (FPZ) is large enough to introduce nonlinear mechanical effects that cannot be neglected. To model the crack growth process in quasi-brittle materials, cohesive crack approach, by the fictitious crack model, is an interesting alternative. This approach was initially applied by [3][4][5] and extensively used in recent numerical applications [6][7][8][9]. The crack propagation modelling of quasi-brittle materials, based on cohesive approach, requires the solution of a nonlinear problem, generally relating the crack opening displacements (COD) to tractions on the crack surfaces.
This paper addressees the analysis of crack growth in quasi-brittle materials using the boundary element method (BEM). BEM is efficient in modelling stress concentration, as only the body's boundary is discretized [10]. Therefore, this numerical technique is well adapted to solve fracture mechanics problems. The BEM formulation applied in this study involves the sub-region technique. Considering this technique, the solid is divided into sub-domains, in which compatibility of displacements and equilibrium of forces and enforced along the interfaces of all sub-domains. As a result, the cohesive cracks are assumed to growth along the body's interfaces. One advantage of this approach concerns the possibility to analyse crack propagation process in structures composed by nonhomogeneous materials.
The proposed formulation is nonlinear because the cohesive tractions along the interfaces are determined according to the crack opening displacement values. The nonlinear problem is solved using classical Newton-Raphson approach, in which the corrections into the cohesive crack tractions are performed by applying a non-equilibrated traction vector, keeping all relevant matrices constant. The results of the proposed formulation are compared with experimental and numerical results to validate and prove its robustness and accuracy.

Fracture mechanics aspects: Cohesive crack model
Linear elastic fracture mechanics is an important approach used to solve many problems in structural engineering, particularly when FPZ surrounding cracks are reduced enough to possibly disregard the nonlinear effects. However, for quasi-brittle materials, the damage zone ahead of the crack tip is large enough to produce nonlinear effects that cannot be ignored. The cohesive crack model accounts for these effects. Using this model, the energy dissipation phenomenon is assumed to occur along a fictitious crack positioned in front of the crack tip, thus reducing the dissipation zone by one dimension. In this study, the energy dissipation process is approximated by a simple softening law, which is assumed to govern the residual resistance of the material along the fictitious crack. This law relates the fictitious crack opening displacement, w ∆ , to residual tensile stresses, σ , acting at the crack surfaces. Several cohesive crack laws that relate cohesive tractions to crack opening displacements have previously been proposed in the literature. Three of these laws are often adopted to perform crack growth analysis in quasi-brittle materials. The simplest law is given by a linear function relating the cohesive tractions to a fictitious crack opening displacement smaller than a critical value, if

BEM integral equations
The boundary element method (BEM) has been widely applied in several engineering fields, such as contact problems, fatigue and fracture mechanics, due to its high precision and robustness in modelling strong stress concentration. Considering a two-dimensional homogeneous elastic domain, Ω , with boundary, Γ , the equilibrium equation, written in terms of displacements, is given by: , , where µ is the shear modulus, υ is the Poisson's ratio, i u are components of the displacement field, and i b are body forces. Using Betti's theorem, the singular integral representation, written in terms of displacements can be obtained, with no body forces, as follows:  [11].
Equation (6) is sufficient to construct the system of algebraic equations to analyse two-dimensional elastic domains. In order to analyse solids containing cracks, the use of only this integral equation to assemble the system of algebraic equations will lead to a singular matrix as both crack surfaces are located along the same geometrical path. Although possible using only the singular integral representation, Eq. (6) requires the definition of a finite gap between the two crack surfaces and a very accurate integral scheme to compute the integral along the quasi singular elements.
The hyper-singular integral, written in terms of traction, is obtained from Eq. (6). First, this equation, written for an internal collocation point, is differentiated in order to obtain the integral representation in terms of strains. Then, using the Hooke's law, the stress integral representation is achieved. Finally, the integral representation of stresses for a boundary collocation point is obtained by carrying out the relevant limits. Then, the Cauchy formula is applied in order to obtain the traction representation as follows: where the terms * kij S and * kij D contain the new kernels computed from * lk P and * lk u , respectively, [11].
Equations (6) and (7) are, as usual, transformed into algebraic relations by dividing the boundary and the interfaces among domains into elements along which displacements and tractions are approximated. Besides that, one has to select a convenient number of collocation points to obtain the algebraic representations. The algebraic equations for boundary nodes are calculated using boundary collocation points either at the element ends, therefore coincident with nodes, or along the element when displacement and traction discontinuities are to be enforced.

Algebraic BEM equations for multi-domain analysis
To simulate the mechanical behaviour of solids composed by multi-domains, the sub-region BEM technique has to be applied. This technique is well known in literature and it has been successfully applied in several researches as [12,13]. In the sub-region BEM approach, the body on analysis is divided into a finite amount of homogeneous subregions interconnected by interfaces. Figure 2 illustrates a simple case of this condition in which a body was divided into two sub-regions.
As previously presented, BEM analyses involving singular and hyper-singular integral representations are performed using Eq. (6) and Eq. (7). When multi-domains are considered, these equations have to be applied at each sub-domain individually. Then, the classical BEM system of algebraic equations is obtained for each sub-region i of the entire solid as follows: where n represents the number of sub-regions used into the discretisation of the entire solid.
The global system of algebraic equations presented in Eq. (9) cannot be solved directly just by imposing the boundary conditions of the problem because along the interfaces neither tractions nor displacements values are known. Therefore it is necessary to enforce the compatibility of displacements and equilibrium of forces along all interfaces. These conditions can be written as follows: The compatibilities conditions, Eq. (10), coupled to the boundary conditions have to be imposed on the global system of equations. By performing a convenient change on the columns of matrices H and G, all known variables are placed at the right hand side of this algebraic system whereas unknown variables, x , are placed at its left hand side, [10]. This system can be presented as follows: Once { } f − is the vector of know boundary values, the system can be solved and the unknowns variables determined. In order to consider the nonlinear analysis of cohesive crack propagation along interfaces, a finite amount of elastic problems must be solved iteratively following the Newton-Raphson approach. At each iteration, the boundary and compatibility conditions of the problem may change. In order to deal these changes, a double storage strategy was adopted. The kernels of the global system presented in Eq. (9) are calculated just once and it is stored until the end of analysis. Then, the final system presented in Eq. (11) is reconstructed at each iteration in order to consider all changes on boundary conditions observed.

Nonlinear solution technique
In the present work, the cohesive crack model is coupled to the BEM sub-region approach in order to simulate crack growth phenomenon in quasi-brittle materials. Into the domain of elastic analyses, compatibility and equilibrium conditions are enforced for displacements and tractions, respectively, along all interfaces of the body. Therefore, in such cases, the elastic mechanical response is achieved independently of the interface solicitation intensity. On the other hand, into the developed formulation, the interfaces are assumed to have a tensile strength limit. Therefore, until this limit, the interfaces mechanical responses follow the Hooke's law as well as the entire domain. When a given point of the solid, positioned at any interface, reaches the tensile strength limit, the material starts the mechanical degradation process. At this point, the interface finishes its linear elastic mechanical behaviour and the cohesive crack approach is assumed.
The cohesive interfaces parameters c t σ and f G are required by cohesive laws to model the mechanical behaviour of interfaces. In order to perform the mechanical analysis of a body composed by multidomains containing cohesive interfaces it is required to apply the load into increments. This procedure is necessary due to the nature of the problem which is nonlinear. Therefore, singular and/or hyper-singular BEM formulations, Eq. (6) and Eq. (7), have to be applied to obtain the increments of displacements and tractions along body's boundaries.
During the incremental analysis, the displacements and tractions are cumulated as follows: At the end of each load increment, the cumulated values of normal tractions at the interface are compared to the tensile material strength at all interface coupled nodes. Since the BEM incremental solution is written with respect to the global coordinate system, the cumulated interface tractions solution must be rotated with respect to the interface normal outward vector before this evaluation. Figure 3 illustrates the rotation procedure which requires the determination of the outward vector to the boundary.
By proceeding the solution of the incremental procedure, in a given load step, the normal traction at a given interface point will become superior to the tensile strength limit. As the interface does not support such intensity of traction, cohesive model starts and the nonlinear procedure has to be applied. The nonlinear problem is solved using Newton-Raphson approach, in which prevision and correction steps are required. The prevision step is always performed assuming the structure in elastic conditions. The corrections are performed considering the cohesive crack model, in which the non-equilibrated traction vector is calculated by the difference between the actual traction value and the traction value given by the adopted cohesive criterion. In order to describe this numerical procedure, the boundary value problem described in Figure 4 is presented, which involves a body represented by two domains divided by a cohesive interface. For the mechanical problem presented in Figure 4, the traction distribution presented in Figure 5 is obtained in the load step in which the tensile strength limit is exceeded. The traction distribution considers the accumulated solution and the traction exceeding are observed on the inferior interface. As the interface is not capable to deal such traction distribution, the traction field along the interface must be corrected in order to achieve a new equilibrated configuration. This configuration is obtained during the correction step of Newton-Raphson scheme. The first cohesive correction is achieved by reapplying the nonequilibrated tractions in the structure and by splitting the sub-regions that belong to the degraded part of the interface. Therefore, the crack opening displacements (COD) at the fictitious crack are obtained and the mechanical degradation process of the corresponding interface is preceded. In order to achieve the new equilibrated condition it is still necessary to consider that, at the fictitious crack or degraded interface part, the tensile strength is not the initial c t σ , as the COD is not null.
The updated value is obtained as a function of the cohesive law adopted. Then, based on the updated tensile resistance value, another nonequilibrated traction must be reapplied on the structure improving the COD at the fictitious crack and, consequently, improving the mechanical degradation. The iterative cohesive process continues until the non-equilibrated tractions become small enough to be considered as negligible according to a stop criterion. It is important to mention that the reapplication of exceeding tractions on the correction step considers that mechanical degraded interface nodes are disconnected and, then, Neumann conditions are prescribed. With respect to the interface strength degradation, the three cohesive laws previously presented are adopted to govern such decrease of resistance during the nonlinear process. After finishing the cohesive corrections, a new load increment is preceded. Again, an elastic previsions is performed and cohesive corrections are required. This procedure is performed until the determination of the structural failure. When the material residual strength of an interface portion becomes null fictitious crack becomes real crack. Then, considering a given load increment in which a portion of the interface is totally degraded, the interface can be divided into three parts: Ci, CZi and Cri as illustrated in Figure 6. The three interface regions are namely as: Ci, connected interface (Non-degraded material); CZi: cohesive zone interface (Fictitious crack) and Cri: Cracked interface (real crack). In addition to the previous discussion, one geometric interpretation of the cohesive crack model is presented in Figure 7 considering linear cohesive law.
According to Figure 7, when the elastic traction is higher than the interface tensile strength limit  iterative procedure is performed until the traction exceeding along all interface nodes become small enough to be neglected according to a stop criterion. In order to verify the convergence of the nonlinear procedure a stop criterion must be defined. The stop criterion adopted in this work evaluates the sum of the exceeding square tractions in the amount of d n degraded interface nodes, at a given cohesive iteration, as shown the following equation:

Applications
In this topic three applications involving cohesive crack growth along interfaces cracks are presented. The first application deals the mechanical analysis of a tensile structure composed by nonhomogeneous materials. The second and third applications concern the mechanical analysis of concrete bended beams where the numerical responses achieved by the implemented BEM formulation are compared with numerical and experimental results available in literature.
For clarity purposes, the term "Singular-BEM" indicates the curves achieved by BEM singular integral equation. On the other hand, the term "Hyper-singular-BEM" is used to indicate the results achieved by BEM hyper-singular integral equation.

Nonhomogeneous tensile structure
The first application deals the mechanical analysis of the nonhomogeneous structure presented in Figure 8. This structure is composed by two materials which have different mechanical properties being connected along one interface. Displacements are prescribed along vertical boundaries, as presented in Figure 8. The cohesive tensile strength of the interface is equal to The nonlinear analyses were performed using 15 load steps and the adopted tolerance for convergence was 4 10 KPa − . In order to solve this problem, each sub-region was discretized into four linear discontinuous boundary elements. The exponential cohesive law was adopted to represent the interface mechanical degradation and both singular and hyper-singular BEM integral equations were applied.
Based on the results presented in Figure 9, the higher stiffness of left region is observed as its displacements are smaller than those verified for the right region. At the end of the analysis, rigid body movement is observed on the right domain whereas the left domain tends to has null strain state, i.e., it returns to its initial configuration before loading.
The cohesive interface tractions behaviour is presented in Figure  10, where its dependency with COD is illustrated. Excellent agreement is observed among the numerical BEM analyses considered and the analytical response of the studied problem. Therefore, the numerical formulation implemented was validated with an analytical approach. Figure 10 shows the undeformed and deformed meshes for steps 5 th , 10 th and 15 th . For clarity purposes, the displacements on this figure were magnified 100 times.

Three point bended beam
The second application of this paper concerns the mechanical analysis of the structure presented in Figure 11, which is a three point bended beam with an initial notch at its middle span. This structure was analysed experimentally by Saleh et al., [14] and numerically by Oliveira et al., [6,7]. The geometry, boundary conditions and material Cohesive interface 2m 2m  properties of the problem are presented in Figure 11. The crack starts growing at the notch tip and it propagates in pure mode I. Therefore, one vertical cohesive interface was considered which starts at the notch tip and finishes at the beam top.
In order to verify the convergence regarding mesh refinement, three different meshes, namely 1, 2 and 3, were considered. The less refined mesh, mesh 1, considers the interface discretized into 15 quadratic boundary elements at each side. At the boundaries with prescribed displacements 20 quadratic boundary elements were used whereas at the complementary boundary 16 linear boundary elements were applied. All elements of this mesh are discontinuous which results a total of 195 nodes. The mesh 2 represents the interface using 15 cubic boundary elements on each side. The boundaries with prescribed displacements were modelled with 24 cubic boundary elements and the remaining boundary with 26 quadratic boundary elements. On this mesh, only discontinuous elements were adopted which results on 307 nodes. Finally, the most refined mesh, mesh 3, contains 40 boundary elements positioned at each side of the interface. The boundaries with prescribed displacements were represented with 24 cubic boundary elements and the complementary boundaries with 26 quadratic boundary elements. In this mesh all boundary elements are discontinuous resulting into 507 nodes. These three meshes are illustrated on Figure 12. In spite of three meshes studied, it was verified that mesh 2 is accurate enough to achieve the convergence results. Therefore, the results presented in this application were obtained using mesh 2. The prescribed displacement at the beam top was applied into 140 load steps and the tolerance for convergence adopted was equal to 5 10 KPa − . The mechanical behaviour of the interface was modelled using the three cohesive laws previously presented. Both singular and hyper-singular integral BEM equations were adopted.
The responses in terms of load versus vertical displacement (at the crack mouth) are presented in Figure 13. The numerical results obtained by the implemented BEM formulations are compared with experimental and numerical responses available in literature. Based on the results presented in Figure 13 one observes that the implemented BEM formulations were capable to reproduce the nonlinear mechanical behaviour introduced by the crack propagation. Both singular and hyper-singular approaches were efficient in modelling the structural behaviour. Moreover, it is important to mention that the structural resistant load was accurately predicted by BEM models. Figure 14 illustrates the deformed mesh obtained by Singular BEM approach and linear cohesive law. For clarity purposes, the displacements in this figure were magnified 100 times. The cohesive tractions and its dependency with COD are presented in Figure 15. The behaviour of tractions for singular and hyper-singular BEM models are compared assuming three cohesive laws adopted. According to this figure, good agreement is observed for both numerical BEM approaches applied. Therefore, the implemented BEM models were validated considering an experimental approach.

Mixed mode crack growth in a bended beam
The last application of this work deals a bended beam containing an initial notch subjected to mixed mode crack growth. This structure was experimentally analysed by [15] and numerically studied by several authors, among them [6,7]. The structural geometry, boundary conditions and material properties are presented in Figure 16. In order to simulate the cohesive crack growth using the implemented BEM models, one interface must be positioned along the crack path observed in the experimental analysis. This path is illustrated in red colour in Figure 16. Two different sub-region compositions were adopted for the mechanical analysis of the structure presented in Figure 16. The first composition divides the entire domain into 2 sub-regions, in which the interface is positioned along the crack growth path. This composition is discretized into 52 cubic boundary elements at the regions in which displacements are prescribed, 10 cubic boundary elements for each interface side and 20 quadratic elements for the complementary structural boundary. The second composition divides the domain into 13 sub-regions which leads to 17 interfaces. The interfaces were discretized with 194 cubic boundary elements; the regions where displacements are prescribed were discretized with 16 cubic boundary elements and the remaining boundary with 18 quadratic elements. In both composition cases, only discontinuous boundary elements were applied. Figure 17 presents both compositions considered in this analysis and its meshes, where meshes 1 and 2 contain 361 and 907 nodes respectively. The prescribed vertical displacement at the beam top was applied into 100 load steps and the considered tolerance for convergence was equal to`4 10 KPa − . The analyses were performed considering three cohesive laws previously presented and singular and hyper-singular BEM formulations. Figure 18 presents the deformed configuration for both compositions adopted considering linear cohesive law. For clarity purposes, the displacements were magnified 100 times. As presented in Figure 18, a new crack appears at the beam bottom. Therefore, this result indicates that micro cracks appeared in this region during the experimental analysis. The load versus displacement curves, for the point of load application, considering mesh 1 is presented in Figure 19. The results for the same variables considering mesh 2 are illustrated in Figure 20.
According to the curves presented in Figures 19 and 20, good agreement is observed among the numerical responses achieved by sub-region BEM approaches implemented and by references [6,7]. The responses achieved by the implemented sub-region BEM approaches are more rigid in comparison with experimental results. It may be explained due to the snap-back phenomenon observed in this

Conclusions
This work addressed the mechanical analysis of crack growth   in quasi-brittle materials using BEM. A short review on nonlinear fracture mechanics is presented as well as cohesive laws. In addition to that, the fictitious crack model, presented in fracture mechanics theory, is coupled to algebraic BEM equations. The bodies analysed are represented using the sub-region technique. The cracks are assumed to growth along the interfaces of the analysed domains. Therefore, the implemented formulation allows the analysis of crack growth in nonhomogeneous media.
Three applications were considered in this work and the numerical responses achieved by BEM are compared against analytical, experimental and numerical solutions available in literature. According to the presented applications, the implemented BEM formulations were capable to represent the interface nonlinear behaviour introduced by the cohesive crack growth. For the second and third applications, the implemented formulations have shown its potential in simulating real crack propagation in concrete beams. However, the snap-back instability was not represented on the numerical results of the last application, once the Newton-Raphson algorithm applied was not adapted to solve such type of numerical instability. Both BEM integral representations provided accurate results for the three cohesive laws considered. Therefore, the crack propagation in homogeneous and nonhomogeneous quasi-brittle media can be consistently analysed by the implemented numerical models.