Formulation and implementation of cohesive fracture models in the symmetric galerkin boundary element method. Study of mode I crack growth

FORMULATION AND IMPLEMENTATION OF COHESIVE FRACTURE MODELS IN THE SYMMETRIC GALERKIN BOUNDARY ELEMENT METHOD. STUDY OF MODE I CRACK GROWTH L. Távara1, V. Manti 1, L. Gray2, F. París1, A. Salvadori2,3 1 Grupo de Elasticidad y Resistencia de Materiales Escuela Técnica Superior de Ingenieros, Universidad de Sevilla Camino de los Descubrimientos s/n, 41092 Sevilla, España. Telf: 954487300 E-mail: ltavar[email protected].es, [email protected], [email protected] 2 Computer Science and Mathematics Division Oak Ridge National Laboratory Oak Ridge, TN 37831-6367, USA e-mail: [email protected] 3 Department of Civil Engineering, University of Brescia, Via Branze 38, 25123 Brescia, Italy e-mail: [email protected] RESUMEN Se desarrolla e implementa, en un código computacional, una formulación simétrica de ecuaciones integrales de contorno para la solución de problemas con presencia de grietas, cuya propagación está controlada por un modelo cohesivo de fractura. Se considera un material homogéneo isótropo y elástico-lineal. El código desarrollado se aplica al ensayo wedge split, comparando los resultados numéricos obtenidos con los datos experimentales disponibles en la bibliografía. Es importante mencionar que el uso del método de los elementos de contorno es bastante atractivo para esta clase de problemas debido a que todas las no linealidades se encuentran localizadas en el contorno de un dominio elástico-lineal. Además el uso del enfoque de Galerkin para la discretización de la formulación simétrica de ecuaciones integrales de contorno asegura la eficiencia, estabilidad y convergencia rápida de la solución numérica. ABSTRACT A symmetric boundary integral formulation is developed and implemented in a computational code for the solution of problems with cracks whose growth is governed by a cohesive fracture model. A homogeneous isotropic and linearelastic material is considered. The computational code has been applied to the wedge split test, and the numerical results with the experimental data presented in the literature. It is important to mention that Boundary Element Methods (BEMs) are very attractive for this class of problems because all the non-linearities are located on the boundary of a linear elastic domain. The symmetric-Galerkin formulation of Boundary Integral Equations (BIEs) ensures efficiency, stability and good convergence properties of the numerical solution. KEYWORDS: Cohesive Zone Model, boundary integral equation, SGBEM. 1. INTRODUCTION The methods employed to simulate crack propagation have been traditionally based on Linear Elastic Fracture Mechanics (LEFM) and have assumed the presence of a crack. This fact made difficult to study crack initiation. Recently, other models have been intensively developed, e.g. the Cohesive Zone Models (CZMs), which assume hypotheses different from those adopted in LEFM, and avoid the presence of a stress singularity at the crack tip. These models are suitable to study both crack initiation and propagation, and also to estimate the fracture energy and the maximum allowable load of a structure. In the present work, the Ortiz  Pandolfi [1] CZM is implemented in a 2D Symmetric Galerkin BEM (SGBEM) code. The original version of this code [2] solved plane elastic problems including several homogenous isotropic linear-elastic materials with traction-free cracks inside a homogenous material. The materials were considered to be perfectly bonded along their interfaces. Constitutive equations of a CZM include a representation of a softening branch, which makes a problem with cohesive cracks strongly non-linear. Thus, the development and implementation of a suitable solution algorithm capable of following the evolution of the cohesive zone (modeling the crack growth) becomes an important issue. An arc-length control combined with a Newton-Raphson algorithm for iterative solution of nonlinear equations is used in the present work. Anales de Mecánica de la Fractura 26, Vol. 2 (2009) 435 The SGBEM has several advantages in comparison with the traditional collocation BEM, e.g., a consistent treatment of subdomain corners or discontinuities of boundary conditions (changing abruptly their value or kind) where traction discontinuities can take place. SGBEM provides the required number of equations without the necessity of additional equations as in the case of the collocational BEM. The SGBEM uses both the strongly singular displacement BIE and the hypersingular traction BIE in such a way that the discretizations of these BIEs leads to a symmetric linear system of algebraic equations, with positive or negative definite diagonal blocks associated to unknown tractions or displacements. 2. SGBEM A short review of the SGBEM is provided in this section. The implementation details of the algorithm employed herein are discussed in Sutradhar et al. [3] and Gray [4], works by Bonnet and co-workers [5, 6] can be recommended as well. 2.1 Boundary Integral Equations The primary boundary integral equation for elasticity is the Somigliana displacement identity: )(),([)( 2 1 )( QuQPTPuP jkjk , (1) 0)](),( dQQtQPU jkj where is the boundary of the domain , P is a point on a smooth part of , and u and denote the displacement and traction vectors, respectively. The integral kernel U(P,Q) is generally taken as the Kelvin fundamental solution for the displacement at Q given a point load at P in an infinite medium. For the SGBEM in general and for a fracture analysis in particular, the Somigliana traction identity is essential. Formally differentiating Eq. (1) with respect to P, applying the Hooke constitutive law, and multiplying by the unit normal vector to at the point P yields: )(),([)( 2 1 )( QuQPSPtP jkjk , (2) 0)](),( dQQtQPT jkj A rigorous derivation of Eq. (2) can be found in the cited references on SGBEM. Notice that, due to the singularity in the kernel functions at r=|P-Q|=0, moving the derivative with respect to P under the integral is, in general, not justified. 2.2 Galerkin discretization scheme In a Galerkin formulation, the displacement and traction integral equations are enforced on average, in the form: B ldPPP ,0)()( B ldPPP ,0)()( (3) where the weight function l(P) is comprised of all shape functions that are equal to one at the node Pl and zero at other nodes. In this way a sufficient number of equations to solve for the boundary unknowns can be obtained. The shape functions themselves are determined by the choice of how the boundary displacements and tractions are interpolated. Note that in the Galerkin discretization of the above BIEs the source and field points P and Q are treated equally, and the weakly singular kernel U and the hypersingular kernel S fulfill the following reciprocity equations: ).P,Q(S)Q,P(S),P,Q(U)Q,P(U jkkjjkkj (4) Thus, if displacements are specified everywhere on the boundary, the displacement equation in (3) leads to a symmetric system of equations for the unknown tractions. Similarly the traction equation in (3) yields a symmetric matrix if tractions are prescribed along the whole boundary. In general, if the displacement equation is employed on the part of the boundary where displacements are specified, while the traction equation is employed on the part of the boundary where tractions are known, then the resulting linear system is symmetric. This follows from the fact that T and T* are adjoint kernels. For standard fracture analysis problems, wherein the boundary conditions on the crack are specified tractions, the symmetry is remarkably simple: the above prescription (writing the traction equation on the crack surface) retains the symmetry, with the proviso that the unknowns on the fracture surface are now the jumps in displacements, and the complementary variables are the sum of the known tractions, see e.g. [4]. At first glance, it might appear that the symmetry would be out of the question for a cohesive fracture model. However, and as shown in [7] for interface cracks and in the present work for cracks in homogenous materials, it is possible (and desirable) to obtain a symmetric formulation. 3. COHESIVE ZONE MODELS Cohesive zone models combine the Strength of Materials formulation for crack initiation with the Fracture Mechanics for crack propagation. Cohesive zone formulations relate displacement discontinuities across the crack, i, to the traction vector, ti, in a zone located ahead of the crack tip [8, 9]. Anales de Mecánica de la Fractura 26, Vol. 2 (2009) 436 A constitutive law relating the cohesive tractions to the relative displacements is required for modelling the behaviour of the material in the process zone. The constitutive law may be formally written as ti = k. i, where k is generally a function of i. A fundamental aspect in the formulation of the constitutive model is the requirement that the energy dissipated at crack propagation must be equal to the fracture toughness, i.e., the following relation must be satisfied: if cii Gdt 0 , (5) where if is the maximum value of a relative displacement 3.1 Ortiz- Pandolfi model The relation between tractions and relative displacements in the cohesive zone proposed in [1], requires the concept of effective opening displacement : 2 2 2 1, (6) Different weights are assigned to the normal opening displacement )( 1 and sliding )( 2 through the parameter . Following [1] the existence of a free energy density per unit undeformed area is postulated. In isothermal conditions, under the assumption of isotropic material, it has the form: ),( q, (7) where q is a suitable collection of internal variables which describe the elastic processes attendant to decohesion. From the first and second laws of thermodynamics, it is possible to show that the cohesive law takes the form: ][gradt. (8) Finally, the evolution of internal variables is governed by a set of kinetics relations of the general form ),( qq f It is assumed that the free energy depends on only through the effective opening displacement. This implies that in (7): 21 nn t2 2 1 2,1 ][ i i grad (10) where ni is the component of the unit vector in the idirection. If no unloading is considered, ( / ) may be taken to be independent of q. In such a case a simple expression for the potential is furnished by Smith and Ferrante´s universal binding law [1]: c c cc ee 11)( , (11) where e 2.71828 is the base of the natural logarithm, c is the maximum cohesive normal traction and c is a characteristic opening displacement. From (11) it is easy to obtain: c c ceet)( (12) It is interesting to specify (10) for pure mode I. The normal tractions are shown in Figure 1 and it reads as follows: c c ceet 1 1 )( (13) 0 2 4 6 1 c 1 0.5 0 0.5 1 1 .5 t1 c Figure 1. Normal tractions across the cohesive surface as a function of 1 with 2=0. 4. CZM AND SGBEM Consider a body of an arbitrary shape B which contains a crack (cohesive or not). The boundary of the body B is composed of the non-crack boundary b and the crack surface c. The crack surface c consists of two coincident surfaces and which represent the upper and lower crack surfaces respectively. As explained in [2] it suffices to discretize the upper crack surface as the crack surfaces are usually symmetrically loaded, i.e . Thus the Somigliana displacement (u-BIE) and traction (t-BIE) identities written for an interior point P take the following form: c c c t c c t b jkjjkjk dQQuQPTQtQPUPu )](),()(),([)( , (14) c jkj dQQuQPT )](),( Anales de Mecánica de la Fractura 26, Vol. 2 (2009) 437 b jkjjkjk dQQuQPSQtQPTPt )](),()(),([)( * . (15) c jkj dQQuQPS )(),( In the Galerkin approximation for the non-crack boundary b, the limit of (14) and (15) is taken as Pu and t, respectively; where u is the portion of the boundary with prescribed displacements uu and t is the portion of the boundary with prescribed tractions tt. As tractions are prescribed on the crack surface , only Equation (15) is written for source points on . c c Discretizing the limit forms of Equations (14) and (15) the following system is obtained in block matrix form: c b cccb bb c b cccb bcbb t t GG G u u HH HH 0, (16) where the subscripts b and c denote the contribution of the non-crack boundary and (upper) crack surface respectively. The vector tb includes known tractions tt and unknown tractions tu. Similarly, ub includes the known displacements uu and unknown displacements ut. Equation (16) can be written in terms of the known and unknown boundary displacement and traction values as: c t u ccctcu tttu utuu c t u ccctcu tctttu ucutuu t t t GGG GG GG u u u HHH HHH HHH 0 0 , (17) where, the subscripts u, t and c represent the terms corresponding to the non-crack boundary with prescribed displacements u, non-crack boundary with prescribed tractions t and the crack surface , respectively. To rearrange Equation (16) into the form [A]{ } = { }, multiply the t-BIEs by -1, and make use of the relation in the crack cohesive zone (which relates the traction along the crack cohesive surface ( = ti) with the crack opening displacement (uc = i) in a specific period of time). In this way we finally arrive at the system written in terms of rates of elastic variables: c x b ii kt c t tctucu tttutu tutuuu c t u ccccctcu tctttu ucutuu tGuH tGuH tGuH u u t kGHHG HHG HHG (18) The final coefficient matrix of this system is symmetric due to the reciprocity relations of the integral kernel tensors as shown in [10] where a similar matrix is obtained for the case of traction free cracks. The only difference is the {3,3} block of the linear system matrix (-Hcc - kGcc). At first sight the (-Hcc - kGcc) block appears nonsymmetric but as the Gcc matrix is a mass matrix this block keeps the desired symmetry. 5. NONLINEAR SOLUTION ALGORITHM The non-linear nature of problems found in structural analysis, as the case of modelling cohesive fracture, often requires the use of sophisticated analysis techniques. The arc-length method is a powerful solution method, allowing for the evolution of the equilibrium states of a problem at various load levels [11]. All these equilibrium states trace the loaddisplacement response of the structure in which the applied load varies proportionally as a function of a unique load parameter called herein . In such a case, for a system with n degrees-of-freedom (DOF), the n+1 unknowns u = (x, ) completely define the problem. Figure 2. Arc-length method with the Newton-Raphson method as iterative scheme. In Figure 2 the variables used in the arc-length method are shown, where the vectors up (p = number of step), in the case of SGBEM, is formed by the unknown displacements and tractions in a converged equilibrium point, while xi (i = number of iteration) is the vector of unknowns in the Newton-Raphson iterations between two converged equilibrium points, before equilibrium is reached. The arc-length scheme has several forms in its discrete formulation; the one used in the present work is the normal-flow algorithm [12], where successive Newton- Raphson iterations converge to the equilibrium solution along a path which is normal (in an asymptotic sense) to the so-called Davidenko flow. The Davidenko flow can u0 u1 (x0, 0 ) u2 1 0 i x0 x1 d1 1 d0, 0 X x0 x1 dx p Anales de Mecánica de la Fractura 26, Vol. 2 (2009) 438 be described by considering a small perturbation to the nonlinear system of equations. The flow chart of the normal flow algorithm is shown in Figure 3, where the following variables are used: Sign[ ] is the Sign function (19a) dI = [Kt]-1{r} (19b) dII = [Kt]-1{fext} (19c) where Kt is the tangential stiffness matrix of the considered structure, {r} is the residual vector and {fext} is the external load and constraints vector. One very important issue of the procedure is the scaling of the known and unknown variables involved in the solution of a nonlinear system of equations. The elements in the final systems should have similar orders of magnitude, so as to aid the performance of the nonlinear numerical solvers. Loop on steps (p = 1, pmax) (Predictor phase) d0 (is assigned only in the first step) sig = Sign [dxp T.dII] d0 = sig . d 0 dx0 = d 0 . dII x0 = u0 + dx0 Loop on iterations (i = 1, imax) (Corrector phase) ]|.[]|[ ]|.[]|[ i xII T i xII i xI T i xII idd dd d dxi = d i . dII + dI xi= xi + dxi convergence test end loop Figure 3. Flow chart of the normal flow arc-length method. 6. NUMERICAL RESULTS In order to verify the capability of the numerical model to reproduce experimental results, the wedge split test for a concrete mix, studied in [13], was modelled by the SGBEM code. The material characteristics are Young´s modulus E = 25200 MPa, Poissons ratio = 0.22, specific fracture energy GIc = 101 J/m2, and the values for the Ortiz model are maximum cohesive stress c = 2.3 MPa and critical opening displacement c = 1.61547x10-5 m. It is important to mention that in this case the parameter is not used, due to the Mode I character of the problem. 30 80 mm 200 200 70 70 60 4 Figure 4. Specimen dimensions of the wedge split test. The specimen dimensions are shown in Figure 4. The numerical results obtained with the SGBEM code are represented by the load  displacement global response curve shown in Figure 5. The solution algorithm was implemented using the full Newton-Raphson method. The prediction phase includes the determination of the arc-length step size at each increment. The first step includes the selection of an appropriate value for the arc length d 0. The arc-length is adjusted from one step to the next using the following simple formula: dp+1= d n . m / n (20) where m is the number of iterations that were required at the previous step and n is the (user specified) desired number of iterations at each step. This procedure allows larger steps to be taken when the solution is converging easily, and forces the solver to take smaller steps when convergence is more difficult. For the present work, n was selected between 3 and 4 so as to achieve the most favourable results in terms of computational time. 0 400 800 1200 1600 2000 2400 2800 3200 0.00E+00 1.00E-04 2.00E-04 3.00E-04 4.00E-04 5.00E-04 Opening Displacement (m) SGBEM Test 1 Test 2 Figure 5.Load Displacement prediction and experimental results [13]. Anales de Mecánica de la Fractura 26, Vol. 2 (2009) 439 7. CONCLUSIONS In the present work a symmetric boundary integral formulation of single-domain problems with cohesive cracks has been proposed. This approach is likely to be suitable for engineering applications involving isotropic materials. The implementation of the formulation in a 2D SGBEM code has been carried out. The introduction of the cohesive zone requires an iterative solution procedure to solve the equations resulting from the boundary integral formulation; the arc-length with the normal flow method has been implemented. The existence and uniqueness of the solution of the equations basically depend on the adopted cohesive law. Non-linear constitutive equations typically present a softening branch and this peculiarity can cause a multiplicity of solutions to the rate problem. As shown by the numerical results presented, the cohesive zone formulation correctly modelled the crack growth behaviour. By adjusting the parameters of the discrete model ( c and c in the case of the Ortiz- Pandolfi model for Mode I), predicting the real behavior of structures should be possible. The present work is a starting point to study the possibilities of different cohesive models presented in literature, as well as different arc-length solver schemes applied in a SGBEM setting. ACKNOWLEDGEMENTS The work was supported by the Junta de Andalucía (Projects of Excellence TEP-1207, TEP 02045 and P08- TEP 04051) and the Spanish Ministry of Education and Science through Projects TRA2005-06764 and TRA2006-08077. The authors thank Dr. A.-V. Phan (University of South Alabama) for the BEM code employed in this work. REFERENCES [1] Ortiz M. and Pandolfi A., Finite-Deformation irreversible cohesive elements for three dimensional crack propagation analysis, International Journal for Numerical Methods in Engineering. 44, 1267-1283, 1999. [2] Phan A. V., Napier J. A. L., Gray L. J., and Kaplan T. Symmetric-Galerkin BEM simulation of fracture with frictional contact. International Journal for Numerical Methods in Engineering, 57(6):835851, 2003. [3] Sutradhar, G. Paulino, Gray L., Symmetric Galerkin Boundary Element Method, Springer, 2008. [4] Gray L. J., Evaluation of singular and hypersingular Galerkin boundary integrals: direct limits and symbolic computation. In V. Sladek and J. Sladek, editors, Singular Integrals in the Boundary Element Method, Advances in Boundary Elements, chapter 2, pages 3384. Computational Mechanics Publishers, 1998. [5] Bonnet M., Boundary Integral Equation Methods for Solids and Fluids. Wiley and Sons, England, 1995. [6] Bonnet M., G. Maier, and C. Polizzotto, Symmetric Galerkin boundary element method. ASME Appl. Mech. Rev., 51(11):669704, 1998. [7] Salvadori, A., A symmetric boundary integral formulation for cohesive interface problems, Computational Mechanics, 32, 381-391, 2003. [8] Dugdale D.S., Yielding of steels sheets containing slits, Journal of Mechanics and Physics of Solids, 8, 100-104, 1960. [9] Hilleborg A. Modeer M. and Petersson P.E., Analysis of a crack formation and crack growth in concrete by fracture mechanics and finite elements, Cement and Concrete Research, 6, 773- 782,1976. [10] Gray L. J. and Paulino G. H., Symmetric Galerkin boundary integral fracture analysis for plane orthotropic elasticity. Computational Mechanics, 20:2633, 1997. [11] Fafard M. and Massicotte B., Geometrical Interpretation of the arc-length method, Computers & Structures, 46, 603-615, 1993. [12] Ragon S.A., Gürdal Z. and Watson L.T., A comparison of three algorithms for tracing nonlinear equilibrium paths of structural systems, International Journal Solid & Structures, 39, 689- 698, 2002. [13] Denarié E. Saouma V.E. Iocco A. Varelas D., Concrete fracture process zone characterization with fiber optics, Journal of Engineering Mechanics, 127:5, 494-502, 2001. Anales de Mecánica de la Fractura 26, Vol. 2 (2009) 440