Boundary Element Method and Tangent Operator Technique applied in the Multi Crack Propagation Modelling of Concrete Structures

Among the usual engineering materials, concrete is widely utilized around world due to its low production cost, geometry adaptability, chemical properties and strength capacities. Moreover, the collapse of engineering structures is caused by the growth of internal discontinuities, such as cracks, which leads to the mechanical failure by fracture. Therefore, the scientific community has applied the fracture mechanics to predict accurately the collapse conditions of concrete. In such material, the mechanical degradation processes, especially due to the crack growth, introduce nonlinear mechanical effects that cannot be disregarded. Consequently, nonlinear models were proposed to represent the concrete fracture, in which the cohesive crack model must be highlighted [1]. This approach was introduced by Barenblatt and Dugdale [2,3] and was extensively used in numerical applications of fracture problems [4-8]. The crack propagation modelling in concrete by the cohesive crack approach requires the solution of a nonlinear system of equations, generally relating the Crack Opening Displacements (COD) to the tractions values along the crack surfaces. The traction values represent the residual material strength at the cohesive zone, which depend on the COD intensity. The residual material strength decrease as the COD values grow.


Introduction
Among the usual engineering materials, concrete is widely utilized around world due to its low production cost, geometry adaptability, chemical properties and strength capacities. Moreover, the collapse of engineering structures is caused by the growth of internal discontinuities, such as cracks, which leads to the mechanical failure by fracture. Therefore, the scientific community has applied the fracture mechanics to predict accurately the collapse conditions of concrete. In such material, the mechanical degradation processes, especially due to the crack growth, introduce nonlinear mechanical effects that cannot be disregarded. Consequently, nonlinear models were proposed to represent the concrete fracture, in which the cohesive crack model must be highlighted [1]. This approach was introduced by Barenblatt and Dugdale [2,3] and was extensively used in numerical applications of fracture problems [4][5][6][7][8]. The crack propagation modelling in concrete by the cohesive crack approach requires the solution of a nonlinear system of equations, generally relating the Crack Opening Displacements (COD) to the tractions values along the crack surfaces. The traction values represent the residual material strength at the cohesive zone, which depend on the COD intensity. The residual material strength decrease as the COD values grow.
In this context, this study addresses the mechanical analysis of multi crack growth in concrete structures using a nonlinear Boundary Element Method (BEM) formulation and an efficient algorithm to solve the nonlinear system of equations based on the Tangent Operator (TO) scheme. The BEM is a robust numerical technique for modelling fracture and crack growth problems as only the boundary surfaces discretization is required [9]. Therefore, this numerical technique is well adapted to solve fracture mechanics problems, as the remeshing procedure becomes a less complex task. Moreover, the stresses singularities present at the crack tip are represented accurately due to the non-requirement of a domain mesh [10]. The numerical formulation applied in this study involves the sub-region BEM technique and the cohesive crack model. The first technique divides the entire domain into sub-regions over which, for elastic problems, compatibility of displacements and equilibrium of forces are enforced along its interfaces. In the cohesive crack problems, the displacement compatibility is no longer verified. Therefore, cohesive laws are required to relate the COD to the tractions values along the crack surfaces. Moreover, cracks are assumed to growth along the sub-region interfaces. One advantage of this approach concerns the possibility to analyse nonlinear fracture problems in nonhomogeneous systems, which is a challenge problem in structural engineering [11].
The nonlinear problem is solved using the TO approach. This procedure assures better convergence and accuracy than the classical Newton approach. Such a technique achieves the solution through initial and corrections steps that account the tangent direction of the global equilibrium trajectory. The tangent search is performed by incorporating the derivative set of cohesive laws into the algebraic BEM equations. This type of operator has been utilized successfully in the literature for dealing many different nonlinear engineering problems. For instance, Botta et al. [12], Oliveira et al. [13] and Leonel et al. [14] applied this operator on the modelling of localization phenomenon, fracture analyses, and contact problems, respectively. The TO is derived for the crack growth analysis in concrete materials, which is the main contribution of this study.
Two applications were performed with the developed nonlinear BEM formulation. The results obtained were compared to numerical and experimental responses available in the literature to illustrate its robustness and accuracy. It was also verified that the TO provides faster convergence and lower computational time consuming for the nonlinear solution when compared to classical Newton method.

The Cohesive Crack Model
The cohesive crack model is an appropriate approach to represent the mechanical degradation processes in quasi-brittle materials, such as the concrete. In this model, the degradation phenomena occur along a fictitious crack positioned ahead the real crack tip. Thus, the material degradation zone is reduced by one dimension. The first studies in which the dissipation zone was reduced either to a curve for 2D problems or to a surface for 3D problems [1][2][3].
In the present study, simple softening laws approximate the residual material strength along the fictitious crack. These laws relate the fictitious crack opening displacement, COD, to the cohesive stresses values. Several cohesive crack laws relating cohesive stresses to the COD have been proposed in the literature. Three of them are often adopted to handle the crack growth analyses in quasi-brittle materials. The simplest law is given by a linear function relating the cohesive stresses to the fictitious COD smaller than the threshold value, COD c . For fictitious crack openings larger than COD c , cohesive stresses are assumed as null. The relations that represent the linear cohesive law are the following: in which E is the Young's modulus and f t the tensile material strength.
Alternatively, the COD and cohesive stresses may be related through the composition of linear functions. Particularly, such law is defined by two linear functions, which is named bi-linear model. This cohesive model is defined by the following equations: For the bi-linear model, the variables f " t ,COD " and COD c are defined as follows: in which G f represents the material fracture energy.
Finally, an exponential expression represents the third cohesive law utilized in this study. Equation (4) provides the analytical expressions for this cohesive model:

Boundary Integral Equations
The BEM has been widely applied to solve engineering problems such as contact, fatigue and crack propagation due to its high accuracy 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, v is the Poisson's ratio, u i are the components of displacements and b i are the body forces. The singular integral representation, written in terms of displacements, is obtained applying the Betti's theorem and the Eq.(5). This integral equation is written, disregarding body forces, as follows: (6) in which p i and u i are the tractions and displacements at the boundary, respectively. The free term c ij is equal to δ ij /2 for smooth boundaries. δ ij is the Kroenecker operator and P * ij and U * ij are the fundamental solutions for tractions and displacements, respectively, written for the source point s and related to the field point f [15]. Equation (6) leads to another important integral equation, known in literature as the hyper singular integral representation. To obtain it, the Eq.(6) has to be differentiated with respect to the directions x and y. Then, the relation among displacements and strains is used in order to determine an integral equation written in terms of strains. Then, a constitutive relation, the generalized Hooke's law for instance, is utilized to obtain an integral equation written in terms of stresses, for a boundary source point. Finally, the Cauchy formula is applied and the hyper singular integral equation is obtained, which is written in terms of tractions as follows: in which S * kij and D * kij contains the new kernels computed from P * ij and U * ij respectively [16,17] and η k is the outward normal vector at the source point.
Equations (6) and (7) can be applied separately to construct the system of algebraic equations for two-dimensional non-cracked domains. However, the application of only one integral equation, either Eq.(6) or Eq. (7), for modelling cracked structures in a single domain leads to the singularities in the system of algebraic BEM equations. Such singularities appear due to the presence of coincident source points at the crack surfaces. Therefore, to analyse structures containing cracks, some special BEM formulations were proposed in the literature. Among them, it is worth mentioning the dual BEM [16,17], the multiregion approach [10] and the dipole BEM formulation [13]. In the present work, the sub-region BEM technique is adopted, in which compatibility of displacements and equilibrium of forces are enforced along the interfaces of the sub-domains. The coupling procedure of sub-region formulation is presented in the next section.
To assemble the system of BEM equations, as usual, Eq.(6) and/ or Eq.(7) are transformed into algebraic relations by discretizing the boundaries and interfaces into elements over which displacements and tractions are approximated by polynomial functions. Besides that, one has to select a convenient number of collocation points to obtain the algebraic representations. The algebraic equations for boundary nodes crack opening displacement (COD) and the crack sliding displacement (CSD). Then, these variables are defined as follows: Where COD (k) is the vector of Crack Opening Displacements whereas CSD (k) is the vector of Crack Sliding Displacements. Each component of these vectors is defined for each pair of coincident nodes of an interface k.
Obviously, when COD (k) =CSD (k) =0, elastic condition is assumed for the interface. Otherwise, the cohesive crack model is incorporated. The cohesive crack propagation starts when at least one component of the vector p η (k) overcomes the material tensile strength f t . In such a case, cohesive behaviour for these tractions follow the stresses values predicted by the cohesive law, which depends on the COD values. Therefore, the problem becomes nonlinear and the tractions at the cohesive regions of the interfaces must satisfy the following constrains: Equation (13) indicates that the normal tractions p η must represent the residual material resistance, which is governed by the cohesive tractions p c (COD (k) ). On the other hand, the second constrain indicates brittle shear behaviour. When a component of the vector COD (k) overcomes the threshold value, COD c , the respective normal traction is assumed as null and the material has no longer a residual resistance to avoid the opening of the crack surfaces. Hence, to consider cohesive mechanical behaviour, the discontinuities defined in Eq. (12) and the constrains presented in Eq. (13) must be incorporated into the system of BEM equations, Eq.(11). Thus: Because the cohesive law was included, Eq.(14) becomes nonlinear due to the dependency of the variables p c and COD.

Nonlinear solution technique using the Tangent Operator approach
To solve the nonlinear problem provided by the cohesive crack model, Eq.(14) has to be rewritten as a function of its variables as follows: Where x B , u η L(k) , COD (k) , u t L(k) , CSD (k) are the variables of the problem, i.e., v.
An increment Δv must be obtained to solve the problem by the TO scheme. Such an increment has to lead to an equilibrium condition and the cohesive constraints have to be attended. Hence, one writes the function Y(v+Δv) as a linear Taylor expansion at the vicinity of a previous equilibrium configuration of the variables v as follows: are calculated considering these boundary collocation points.
After determining the displacement and traction fields at the boundaries, internal values for displacements, stresses and strains can be achieved. Internal displacements are determined using the integral Eq.(6) with the source point located inner the boundary. In this case, the free term c ij becomes δ ij . On the hand, the stress field at internal nodes is obtained through the integral stress representation [15].

Algebraic BEM Equations
One interesting approach to address crack growth problems through BEM is the sub-region BEM technique. In this approach, the entire domain is divided into sub-domains over which Eq.(6) and/ or Eq.(7) are applied. Then, the compatibility of displacements and equilibrium of forces are enforced along the sub-domain interfaces. As a result, the cracks are assumed to growth along such interfaces.
The crack surfaces, i.e. the interfaces, are discretized by oppositely oriented boundary elements. To analyse the crack growth in concrete, the cohesive crack model is incorporated into the interface elements. Thus, the cohesive law governs the COD and stresses values along the cohesive crack and, consequently, the nonlinear mechanical behaviour is introduced.
The classical set of algebraic BEM equations is written for a given multi-region domain taking into account the location of the collocation points, i.e. the source points. These points may be located at the left (L) or right (R) opposite boundary elements of a given interface k, or even at the external boundary (B). Hence, for a domain composed by NI interfaces, the BEM algebraic equations is presented as follows: , p η s(k) and p t s(k) as the normal and tangent components of displacements and tractions at a given side S, either L or R, of an interface k, respectively. Over the algebraic equations presented in Eq.(8), one imposes the boundary conditions at the external boundary, as usual in BEM, which leads to the following: in which, H η , H t and G η ,G t are the sub-matrices of H and G that multiply the normal and tangent displacements and tractions after the rotation, respectively. x B are the boundary unknown values, A B is composed by the columns from H B and G B , and F B contains the contribution of prescribed values at the boundary. To solve such problem, the equilibrium condition must be enforced along the interface boundary elements. Then, over such elements, the following conditions are imposed: Therefore, based on Eq.(10), one rewrites Eq.(9) as follows: The displacements calculated at the interfaces are used to define the (16) When the elastic condition is violated, the elastic prediction, i.e. COD (k) =CSD (k) =0, provides normal tractions higher than the tensile material strength, p n >f t =p c (COD=0). Therefore, the solution with null discontinuities is no longer acceptable. The difference between p η and p c (COD) is the normal traction exceeding vector Δp η exc , which must be reapplied along the cohesive interface regions to satisfy the constrains in Eq. (13). However, after reapplying it, the equilibrium condition will be violated, i.e, Y(v)=ΔF≠0. The vector ΔF represents the nonequilibrated force vector, which appears after the imposition of the cohesive constrains.
The solution of the problem requires Y(v+Δv)=0. Therefore, from is computed by the following differentiation: From Eq. (17), an approximated solution for the increment Δv, which leads to the solution of the problem, can be calculated as follows: In the system presented in Eq. (18), the matrix that multiplies {-ΔF} is the so called the Tangent Operator(TO). This matrix includes information about the derivative of the cohesive laws. Therefore, when incorporated into the BEM equations, it makes the search for the nonlinear solution along the tangent direction of the global equilibrium trajectory.
When linear or partially cohesive laws are adopted, the term ∂p c /∂COD (k) becomes constant. Consequently, the TO also becomes constant as well. Therefore, the solution Δv for a given load step can be obtained with only one iteration. On the other hand, when nonlinear cohesive laws are involved, such as an exponential one, the solution Δv must be obtained iteratively as the derivative terms ∂p C /∂COD depends on the actual openings COD (k) . At the first iteration, the increment Δv (1) must be computed form Eq. (18). Therefore, the initial traction exceeding vector Δp η exc(1) must be accounted in ΔF and the tangent terms ∂p C /∂COD are computed assuming COD=0. After the first increment, non-null COD are obtained. However, because the cohesive law is no longer linear, a new traction exceeding vector Δp η exc(2) is observed and must be reapplied at the cohesive regions. Therefore, the second iterative increment Δv (2) is computed from Eq. (18) in which the tangent terms ∂p C /∂COD (k) were updated with the new values of COD. The stop criterion for such iterative procedure was based on the norm of the traction exceeding vector, which must be smaller than a prescribed tolerance for achieving the converged solution.
As the boundary loads are prescribed incrementally, the TO procedure must be performed at each incremental step.

Applications
Two applications of fracture in concrete specimens are presented. The first concerns the three-point bended beam whereas the second addresses the four-point mixed mode multiple cohesive crack growth.

Three-point bended-notched beam
The first application of this study concerns a mode I fracture test, which was analysed experimentally by Saleh [18] and numerically by Oliveira et al. [13]. The specimen consists of a three point bended concrete beam with an initial notch. To define the three cohesive laws, two parameters were considered: the tensile material strength f t =3 MPa and the material fracture energy G t =75 N/m. Figure 1 illustrates the structure analysed, in which the geometric data and the material elastic proprieties are also presented.
The structure was divided into two sub-regions, with a vertical interface, in which the main crack propagates in mode I. The boundary mesh adopted is composed by 55 cubic and 29 quadratic discontinuous boundary elements with 307 source points. The displacements prescribed were imposed into 100 load steps and the tolerance for convergence was assumed as 10 -6 MPa. Figure 2 shows the initial and the final deformed structural configuration, with displacements magnified 100 times.
The displacements fields in the last load step are presented in Figure 3. It is worth mentioning the physical discontinuity introduced by the crack.
During the numerical analysis, the interface tractions were monitored and its cohesive behaviour is illustrated in Figure 4 (left) for linear, bilinear and exponential cohesive laws. As illustrated in this figure, the cohesive crack becomes a real crack in the end of the incremental load process. The nonlinear structural behaviour was also monitored through the load-displacement curves. The load was calculated as the equivalent force at the upper beam surface.  The displacement values were determined by averaging the vertical displacements at the crack mouth. Such nonlinear responses are illustrated in Figure 4 (right), in which experimental and numerical curves obtained by previous references are also illustrated.
The developed BEM formulation was capable to represent accurately the nonlinear mechanical behaviour of the cracked specimen, as illustrated in Figure 4. Moreover, the BEM model with the three cohesive laws achieved the resistant structural load with considerable accuracy if compared to the experimental value. The softening part of the numerical curves reproduce quite well the experimental results.
The convergence path obtained by the TO and the classical Newton approaches are presented in Figure 5. As presented in this figure, the TO approach makes the search for the equilibrium configuration by the tangent direction of the global equilibrium trajectory. Then, it is quite efficient in comparison with the classical approach.
To quantify the efficiency of the TO approch and the classical Newton scheme, the number of iteration required during the analyses were monitorated and presented in Table 1. The improvement in terms of the required iterations to achieve the convergence is remarkable with the TO. Such reduction reaches, for the best case, 92% and 77%, for the worst case. Then, the TO approach is recommended.

Four-point mixed mode multiple crack growth
The second application addresses the mechanical analysis of a mixed mode multiple cohesive crack growth case in a concrete specimen. The specimen is subjected to a four-point shear beam test, which as analysed numerically by Carpinteri [19] with the Finite Element Method (FEM) and by Saleh et al. [6] with the BEM. Both authors adopted the linear cohesive law to describe the residual material strength along the cohesive crack. The crack paths determined in Saleh et al. and Carpinteri [6,19] were adopted as interfaces in the model proposed in this study. Therefore, the entire domain was split into three sub-regions as presented in Figure 6, which also shows geometric data and the material elastic properties considered. In the analyses performed by the proposed BEM model, the cohesive laws linear, bilinear and exponential were adopted. The tensile material strength was defined as f t =2 MPa and the material fracture energy was assumed equal to G t =100 N/m. The adopted boundary element mesh was composed by 72 cubic and 30 quadratic discontinuous boundary elements with 378 source points. The displacements prescribed and the stop tolerance criterions were assumed as the same of the previous application. Figure 7 illustrates the initial and the final deformed configuration of the structure, with displacements magnified 500 times.
The displacements fields in the last load step are presented in Figure  8. The colour scale graphs show the physical discontinuity introduced by the cracks.
The load was considered as the equivalent force at the beam upper      face. The displacements values were considered as the prescribed ones. Figure 9 (left) shows the cohesive traction responses and Figure  9 (right) shows the load-displacement curves. Good agreement is observed among the numerical results considered. The developed BEM formulation was capable to represent the nonlinear structural behaviour for this application. As authors Saleh et al. and Carpinteri [6,19] used the linear cohesive law, the results obtained with such a law presented the best agreement. However, the three cohesive laws applied given good accuracy in terms of maximum resistant load.
The numerical efficiency was compared between the TO approach and the classical Newton method. The efficiency in terms of the required amount of iterations is presented in Table 2. Strong improvement is achieved with the TO, 93% for the best case and 84% for the worst case.

Conclusions
A BEM formulation has been presented for modelling the fracture of concrete structures using the cohesive crack model. For the considered applications, the results obtained showed strong agreement with numerical and experimental responses available in the literature. Furthermore, the TO scheme has demonstrated its efficiency as it required an amount of iterations considerable lesser than the classical Newton approach. Therefore, such an operator may lead to remarkable gains in terms of computational time consuming when applied in most complex engineering problems involving multi-fractured structures.