The fracture failure of three-dimensional concrete structures subjected to concentrated loadings using the Boundary Element Method

: This study applies the Boundary Element Method (BEM) for the fracture failure modelling of three- dimensional concrete structures subjected to concentrated boundary conditions. The non-requirement of domain mesh by the BEM enables high accuracy in the domain fields assessment in addition to less complex remeshing procedures during crack propagation. However, concentrated boundary conditions often occur in fracture mechanics. The Lagrangian version of the BEM enforces such boundary conditions approximately by small length elements, which lead to numerical instabilities or even inaccurate problem representation. This study proposes a formulation for representing properly concentrated boundary conditions within the Lagrangian BEM framework. Nonlinear fracture mechanics describes the material failure processes herein. The classical cohesive crack approach governs the nonlinear energy dissipation processes, in which constant and tangent operators solve the resulting nonlinear system. Three applications demonstrate the accuracy of the proposed formulation, in which the BEM responses are compared against numerical and experimental results available in the literature.


INTRODUCTION
The accurate prediction of material failure processes is a concern for structural engineers. The knowledge and effective modelling of such processes enable the design of efficient and safe structural elements. Because cracks explain consistently several mechanical-material failures, fracture mechanics theories emerged and demonstrated its robustness in the past decades [1], [2]. The material discontinuity associated to the crack leads to the stresses concentrations at its tip, which is the preferential triggering failure point. The Linear Elastic Fracture Mechanics (LEFM) predicts correctly the mechanical behaviour when the fracture process zone (FPZ) surrounding the crack tip is sufficiently small. In this case, Stress Intensity Factors (SIFs) govern the collapse process due to the stresses singularities at the crack tip. On the other hand, large FPZ leads to relevant nonlinear effects associated to the material damage. The cohesive crack approach represents accurately these nonlinear effects, in which cohesive stresses indicate the residual material resistance at the FPZ. Besides, cohesive laws associate the cohesive stresses to the crack opening displacements. It is worth mentioning that the mutual dependence among cohesive stresses and crack opening displacements values triggers the nonlinear problem, in which efficient iterative schemes allow the accurate solution [3].
Analytical approaches handled initially fracture problems, particularly within LEFM framework. However, the complexities upon geometry and boundary conditions limit the proper solution to few cases. Then, robust solution strategies in this domain require numerical methods, such as the Finite Element Method (FEM). However, its domain mesh leads to the complex remeshing schemes during the crack propagation modelling. Besides, the FEM requires fine meshes surrounding the crack tip for the accurate representation of stresses concentrations. Despite these limitations, some FEM approaches addressed fracture problems. It is worth mentioning Bittencourt et al. [4] and Bouchard et al. [5] in two-dimensional and Bremberg and Dhondt [6], Nejati et al. [7] and Wowk et al. [8] in three-dimensional LEFM modelling, among many others. Besides, Galvez et al. [9] and Dimitri et al. [10] addressed two-dimensional cohesive fracture problems whereas Gaedicke et al. [11], Chen et al. [12] and Skar and Poulsen [13] focused on specific aspects of the three-dimensional cohesive approaches.
The above-mentioned FEM limitations can be avoided by using adequate functions into the mechanical fields description, particularly the displacements. These special functions provide a priori knowledge about the space solutions, which enable accurate mechanical representation. This strategy triggers the eXtended Finite Element Method (XFEM). Then, the SIFs can be assessed accurately within relatively coarse meshes and the remeshing process efforts can be minimised or even avoided. Some contributions of XFEM formulations in LEFM are available in [14]- [16], for instance. In the nonlinear fracture mechanics field, it is worth mentioning the developments provided by Duarte et al. [17], Moes and Belytschko [18], Mohammadnejad and Khoei [19], Jaśkowiec et al. [20] and Ferte et al. [21], among others.
On the other hand, the Boundary Element Method (BEM) demonstrates excellent performance in fracture problems, especially in the three-dimensional context as demonstrated by Cordeiro and Leonel [22], [23] and Mi and Aliabadi [24], [25]. The BEM does not require domain mesh, which enables accurate assessment of internal mechanical fields and non-complex remeshing procedures during the crack growth. Several BEM formulations have been proposed in the literature for the solution of fracture problems, such as single-domain technique [26], multi-domain technique [27], dipoles approach [28] and cells with embedded discontinuities [29], for instance. Nevertheless, the robust BEM approach in fracture mechanics is the dual BEM (DBEM) formulation [24], [30]. The DBEM utilises two different integral equations: displacement or singular boundary integral equation and traction or hypersingular boundary integral equation. This approach demonstrated its accuracy in the literature, in which two-dimensional [31], [32] and three-dimensional problems [33]- [36] have been addressed properly. However, the cohesive crack growth modelling in three-dimensional context by BEM has not been explored totally in the literature, which inspired the developments presented herein.
Enrichment techniques have been explored vastly in the context of XFEM. However, its application in the coupling of BEM and fracture mechanics is quite recent. Simpson and Trevelyan [37], [38] were the pioneers in this type of enrichment within BEM framework, which improved the displacement solutions near the crack tip. Besides, extensions of this type of approach in LEFM context are available for two-dimensional isotropic [39], [40] and anisotropic [41] materials. Despite these advances, three-dimensional XBEM formulations for crack propagation and fracture problems have been explored marginally in the literature especially for the accurate representation of boundary conditions in nonlinear fracture mechanics field.
This study proposes a BEM formulation for the fracture failure modelling of three-dimensional structures composed of quasi-brittle materials and subjected to concentrated boundary conditions. The classical cohesive crack approach governs the nonlinear problems, which enables the accurate fracture and crack propagation modelling in quasi-brittle materials. Two different nonlinear solution techniques solve the nonlinear problem: the constant operator and the tangent operator. The constant operator keeps constant all relevant influence matrices during the iterative procedure. On the other hand, the tangent operator incorporates the cohesive law derivatives, which improves the convergence rate. The performance of each nonlinear solution scheme and its accuracy have been tested in two applications. Despite being an ordinary task for FEM-like models, the enforcement of concentrated boundary conditions in the context of elastostatics is not a standard operation within the Lagrangian version of BEM. This particular boundary condition type can be enforced approximately by small length or even narrow boundary elements. Nevertheless, because of the singular nature of the BEM fundamental solutions, collocation points closely positioned lead to the instabilities into the algebraic system of equations and inaccurate problem representation. It is worth stressing that the topic has been explored marginally in the literature in conventional elastostatic problems. Qian [42] utilised the Betti's reciprocity theorem and obtained a boundary integration equation with incorporated concentrated forces. Wang and Wang [43] obtained similar developments accounting for spline approximation in the boundary integration equation. Recently, Zhou et al. [44] proposed singular logarithmic shape functions for approximating tractions and displacements along the element domain. In all cases, the validity and accuracy of the proposed methods require further improvements, which also motivated the developments presented herein. It is worth stressing that concentrated forces can be handled straightforwardly within BEM framework when such forces are at the material domain. Kamiya and Sawaki [45] show the solution for this case in plates. However, the proper modelling of concentrated forces and concentrated supports at the body's boundary is still a challenge nowadays. This problem is solved elegantly herein by expanding/enriching the traction field along the boundary elements subjected to the concentrated boundary conditions. This scheme avoids singularities and enables the adequate representation of concentrated loads and supports. Three applications demonstrate the accuracy of the proposed formulation. The responses achieved by the proposed formulation are compared against numerical and experimental results available in the literature.

Boundary Integral Equations
Integral equations describe the mechanical behaviour of solids within BEM framework. Particularly, the DBEM formulation accurately represents this behaviour in cracked bodies. The DBEM requires two different integral equations because collocation points share the same position along the crack surface. The weighted residual method, the BEM fundamental solutions and elasticity equations enable the displacement or singular Boundary Integral Equation (BIE), which for a three-dimensional domain composed of boundary and nil body forces is as follows [25]: in which and are the source (collocation) and field points, respectively. The terms * and * indicate traction and displacement Kelvin fundamental solutions. is the free term, which is equal to 2 for smooth boundaries. is the Kronecker delta. and represent the displacements and tractions along the boundary. The fundamental solutions presented in Equation 1 are: indicates the components of the distance vector and = ‖ ‖ = � − � its norm; represents the components of the normal outward vector; and are the shear modulus and Poisson ratio, respectively. Because of the singularities into the fundamental solutions, the integral on the left-hand side of Equation 1 requires the evaluation in the Cauchy Principal Value sense.
The Equation 1 and classical elasticity expressions enable the well-known traction or hypersingular BIE. This BIE can be obtained by the differentiation of the displacement BIE in terms of the field point position. Then, one applies the relation between displacements and strains, the Hooke's law, and the Cauchy equilibrium formula, which lead to the following [24]: where * and * are the derivative kernels associated to * and * , respectively. Then, the singularities in these kernels are of order * = 1 3 and * = 1 2 . Because of the hyper singularity on kernel * , this integral term requires the evaluation into Hadamard Finite Part sense. For sake of completeness, the hyper singular kernels are as follows [24], [25]:  Figure 1 illustrates the collocation strategy of DBEM at the crack surfaces. Then, the displacement BIE discretises the external boundary and one crack surface whereas the traction BIE discretises the opposite crack surface. In the DBEM, the displacement BIE is written for collocations at the upper crack surface + as follows:

The DBEM collocation strategy and numerical assembling process
Complementarily, the traction BIE is written at the lower crack surface − as follows: It is worth mentioning that Equations 5 and 6 contain additional free terms, which account for the collocations at the same position. Besides, + and − at upper and lower crack surfaces, respectively, share the same position and have opposite normal outward vectors. The geometry and mechanical fields on Equations 1, 5 and 6 have been approximated by standard quadrilateral isoparametric linear Lagrangian elements. The boundary mesh can be composed of continuous, semi-continuous and discontinuous boundary elements. These element-types have been utilised along the external boundary according to the geometry and boundary conditions continuities. Solely discontinuous boundary elements have been utilised along the crack surfaces as classically required by the dual BEM approach. Besides, the Gauss-Legendre numerical scheme handles the integration over non-singular boundary elements, in which sub-element technique improves its accuracy [23]. The singularity subtraction technique proposed by Guiggiani and Gigante [46] and Guiggiani et al. [47] regularizes the kernels over singular-elements and enables the accurate assessment of singular and hypersingular BIE. It is worth mentioning that Cordeiro and Leonel [23] presents the regularization expressions for hypersingular kernels, which have been utilised herein.
The integral equations Equations 1, 5 and 6 lead to the classical algebraic representation of BEM, i.e., = . Thus, and are the DBEM influence matrices whereas the vectors and indicate the displacement and tractions at the collocation points, respectively. contains the integration of kernels * and * and contains the integration of kernels * and * . This linear system of equations can be solved by enforcing the boundary conditions, as usual in BEM. Then, all unknowns are moved to the left-hand side, and all knows are moved to the right-hand side of this equality. Hence, the final linear system of equations = can be solved and the unknowns at the entire boundary, , determined.

Concentrated force formulation
Structural modelling frequently idealises tractions over small surfaces as concentrated forces. This idealisation has been accepted fully in solid mechanics domain. However, the accurate representation of concentrated forces is a complex task for BEM, which applies tractions at the boundary. This condition has been enforced approximately by small length elements, which lead to the collocation points closely positioned in the boundary mesh. Because of the singular nature of the BEM fundamental solutions, this condition triggers numerical instabilities upon the algebraic system of equations, especially in three-dimensional problems.
This study proposes a simple and effective scheme for the accurate representation of concentrated forces. The proposed formulation utilises the Dirac delta function for improving the approximation over the traction field. Therefore, the traction field over the boundary can be rewritten as follows: where indicates the concentrated force along the j direction applied at the point of the boundary mesh. indicates the shape functions. It is worth emphasizing that → ∞ at , which characterizes the concentrated force. Besides, the formulation enables the accurate representation of concentrated forces applied at different positions along the boundary. ( − ) indicates the Dirac delta function.
Then, the traction form in Equation 7 enables the updating of the right-hand side of Equation 1 as follows: It is worth mentioning that the first integral term on the last equation is part of the classical BEM formulation. On the other hand, the second term can be further modified by applying the Dirac delta sifting property, which transforms this integral term in a simple evaluation of the fundamental solution at points and , as follows: Therefore, the improvement term becomes a simple multiplication, which does not require additional numerical integration procedures. Because the improvement terms are known in this numerical procedure, such terms lead to an independent vector. This vector simply adds to the vector of known values at the boundary. Therefore, the algebraic system of equations can be updated as follows: where and + indicate the vector arising from the multiplication of the fundamental solution * and the concentrated forces whereas − results from the multiplication of the fundamental solution * and the concentrated forces at the hypersingular crack surface. The superscripts + and -indicate the mechanical variables along the upper and lower crack surfaces, respectively.

Concentrated supports formulation
Concentrated support conditions have been utilised largely in solid mechanics problems. Such boundary condition idealises supports on small surfaces as concentrated in a single point. Because of the singularity characteristics of the fundamental solutions previously mentioned, the accurate representation of this boundary condition is not a standard task for BEM.
Concentrated supports trigger concentrated forces, which are the supports reactions. Therefore, the formulation presented on section 3.1 can be utilised for the accurate modelling of this boundary condition. Therefore, the concentrated forces described in Equation 7 and Equation 8 represent the support reactions herein, which are unknown in the problem. Because additional unknowns have been included in the problem (the supports reactions itself), this scheme requires additional compatibility conditions. These additional compatibilities are the polynomial approximations for displacements into the boundary element containing the support, which can be expressed by: where 1 and 2 are the dimensionless coordinates associated to the concentrated support position.
It is worth mentioning that the strategy does not modify the conventional influence BEM matrices. It requires a slightly modification on the algebraic system of equation, which accommodates new lines and columns associated to the concentrated support conditions. Thus, the algebraic system of equations in this scheme is as follows: in which results from the columns change procedure between and matrices. contains the influence terms associated to the unknowns at the boundary. represents the unknowns at the boundary whereas indicates the influence vector associated to the known/prescribed values at the boundary. indicates the support reaction values. The sub-matrix is composed of shape function values evaluated from Equation 11, i.e., = ( , ) ( 1 , 2 ), in which ( , ) is the element incidence associated to the global degree of freedom in the element . ̄ describes the concentrated displacement values prescribed at the boundary, which assumes nil values in the case of static supports. It is worth remarking that has similar meaning of in Equation 7.
contains the evaluation of the fundamental solution at source point and support point, which is as follows: where only the column , associated to the direction of the prescribed displacement at the point , has been accounted into the assembling system. The fundamental solutions * and * have been utilised according to the BIE associated to the source point.
Both concentrated boundary conditions have been prescribed into global coordinates, i.e., coordinates x, y and z. Based on this information, a simple Newton-Raphson scheme achieves the dimensionless coordinates 1 and 2 within the boundary element domain. For sake of completeness, Appendix 1 presents a flowchart, in which this procedure is described clearly.

NONLINEAR FRACTURE MECHANICS
The LEFM represents properly the mechanical behaviour of cracked solids in which the FPZ surrounding the crack tip is small in comparison to the crack dimensions. However, various materials do not follow this pattern. In quasibrittle materials, for instance, the FPZ is sufficiently large and triggers nonlinear effects that cannot be disregarded. The cohesive crack model represents accurately the residual material strength and the energy dissipation phenomena in the FPZ of quasi-brittle materials. Then, cohesive tractions close the crack surfaces and represent the residual material strength, as illustrated in Figure 2. Besides, the cohesive traction values have been associated to the crack opening displacement (COD) by cohesive laws. Several cohesive laws have been proposed in the literature to govern the cohesive traction behaviour [48]- [50]. Three of them have been applied often in this problem: linear, bi-linear and exponential, as illustrated in Figure 3. The linear cohesive law associates and by a linear function. This law requires the material tensile strength and the threshold opening displacement values. Then, this law is as follows: The bi-linear function leads to the bi-linear cohesive law, which is as follows: where the additional parameters * , * and have been presented by Petersson [51] accounting for concrete: * = 3 * = 0.8 = 3.6 (16) where indicates the fracture energy.
Finally, the exponential cohesive law associates and as follows: The classical cohesive laws represent the material nonlinearity in the modelling. Thus, the cohesive traction values depend on the COD and vice-versa. This nonlinear problem has been solved in the present study by the Newton-Raphson like schemes, in which the DBEM formulation provides the mechanical fields' values. The nonlinear solution technique is presented in the following.

Algebraic description
The nonlinear formulation utilises the BEM integral equations previously presented. Thus, this BEM approach provides the displacements and tractions values along the FPZ, which is named herein as cohesive interfaces. In addition, the proposed boundary fields improvement schemes enable the proper modelling of concentrated boundary conditions, which is an element of novelty of this study.
The algebraic formulation description starts by splitting the classical influence BEM matrices. Thus, such matrices account for the collocation point positions, which are either external boundary (b), left (L) or right (R) cohesive interfaces. Then, the classical algebraic system of equations can be rewritten as follows: The mechanical fields along the FPZ can be described in terms of the local coordinate system illustrated in Figure 4. Such operation enables the straightforward enforcement of the crack displacement discontinuities, i.e. the crack opening displacements. Thus, the variables at the right and left crack surfaces can be rotated accordingly, which leads to the following: where , , , , , , , , , , and indicate the product between the rotation matrix and the corresponding influence matrices. and indicate tangential crack surface directions whereas refers to the normal crack surface direction. enables the introduction of the crack displacements discontinuities: crack opening displacement (COD), crack sliding displacement (CSD) and crack tearing displacement (CTD). Such variables account for the local coordinate system defined in Figure 4. Besides, it has been defined as follows:

Equation 20 modifies Equation 19
, which eliminates the left crack displacements in the formulation. Besides, the latter equation can be further modified by enforcing the equilibrium conditions along the crack surfaces, i.e., , and . In addition, the cohesive approach predicts as a function of COD. Finally, one enforces the boundary conditions at the external boundary. All those algebraic manipulations lead to the following: in which matrix multiplies the unknown fields at the external boundary, obtained from the columns change procedure of BEM, i.e., by enforcing the boundary conditions. Besides, results from the product between the known values at the boundary and the respective influence values.
The cohesive law governs ( ), which describes the material nonlinearity at the FPZ during the crack propagation. Because Equation 21 is nonlinear, it requires an incremental-iterative nonlinear solution technique based on trial and correction steps. In the present study, these solutions have been proposed in the context of Newton-Raphson like schemes. Thus, this approach assumes an elastic prediction as trial step. Afterwards, the correction steps follow the cohesive law threshold values. Each step leads to the increment values over each mechanical variable in the analysis. Such increments must be added, which provide its accumulated values. Obviously, ( ) updates at each correction step according to the adopted cohesive law. The tractions values at the FPZ predicted in the elastic step may differ from the limit values assumed by the cohesive law. This difference may also appear during the iterations of the correction steps. This difference is named as exceeding traction . Then, the ( ) updates as follows: The nonlinear solution process converges when the norm of is smaller than a prescribed tolerance. Then, the non-equilibrated traction vector is sufficiently small. In the present study, the correction steps, and consequently the search for the equilibrium configuration at each load step, have been solved by two different techniques: Constant Operator (CO) and Tangent Operator (TO). In both, the solution is incremental. However, CO keeps constant all relevant influence matrices during the correction steps. On the other hand, TO updates the search for the equilibrium configuration accounting for the tangent direction provided by the cohesive law. This study assesses the performance of each nonlinear solution technique in three-dimensional cohesive crack growth, which is an element of novelty.

Solution Technique by Constant Operator (CO)
The CO is a nonlinear solution technique, in which all relevant influence matrices are kept constant during the iterative process. This approach utilises the initial body stiffness and the search for the equilibrium configuration reapplies the exceeding traction within the initial stiffness [27]. In this approach, the first step refers to the elastic prediction, which provides the increments on the mechanical fields. Then, the exceeding traction i is evaluated from Equation 22. Thus, for a given correction step i, the exceeding traction leads to the following increments in the mechanical fields: (23) where the increments on the mechanical fields are calculated and added to the responses from the previous iterations. Besides, the search for = 0 occurs simultaneously for all collocations at the cohesive interface. Then, the convergence occurs when � � is smaller than a prescribed tolerance, which indicates that the non-equilibrated traction vector is sufficiently small.

Solution Technique by Tangent Operator (TO)
The solution technique via TO requires the definition of a residue from the algebraic equations presented in Equation 21 [27]. Such residue accounts for the problem variables as follows: where indicates the vector of mechanical fields. It is worth mentioning that the equilibrium configuration leads to the nil residue value. Then, in the (i-1)-th iteration, the coupling of Equation 22 in Equation 24 enables the definition of the residue variation as follows: The TO determines the increment set leading to the nil residue ( + )value. For this purpose, one assumes the residue function as continuous. Therefore, this function can be expanded in a linear Taylor's expansion around . Finally, one truncates this expansion to the first term, which provides the corrections as follows: where represents the cohesive law derivative, which modifies the BEM influence matrices. Then, the influence matrices associated to such derivative compose the tangent operator. Therefore, the TO incorporates the mechanical degradation properties into the algebraic system of equations, which enables the efficient search for the nonlinear solution. Then, the equilibrium trajectory follows the tangent direction of the global equilibrium path, which names this solution technique. It is worth mentioning that linear functions often describe the cohesive laws, such as observed in linear and bi-linear cohesive laws, see Figure 3. In such cases, becomes constant and the convergence can be achieved by solely one iteration. Then, the TO solution scheme tends to be numerically efficient when compared to the classical procedures. In this scheme, the iterative process stops when � � is smaller than a prescribed tolerance.

APPLICATIONS
Three applications demonstrate the accuracy of the proposed formulations. The first application handles a LEFM problem in pure mode I, in which the BEM predictions have been compared to the responses provided by an analytical solution and FEM/ANSYS. The second application illustrates the cohesive cracks propagation modelling in the fourpoint bending test. This application demonstrates the robustness of the proposed formulation in the multiple cohesive crack growth modelling. The last application presents the nonlinear fracture analysis of the Brokenshire torsion test, which is subjected to complex crack propagation in mixed modes I-II-III.

Pure mode I problem
The first application handles the mechanical analysis of the prismatic cracked solid illustrated in Figure 5. This structure is subjected to concentrated boundary conditions (forces) at its top and bottom surfaces.  Figure 5. The meshes A, B, C and D contain 543, 1,482, 2,887 and 4,755 linear quadrilateral boundary elements, which lead to 1,156, 2,972, 5,636 and 9,148 collocation points, respectively. The FEM model utilises a mesh with 609,508 nodes and 138,240 hexahedral quadratic elements (SOLID 186). It is worth mentioning that a previous convergence analysis enabled the FE mesh. For sake of completeness, an additional analysis accounts for uniform traction applied at the structure top instead of concentrated force. The results of this modelling hypothesis demonstrate the importance of handling accurately the prescribed boundary conditions. In such modelling, concentrated supports avoid rigid body motion, similarly as the previously mentioned cases. Gross et al. [52] and Tada et al. [53] present the analytical plane solution assuming LEFM behaviour, in which the crack mouth opening displacement (CMOD) can be determined as follows: for plane strain and ′ = for plane stress (27) where is the applied normal traction, is the crack length, is the solid width, is the Young's modulus and is the Poisson ratio. In this application, = 2.5 and = 5.0. Then, the analytical CMOD value is 0.0448 (displacement unity) for plane strain case. It is worth stressing that the analytical solution assumes the application of uniform tractions at the structure ends, as illustrated in Figure 6. Therefore, this application demonstrates the accuracy of BEM in fracture problems. Additionally, the Saint-Venant principle enables the performance assessment between FEM and BEM within concentrated boundary conditions framework. The concentrated force scheme has been applied initially in this analysis. This approach leads to consistent displacement behaviour as illustrated in Figure 7, which accounts for the displacement along 2 direction for mesh D. It is worth emphasizing the displacement discontinuity caused by the notch in addition to the good agreement with FEM/ANSYS predictions.   Figure 8 illustrates the CMOD responses at the crack front points as a function of the boundary mesh refinement considering the concentrated force formulation (conc) and the uniform traction application (dist). Besides, Figure 9 presents the convergence behaviour for CMOD at the point (0.0; 25.0; 2.5) accounting for the percentual error between the BEM and analytical responses approaches for the concentrated forces at the structural ends (Conc) and for the uniform traction along these faces (Dist). The proposed scheme provided good agreement with the reference solution. As expected, the mesh D provided the best solutions, in which the displacements at the point (0.0; 25.0; 2.5) obtained by mesh D differs from the reference solution by 1.50%. On the other hand, the FEM response at the same point differs from the analytical solution by 12.37%, which indicates the superior performance of BEM in the problem. Because the analytical solution considers the remote uniform traction, there is an excellent agreement between the analytical predictions and the BEM modelling assuming uniform traction. Thus, these results illustrate the accuracy of BEM in fracture mechanics problems in addition to the superior performance of BEM in comparison with FEM when concentrated boundary conditions are observed.    For sake of completeness, the concentrated supports formulation has been applied also in this application to attest its accuracy. Then, prescribed concentrated displacements are applied instead of concentrated forces in the following results. Moreover, to keep the comparison reference, the concentrated displacements have been enforced at the same positions of the concentrated forces modelling previously presented. Besides, the prescribed displacement value refers to the displacement evaluated in the respective concentrated force modelling. Therefore, the reaction force achieved by the concentrated displacement formulation must agree with the prescribed concentrated force value. Table 1 presents the reaction force values obtained by each BEM mesh and the correspondent prescribed displacement. Good agreement is observed once the prescribed force value is 25. Additionally, Figure 11 presents the CMOD responses for concentrated force (force) and concentrated supports (displ) conditions. Excellent agreement is observed for the correspondent meshes, which illustrates the accuracy of the concentrated displacement formulation.  The results presented for the BEM models in the last figures demonstrate the accuracy of the proposed boundary conditions representation scheme. Good agreement has been observed as the number of degrees of freedom increase, as expected. Moreover, smooth convergence behaviour was observed, which demonstrates the numerical stability of these formulations. Besides, it is remarkable the superior performance of BEM in comparison to FEM in this problem, in which BEM achieved accurate responses with less degrees of freedom.

Four-point bending problem. Multiple cohesive crack growth
This application deals with the cohesive crack growth modelling of the specimen illustrated in Figure 12. The cracks positions and boundary conditions trigger the multiple crack propagation in mixed fracture mode I-II. Bocca et al. [54] analysed this specimen experimentally whereas Saleh and Aliabadi [55] performed the numerical analysis using a two-dimensional BEM approach. The specimen geometry is parallelepiped with length of 0.8m, height of 0.2m, thickness of 0.2m and notches depth of 0.04m. The material properties are, according to Bocca et al. [54]: Young's modulus = 27 , Poisson ratio = 0.1, tensile material strength = 2 and fracture energy = 100 .  Figure 13 illustrates the final boundary element mesh in this analysis, which contains 2,347 collocation points and 1,379 quadrilateral linear boundary elements. Besides, Figure 14 presents the discretisation along the cohesive interfaces for the last load step, in which each of them is composed of 320 collocation points and 80 quadrilateral linear boundary elements. This figure also illustrates the crack surfaces geometry determined by the experimental analysis [54] and the two-dimensional numerical modelling [55], in which good agreement has been observed among the predictions of the proposed formulation and that provided by [54].  -Crack geometry. present study, experimental [54] and numerical [55] The proposed boundary conditions representation schemes (conc) handle the boundary conditions in this application. Then, the prescribed displacements have been enforced in concentrated form by the concentrated supports approach. For sake of simplicity, concentrated supports have been applied at the centre of each boundary element positioned along the supports' regions indicated in Figure 12. Displacements along directions 1 and 2 are nil in the following ranges 1 = 0.32; 2 = 0; 0 ≤ 3 ≤ 0.2 and 1 = 0.8; 2 = 0; 0 ≤ 3 ≤ 0.2 . Besides, the points (0.32; 0.0; 0.0)m and (0.8; 0.0; 0.0)m have concentrated supports along 3 direction, which avoid rigid body motion. For the sake of comparison, the standard boundary conditions application (dist) has been enforced along the entire supports' regions indicated in Figure 12. The non-nil displacements have been applied within 16 increments, in which the tolerance for convergence is 10 −3 , in terms of � �. In addition, the cohesive laws linear, bi-linear and exponential represent the residual material strength at the FPZ in this analysis. The CO and TO solve the nonlinear equations in this application. It is worth mentioning that Saleh and Aliabadi [55] applied solely the linear cohesive law in their numerical modelling. Figure 15 illustrates the performance of the proposed nonlinear BEM formulations in this application. The dimensionless load and deflection values have been utilised for this purpose. The dimensionless load relates the reaction forces 1 and 2 at the supports to the prescribed displacement to geometrical and material parameters as follows: 1 2 adim T P P F f dt   (28) in which is the tensile material strength and and are, respectively, the specimen's height and the thickness. Besides, the dimensionless displacement corresponds to the division between 6.10 4 times the applied displacement and the specimen's height. As expected, the numerical responses demonstrate good agreement at the elastic part of the curve. Besides, the BEM results strongly agree to the reference [55] for the linear cohesive law, once the reference applied solely such cohesive law. In addition, the peak loads predicted by BEM through bi-linear and exponential cohesive laws are slightly smaller than the reference value, as expected. However, the latter cohesive law leads to the smooth behaviour during the softening part. The BEM solutions evaluated by CO and TO are in excellent agreement among them. Such behaviour indicates the accuracy of the nonlinear solution schemes into the determination of the equilibrium configuration at each load step. Finally, the application of boundary conditions in conventional form (dist) leads to rigid response, as expected. Because the experimental and numerical references account concentrated boundary conditions, the conventional scheme does not represent accurately the boundary conditions in this problem. Then, the application of boundary conditions using the proposed formulation (conc) enables accurate responses in comparison to the conventional scheme. Table 2 presents the number of iterations required by each nonlinear solution technique in the search for the equilibrium configuration. Then, this table enables the efficiency analysis by the ratio TO/CO. The TO leads to reductions superior to 82% in this application, which demonstrates its efficiency. For sake of clarity, Figure 16 illustrates the displacement values along the 2 direction and the deformed shape in the last load step 1,000 times magnified. The linear cohesive law and the TO nonlinear solution schemes provided these values. It is worth remarking the strong discontinuities into the displacement fields introduced by the cracks.

Brokenshire torsion test: complex mixed mode crack propagation path
The Brokenshire test refers to a torsion test in a notched prismatic specimen composed of concrete. The experimental details and results are in Jefferson et al. [56]. Besides, many authors [21], [48], [49] have analysed this problem via FEM/XFEM approaches because it requires a full three-dimensional formulation. Nevertheless, by the best knowledge of the authors, this problem has not been solved by BEM until the present. Figure 18 illustrates the geometry and boundary conditions. The prismatic solid has length of 0.4m, width of 0.1m and height of 0.1m. There is a notch inclined by 45º at the centre, which depth is 0.05m. Besides, the steel frame is placed at both specimen ends, where the boundary conditions have been applied. It is worth mentioning that the steel frame has been accounted in the modelling and represented by the classical sub region BEM technique. The material properties [56] are: Young's modulus of concrete = 35 , Poisson ratio of concrete = 0.2, Young's modulus of steel frame = 200 , Poisson ratio of steel frame = 0.3, energy fracture = 120 and tensile material strength = 1.6 . The supports positioned along 3 = 0 restrict displacements along directions 1 , 2 and 3 . The other supports in this modelling restrict the displacement solely along the 2 direction. Besides, the prescribed displacement is 0.8mm along 2 direction, which has been applied within 16 load steps. The boundary mesh is composed of 7,430 collocation points and 2,592 quadrilateral boundary elements with linear approximation as illustrated in Figure 19. The CO and TO achieve the equilibrium configuration in the nonlinear problem, in which linear, bi-linear and cohesive laws represent the material residual resistance at the FPZ. The tolerance for convergence is 10 −3 , in terms of � �. Because of the accurate results achieved in the previous application, solely the proposed BEM formulation is applied herein. Besides, the experimental scheme described in Jefferson et al. [56] accounts for concentrated boundary conditions.  The numerical responses obtained by the proposed BEM formulation are compared to experimental results presented by Jefferson et al. [56], particularly, the pair values of load and crack mouth opening displacement (CMOD) between the points A and B illustrated in Figure 18. Figure 20 illustrates this comparison, in which one observes good agreement among numerical and experimental responses. It is worth emphasizing that the application involves a complex mixed mode fracture problem, in which modes I and III are preponderant. Despite the good agreement observed among numerical and experimental results, the differences observed in Figure 20 might reduce if the material strength is penalized accordingly for mode III solicitation, once the classical cohesive model utilised herein penalizes solely the material strength along the perpendicular direction to the crack surfaces.
The proposed formulation provides accurate result for the peak load value, especially through the linear cohesive law. Besides, all cohesive laws lead to the adequate modelling of the softening part. In addition, the CO and TO responses are in excellent agreement, which demonstrates the numerical stability of the proposed nonlinear solution techniques.  Table 3 presents the number of iterations required by CO and TO in the search for the equilibrium configuration in this application. Analogously to the previous application, TO demonstrates excellent performance. In addition to the accuracy, TO lead to reductions superior to 58% in the present analysis.  Figure 21 presents the evolution of the mechanical degradation processes at the FPZ. In this figure, zero indicates free traction crack surfaces and one represents continuous material. The results in this figure account for the linear cohesive law and CO solution technique. It is worth remarking the significant increase of the FPZ size between increments 8 and 12, in which the softening behaviour start. The crack path predicted by the BEM approach has been compared to the experimental response of Jefferson et al. [56] and the XFEM results obtained by Kaczmarczyk et al. [57] in Figure 22. This figure illustrates good agreement between the numerical approaches. Besides, one observes good agreement among numerical and experimental results in the beginning of the propagation path and slight differences near the specimen collapse.

CONCLUSIONS
This study handled the mechanical representation and crack propagation modelling of three-dimensional cracked quasi-brittle materials by the BEM. In addition, this study presented a simple and effective scheme for representing accurately concentrated boundary conditions within BEM framework. Despite being a standard task for FEM-like models, concentrated boundary conditions are challenge for BEM, which applies boundary conditions distributed over elements. Concentrated boundary conditions have been applied approximately by small size boundary elements by the Lagrangian version of BEM, which requires fine meshes. Besides, singularities may appear in the final algebraic system of equations because of the singular nature of the BEM fundamental solutions. Therefore, the proposed formulation contributes for the accurate representation of concentrated boundary conditions within BEM. The solutions obtained by the proposed formulation are accurate than the responses provided by the standard BEM approach for the same mesh refinement level, as illustrated in considered applications. Besides, the proposed formulation led to good performance in comparison to experimental tests, as illustrated in all applications.
Further extensions of the proposed schemes can consider three-dimensional isogeometric BEM approaches, in which the representation of concentrated boundary conditions is still more complex than the observed in classical Lagrangian approaches.