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

Falha por fratura de estruturas em concreto submetidas a carregamentos concentrados pelo Método dos Elementos de Contorno

Abstract

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.

Keywords:
concentrated boundary conditions; nonlinear fracture mechanics; enrichment functions; three-dimensional fracture modelling

Resumo

Este estudo utiliza o Método dos Elementos de Contorno (MEC) para a modelagem da falha por fratura de estruturas tridimensionais em concreto submetidas a condições de contorno concentradas. A ausência de malha de domínio pelo MEC possibilita elevada precisão na avaliação dos campos mecânicos internos além de procedimentos menos complexos durante a propagação da fissura. No entanto, condições de contorno concentradas são frequentemente encontradas em mecânica da fratura. A versão Lagrangiana do MEC impõe tais condições de contorno aproximadamente por meio de elementos pequenos, os quais levam a instabilidades numéricas ou mesmo a representação inadequada do problema. Este estudo propõe uma formulação para a representação adequada de condições de contorno concentradas dentro da abordagem Lagrangiana do MEC. A mecânica da fratura não linear descreve os processos de falha material. A abordagem coesiva governa os processos não lineares de dissipação de energia, sendo os operadores constante e tangente responsáveis por solucionar o sistema não linear resultante. Três aplicações demonstram a precisão da formulação proposta, na qual as respostas do MEC são comparadas às soluções numéricas e experimentais disponíveis na literatura.

Palavras-chave:
condições de contorno concentradas; mecânica da fratura não linear; funções de enriquecimento; Modelagem da fratura tridimensional

1 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 [¹], [²]. 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 [³].

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. [⁴] and Bouchard et al. [⁵] in two-dimensional and Bremberg and Dhondt [⁶], Nejati et al. [⁷] and Wowk et al. [⁸] in three-dimensional LEFM modelling, among many others. Besides, Galvez et al. [⁹] and Dimitri et al. [¹⁰] addressed two-dimensional cohesive fracture problems whereas Gaedicke et al. [¹¹], Chen et al. [¹²] and Skar and Poulsen [¹³] 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 [¹⁴]–[¹⁶], for instance. In the nonlinear fracture mechanics field, it is worth mentioning the developments provided by Duarte et al. [¹⁷], Moes and Belytschko [¹⁸], Mohammadnejad and Khoei [¹⁹], Jaśkowiec et al. [²⁰] and Ferte et al. [²¹], 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 [²²], [²³] and Mi and Aliabadi [²⁴], [²⁵]. 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 [²⁶], multi-domain technique [²⁷], dipoles approach [²⁸] and cells with embedded discontinuities [²⁹], for instance. Nevertheless, the robust BEM approach in fracture mechanics is the dual BEM (DBEM) formulation [²⁴], [³⁰]. 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 [³¹], [³²] and three-dimensional problems [³³]–[³⁶] 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 [³⁷], [³⁸] 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 [³⁹], [⁴⁰] and anisotropic [⁴¹] 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 [⁴²] utilised the Betti’s reciprocity theorem and obtained a boundary integration equation with incorporated concentrated forces. Wang and Wang [⁴³] obtained similar developments accounting for spline approximation in the boundary integration equation. Recently, Zhou et al. [⁴⁴] 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 [⁴⁵] 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.

2 THE BOUNDARY ELEMENT METHOD

2.1 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 $\Omega$ composed of boundary $\Gamma$and nil body forces is as follows [²⁵]:

$c_{ij}\left({\mathit{x}}^s\right)u_j\left({\mathit{x}}^s\right)+{\int }_{\Gamma }{P}_{ij}^{\mathrm{*}}({\mathit{x}}^s,{\mathit{x}}^f)u_j\left({\mathit{x}}^f\right)d\Gamma ={\int }_{\Gamma }{U}_{ij}^{\mathrm{*}}({\mathit{x}}^s,{\mathit{x}}^f)p_j\left({\mathit{x}}^f\right)d\Gamma$ (1)

in which ${\mathit{x}}^s$ and ${\mathit{x}}^f$ are the source (collocation) and field points, respectively. The terms $P^{*}$ and $U^{*}$ indicate traction and displacement Kelvin fundamental solutions. $c_{ij}$ is the free term, which is equal to $\frac{{\delta }_{ij}}{2}$ for smooth boundaries. ${\delta }_{ij}$ is the Kronecker delta. $u$ and $p$ represent the displacements and tractions along the boundary. The fundamental solutions presented in Equation 1 are:

$U_{ij}^{\mathrm{*}}\left({\mathit{x}}^s,{\mathit{x}}^f\right)=\frac{1}{16\pi \mu \left(1-\upsilon \right)r}\left[\left(3-4\upsilon \right){\delta }_{ij}+r_{,i}{r}_{,j}\right]T_{ij}^{\mathrm{*}}\left({\mathit{x}}^s,{\mathit{x}}^f\right)=\frac{1}{8\pi (1-\upsilon )r^2}\left\{\frac{\partial r}{\partial \eta }\left[\right(1-2\upsilon ){\delta }_{ij}+3r_{,i}{r}_{,j}]-(1-2\upsilon )(r_{,i}{\eta }_j+r_{,j}{\eta }_i)\right\}$ (2)

where $r_i$ indicates the components of the distance vector and $r=‖r‖=‖{\mathit{x}}^f-{\mathit{x}}^s‖$ its norm; ${\eta }_i$ represents the components of the normal outward vector; $\mu$ and $\upsilon$ 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 [²⁴]:

$\frac{1}{2}{p}_j\left({\mathit{x}}^s\right)+{\eta }_i\left({\mathit{x}}^s\right){\int }_{\Gamma }{S}_{kij}^{\mathrm{*}}({\mathit{x}}^s,{\mathit{x}}^f)u_k\left({\mathit{x}}^f\right)d\Gamma ={\eta }_i\left({\mathit{x}}^s\right){\int }_{\Gamma }{D}_{kij}^{\mathrm{*}}({\mathit{x}}^s,{\mathit{x}}^f)p_k\left({\mathit{x}}^f\right)d\Gamma$ (3)

where $S^{*}$ and $D^{*}$ are the derivative kernels associated to $P^{*}$ and $U^{*}$, respectively. Then, the singularities in these kernels are of order $S^{*}=\frac{1}{r^3}$ and $D^{*}=\frac{1}{r^2}$. Because of the hyper singularity on kernel $S^{*}$, this integral term requires the evaluation into Hadamard Finite Part sense. For sake of completeness, the hyper singular kernels are as follows [²⁴], [²⁵]:

$D_{k\mathcal{l}j}^{\mathrm{*}}\left(\mathit{x},\widehat{\mathit{x}}\right)=\frac{1}{8\pi \left(1-\nu \right)r^2}\left[\left(1-2\nu \right)\left({\delta }_{k\mathcal{l}}{r}_{,j}+{\delta }_{jk}{r}_{,\mathcal{l}}-{\delta }_{\mathcal{l}j}{r}_{,k}\right)+3\left(r_{,\mathcal{l}}{r}_{,j}{r}_{,k}\right)\right]$

$S_{k\mathcal{l}j}^{\mathrm{*}}(\mathit{x},\widehat{\mathit{x}})=\frac{\mu }{4\pi (1-\nu )r^3}\left\{3\frac{\partial r}{\partial n}\right[(1-2\nu ){\delta }_{\mathcal{l}j}{r}_{,k}+\nu ({\delta }_{\mathcal{l}k}{r}_{,j}+{\delta }_{jk}{r}_{,\mathcal{l}})-5r_{,\mathcal{l}}{r}_{,j}{r}_{,k}]$ $+3\nu ({\eta }_{\mathcal{l}}{r}_{,j}{r}_{,k}+{\eta }_j{r}_{,\mathcal{l}}{r}_{,k})+(1-2\nu )(3{\eta }_k{r}_{,\mathcal{l}}{r}_{,j}+{\eta }_j{\delta }_{\mathcal{l}k}+{\eta }_{\mathcal{l}}{\delta }_{jk})-(1-4\nu ){\eta }_k{\delta }_{\mathcal{l}j}\}$ (4)

2.2 The DBEM collocation strategy and numerical assembling process

Collocation points at the same geometrical position lead to the singularities in the algebraic system of equations if Equation 1 or Equation 3 separately discretise the boundary. The singularities in the fundamental solutions explain such behaviour. Nevertheless, this problem can be avoided if Equation 1 and Equation 3 are coupled properly, which triggers the DBEM approach. 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.

Figure 1
- Discretisation in the DBEM

In the DBEM, the displacement BIE is written for collocations at the upper crack surface ${\Gamma }_c^{+}$ as follows:

$c_{ij}\left({\mathit{x}}^{s+}\right)u_j\left({\mathit{x}}^{s+}\right)+c_{ij}\left({\mathit{x}}^{s-}\right)u_j\left({\mathit{x}}^{s-}\right)+{\int }_{\Gamma }{P}_{ij}^{\mathrm{*}}({\mathit{x}}^{s+},{\mathit{x}}^f)u_j\left({\mathit{x}}^f\right)d\Gamma ={\int }_{\Gamma }{U}_{ij}^{\mathrm{*}}({\mathit{x}}^{s+},{\mathit{x}}^f)p_j\left({\mathit{x}}^f\right)d\Gamma$ (5)

Complementarily, the traction BIE is written at the lower crack surface ${\Gamma }_c^{-}$as follows:

$\frac{1}{2}{p}_j\left({\mathit{x}}^{s-}\right)-\frac{1}{2}{p}_j\left({\mathit{x}}^{s+}\right)+{\eta }_i\left({\mathit{x}}^{s-}\right){\int }_{\Gamma }{S}_{kij}^{\mathrm{*}}({\mathit{x}}^{s-},{\mathit{x}}^f)u_k\left({\mathit{x}}^f\right)d\Gamma ={\eta }_i\left({\mathit{x}}^{s-}\right){\int }_{\Gamma }{D}_{kij}^{\mathrm{*}}({\mathit{x}}^{s-},{\mathit{x}}^f)p_k\left({\mathit{x}}^f\right)d\Gamma$ (6)

It is worth mentioning that Equations 5 and 6 contain additional free terms, which account for the collocations at the same position. Besides, ${\mathit{x}}^{s+}$ and ${\mathit{x}}^{s-}$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 [²³]. The singularity subtraction technique proposed by Guiggiani and Gigante [⁴⁶] and Guiggiani et al. [⁴⁷] regularizes the kernels over singular-elements and enables the accurate assessment of singular and hypersingular BIE. It is worth mentioning that Cordeiro and Leonel [²³] 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., $\mathit{H}\mathit{u}=\mathit{G}\mathit{p}$. Thus, $\mathit{H}$ and $\mathit{G}$ are the DBEM influence matrices whereas the vectors $\mathit{u}$ and $\mathit{p}$ indicate the displacement and tractions at the collocation points, respectively. $\mathit{H}$ contains the integration of kernels $S^{*}$and $P^{*}$and $\mathit{G}$ contains the integration of kernels $D^{*}$ and $U^{*}$. 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 $\mathit{A}\mathit{x}=\mathit{b}$ can be solved and the unknowns at the entire boundary, $\mathit{x}$, determined.

3 THE REPRESENTATION OF CONCENTRATED BOUNDARY CONDITIONS ALONG THE BOUNDARY

3.1 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:

$p_j({\xi }_1,{\xi }_2)=N^{\alpha }({\xi }_1,{\xi }_2)p_j^{\alpha }+\Delta (\mathit{x}-{\mathit{x}}^l)F_j^l$ (7)

where $F_j^l$ indicates the concentrated force along the j direction applied at the $l$ point of the boundary mesh. $N$indicates the shape functions. It is worth emphasizing that $p\to \infty$ at ${\mathit{x}}^l$, which characterizes the concentrated force. Besides, the formulation enables the accurate representation of $NC$ concentrated forces applied at different positions along the boundary. $\Delta (\mathit{x}-{\mathit{x}}^l)$ 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:

${\int }_{\Gamma }{U}_{ij}^{\mathrm{*}}{p}_j\left({\mathit{x}}^f\right)d\Gamma ={\int }_{\Gamma }{U}_{ij}^{\mathrm{*}}\left[N_{\alpha }\right({\xi }_1,{\xi }_2)p_j^{\alpha }+\Delta (\mathit{x}-{\mathit{x}}^l\left)F_j^l\right]d\Gamma ={\int }_{\Gamma }{U}_{ij}^{\mathrm{*}}{N}_{\alpha }({\xi }_1,{\xi }_2)p_j^{\alpha }d\Gamma +{\int }_{\Gamma }{U}_{ij}^{\mathrm{*}}\Delta (\mathit{x}-{\mathit{x}}^l)F_j^{l}d\Gamma$ (8)

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 ${\mathit{x}}^s$ and ${\mathit{x}}^l$, as follows:

${\int }_{\Gamma }{U}_{ij}^{\mathrm{*}}\Delta (\mathit{x}-{\mathit{x}}^l)F_j^{l}d\Gamma =F_j^l{\int }_{\Gamma }{U}_{ij}^{\mathrm{*}}\Delta (\mathit{x}-{\mathit{x}}^l)d\Gamma =U_{ij}^{\mathrm{*}}({\mathit{x}}^s,{\mathit{x}}^l)F_j^l$ (9)

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:

$\left[\begin{array}{c}\mathit{H}\\ {\mathit{H}}^{+}\\ {\mathit{H}}^{-}\end{array}\right]\left\{\begin{array}{c}\mathit{u}\\ {\mathit{u}}^{+}\\ {\mathit{u}}^{-}\end{array}\right\}=\left[\begin{array}{c}\mathit{G}\\ {\mathit{G}}^{+}\\ {\mathit{G}}^{-}\end{array}\right]\left\{\begin{array}{c}\mathit{p}\\ {\mathit{p}}^{+}\\ {\mathit{p}}^{-}\end{array}\right\}+\left\{\begin{array}{c}{\mathit{b}}^X\\ {\mathit{b}}^{X^{+}}\\ {\mathit{b}}^{X^{-}}\end{array}\right\}$ (10)

where ${\mathit{b}}^X$ and ${\mathit{b}}^{X+}$ indicate the vector arising from the multiplication of the fundamental solution $U^{*}$ and the concentrated forces whereas ${\mathit{b}}^{X-}$ results from the multiplication of the fundamental solution $D^{*}$ and the concentrated forces at the hypersingular crack surface. The superscripts + and – indicate the mechanical variables along the upper and lower crack surfaces, respectively.

3.2 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:

${\bar{u}}_j^l=N^{\alpha }({\xi }_1^l,{\xi }_2^l)u_j^{\alpha }$ (11)

where ${\xi }_1$ and ${\xi }_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:

$\left[\begin{array}{cc}\mathit{A}& -{\mathit{G}}^{\mathit{S}}\\ {\mathit{A}}^{\mathit{X}}& 0\end{array}\right]\left\{\begin{array}{c}\mathit{x}\\ \mathit{R}\end{array}\right\}=\left\{\begin{array}{c}\mathit{b}\\ \bar{\mathit{u}}\end{array}\right\}$ (12)

in which $\mathit{A}$ results from the columns change procedure between $\mathit{H}$ and $\mathit{G}$ matrices. $A$ contains the influence terms associated to the unknowns at the boundary. $\mathit{x}$ represents the unknowns at the boundary whereas $b$ indicates the influence vector associated to the known/prescribed values at the boundary. $R$indicates the support reaction values. The sub-matrix ${\mathit{A}}^X$ is composed of shape function values evaluated from Equation 11, i.e., $A_{lj}^X=N_{\alpha \left(j,e\right)}\left({\xi }_1^l,{\xi }_2^l\right)$, in which $\alpha \left(j,e\right)$ is the element incidence associated to the global degree of freedom $j$ in the element $e$. $\bar{\mathit{u}}$ describes the concentrated displacement values prescribed at the boundary, which assumes nil values in the case of static supports. It is worth remarking that $\mathit{R}$ has similar meaning of $F$ in Equation 7. ${\mathit{G}}^{\mathit{S}}$ contains the evaluation of the fundamental solution at source point and support point, which is as follows:

${\mathit{G}}_{il}^A=\left[\begin{array}{c}{U}_{1k}^{\mathrm{*}}\left(x_i^s,x^l\right)\\ U_{2k}^{\mathrm{*}}\left(x_i^s,x^l\right)\\ U_{3k}^{\mathrm{*}}\left(x_i^s,x^l\right)\end{array}\right]$ (13)

where only the column $k$, associated to the direction of the prescribed displacement at the point $x^l$, has been accounted into the assembling system. The fundamental solutions $U^{*}$ and $D^{*}$ 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 ${\xi }_1$ and ${\xi }_2$ within the boundary element domain. For sake of completeness, ^{Appendix 1} presents a flowchart, in which this procedure is described clearly.

4 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 quasi-brittle 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 $p_{\eta }$ 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.

Figure 2
- Cohesive crack approach. Cohesive tractions.

Several cohesive laws have been proposed in the literature to govern the cohesive traction behaviour [⁴⁸]–[⁵⁰]. Three of them have been applied often in this problem: linear, bi-linear and exponential, as illustrated in Figure 3.

Figure 3
- Cohesive laws: (a) linear, (b) bi-linear and (c) exponential

The linear cohesive law associates $p_{\eta }$ and $COD$ by a linear function. This law requires the material tensile strength $f_T$ and the threshold opening displacement $CO{D}_c$ values. Then, this law is as follows:

$p_{\eta }\left(COD\right)=\left\{\begin{array}{cc}{f}_T\left(1-\frac{COD}{CO{D}_C}\right),& \mathrm{}if\mathrm{}0\le COD\le CO{D}_c\\ 0& if\mathrm{}COD>CO{D}_c\end{array}\right.$ (14)

The bi-linear function leads to the bi-linear cohesive law, which is as follows:

$p_{\eta }\left(COD\right)=\left\{\begin{array}{cc}{f}_T-\left(\frac{f_T-f_T^{\mathrm{*}}}{CO{D}^{\mathrm{*}}}\right)COD& if\mathrm{}0\le COD\le CO{D}^{\mathrm{*}}\\ \frac{f_T^{\mathrm{*}}COD}{CO{D}^{\mathrm{*}}-CO{D}_c}+f_T^{\mathrm{*}}\left(1-\frac{CO{D}^{\mathrm{*}}}{CO{D}^{\mathrm{*}}-CO{D}_c}\right)& if\mathrm{}CO{D}^{\mathrm{*}}\le COD\le CO{D}_c\\ 0& if\mathrm{}COD>CO{D}_c\end{array}\right.$ (15)

where the additional parameters $f_T^{*}$, $CO{D}^{*}$ and $CO{D}_c$ have been presented by Petersson [⁵¹] accounting for concrete:

$f_T^{\mathrm{*}}=\frac{f_T}{3}CO{D}^{\mathrm{*}}=0.8\frac{G_{\sigma }}{f_T}CO{D}_c=3.6\frac{G_{\sigma }}{f_T}$ (16)

where $G_{\sigma }$ indicates the fracture energy.

Finally, the exponential cohesive law associates $p_{\eta }$ and $COD$ as follows:

$p_{\eta }\left(COD\right)=f_T\exp \left(-\frac{f_T}{G_{\sigma }}COD\right)$ (17)

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.

5 NONLINEAR SOLUTION TECHNIQUE

5.1 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:

${\mathit{H}}^b{\mathit{u}}^b+{\mathit{H}}^R{\mathit{u}}^R+{\mathit{H}}^L{\mathit{u}}^L={\mathit{G}}^b{\mathit{p}}^b+{\mathit{G}}^R{\mathit{p}}^R+{\mathit{G}}^L{\mathit{p}}^L$ (18)

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:

${\mathit{H}}^b{\mathit{u}}^b+{\mathit{H}}_{\eta }^R{\mathit{u}}_{\eta }^R+{\mathit{H}}_l^R{\mathit{u}}_l^R+{\mathit{H}}_t^R{\mathit{u}}_t^R+{\mathit{H}}_{\eta }^L{\mathit{u}}_{\eta }^L+{\mathit{H}}_l^L{\mathit{u}}_l^L+{\mathit{H}}_t^L{\mathit{u}}_t^L={\mathit{G}}^b{\mathit{p}}^b+{\mathit{G}}_{\eta }^R{\mathit{p}}_{\eta }^R+{\mathit{G}}_l^R{\mathit{p}}_l^R+{\mathit{G}}_t^R{\mathit{p}}_t^R+{\mathit{G}}_{\eta }^L{\mathit{p}}_{\eta }^L+{\mathit{G}}_l^L{\mathit{p}}_l^L+{\mathit{G}}_t^L{\mathit{p}}_t^L$ (19)

where ${\mathit{H}}_{\eta }^R$, ${\mathit{H}}_l^R$, ${\mathit{H}}_t^R$, ${\mathit{H}}_{\eta }^L$, ${\mathit{H}}_l^L$, ${\mathit{H}}_t^L$, ${\mathit{G}}_{\eta }^R$, ${\mathit{G}}_l^R$, ${\mathit{G}}_t^R$, ${\mathit{G}}_{\eta }^L$, ${\mathit{G}}_l^L$ and ${\mathit{G}}_t^L$ indicate the product between the rotation matrix and the corresponding influence matrices. $l$ and $t$ indicate tangential crack surface directions whereas $\eta$ refers to the normal crack surface direction.

Figure 4
- Local coordinate system for nonlinear fracture analysis

Equation 19 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:

$COD=-u_{\eta }^L-u_{\eta }^{R}CSD=u_l^R+u_l^{L}CTD=u_t^R-u_t^L$ (20)

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.,$p_{\eta }^R=p_{\eta }^L,p_l^R=p_l^L\text{and}{p}_t^R=p_t^L$. In addition, the cohesive approach predicts $p_{\eta }$ as a function of COD. Finally, one enforces the boundary conditions at the external boundary. All those algebraic manipulations lead to the following:

${\mathit{A}}_b{\mathit{x}}_b+\left({\mathit{H}}_{\eta }^R-{\mathit{H}}_{\eta }^L\right){\mathit{u}}_{\eta }^R+\left({\mathit{H}}_l^R-{\mathit{H}}_l^L\right){\mathit{u}}_l^R+\left({\mathit{H}}_t^R+{\mathit{H}}_t^L\right){\mathit{u}}_t^R-{\mathit{H}}_{\eta }^L\mathit{C}\mathit{O}\mathit{D}+{\mathit{H}}_l^L\mathit{C}\mathit{S}\mathit{D}-{\mathit{H}}_t^L\mathit{C}\mathit{T}\mathit{D}={\mathit{f}}_b+\left({\mathit{G}}_l^R+{\mathit{G}}_l^L\right){\mathit{p}}_l+({\mathit{G}}_t^R+{\mathit{G}}_t^L){\mathit{p}}_t+({\mathit{G}}_{\eta }^R+{\mathit{G}}_{\eta }^L){\mathit{p}}_{\eta }\left(\mathit{C}\mathit{O}\mathit{D}\right)$ (21)

in which matrix ${\mathit{A}}_b$multiplies the unknown fields ${\mathit{x}}_b$ at the external boundary, obtained from the columns change procedure of BEM, i.e., by enforcing the boundary conditions. Besides, ${\mathit{f}}_b$ results from the product between the known values at the boundary and the respective influence values.

The cohesive law governs$p_{\eta }\left(COD\right)$, 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, $p_{\eta }\left(COD\right)$ 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 $\Delta p_{\eta }^{exc}$. Then, the $p_{\eta }\left(COD\right)$ updates as follows:

$p_{\eta }\left(COD\right)=p_{\eta }+\Delta p_{\eta }^{exc}$ (22)

The nonlinear solution process converges when the norm of $\Delta p_{\eta }^{exc}$ 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.

5.2 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 [²⁷]. In this approach, the first step refers to the elastic prediction, which provides the increments $\delta$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:

$$\begin{array}{l}{A}_b\delta x_b{}^{(i)}+\left(H_{\eta }^R-H_{\eta }^L\right)\delta u_{\eta }^{R^R}{}^{(i)}+\left(H_l^R-H_l^L\right)\delta u_l^{R^{(i)}}\\ +\left(H_t^R+H_t^L\right)\delta u_t^{R^{(i)}}-H_{\eta }^L\delta CO{D}^{(i)}-H_l^L\delta CS{D}^{(i)}+H_t^L\delta CT{D}^{(i)}=\delta f_b+\left(G_l^R+G_l^L\right)\text{Δ}{p}_l^{exc(i)}+\left(G_t^R+G_t^L\right)\text{Δ}{p}_t^{exc(i)}\\ +\left(G_{d\eta }^f+G_{e\eta }^f\right)\text{Δ}{p}_{\eta }^{exc(i)}\end{array}$$ (23)

where the increments on the mechanical fields are calculated and added to the responses from the previous iterations. Besides, the search for $\Delta p_{\eta }^{exc}=0$ occurs simultaneously for all collocations at the cohesive interface. Then, the convergence occurs when $‖\Delta {\mathit{p}}_{\eta }^{exc}‖$ is smaller than a prescribed tolerance, which indicates that the non-equilibrated traction vector is sufficiently small.

5.3 Solution Technique by Tangent Operator (TO)

The solution technique via TO requires the definition of a residue $\mathit{R}$ from the algebraic equations presented in Equation 21 [²⁷]. Such residue accounts for the problem variables as follows:

$\mathit{R}\left(\mathit{X}\right)={\mathit{A}}_b{\mathit{x}}_b+({\mathit{H}}_{\eta }^R-{\mathit{H}}_{\eta }^L){\mathit{u}}_{\eta }^R+({\mathit{H}}_l^R-{\mathit{H}}_l^L){\mathit{u}}_l^R+({\mathit{H}}_t^R+{\mathit{H}}_t^L){\mathit{u}}_t^R-{\mathit{H}}_{\eta }^L\mathit{C}\mathit{O}\mathit{D}+{\mathit{H}}_l^L\mathit{C}\mathit{S}\mathit{D}-{\mathit{H}}_t^L\mathit{C}\mathit{T}\mathit{D}-{\mathit{f}}_b-({\mathit{G}}_l^R+{\mathit{G}}_l^L){\mathit{p}}_l-({\mathit{G}}_t^R+{\mathit{G}}_t^L){\mathit{p}}_t-({\mathit{G}}_{\eta }^R+{\mathit{G}}_{\eta }^L){\mathit{p}}_{\eta }\left(\mathit{C}\mathit{O}\mathit{D}\right)$ (24)

where $\mathit{X}$ 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:

$\mathit{R}\left({\mathit{X}}^{(i-1)}\right)=-({\mathit{G}}_{\eta }^R+{\mathit{G}}_{\eta }^L)\Delta {{\mathit{p}}_{\eta }^{exc}}^{\left(i\right)}$ (25)

The TO determines the increment set $\delta \mathit{X}$ leading to the nil residue $\mathit{R}\left(\mathit{X}+\delta \mathit{X}\right)$value. For this purpose, one assumes the residue function as continuous. Therefore, this function can be expanded in a linear Taylor’s expansion around $\mathit{X}$. Finally, one truncates this expansion to the first term, which provides the corrections $\delta \mathit{X}$ as follows:

$\mathit{R}(\mathit{X}+\delta \mathit{X})=0\Rightarrow \mathit{R}\left({\mathit{X}}^{(i-1)}\right)+\frac{\partial \mathit{R}}{\partial \mathit{X}}\delta {\mathit{X}}^{\left(i\right)}=0\Rightarrow {\mathit{A}}_b\delta {\mathit{x}}_b^{\left(i\right)}+({\mathit{H}}_{\eta }^R-{\mathit{H}}_{\eta }^L)\delta {\mathit{u}}_{\eta }^R+({\mathit{H}}_l^R-{\mathit{H}}_l^L)\delta {\mathit{u}}_l^R+({\mathit{H}}_t^R+{\mathit{H}}_t^L)\delta {\mathit{u}}_t^R-({\mathit{H}}_{\eta }^L-({\mathit{G}}_{\eta }^R+{\mathit{G}}_{\eta }^L\left)\frac{\partial p_{\eta }}{\partial COD}\right)\delta \mathit{C}\mathit{O}\mathit{D}+{\mathit{H}}_l^L\delta \mathit{C}\mathit{S}\mathit{D}-{\mathit{H}}_t^L\delta \mathit{C}\mathit{T}\mathit{D}=\delta {\mathit{f}}_b+({\mathit{G}}_l^R+{\mathit{G}}_l^L)\delta {\mathit{p}}_l+({\mathit{G}}_t^R+{\mathit{G}}_t^L)\delta {\mathit{p}}_t+({\mathit{G}}_{\eta }^R+{\mathit{G}}_{\eta }^L)\Delta {\mathit{p}}_{\eta }^{exc\left(i\right)}$ (26)

where $\frac{\partial p_{\eta }}{\partial COD}$ 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, $\frac{\partial p_{\eta }}{\partial COD}$ 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 $‖\Delta {\mathit{p}}_{\eta }^{exc}‖$ is smaller than a prescribed tolerance.

6 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 four-point 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.

6.1 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. The numerical responses provided by the BEM formulations are compared against the analytical solution presented in Gross et al. [⁵²] and Tada et al. [⁵³] for LEFM condition, in addition to the responses of an equivalent model in FEM/ANSYS. Concentrated supports at (5.0; 25.0; 1.25), (5.0; 25.0; 2.5) and (5.0; 25.0; 3.75) avoid rigid body motion in this application. The material properties are: Young’s modulus 1000 (stress unity) and Poisson ratio 0.3. Furthermore, four different meshes have been utilised in this mechanical modelling as illustrated in 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.

Figure 5
- Notched prismatic solid and boundary conditions. Dimensions in length unity and Force in force unity. Geometry, DBEM and FEM meshes

Gross et al. [⁵²] and Tada et al. [⁵³] present the analytical plane solution assuming LEFM behaviour, in which the crack mouth opening displacement (CMOD) can be determined as follows:

$CMOD=\frac{4\sigma a}{E\mathrm{\text{'}}}{V}_1\left(\frac{a}{b}\right)V_1\left(\frac{a}{b}\right)=\frac{1.46+3.42\left[1-\mathit{cos}\left(\frac{\pi a}{2b}\right)\right]}{{\mathit{cos}}^2\left(\frac{\pi a}{2b}\right)}E\mathrm{\text{'}}=\frac{E}{1-{\nu }^2}\text{for plane strain and}E\mathrm{\text{'}}=E\text{for plane stress}$ (27)

where $\sigma$ is the applied normal traction, $a$ is the crack length, $b$ is the solid width, $E$ is the Young’s modulus and $\nu$ is the Poisson ratio. In this application, $a=2.5$ and $b=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.

Figure 6
– Analytical solution geometry and boundary conditions [⁵²], [⁵³]

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 $x_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 7
– Deformed shape and displacements along $x_2$ direction: Mesh D and FEM/ANSYS^®

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.

Figure 8
- CMOD (displacement unity) results for each mesh and reference: concentrated force application

Figure 9
- CMOD percentual error (0.0; 25.0; 2.5).

Figure 10 illustrates the Von-Mises stress values evaluated by FEM/ANSYS and BEM for an internal line at the following position: $x_1=3.75$, $0\le x_2\le 50.0$ and $x_3=2.5$. As expected, stresses concentrations have been observed along the structure ends, in which the BEM represented properly the mechanical effects triggered by the concentrated forces. Besides, one emphasizes the complex stress behaviour at the coordinate $x_2=25.0$ caused the crack tip, which has been identified properly by the BEM.

Figure 10
– Von-Mises stress values for internal line $x_1=3.75$, $0\le x_2\le 50.0$, $x_3=2.5$

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.

Table 1
– Prescribed displacements and reaction values for concentrated supports case

Mesh	Prescribed displacement	Reaction value
A	0.03037	25.00155
B	0.03559	25.00233
C	0.04008	25.00273
D	0.04434	25.00297

Figure 11
– Comparison between CMOD for concentrated force and concentrated displacement at structural ends.

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.

6.2 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. [⁵⁴] analysed this specimen experimentally whereas Saleh and Aliabadi [⁵⁵] 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. [⁵⁴]: Young’s modulus$E=27GPa$, Poisson ratio $\nu =0.1$, tensile material strength $f_T=2MPa$ and fracture energy $G_{\sigma }=100\frac{N}{m}$.

Figure 12
- Four-point bending test with two initial notches

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 [⁵⁴] and the two-dimensional numerical modelling [⁵⁵], in which good agreement has been observed among the predictions of the proposed formulation and that provided by [⁵⁴].

Figure 13
- Boundary mesh

Figure 14
– Crack geometry. present study, experimental [⁵⁴] and numerical [⁵⁵]

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 $x_1$ and $x_2$ are nil in the following ranges $x_1=0.32;x_2=0;0\le x_3\le 0.2m$ and $x_1=0.8;x_2=0;0\le x_3\le 0.2m$. Besides, the points (0.32; 0.0; 0.0)m and (0.8; 0.0; 0.0)m have concentrated supports along $x_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}kPa$, in terms of $‖\Delta {\mathit{p}}_{\eta }^{exc}‖$. 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 [⁵⁵] 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 $P_1$ and $P_2$ at the supports to the prescribed displacement to geometrical and material parameters as follows:

$$F_{\text{adim}}=\frac{P_1+P_2}{f_{T}dt}$$ (28)

in which $f_T$ is the tensile material strength and $d$ and $t$ 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.

Figure 15
- Dimensionless load x deflection curves

As expected, the numerical responses demonstrate good agreement at the elastic part of the curve. Besides, the BEM results strongly agree to the reference [⁵⁵] 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.

Table 2
- Comparative of iterations. Required number of iterations

Cohesive Law	Concentrated supports			Distributed supports
	Constant Operator	Tangent Operator	Reduction	Constant Operator	Tangent Operator		Reduction
Linear	366	42	88.5%	419	49		88.3%
Bi-linear	282	41	85.5%	270	45		83.3%
Exponential	323	53	83.6%	336	61		81.8%

For sake of clarity, Figure 16 illustrates the displacement values along the $x_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.

Figure 16
- Deformed shape and displacements at $x_2$ direction

Finally, Figure 17 illustrates the evolution of the mechanical degradation processes at the FPZ. The numbers in this figure indicate the residual resistance at the FPZ. Then, zero indicates free traction crack surfaces and one refers to the continuous material. The results in this figure account for the bilinear cohesive law and CO solution techniques at load increments 14 and 15. It is worth remarking the significant mechanical degradation observed between these load increments.

Figure 17
- Residual material resistance at the FPZ.

6.3 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. [⁵⁶]. Besides, many authors [²¹], [⁴⁸], [⁴⁹] 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 [⁵⁶] are: Young’s modulus of concrete $E_c=35GPa$, Poisson ratio of concrete ${\nu }_c=0.2$, Young’s modulus of steel frame $E_s=200GPa$, Poisson ratio of steel frame ${\nu }_s=0.3$, energy fracture $G_{\sigma }=120\frac{N}{m}$ and tensile material strength $f_T=1.6MPa$.

Figure 18
- Brokenshire torsion test: prismatic solid.

The supports positioned along $x_3=0$ restrict displacements along directions$x_1$, $x_2$ and $x_3$. The other supports in this modelling restrict the displacement solely along the $x_2$ direction. Besides, the prescribed displacement is 0.8mm along $x_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}kPa$, in terms of $‖\Delta p_{\eta }^{exc}‖$. 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. [⁵⁶] accounts for concentrated boundary conditions.

Figure 19
- Initial boundary element mesh. (a) Complete model (b) concrete specimen (c) Steel frame

The numerical responses obtained by the proposed BEM formulation are compared to experimental results presented by Jefferson et al. [⁵⁶], 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.

Figure 20
- Load x CMOD curves

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.

Table 3
- Comparative of iterations. Required number of iterations

Cohesive Law	Constant Operator	Tangent Operator	Reduction
Linear	175	68	61.1%
Bi-linear	192	68	64.6%
Exponential	189	79	58.2%

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. [⁵⁶] and the XFEM results obtained by Kaczmarczyk et al. [⁵⁷] 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.

Figure 21
- Residual material resistance at the FPZ

Figure 22
- Comparison of crack propagation paths obtained by different approaches and the proposed BEM formulation

7 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.

APPENDIX 1 SEARCH FOR DIMENSIONLESS COORDINATES

Because concentrated boundary conditions have been applied in global coordinates, the search for the correspondent element dimensionless coordinates can be performed by a Newton-Raphson scheme. The flowchart below describes the procedure utilised in the proposed formulation.

ACKNOWLEDGEMENTS

This research project was sponsored by grant #2019/18795-6, São Paulo Research Foundation (FAPESP).

• Financial support: Fundação de Amparo e Pesquisa do Estado de São Paulo - FAPESP.
• Data Availability: Data associated with this paper is not available.
• How to cite: M. Rocha, and E. D. Leonel, “The fracture failure of three-dimensional concrete structures subjected to concentrated loadings using the Boundary Element Method,” Rev. IBRACON Estrut. Mater., vol. 15, no. 5, e15504, 2022, https://doi.org/10.1590/S1983-41952022000500004

REFERENCES

• ¹
  B. Cotterell, "The past, present, and future of fracture mechanics," Eng. Fract. Mech., vol. 69, no. 5, pp. 533–553, 2002, http://dx.doi.org/10.1016/S0013-7944(01)00101-1
  » http://dx.doi.org/10.1016/S0013-7944(01)00101-1
• ²
  L. Marsavina, F. Berto, R. Negru, D. Serban, and E. Linul, "An engineering approach to predict mixed mode fracture of PUR foams based on ASED and micromechanical modelling," Theor. Appl. Fract. Mech., vol. 91, pp. 148–154, 2017, http://dx.doi.org/10.1016/j.tafmec.2017.06.008
  » http://dx.doi.org/10.1016/j.tafmec.2017.06.008
• ³
  F. Yan, X. Feng, J. Lv, P. Pan, and S. Li, "Continuous-discontinuous cellular automaton method for cohesive crack growth in rock," Eng. Fract. Mech., vol. 188, pp. 361–380, 2018, http://dx.doi.org/10.1016/j.engfracmech.2017.09.007
  » http://dx.doi.org/10.1016/j.engfracmech.2017.09.007
• ⁴
  T. Bittencourt, P. Wawrzynek, A. Ingraffea, and J. Sousa, "Quasi-automatic simulation of crack propagation for 2D LEFM problems," Eng. Fract. Mech., vol. 55, no. 2, pp. 321–334, 1996, http://dx.doi.org/10.1016/0013-7944(95)00247-2
  » http://dx.doi.org/10.1016/0013-7944(95)00247-2
• ⁵
  P. Bouchard, F. Bay, and Y. Chastel, "Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria," Comput. Methods Appl. Mech. Eng., vol. 192, no. 35-36, pp. 3887–3908, 2003, http://dx.doi.org/10.1016/s0045-7825(03)00391-8
  » http://dx.doi.org/10.1016/s0045-7825(03)00391-8
• ⁶
  D. Bremberg and G. Dhondt, "Automatic crack-insertion for arbitrary crack growth," Eng. Fract. Mech., vol. 75, no. 3-4, pp. 404–416, 2008, http://dx.doi.org/10.1016/j.engfracmech.2007.01.003
  » http://dx.doi.org/10.1016/j.engfracmech.2007.01.003
• ⁷
  M. Nejati, A. Paluszny, and R. Zimmerman, "On the use of quarter-point tetrahedral finite elements in linear elastic fracture mechanics," Eng. Fract. Mech., vol. 144, pp. 194–221, 2015, http://dx.doi.org/10.1016/j.engfracmech.2015.06.055
  » http://dx.doi.org/10.1016/j.engfracmech.2015.06.055
• ⁸
  D. Wowk, L. Alousis, and P. Underhill, "An adaptive remeshing technique for predicting the growth of irregular crack fronts using p-version finite element analysis," Eng. Fract. Mech., vol. 207, pp. 36–47, 2019, http://dx.doi.org/10.1016/j.engfracmech.2018.12.002
  » http://dx.doi.org/10.1016/j.engfracmech.2018.12.002
• ⁹
  J. C. Gálvez, J. Planas, J. M. Sancho, E. Reyes, D. A. Cendón, and M. J. Casati, "An embedded cohesive crack model for finite element analysis of quasi-brittle materials," Eng. Fract. Mech., vol. 109, pp. 369–386, 2013.
• ¹⁰
  R. Dimitri, P. Cornetti, V. Mantič, M. Trullo, and L. De Lorenzis, "Mode-I DEBONDING of a double cantilever BEAM: A comparison BETWEEN cohesive CRACK modeling and Finite Fracture Mechanics," Int. J. Solids Struct., vol. 124, pp. 57–72, 2017.
• ¹¹
  C. Gaedicke, J. Roesler, and F. Evangelista, “Three-dimensional cohesive crack model prediction of the flexural capacity of concrete slabs on soil,” Eng. Fract. Mech., vol. 94, pp. 1–12. May 2012. Accessed: Sep. 15, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0013794412001749
  » https://www.sciencedirect.com/science/article/pii/S0013794412001749
• ¹²
  X. Chen, X. Deng, and M. Sutton, "Simulation of stable tearing crack growth events using the cohesive zone model approach," Eng. Fract. Mech., vol. 99, pp. 223–238, 2013, http://dx.doi.org/10.1016/j.engfracmech.2012.12.017
  » http://dx.doi.org/10.1016/j.engfracmech.2012.12.017
• ¹³
  A. Skar and P. Poulsen, "3-D cohesive finite element model for application in structural analysis of heavy duty composite pavements," Constr. Build. Mater., vol. 101, pp. 417–431, 2015, http://dx.doi.org/10.1016/j.conbuildmat.2015.10.052
  » http://dx.doi.org/10.1016/j.conbuildmat.2015.10.052
• ¹⁴
  K. Agathos, E. Chatzi, and S. Bordas, "Multiple crack detection in 3D using a stable XFEM and global optimization," Comput. Mech., vol. 62, no. 4, pp. 835–852, 2018, http://dx.doi.org/10.1007/s00466-017-1532-y
  » http://dx.doi.org/10.1007/s00466-017-1532-y
• ¹⁵
  K. Agathos, S. Bordas, and E. Chatzi, "Improving the conditioning of XFEM/GFEM for fracture mechanics problems through enrichment quasi-orthogonalization," Comput. Methods Appl. Mech. Eng., vol. 346, pp. 1051–1073, 2019., http://dx.doi.org/10.1016/j.cma.2018.08.007
  » http://dx.doi.org/10.1016/j.cma.2018.08.007
• ¹⁶
  R. Tian, L. Wen, and L. Wang, "Three-dimensional improved XFEM (IXFEM) for static crack problems," Comput. Methods Appl. Mech. Eng., vol. 343, pp. 339–367, 2019, http://dx.doi.org/10.1016/j.cma.2018.08.029
  » http://dx.doi.org/10.1016/j.cma.2018.08.029
• ¹⁷
  C. Duarte, O. Hamzeh, T. Liszka, and W. Tworzydlo, "A generalized finite element method for the simulation of three-dimensional dynamic crack propagation," Comput. Methods Appl. Mech. Eng., vol. 190, no. 15-17, pp. 2227–2262, 2001, http://dx.doi.org/10.1016/s0045-7825(00)00233-4
  » http://dx.doi.org/10.1016/s0045-7825(00)00233-4
• ¹⁸
  N. Moës and T. Belytschko, "Extended finite element method for cohesive crack growth," Eng. Fract. Mech., vol. 69, no. 7, pp. 813–833, 2002, http://dx.doi.org/10.1016/s0013-7944(01)00128-x
  » http://dx.doi.org/10.1016/s0013-7944(01)00128-x
• ¹⁹
  T. Mohammadnejad and A. Khoei, "An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model," Finite Elem. Anal. Des., vol. 73, pp. 77–95, 2013, http://dx.doi.org/10.1016/j.finel.2013.05.005
  » http://dx.doi.org/10.1016/j.finel.2013.05.005
• ²⁰
  J. Jaśkowiec and F. van der Meer, "A consistent iterative scheme for 2D and 3D cohesive crack analysis in XFEM," Comput. Struc., vol. 136, pp. 98–107, 2014, http://dx.doi.org/10.1016/j.compstruc.2014.01.029
  » http://dx.doi.org/10.1016/j.compstruc.2014.01.029
• ²¹
  G. Ferté, P. Massin, and N. Moës, "3D crack propagation with cohesive elements in the extended finite element method," Comput. Methods Appl. Mech. Eng., vol. 300, pp. 347–374, 2016, http://dx.doi.org/10.1016/j.cma.2015.11.018
  » http://dx.doi.org/10.1016/j.cma.2015.11.018
• ²²
  S. Cordeiro and E. Leonel, "Mechanical modelling of three-dimensional cracked structural components using the isogeometric dual boundary element method," Appl. Math. Model., vol. 63, pp. 415–444, 2018, http://dx.doi.org/10.1016/j.apm.2018.06.042
  » http://dx.doi.org/10.1016/j.apm.2018.06.042
• ²³
  S. Cordeiro and E. Leonel, "An improved computational framework based on the dual boundary element method for three-dimensional mixed-mode crack propagation analyses," Adv. Eng. Softw., vol. 135, pp. 102689, 2019, http://dx.doi.org/10.1016/j.advengsoft.2019.102689
  » http://dx.doi.org/10.1016/j.advengsoft.2019.102689
• ²⁴
  Y. Mi and M. Aliabadi, "Dual boundary element method for three-dimensional fracture mechanics analysis," Eng. Anal. Bound. Elem., vol. 10, no. 2, pp. 161–171, 1992, http://dx.doi.org/10.1016/0955-7997(92)90047-b
  » http://dx.doi.org/10.1016/0955-7997(92)90047-b
• ²⁵
  Y. Mi and M. Aliabadi, "Three-dimensional crack growth simulation using BEM," Comput. Struc., vol. 52, no. 5, pp. 871–878, 1994, http://dx.doi.org/10.1016/0045-7949(94)90072-8
  » http://dx.doi.org/10.1016/0045-7949(94)90072-8
• ²⁶
  B. Yang and K. Ravi-Chandar “A single-domain dual-boundary-element formulation incorporating a cohesive zone model for elastostatic cracks,” in Recent Advances in Fracture Mechanics, W. G. Knauss, R. A. Schapery, Eds., Dordrecht: Springer, 1998, pp. 115–144.
• ²⁷
  S. Cordeiro and E. Leonel, "Cohesive crack propagation modelling in wood structures using BEM and the Tangent Operator Technique," Eng. Anal. Bound. Elem., vol. 64, pp. 111–121, 2016, http://dx.doi.org/10.1016/j.enganabound.2015.11.013
  » http://dx.doi.org/10.1016/j.enganabound.2015.11.013
• ²⁸
  L. Almeida, E. Lima Jr., and J. Barbirato, "Cohesive crack propagation analysis using a dipole BEM formulation with tangent operator," Theor. Appl. Fract. Mech., vol. 109, pp. 102765, 2020, http://dx.doi.org/10.1016/j.tafmec.2020.102765
  » http://dx.doi.org/10.1016/j.tafmec.2020.102765
• ²⁹
  R. Peixoto, G. Ribeiro, and R. Pitangueira, "A boundary element method formulation for quasi-brittle material fracture analysis using the continuum strong discontinuity approach," Eng. Fract. Mech., vol. 202, pp. 47–74, 2018, http://dx.doi.org/10.1016/j.engfracmech.2018.09.012
  » http://dx.doi.org/10.1016/j.engfracmech.2018.09.012
• ³⁰
  H. Hong and J. Chen, "Derivations of Integral Equations of Elasticity," J. Eng. Mech., vol. 114, no. 6, pp. 1028–1044, 1988, http://dx.doi.org/10.1061/(asce)0733-9399(1988)114:6(1028)
  » http://dx.doi.org/10.1061/(asce)0733-9399(1988)114:6(1028)
• ³¹
  R. Price and J. Trevelyan, "Boundary element simulation of fatigue crack growth in multi-site damage," Eng. Anal. Bound. Elem., vol. 43, pp. 67–75, 2014, http://dx.doi.org/10.1016/j.enganabound.2014.03.002
  » http://dx.doi.org/10.1016/j.enganabound.2014.03.002
• ³²
  H. C. Andrade and E. Leonel, "The multiple fatigue crack propagation modelling in nonhomogeneous structures using the DBEM," Eng. Anal. Bound. Elem., vol. 98, pp. 296–309, 2019, http://dx.doi.org/10.1016/j.enganabound.2018.10.026
  » http://dx.doi.org/10.1016/j.enganabound.2018.10.026
• ³³
  I. Benedetti, M. Aliabadi, and G. Davì, "A fast 3D dual boundary element method based on hierarchical matrices," Int. J. Solids Struct., vol. 45, no. 7-8, pp. 2355–2376, 2008, http://dx.doi.org/10.1016/j.ijsolstr.2007.11.018
  » http://dx.doi.org/10.1016/j.ijsolstr.2007.11.018
• ³⁴
  R. Citarella and F. Buchholz, "Comparison of crack growth simulation by DBEM and FEM for SEN-specimens undergoing torsion or bending loading," Eng. Fract. Mech., vol. 75, no. 3-4, pp. 489–509, 2008, http://dx.doi.org/10.1016/j.engfracmech.2007.03.039
  » http://dx.doi.org/10.1016/j.engfracmech.2007.03.039
• ³⁵
  P. Carlone, R. Citarella, M. Sonne, and J. Hattel, "Multiple crack growth prediction in AA2024-T3 friction stir welded joints, including manufacturing effects," Int. J. Fatigue, vol. 90, pp. 69–77, 2016, http://dx.doi.org/10.1016/j.ijfatigue.2016.04.004
  » http://dx.doi.org/10.1016/j.ijfatigue.2016.04.004
• ³⁶
  R. Citarella, P. Carlone, R. Sepe, and M. Lepore, "DBEM crack propagation in friction stir welded aluminum joints," Adv. Eng. Softw., vol. 101, pp. 50–59, 2016, http://dx.doi.org/10.1016/j.advengsoft.2015.12.002
  » http://dx.doi.org/10.1016/j.advengsoft.2015.12.002
• ³⁷
  R. Simpson and J. Trevelyan, "A partition of unity enriched dual boundary element method for accurate computations in fracture mechanics," Comput. Methods Appl. Mech. Eng., vol. 200, no. 1-4, pp. 1–10, 2011, http://dx.doi.org/10.1016/j.cma.2010.06.015
  » http://dx.doi.org/10.1016/j.cma.2010.06.015
• ³⁸
  R. Simpson and J. Trevelyan, "Evaluation of J1 and J2 integrals for curved cracks using an enriched boundary element method," Eng. Fract. Mech., vol. 78, no. 4, pp. 623–637, 2011, http://dx.doi.org/10.1016/j.engfracmech.2010.12.006
  » http://dx.doi.org/10.1016/j.engfracmech.2010.12.006
• ³⁹
  X. Peng, E. Atroshchenko, P. Kerfriden, and S. Bordas, "Linear elastic fracture simulation directly from CAD: 2D NURBS-based implementation and role of tip enrichment," Int. J. Fract., vol. 204, no. 1, pp. 55–78, 2016, http://dx.doi.org/10.1007/s10704-016-0153-3
  » http://dx.doi.org/10.1007/s10704-016-0153-3
• ⁴⁰
  H. Andrade and E. Leonel, "An enriched dual boundary element method formulation for linear elastic crack propagation," Eng. Anal. Bound. Elem., vol. 121, pp. 158–179, 2020, http://dx.doi.org/10.1016/j.enganabound.2020.09.007
  » http://dx.doi.org/10.1016/j.enganabound.2020.09.007
• ⁴¹
  G. Hattori, I. Alatawi, and J. Trevelyan, "An extended boundary element method formulation for the direct calculation of the stress intensity factors in fully anisotropic materials," Int. J. Numer. Methods Eng., vol. 109, no. 7, pp. 965–981, 2016, http://dx.doi.org/10.1002/nme.5311
  » http://dx.doi.org/10.1002/nme.5311
• ⁴²
  J. Qian, "Boundary element method with concentrated force," Eng. Fract. Mech., vol. 7, no. 2, pp. 140–144, 1990.
• ⁴³
  Z. Wang and Y. Wang, "Spline boundary element method for the analysis of elastic body subjected to concentrated forces," Shanghai J Mec, vol. 15, no. 1, pp. 55–60, 1994.
• ⁴⁴
  W. Zhou, Q. Yue, Q. Wang, Y. Feng, and X. Chang, "The boundary element method for elasticity problems with concentrated loads based on displacement singular elements," Eng. Anal. Bound. Elem., vol. 99, pp. 195–205, 2019., http://dx.doi.org/10.1016/j.enganabound.2018.11.016
  » http://dx.doi.org/10.1016/j.enganabound.2018.11.016
• ⁴⁵
  N. Kamiya and Y. Sawaki, “An efficient BEM for some inhomogeneous and nonlinear problems,” in Proc. 7th Int . Conf., Italy, 1985, pp. 13-59.
• ⁴⁶
  M. Guiggiani and A. Gigante, "A general algorithm for multidimensional cauchy principal value integrals in the boundary element method," J. Appl. Mech., vol. 57, no. 4, pp. 906–915, 1990, http://dx.doi.org/10.1115/1.2897660
  » http://dx.doi.org/10.1115/1.2897660
• ⁴⁷
  M. Guiggiani, G. Krishnasamy, T. Rudolphi, and F. Rizzo, "A General Algorithm for the Numerical Solution of Hypersingular Boundary Integral Equations," J. Appl. Mech., vol. 59, no. 3, pp. 604–614, 1992, http://dx.doi.org/10.1115/1.2893766
  » http://dx.doi.org/10.1115/1.2893766
• ⁴⁸
  E. Lorentz, "A nonlocal damage model for plain concrete consistent with cohesive fracture," Int. J. Fract., vol. 207, no. 2, pp. 123–159, 2017, http://dx.doi.org/10.1007/s10704-017-0225-z
  » http://dx.doi.org/10.1007/s10704-017-0225-z
• ⁴⁹
  I. Zreid and M. Kaliske, "A gradient enhanced plasticity–damage microplane model for concrete," Comput. Mech., vol. 62, no. 5, pp. 1239–1257, 2018, http://dx.doi.org/10.1007/s00466-018-1561-1
  » http://dx.doi.org/10.1007/s00466-018-1561-1
• ⁵⁰
  A. Hillerborg, M. Modéer, and P. Petersson, "Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements," Cement Concr. Res., vol. 6, no. 6, pp. 773–781, 1976, http://dx.doi.org/10.1016/0008-8846(76)90007-7
  » http://dx.doi.org/10.1016/0008-8846(76)90007-7
• ⁵¹
  P. Petersson, “Crack growth and development of fracture zones in plain concrete and similar materials,” Ph. D. thesis, Lund Inst. Tech. Sweden, 1981.
• ⁵²
  B. Gross, E. Roberts, and J. E. Srawley, "Elastic displacements for various edge-cracked plate specimens," Int. J. Fract. Mech., vol. 4, no. 3, pp. 267–276, 1968, http://dx.doi.org/10.1007/BF00185262
  » http://dx.doi.org/10.1007/BF00185262
• ⁵³
  H. Tada, P. C. Paris, and G. R. Irwin, Stress Analysis of Cracks Handbook United States: American Society of Mechanical Engineers, 2001.
• ⁵⁴
  P. Bocca, A. Carpinteri, and S. Valente, "Size effects in the mixed mode crack propagation: Softening and snap-back analysis," Eng. Fract. Mech., vol. 35, no. 1-3, pp. 159–170, 1990, http://dx.doi.org/10.1016/0013-7944(90)90193-k
  » http://dx.doi.org/10.1016/0013-7944(90)90193-k
• ⁵⁵
  A. Saleh and M. Aliabadi, "Crack growth analysis in concrete using boundary element method," Eng. Fract. Mech., vol. 51, no. 4, pp. 533–545, 1995, http://dx.doi.org/10.1016/0013-7944(94)00301-w
  » http://dx.doi.org/10.1016/0013-7944(94)00301-w
• ⁵⁶
  A. Jefferson, B. Barr, T. Bennett, and S. Hee, "Three dimensional finite element simulations of fracture tests using the Craft concrete model," Comput. Concr., vol. 1, no. 3, pp. 261–284, 2004, http://dx.doi.org/10.12989/cac.2004.1.3.261
  » http://dx.doi.org/10.12989/cac.2004.1.3.261
• ⁵⁷
  Ł. Kaczmarczyk, M. Nezhad, and C. Pearce, "Three-dimensional brittle fracture: configurational-force-driven crack propagation," Int. J. Numer. Methods Eng., vol. 97, no. 7, pp. 531–550, 2013, http://dx.doi.org/10.1002/nme.4603
  » http://dx.doi.org/10.1002/nme.4603

Edited by

Editors: Samir Maghous, Guilherme Aris Parsekian.

Publication Dates

• Publication in this collection
  09 Mar 2022
• Date of issue
  2022

History

• Received
  15 Sept 2021
• Accepted
  05 Feb 2022