Cohesive crack growth modelling based on an alternative nonlinear BEM formulation

Engineering Fracture Mechanics 111 (2013) 86–97

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Cohesive crack growth modelling based on an alternative nonlinear BEM formulation Hugo Luiz Oliveira, Edson Denner Leonel ⇑ University of São Paulo, School of Engineering of São Carlos, Department of Structural Engineering, Av. Trabalhador São Carlense, 400, 13566-590 São Carlos, SP, Brazil

a r t i c l e

i n f o

Article history: Received 1 April 2013 Received in revised form 27 August 2013 Accepted 9 September 2013

Keywords: Cohesive crack growth Dipoles of stresses Initial stress field Quasi-brittle fracture Boundary element method

a b s t r a c t Fracture mechanics and crack propagation problems have been widely studied by the scientific community in recent years, because crack growth phenomenon can explain the failure of structures. In order to model accurately the structural behaviour of complex engineering structures numerical techniques are required. In this regard, the boundary element method (BEM) has been widely used to solve complex engineering problems, especially those where its mesh dimension reduction includes advantages on the modelling. This paper addresses the analysis of crack propagation in quasi-brittle materials using an alternative BEM formulation. In this type of problem, the damaged zone ahead the crack tip is modelled based on the fictitious crack model. Therefore, the residual resistance of the damaged zone is represented by cohesive stresses, which tends to close the crack faces. The alternative BEM formulation proposed aims at modelling the cohesive stresses using the domain term of the direct integral representation. This term is degenerated, in order to be non null only at the fictitious crack path. As a result of the domain term manipulation, appears a dipole of stresses, which will govern the cohesive stresses. It leads a nonlinear formulation, where the cohesive stresses are determined according the crack opening displacement. Linear, bi-linear and exponential cohesive laws are used to model the residual resistance of the cohesive zone. The results achieved by the proposed formulation are compared with experimental and numerical ones available in literature, in order to validate and prove its robustness and accuracy. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The collapse of structures caused by discontinuities’ growth can only be accurately analysed since all discontinuities be properly modelled by the theory adopted. In this regard, fracture mechanics is recognized as a robust and accurate theory for modelling crack growth phenomenon and consequently the collapse of structures due to crack (discontinuity) propagation. This subject has received large attention by the scientific community in recent years [1,2]. Especially in the context of the development of numerical formulations and mathematical techniques capable at modelling crack growth processes in complexes systems and bodies. The linear elastic fracture mechanics is accurately applied when the size of the fracture process zone (FPZ) in front of the crack tip is small compared to the crack size or even the specimen size. This case is observed in brittle bodies. When quasibrittle materials are considered, other models must be used to take into account the mechanical behaviour of the FPZ, because its length is not small enough. The cohesive crack model is the simplest of such nonlinear fracture mechanics models. ⇑ Corresponding author. Tel.: +55 16 3373 8211; fax: +55 16 3373 9482. E-mail address: [email protected] (E.D. Leonel). 0013-7944/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engfracmech.2013.09.003

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

87

Nomenclature FBZ Dw Dwc

fracture process zone crack opening displacement threshold value of crack opening displacement r cohesive stresses rct material tensile resistance hp crack propagation angle X0 initial stress region ulk ; plk and elkj fundamental solutions for displacements, tractions and strains bi components of body forces E Young’s modulus H; G corresponding matrices to take into account displacement and traction effects l shear modulus gk cosines of normal direction hp crack propagation angle t Poisson’s ratio K and K 0 algebraic matrices that include the integral kernels due to dipoles Q vector describing dipoles q dipoles of stresses

Considering the cohesive crack model, the propagation is governed by a stress–displacement relation across the crack faces near the tip. This model was introduced in the early sixties for ductile materials by Barenblatt [3] and Dugdale [4]. In the late seventies, Hillerborg et al. [5] introduced the concept of fracture energy into the cohesive crack model and proposed a number of traction–displacement relationships for concrete. The cohesive crack model has been included in several numerical methods as boundary element method (BEM), mesh less methods and finite element method (FEM). Several algorithms have been proposed in order to solve the nonlinearities coming from the FPZ [6,7]. When cohesive crack growth is simulated, one of the most significant difficulties that arise is the requirement of mesh’s conformity to crack geometry. It needs remeshing procedures during crack propagation, leading a high computational cost for large and complexes geometries and boundary conditions. In this regard, FEM is not a numerical technique completely adapted to solve this problem. Due to its domain mesh, the remeshing procedures are not a simple task. Moreover, FEM requires a very fine mesh in order to represent the stress singularities at the crack tip, leading the construction of large stiffness’ matrix and, consequently, low performance. Despite these difficulties, there are some works in literature associating the cohesive crack growth problem to FEM, as presented by [8,9]. On the other hand, BEM has been extensively used to solve many problems in engineering, especially those problems related to fracture mechanics domain. In such problems, this numerical method requires only boundary and crack surface discretizations. Therefore, remeshing difficulties are considerable avoided. Moreover, compared to domain mesh approaches, BEM is more efficient, particularly in solving mixed mode crack propagation problems, due to its efficiency in stress concentration modelling and mesh reduction aspects [10–14]. In addition to that, it is important to emphasize that crack propagation modelling in quasi-brittle materials, based on cohesive approach, leads the solution of a nonlinear problem, generally relating the crack opening displacement to tractions on the crack surfaces. BEM has demonstrated to be a very efficient and robust technique in solving nonlinear problems as presented by [6,15,16]. Therefore, BEM is recommended for modelling cohesive crack growth. BEM has been used by many researchers in order to handle cohesive crack growth modelling. Among them, it is worth to mention [17], where the dual BEM formulation was used to simulate cohesive mixed mode crack growth in plane structures composed by concrete [18] used the single domain dual boundary element formulation as a boundary element approach. In this formulation, the cohesive zone is incorporated into the formulation, resulting in a nonlinear problem. The multi-domain BEM and the symmetric Galerkin BEM were used by [19,20], respectively, for analyses of cohesive crack growth. For the same problem, [21] used multi-zone BEM. During decades, numerical methods have been used and regarded as suitable tools to predict the fracture and failure of engineering structures. Due to the large advances in the field of science and engineering, the need to analyse larger and more complex structures in numerical form has become a challenger issue. Therefore, the development of numerical tools and more efficient algorithms regarding the computational cost and accuracy of results is quite essential. In this context, this paper presents an alternative BEM formulation based on dipoles of stresses applied to model cohesive crack propagation. The proposed formulation aims at modelling the cohesive stresses, on the fictitious crack, using the domain term of the direct BEM integral representation. This term is modified, in order to be non null only at the fictitious crack path. As a result of the domain term manipulation, appears a dipole of stresses, which will govern the cohesive stresses. Therefore, the dipoles are defined in order to introduce corrective stresses into domain according the cohesive law adopted. This formulation solves

88

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

the crack growth problem using three algebraic equations per source point positioned at the crack path. By comparing it with the classical dual BEM, that uses four algebraic equations per crack source point, it is a great advantage. The proposed formulation is nonlinear, since the cohesive stresses are governed by the crack opening displacement values (cohesive law). The nonlinear problem is solved using classical Newton–Raphson scheme, considering prevision and correction steps, in which the corrections into the cohesive crack stresses are performed by applying a non-equilibrated stress vector, keeping constant all relevant matrices. The cohesive stresses were modelled using three cohesive laws: linear, bi-linear and exponential. The results achieved by the proposed formulation are compared with experimental and numerical ones, in order to validate and prove its robustness and accuracy. 2. Fracture mechanics aspects 2.1. Cohesive crack model The linear elastic fracture mechanics is an important approach to solve many problems in structural engineering. Particularly those where the FPZ surrounding cracks is reduced enough to be possible to neglect its nonlinear effects. However, for quasi-brittle materials, the damaged zone ahead the crack tip is large enough to produce nonlinear effects that cannot be neglected. The cohesive crack model is an appropriate approach to take into account these effects. Considering this model, the energy dissipation phenomenon is assumed to occur along a fictitious crack positioned in front of the crack tip, reducing, therefore, by one the dimension of the dissipation zone. In this work, the energy dissipation process is approximated by a simple softening law, which is assumed to govern the residual resistance of the material along the fictitious crack. This law relates the fictitious crack opening displacement, Dw, to residual tensile stresses, r, acting at the crack surfaces. Several cohesive crack laws relating cohesive stresses and crack opening displacement have already been proposed in literature. Three of them are often adopted to carry out crack growth analysis in quasi-brittle materials. The simplest law is given by a linear function relating the cohesive stresses to fictitious crack opening displacement smaller than a critical value, Dwc. For fictitious crack openings larger than Dwc, cohesive stresses are assumed equal to zero, Fig. 1(a). The equations that represent the linear cohesive law are the following:

r ¼ Ee   if e 6 ec rt ðDwÞ ¼ rct 1  DDwwc if 0 6 Dw 6 Dwc rt ðDwÞ ¼ 0 if Dw > Dwc

ð1Þ

in which rct indicates the material tensile resistance. An alternative model relating cohesive stresses and fictitious crack opening displacement is the bi-linear model, Fig. 1(b). This criterion is mathematically represented by the following equations:

r ¼ Ee  c 00  rt ðDwÞ ¼ rct  rDt wr00 t Dw 00



rt ðDwÞ ¼ Dwr00tDDwwc þ r00t 1  DwD00 wD00 wc rt ðDwÞ ¼ 0



if e 6 ec if 0 6 Dw 6 Dw00 if Dw00 6 Dw 6 Dwc if Dw > Dwc

Considering bi-linear model, the variables

r00t ¼

rct 3

Dw00 ¼

0:8Gf

rct

Dwc ¼

ð2Þ

r00t ; Dw00 and Dwc are defined as:

3:6Gf

rct

ð3Þ

where Gf represents the material fracture energy. The third cohesive crack model considered in this work is represented by an exponential law, Fig. 1(c). Eq. (4) gives the analytical expressions for this cohesive law:

Fig. 1. Cohesive models. (a) Linear model, (b) bi-linear model and (c) exponential model.

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

r ¼ Ee

if

rc Gt Dw

rt ðDwÞ ¼ rct e

e 6 ec

ð4Þ

if Dw > 0

f

89

2.2. Crack growth scheme Considering cohesive crack growth analysis, one major task is to determine when the cracks grow. In this context, the real stress state at the crack tip has to be compared with an ultimate stress state, given by an adopted criterion. In this work, the ultimate stresses are calculated using the Rankine’s model, which has also been successfully used by [6,17]. In order to achieve accurately the stress state at the tip, a polynomial extrapolation process has been adopted to take into account the continuity of stresses along all points surrounding the crack tip. Several internal points are defined ahead of crack tip, as shown in Fig. 2. The number of semi-circles and internal points into each one can be chosen according the desired accuracy. The polynomial degree, which will approximate the stress state, is chosen according the number of the semi-circles. Considering the case illustrated in Fig. 2, the extrapolation process is performed using third degree polynomial. The extrapolation is performed for each radial line. Thus, the real stress state at the tip is obtained by averaging all extrapolated values. The direction of the new crack length increment, hp, is given by the direction perpendicular to the maximum circumferential tensile stress. According this criterion, the crack growth direction is given by:

 hp ¼ ArcTan



sxy

ð5Þ

r1  ry

where r1 is the maximum tensile stress value. The cracks grow when the maximum tensile stress at the crack tip is larger than the ultimate tensile value given by Rankine’s criterion. The crack length increment is determined adjusting the size of the element in such a way that the stress state at the new fictitious crack tip be equal to ultimate tensile limit. This length is determined by minimizing the difference between these two stresses states. In this work, this simple minimization problem is solved using the optimization algorithm Golden Section. 3. Bem integral equations Considering an elastic body with domain X and boundary C subjected to an initial stress field r0jk acting into the domain, the following integral representation can be written, as presented by [22].

cilk uik þ

Z C

plk uk dC ¼

Z C

ulk pk dC þ

Z X0

r0jk eljk dX0

ð6Þ

where ulk , plk and elkj are the well known fundamental solutions for displacements, tractions and strains, respectively [22]. X0 is the initial stress region, which for the problem studied represents the FPZ. The initial stress region is a wide thin region C C C C limited by a boundary C , where C ¼ C1 [ C2 as presented in Fig. 3.

Fig. 2. Distribution of internal points ahead the crack tip.

90

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

Fig. 3. Initial stress region.

In order to develop the cohesive fracture formulation, the initial stress integral shown in Eq. (6) has to be conveniently manipulated. Initially, it can be integrated over the initial stress region leading to:

Z X0

r0jk eljk dX0 ¼

Z

C

C

C

ulj r0jk gk dC 

Z X0

ulj r0jk;k dX0

ð7Þ

C

in which gk is the normal vector to boundary C and the term r0jk gk represents the tractions at the initial stress region boundaries, p01 j . Therefore, the first integral on right hand side of Eq. (7) can be rewritten as:

Z

C

C

C

ulj ðf ; SÞp01 j ðSÞdC ¼

Z

C

þ ulj ðf ; Sþ Þp01 j ðS ÞdC1 þ

C

C1

Z

C

 ulj ðf ; S Þp01 j ðS ÞdC2

C

C2

C

ð8Þ

C

where S+ and S are field points positioned at boundaries C1 and C2 , respectively. Assuming that the width of initial stress region is too small compared to its length; thus, the kernels defined at S+ and S can be rewritten for the middle path of the region, S. Therefore, these kernels can be redefined using Taylor series expansion as follows:

ulj ðf ; Sþ Þ ¼ ulj ðf ; SÞ þ ulj ðf ; S Þ ¼ ulj ðf ; SÞ 

@ulj ðf ;SÞ @x1 @ulj ðf ;SÞ @x1

a

ð9Þ

a

Considering the expansion terms presented by Eq. (9), it can be used to rewrite Eq. (8) as:

Z C

C

ulj ðf ; SÞp01 j ðSÞd

C

C ¼

"

Z

ulj ðf ; SÞ

C C1

C

þ

#

@ulj ðf ; SÞ

a

@x1

þ p01 j ðS Þd

C 1

C þ

Z C C2

" ulj ðf ; SÞ



@ulj ðf ; SÞ @x1

#

C

 a p01 j ðS ÞdC2

ð10Þ

C

Assuming that dC1 ¼ dC and dC2 ¼ dC, the equation above can be rewritten as:

Z

C

C

C

ulj ðf ; SÞp01 j ðSÞdC ¼

Z

@ulj ðf ; SÞ

C

@x1

p01 j ðSÞ2adC

ð11Þ

The domain integral presented on right hand side of Eq. (7), which includes stress derivatives, must also be treated. This C C integral term can be transformed into integrals written on boundaries C1 and C2 . In this regard, the width of the initial stress region is assumed as thin enough when compared to its length. Consequently, the stress’ variation along direction x1 is null. Therefore:

r0jk;k ¼

@ r0jk @xk

¼

@ r0jk @xm @ r0jk @x2 @ r0jk @ðr0jk tk Þ ¼ ¼ tk ¼ @x2 @xm @xk @x2 @xk @x2

ð12Þ

in which tk contains the outward normal components along x2 . Based on these assumptions, the domain integral presented on right hand side of Eq. (7) can be rewritten as:



Z X0

ulj r0jk;k dX0 ¼ 

Z X0

ulj

@ðr0jk tk Þ @x2

dX0 ¼ 

Z x2

ulj

@ðr0jk t k Þ @x2

2adx2

ð13Þ

Eq. (13) can be integrated over direction x2 leading to:



Z x2

ulj

@ðr0jk t k Þ @x2

x

1

2adx2 ¼ ½ulj ðr0jk tk Þ2ax22 þ 2

Z x2

@ulj @x2

ðr0jk t k Þ2adx2

ð14Þ

91

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

As the width of the initial stress zone is small, the first term on right hand side of Eq. (14) becomes null. Considering that dxC2 ¼ dC and r0jk tk ¼ p02 j , Eq. (14) can be rewritten as:



Z

X0

Z

ulj r0jk;k dX0 ¼

@ulj @x2

C

p02 j 2adC

ð15Þ

where p02 j indicates the tractions aligned to direction x2 . Therefore, based on results presented in Eqs. (11) and (15), the domain integral shown in Eq. (6) can be transformed as:

Z X0

r0jk eljk dX0 ¼

Z

@ulj ðf ; SÞ @x1

C

p01 j ðSÞ2adC þ

Z C

@ulj @x2

p02 j 2adC ¼

Z C

@ulj @xk

p0k j 2adC

ð16Þ

In order to avoid the local aspect of this deduction, Eq. (16) can be finally rewritten as:

Z X0

r0jk eljk dX0 ¼

Z

@ulj

C

@xk

p0k j 2adC ¼

Z C

@ulj @xm 0k p 2adC @xm @xk j

ð17Þ

Therefore, a domain integral defined over any thin zone was transformed into a line integral. Based on this mathematical manipulation, stresses and displacements analyses for domains where nonlinear behaviours are assumed inside particular narrow regions can be performed. For instance, crack growth analysis in which cohesive stresses are assumed over a finite strip. Although a numerical algorithm based on the definition of a thin finite process zone could be derived, the nonlinear zone must be wide enough to guarantee initial stress finite values. When the thickness of the thin region goes to zero, infinite initial stresses are required due to the nature of the problem. Thus, in order to write a proper integral term for which the strip thickness is assumed as zero, a new tensor, denominated dipole, q must be defined as follows:

@xm 0k p 2a ¼ qm j @xk j

ð18Þ

This new variable, which is represented in Fig. 4, is given by finite values when the initial stresses goes to infinity. Then, considering this new variable and the discussion presented above, Eq. (6) can be redefined as:

cilk uik þ

Z C

plk uk dC ¼

Z C

ulk pk dC þ

Z C

Gljm qm j dC

ð19Þ

where Gljm ¼ @ulj =@xm . This new kernel is defined by:

Gljm ¼ ulj;m ¼

1 fð3  4tÞr;m dlj  r ;j dlm  r ;l djm þ 2r ;j r;l r ;m g 8pGð1  tÞr

ð20Þ

In order to derive the stress integral representation, similar developments can be followed. In this regard, Eq. (19) has to be derived and the Hooke´s law applied. This procedure leads to:

rim þ

Z

Simk ðf ; SÞuk ðSÞdC ¼

C

Z

Dimk ðf ; SÞpk ðSÞdC þ C

Z C

Gljim ðf ; SÞqlj ðSÞdC þ g ljim ðrlj ðpÞÞ

ð21Þ

in which g ljim ðrlj ðpÞÞ is a free term that appears due to the singularity of the problem. The kernel Gljim is given by:

Gljim ¼

1 4pð1  tÞr 2



ð1  2tÞðdij dlm þ djm dil  djl dim Þ  2ð1  2tÞðr ;m r ;l dij þ r ;i r;l dmj  r ;j r ;l dim Þ þ2ðr ;j r ;m dil þ r ;j r ;i dml þ r ;m r ;i djl Þ  8r ;i r;j r ;m r ;l

Fig. 4. Dipoles representation (a) opening and (b) sliding.



ð22Þ

92

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

Other parameters widely important for modelling fracture mechanics problems are the crack opening displacements. Regarding the proposed formulation, it can be obtained from Eq. (19), which has to be applied for collocation points taken over the dipole’s line, where discontinuities of displacements and stresses will appear. This equation has to be written for points symmetrically positioned at S+ and S, as indicated in Fig. 3. As a result, different displacement representations are achieved for these points. By subtracting them, the crack opening displacements are obtained as:



Dw1 Dw2

"

 ¼

12t 2Gð1tÞ

0

0

1 G

#

q11 q12

 ð23Þ

It is worth to stress that the presence of dipoles leads to a displacement discontinuity. On Eq. (23), Dw1 represents the crack opening displacement (mode I) and Dw2 the crack sliding displacement (mode II). The third dipole component does not result displacement discontinuities. Following the same steps performed above, the stress equation can be written for points located along C. After performing properly the limits, in order to take the point to the initial stress region, the free term g ljim ðrlj ðpÞÞ shown in Eq. (21) is achieved. It will be r0jk for source points positioned at the dipole’s line and null elsewhere. Although a dipole tensor is obtained, it is not hard to verify that the component q22 can be written in terms of q11 . Assuming that the normal stress component is continuous along C and a jump equal to r0jk is observed in the direction perpendicular to the crack path, the following value for q22 can be derived:

q22 ¼

t 1t

q11

ð24Þ

4. Algebraic equations In previous section, the integral representations with additional terms for modelling stresses and displacements discontinuities were derived. The boundary discretization leads to the known H and G matrices, which take into account boundary displacements and tractions. The remaining term can also be written into its algebraic form if the variable qm j is assumed to be approximate by standard shape functions along the crack line. After discretizing the whole crack path by elements and performing properly the integrals over them, Eqs. (19) and (21) assume the following algebraic representations:

HU ¼ GP þ KQ

ð25Þ

r þ H 0 U ¼ G0 P þ K 0 Q

ð26Þ

in which K and K0 include the integral kernels due to dipoles Q. The boundary conditions can be applied on Eq. (25). As usual in BEM formulations, all unknown boundary values are stored in a vector X and all known in a vector F. Therefore:

AX ¼ BF þ KQ ) X ¼ M þ RQ

ð27Þ 1

where matrices A and B result from columns change between H and G matrices. M and R are defined as: M = A BF and R = A1K. Similarly, the boundary conditions can be applied on Eq. (26). Performing this step together with the result presented by Eq. (27) one obtains:

r þ A0 X ¼ B0 F þ K 0 Q ) r ¼ N þ SQ

ð28Þ

in which A0 and B0 result from columns change between H0 and G0 matrices. The terms N and S are equal to N = B0 F  A0 M and S = K0  A0 R, respectively. In this paper, the integrals that lead the algebraic equations Eqs. (27) and (28) were evaluated by Gauss–Legendre numerical scheme accomplished with a sub-element technique. Based on these procedures, Eqs. (19) and (21) are transformed into algebraic representations with very low integration error. It is worth to mention that due to the singularity present in the kernel Gljim , only discontinuous elements have to be used at the crack path. The cohesive crack propagation problem can be solved using Eqs. (27) and (28). The first one leads the solution of boundary values whereas the second equation allows the determination of stresses and displacements at crack surfaces. It is worth to emphasize that boundary values are determined since dipoles’ values in equilibrium configuration are achieved. The dipoles are calculated in the context of nonlinear solutions, which are performed in incremental form. In this paper, the nonlinear problem is solved considering a Newton–Raphson scheme, thus, prevision and correction steps are used. In this regard, the exceeding stresses are reapplied on the structure and a dipoles variation is determined, as presented in Eq. (28). As a result of this dipoles variation, the crack opening displacements variation is determined, Eq. (23), which leads a new stress state of equilibrium. This last one is compared with the stresses due to the external loading. The differences between these two stresses’ state are reapplied until the residual stresses vector norm be smaller than a specified tolerance, indicating the convergence. This classical procedure is known as constant operator, since all relevant matrices are kept constant during the iterative process.

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

93

The tolerance to stop the iterative process within an increment of load is applied on the variation of the crack opening displacement corrections, i.e., wi  wi1 6 tolerance. Moreover, it is worth to remark that the total crack opening displacement is always required to compute the local tangent operator for the next iteration. The BEM formulation proposed in this paper presents some advantages when compared with classical BEM approaches as dual BEM and multi-domain technique, which are techniques widely known in literature. The dual BEM is one of the most popular BEM formulations used to analyse random crack growth [14,17]. In this formulation, the singular integral representation (displacement integral equation) is used to determine the algebraic representation related to collocation points defined along one crack surface, whereas the hyper-singular integral representation (traction integral equation) is used to obtain the algebraic representation related to collocation points located along the opposite crack surface. For all external boundaries, the singular representation applied to the boundary collocation points is sufficient to obtain the required algebraic relations. Therefore, dual BEM requires four integral equations (two displacements and two tractions) per source node along the crack path. By comparing it against the proposed formulation, it can be observed that dipoles’ approach requires only three integral equations per source node along the crack path (only those equations related to stresses correction). Therefore, it makes the dipoles’ formulation faster and efficient in the context of computational effort. In addition to that, the nonlinear problem is solved by dual BEM using algebraic equations written over all collocation points whereas dipoles’ formulation solves the same problem using Eq. (28), which relates the dipoles to stresses on the crack line. Then, by comparing it against dual BEM another advantage is observed. Regarding the dipoles’ formulation, the nonlinear problem is solved using an algebraic system smaller than the system used by dual BEM, leading a faster and accurate numerical solution. Considering multi-domain technique, cracks are assumed to occur among the multi-domains interfaces. When linear elastic mechanical behaviour is observed, equilibrium of forces and compatibility of displacements are enforced along all interfaces. The crack growth case is simulated by modifying those conditions, in order to take into account the mechanical influence of the cracks (discontinuities) [16,23]. Therefore, the multi-domain technique is not capable to represent accurately the random crack propagation phenomenon. Based on this technique, the crack path must be known a priori, in order to allow the definition of convenient number of domains, its locations and the mesh among all interfaces. It is important to stress that those inconvenient are not observed if dipoles’ formulation be adopted. Moreover, multi-domain technique requires eight algebraic equations per point along the crack path to assess the structural mechanical behaviour whereas dipoles’ formulation requires three. Besides that, the nonlinear problem is solved by multi-domain technique using all algebraic equations evaluated over the boundaries of all domains. On the other hand, dipoles’ formulation solves the same problem using Eq. (28), which relates dipoles and stresses only at the cohesive crack. Therefore, dipoles’ approach solves the nonlinear problem using a smaller system of equations when compared with multi-domain approach and dual BEM. 5. Applications In order to illustrate the applicability and robustness of the proposed nonlinear BEM formulation, numerical analyses were performed considering crack propagation modelling in plane structures composed by concrete, which is a quasi-brittle material. Three applications are considered. The first structure analysed has an analytical response, which was used to validate the formulation. The second and third applications deal with bending structures, which were analysed experimentally and numerically by other researches. 5.1. Example 1 The first example deals with the structure presented in Fig. 5, which is a plane structure clamped at its left side and subjected to a prescribed displacement at its right boundary. The prescribed displacement leads the appearance of a predom-

Fig. 5. Structure analysed.

94

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

Fig. 6. Load displacement curves.

inant uniaxial stress distribution into the structure. One initial crack is assumed to appears, at the middle span, when the imposed displacement produces a stress level higher than the material tensile resistance. The following material properties were adopted: tensile strength rct ¼ 3:0 MPa, Young’s modulus E ¼ 30000 MPa, Poisson ratio t = 0.0 and fracture energy Gf ¼ 15 kN=m. In order to model the crack propagation in this structure, three cohesive laws were adopted: linear, bi-linear and exponential. The total prescribed displacement, 0.020 m, was applied into 100 increments and the convergence was verified with a tolerance of 104, based on the norm of non-equilibrated stress vector. It is worth to mention that this structure was chosen because it has an analytical solution available according the Theory of Elasticity. Therefore, the numerical results achieved by the proposed formulation will be compared with an analytical solution. The results obtained by the proposed formulation, in terms of load x displacement, are compared with those given by analytical solutions considering three cohesive laws cited in Section 2.1. These curves are shown in Fig. 6, where the displacement was monitored at the middle right boundary. According the figure above, excellent agreement is observed among analytical and numerical responses. The proposed formulation was capable to represent, with very low errors, the predicted analytical results. Therefore, the accuracy of the proposed formulation was tested and approved. It is also important to cite that, when the cohesive stresses become null, the structure observes a rigid body movement. The proposed formulation was capable to reproduce this mechanical behaviour, which is impossible to simulate with FEM. 5.2. Example 2 The three point bending test considered in this example is shown in Fig. 7. The geometry is given by its length, 800 mm, height, 200 mm, and a central notch of 50 mm deep. The experimental results are shown in [17], from where the following properties were obtained: tensile strength rct ¼ 3:0 MPa; Young’s modulus E ¼ 30000 MPa; Poisson ratio t = 0.15 and fracture energy Gf ¼ 75 N=m. In order to model the nonlinear crack propagation in this structure, three cohesive laws were adopted: linear, bi-linear and exponential. The load was applied into 24 increments and the convergence has been verified with a tolerance of 105, based on the norm of non-equilibrated stress vector. The load x displacement curves achieved in this analysis are shown in Fig. 8, where the results obtained by the proposed nonlinear formulation are compared with experimental [17] and numerical responses, dual BEM [6] and multi-domain technique [23], both available in literature.

Fig. 7. Three point bending specimen.

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

95

Fig. 8. Load x displacement curve.

Fig. 9. Crack growth path observed.

According this figure, good agreement is observed between experimental results and those achieved by the proposed formulation. Moreover, the results obtained by dipole’ formulation are equivalent with those found by dual BEM formulation, indicating its robustness. The crack growth path observed in this application is shown in Fig. 9. According this figure, it can be observed that the crack grows in mode I. 5.3. Example 3 The four point bending beam considered in this example is shown in Fig. 10. The structural geometry is given by its length of 675 mm, height of 150 mm and central notch of 75 mm deep. The material properties were taken from [24], who have performed the experimental test: tensile strength rct ¼ 3:0 MPa; Young’s modulus E ¼ 37:000 MPa; Poisson ratio t = 0.20 and fracture energy Gf ¼ 69 N=m. Considering the present analysis, three cohesive laws were adopted: linear, bi-linear and exponential. The load was applied, for all cases, into 24 increments and the adopted tolerance within each increment was equal to 104, based on the norm of non-equilibrated stresses. The load x displacement curves achieved in this application are presented in Fig. 11, where the results obtained by the proposed formulation are compared against experimental [24] and numerical results, dual BEM [6] and FEM [25]. According this figure, it can be observed that the proposed formulation is capable to represent the structural nonlinear behaviour introduced by FPZ. However, it leads to more rigid responses when compared with experimental and numerical results. It may be explained due to snap-back behaviour observed on experimental test. This nonlinear behaviour can be

96

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

Fig. 10. Four point bending beam. Dimensions in mm.

Fig. 11. Load x displacement curves.

Fig. 12. Crack growth path during the loading process.

H.L. Oliveira, E.D. Leonel / Engineering Fracture Mechanics 111 (2013) 86–97

97

modelled using arc-length algorithm. However, this solution technique is not the focus of this paper. Similar numerical instabilities were observed on results presented by [25], where a brittle behaviour was observed after the peak load. Fig. 12 illustrates the crack growth path observed during the beam loading. According this figure it can be observed that the crack grew in mixed mode. The crack path achieved by dipoles’ formulation was compared with the experimental response given by [24] and with the numerical response obtained by XFEM model [26]. It can be observed good agreement among those predicted paths. 6. Conclusions In this paper, the crack growth process in quasi-brittle materials was studied adopting the cohesive crack approach. This complex structural problem can be modelled solving a nonlinear system of equations, which appears due to the dependency between crack opening displacement and cohesive stresses along the crack path. In order to simulate this nonlinear structural problem, BEM has shown to be an accurate and efficient numerical technique. An alternative BEM formulation based on initial stresses field, which leads the appearance of dipoles of stresses, is proposed in order to simulate the nonlinear effects caused by damaged zones present in narrow regions into the domains. The dipoles are determined based on the domain term of the direct BEM integral representation, which is assumed to be non null only at the fictitious crack path. It is worth to emphasize that the proposed formulation requires three algebraic equations to describe the crack mechanical behaviour whereas the classical dual BEM needs four. It leads to a faster computational performance of the proposed formulation when compared against classical approaches. The nonlinear problem was solved using a classical Newton–Raphson scheme, where all relevant matrices are kept constant along the iterative process and the equilibrium configuration is achieved by reapplying the non-equilibrated stresses. The proposed formulation was tested in three applications. The first has an analytical response whereas the second and third have experimental and numerical responses. These responses were used to validate and prove the accuracy and robustness of the proposed formulation. In both applications, the proposed formulation presented good performance and agreement with known results. Acknowledgement Sponsorship of this research project by the São Paulo State Foundation for Research (FAPESP), Project Number 2011/ 07771-7 is greatly appreciated. References [1] Budyn E, Zi G, Moes N, Belytschko T. A method for multiple crack growth in brittle materials without remeshing. Int J Numer Meth Eng 2004;61:1741–70. [2] Oliveira HL, Leonel ED. Dual BEM formulation applied to analysis of multiple crack propagation. Key Eng Mater 2013;560:99–106. [3] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962;7:55–129. [4] Dugdale DS. Yelding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4. [5] Hillerborg A, Modeer M, Peterson PE. Analysis of crack formation and crack growth in concrete by mean of failure mechanics and finite elements. Cem Concr Res 1976;6:773–82. [6] Leonel ED, Venturini WS. Non-linear boundary element formulation with tangent operator to analyse crack propagation in quasi-brittle materials. Eng Anal Boundary Elem 2010;34:122–9. [7] Ameen M, Raghuprasad B. A hibrid technique of modeling of cracks using displacement discontinuity and direct boundary element method. Int J Fract 1994;67:335–43. [8] Carpinteri A. Post-peak and post-bifurcation analysis of cohesive crack propagation. Eng Fract Mech 1989;32:265–78. [9] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46:131–50. [10] Leonel ED, Venturini WS, Chateauneuf A. A BEM model applied to failure analysis of multi-fractured structures. Eng Fail Anal 2011;18:1538–49. [11] Leonel ED, Chateauneuf A, Venturini WS. Probabilistic crack growth analyses using a boundary element model: applications in linear elastic fracture and fatigue problems. Eng Anal Bound Elem 2012;36:944–59. [12] Leonel ED, Venturini WS. Dual boundary element formulation applied to analysis of multi-fractured domains. Eng Anal Bound Elem 2010;34:1092–999. [13] Yan X. Automated simulation of fatigue crack propagation for two-dimensional linear elastic fracture mechanics problems by boundary element method. Eng Fract Mech 2007;74:2225–46. [14] Leonel ED, Venturini WS. Multiple random crack propagation using a boundary element formulation. Eng Fract Mech 2011;78:1077–90. [15] Kumar V, Mukherjee S. Boundary-integral equation formulation for time-dependent inelastic deformation in metals. Int Mech Sci 1977;19(12):713–24. [16] Leonel ED, Venturini WS. Non-linear boundary element formulation applied to contact analysis using tangent operator. Eng Anal Bound Elem 2011;35:1237–47. [17] Saleh AL, Aliabadi MH. Crack-growth analysis in concrete using boundary element method. Eng Fract Mech 1995;51(4):533–45. [18] Yang B, Ravi-Chandar K. A single-domain dual-boundary-element formulation incorporating a cohesive zone model for elastostatic cracks. Int J Frac 1998;93:115–44. [19] Cen, Z; Maier, G. Bifurications and instabilities in fracture of cohesive-softening structures: a boundary element analysis. Fract Eng Mater 1992. [20] Maier G, Novati G, Cen Z. Symmetric galerkin boundary element method for quasi-brittle-fracture and frictional contact problems. Compd Mech 1993;13:74–89. [21] Chen T, Wang B, Cen Z, Wu Z. A symmetric galerkin multi-zone boundary element method for cohesive crack growth. Eng Frac Mech 1999;63:591–609. [22] Brebbia CA, Dominguez J. Boundary elements: an introductory course. Second edition. Southampton: WIT Press; 1992. [23] Ferreira MDC, Venturini WS. Inverse analysis for two-dimensional structures using the boundary element method. Adv Eng Softw 2010;41:1061–72. [24] Galvez JC, Elices M, Guinea GV, Planas J. Mixed mode fracture of concrete under proportional and nonproportional loading. Int J Fract 1998;94:267–84. [25] Cendon DA, Galvez JC, Elices M, Planas J. Modelling the fracture of concrete under mixed loading, Int. J Fract 2000;103:293–310. [26] Abbas S, Alizada A, Fries T-P. Model-independet approaches for the XFEM in fracture mechanics. Int J Numer Meth Eng 2010;1. 0–10.