A discrete model for cyclic mode I loading

Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 www.elsevier.com/locate/cam A discrete model for cyclic mode I loading K. De Profta;∗ , L.J. Sluysb , W.P. De Wildea a Department of Mechanics of Materials and Constructions, Free University of Brussels, Pleinlaan 2, Brussel 1050, Belgium b Delft University of Technology, Stevinweg 1, Delft 2628 CN, The Netherlands Received 12 September 2002; received in revised form 19 May 2003 Abstract The cyclic behaviour of a double-edge notched specimen loaded in tension is studied. Cracks in the material are modelled by displacement discontinuities that can propagate during computation. Within these discontinuities, a cohesive zone model is used. The model assumes an additive split of the inelastic jump into a recoverable and an unrecoverable part. The in uence of model parameters and discretisation is studied and the results have been compared with experimental data. c 2003 Elsevier B.V. All rights reserved.
Keywords: Damage-plasticity; Cohesive zone; Partition of unity 1. Introduction The numerical simulation of fracture of quasi-brittle materials is a challenging subject for many years now. Di erent kinds of models were developed over the years, all trying the simulate fracture independent of the discretisation. In the continuum approach, nonlocal [5] and gradient [4] models have shown to be an e ective tool to overcome mesh dependency. These models need a very ne mesh within the fracture zone, in order to capture the steep deformation gradient. Another way of regularising the problem is to use discrete methods like the cohesive surface methodology or embedded discontinuities. The rst introduces cohesive surfaces at inter element boundaries, allowing cracks to initiate and propagate [8]. The cohesive zones are present from the beginning of the calculation. Since the crack path is not known a priori, the mesh must be re ned in order to have enough possible crack trajectories. The embedded discontinuities [6] allow cracks (discontinuities) to initiate during computation and to run through elements. ∗ Corresponding author. E-mail address: [email protected] (K. De Proft). c 2003 Elsevier B.V. All rights reserved. 0377-0427/$ - see front matter  doi:10.1016/j.cam.2003.05.007 136 K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 In this paper, the cracks are modelled as a displacement discontinuity (strong discontinuity). Using the partition of unity concept, nodes are locally enhanced by adding extra degrees of freedom with the Heaviside function as enhanced basis [7]. In the rst section, the kinematics of the displacement eld and the governing equations are derived for two nonintersecting discontinuities. A combined damage-plastic model is introduced as cohesive zone model in Section 3. Numerical simulation of a double edge-notched specimen subjected to monotonic and cyclic loading is discussed in Section 4. 2. Cohesive zone model based on partitions of unity 2.1. Kinematics of a displacement jump Consider a body crossed by a discontinuity d as shown in Fig. 1. The discontinuity splits the body into two parts, + and − . The normal to the discontinuity is pointed towards + . The displacement eld within the body can be split up into a continuous and a discontinuous part u = û + Hd ũ; (1) where û and ũ are continuous elds and Hd is the Heaviside function centred at the discontinuity d . The in nitesimal strain eld can be derived by taking the symmetric gradient of the displacement eld ” = ∇s û + Hd ∇s ũ + d (ũ ⊗ n)s ; (2) where d is the Dirac delta distribution centred at the discontinuity. The presence of this Dirac delta distribution makes the strain eld at the discontinuity unbounded. Fig. 1. Body crossed by a discontinuity d . K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 137 The kinematic description of the displacement eld crossed by a discontinuity can be easily generalised to the problem of a volume crossed by m non-intersecting discontinuities u = û + m
Hi u˜i : (3) i=1 The strain eld can be derived in the same way as for the one discontinuity case. For intersecting discontinuities, additional terms are needed to describe the intersection. 2.2. Partition of unity concept Finite element shape functions, Ni , form partitions of unity since n
Ni = 1; (4) i=1 where n is the number of discrete points. Duarte and Oden [2] showed that a displacement eld, can be interpolated making use of these partitions of unity   k n

bij j  ; N i  ai + u= eld, such as the (5) j=1 i=1 where ai and bij are discrete values and j is an enhanced basis with k basis terms. When compared with ‘classical’ nite elements, it can be seen that with partition of unity the regular interpolation of the displacement eld can be enriched with enhanced terms. When the Heaviside function is chosen as enhanced basis, the displacement eld can be interpolated as u = Na + m
Hi Nbi ; (6) i=1 where a are the regular degrees of freedom and bi are the enhanced degrees of freedom related to crack i. Comparing with Eq. (3), it can be seen that û = Na; (7) u˜i = Nbi : (8) 2.3. Governing nite element equations The weak form of the virtual work equation without body forces reads:   s ∇  :  d =  · t dS; (9) S where  is taken from the set of admissible displacement variations and S is the outer surface where external tractions are applied. Using a Galerkin approach, the admissible displacement variations can be decomposed in the same manner as the actual displacement eld. Inserting the expression for 138 K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 multiple (in this case m) nonintersecting discontinuities, the virtual work equation can be rewritten in  m  m 

s s i (˜i ⊗ ni )s :  d Hi ∇ ˜i :  d + ∇ ˆ :  d + i=1 i=1 =  ˆ · t dS + S m 
˜i · t dS: (10) S i=1 Taking rst variations of ˆ (˜i =0) and then variations of ˜i (ˆ = 0), a set of m + 1 equations is obtained.   s ∇ ˆ :  = ˆ · t dS; (11) S  Hi ∇s ˜i :  d +  ˜i · ti di = 0 (i = 1; : : : ; m); (12) i where ti are the traction forces working at the discontinuity i . To obtain Eq. (11), the enhanced displacement eld, ũ, is assumed to be zero where essential boundary conditions are imposed. The discretised equations are obtained by inserting the discretised forms for the displacement and the strain eld.   T B  d = N T · t dS; (13) S  T H i B  d +  N T · ti di = 0 (i = 1; : : : ; m); (14) i In these equations, the continuum and the discontinuity response are completely split up. The continuum is assumed to remain elastic during the complete computation, so the stress rate is de ned as ˙ = C e ”˙ = C e (B ȧ + Hi B ḃi ); (15) where C e is the elastic material tensor. The traction rate is de ned as T˙i = DN ḃi ; (16) where D is the material tangent for the discontinuity. To obtain the linearized governing equations, the stress and traction rate are inserted in Eq. (13). After some mathematical manipulations, the linearized set of equations is given by     int; t    fa  Kaa Kab1 : : : Kabm  da   faext; t+dt                       int; t        Kb1 a Kb1 b1 : : : Kb1 bm   db1   0   fb1     (17) = −  .  .. .. ..   ..   ..  ; ..  ..         . . . . . .                            int; t dbm Kbm a Kbm b1 : : : Kbm bm 0 fbm K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 139 where  BTC eB d ; (18) Kabj =  Hj B T C e B d ; (19) K bj a =  Hj B T C e B d ; (20) K bi bj =  Hi Hj B T C e B d ; (21) K bj bj =  Kaa = T e Hj B C B d +  N T DN dj ; (22) j It can be seen from previous equations that when C e and D are symmetric, the global system of equations remains symmetric. Eq. (17) is written on element level. It is assumed that only one discontinuity j crosses the element but the element is in uenced by discontinuity i. 2.4. Numerical implementation The cohesive zone model based on partitions of unity allows discontinuities to propagate independent of a nite element discretization. For the propagation of a discontinuity, a criterion is evaluated in the elements in front of the crack tips. For mode I cracking, a discontinuity can be initiated perpendicular to the maximal principal stress direction when the tensile strength in the element is exceeded. More general, a yield surface can be used to decide whether the discontinuity should propagate or not. In this case, a bifurcation analysis may be necessary to compute the propagation direction. Furthermore, discontinuities can only grow at the end of a time step, so that no discontinuities are introduced in nonequilibrium states. As was stated in [7], this procedure ensures quadratic convergence of the Newton–Raphson process. Once a discontinuity has propagated, the nodes of the crossed element are enhanced. The support of the discontinuity tip remains restrained, so that the displacement jump at the tip is zero. During computation, only one discontinuity is allowed to cross an element. When a new discontinuity enters an element, which is already crossed by a discontinuity, the discontinuity is arrested. When two discontinuities reach the same elastic element, the discontinuities can link up. In order to assure an adequate integration of crossed elements, the integration scheme is changed in these elements. The integration scheme, proposed in [7], is used. 3. Cohesive zone model The starting point for the cohesive zone model used in this paper is the Rankine yield surface. When n and t represent the normal and tangential direction to the discontinuity, respectively, the 140 K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 Rankine yield surface can be written as f = T n − ft ; (23) where Tn is the normal component of the traction vector and ft is the tensile strength of the material. For modelling softening materials, the tensile strength is allowed to decrease during loading. The decrease is driven by the maximal plastic deformation reached.   ft0 pl ft = ft0 exp − (24) u : Gf n For interface elements, the plastic deformation and the updated stress can be found using a classical elasto-plastic computation. Because interface elements are present from the beginning of the computation, they also develop an elastic part. This elastic part should be minimised which is achieved by using a very high dummy elastic sti ness. For the methodology used here, discontinuities are introduced only when the tensile strength of the material is exceeded. As a consequence the deformation in the discontinuity is completely inelastic. As was shown in [3], the elasto-plastic model should evolve to a rigid plastic computation, where no elastic part exists. For the Rankine yield surface, a rigid plastic analysis is relative simple. The new position of the yield surface can be computed assuming that the jump is completely plastic. The new value for the normal traction is immediately found.   ft0 jump Tn = ft0 exp − (25) u Gf n and the tangential sti ness can be computed as   ft0 jump jump f2 u̇ n : un Ṫ n = − t0 exp − Gf Gf (26) For the Rankine yield surface, these relations can also be obtained starting from the elasto-plastic rate equations. For an associative plastic material, the elasto-plastic rate equations are given by   C e nnT C e e u̇; (27) Ṫ = C − T e n C n+h where n = @f=@T, C e is the elastic sti ness matrix and h is the softening modulus. For the Rankine yield surface, nT = (1; 0) so Eq. (27) simpli es to   E2 u̇ n Ṫ n = E − (28) E+h or  h u̇ n Ṫ n = 1 + h=E  The elastic deformations becomes negligible when elastic sti ness goes to in nity,   h u̇ n = hu̇ n Ṫ n = lim E →∞ 1 + h=E (29) (30) K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 141 and   2 ft0 ft0 h=− exp − un : Gf Gf (31) Eq. (30) is similar to the tangential sti ness found with the rigid plastic model and integration of Eq. (30) results in Eq. (25). For the modelling of cyclic behaviour, also the unloading behaviour should be studied. Assume that the inelastic deformation can be decomposed in a recoverable—damage—part and an unrecoverable— plastic—part. The index n is skipped for clarity: ujump = upl + ud = (1 − )ujump + ujump ; where is a model parameter. A fully damage model is found for model is obtained when → 0. The unloading is then obtained by   ft0 ft0 max Tn = max exp − (un − (1 − )unmax ) u un Gf n (32) → 1, while a fully plastic and the unloading sti ness,   ft0 ft0 max u̇ n : Ṫ n = max exp − un un Gf (33) (34) For the tangential direction, elastic behaviour is assumed. 4. Numerical examples The proposed model is used to simulate the tensile behaviour of a double edge-notched stone specimen. The geometry of the specimen is given in Fig. 2. In the experiments, the lower and upper boundary are glued to loading platens, so that no rotation of the boundaries is possible. The load is applied under displacement control. The opening displacement of both notches is averaged and is used as control signal for the displacement control. The numerical simulations are performed under load control with arc length control. For the arc length, the averaged opening displacement of the notches is used as control signal. At the upper boundary all degrees of freedom are constrained while at the lower boundary, no displacements Fig. 2. Geometry of the double edge notched specimen. 142 K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 372 elements 660 elements (a) (b) (c) Fig. 3. Mesh with (a) 660 elements and (b) 372 elements and (c) obtained crack paths. along the x-axis are allowed. The rotation of the specimen is prevented by expressing the symmetry condition along the y-axis. The steel loading plate for the upper boundary is also modelled. For the simulations, a quadratic triangular element was used. The behaviour of the continuum is assumed to remain elastic. The elastic material properties are: Young’s modulus E = 40 GPa and a Poisson’s ratio  = 0:29. The behaviour of the discontinuity is described by the Rankine model where a tensile strength ft0 = 6:5 MPa and a fracture energy Gf = 0:05 N=mm is chosen. Two meshes with 372 elements and 660 elements are analysed. Fig. 3 shows used meshes and the obtained crack path for the two meshes. As can be seen the crack paths for di erent meshes coincide remarkably well. The numerically simulated crack path also shows good agreement with the experimentally observed crack path. Apart from local information, i.e. the monitoring of the crack path, also global information is recorded. The averaged opening of the notches d, is plotted against the applied load in Fig. 4. Comparing the obtained load deformation curve for the mesh with 372 elements with the curve for the mesh with 660 elements shows that there is only a small di erence near the end of the post peak branch. The predicted peak load is for both simulations identical. Although small di erences, the cohesive zone model based on partition of unity seems to simulate fracture independent of the discretization. Al further simulations are carried out with the coarser mesh. Fig. 5 shows the load–deformation curve of a specimen under cyclic loading. The parameter controlling the amount of recoverable deformation is varied from 0 to 1. Fig. 5a shows the load– deformation curve for = 0:7. A fully plastic model (Fig. 5b) is obtained for = 0:0 while the fully damage model (Fig. 5c) is found for = 1:0. Finally, model parameters are tted to give a good correspondence with experiments [1]. The following set is chosen: ft0 = 6:5 MPa, Gf = 0:05 N=mm and = 0:7. The resulting load deformation curve is given in Fig. 6. The value of shows that the behaviour of the stone is more damage related than plasticity related. Although for a good representation of unloading behaviour, the contribution of combined plasticity and damage is crucial. K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 143 372 elements 660 elements 2000 P (N) 1500 1000 500 0 0 0.005 0.01 0.015 d (mm) 2000 2000 1500 1500 P(N) P(N) Fig. 4. Load–deformation curve for DEN tensile test. 1000 1000 500 500 0 (a) 0 0 0.005 0.01 0.015 d (mm) (b) 0 0.005 0.01 0.015 d (mm) 2000 P(N) 1500 1000 500 0 (c) 0 0.005 Fig. 5. Load–deformation curve for (a) 0.01 0.015 d (mm) = 0:7, (b) = 0:0 and (c) = 1:0. 5. Conclusions The numerical simulations of a double-edge notched specimen under tensile loading are presented. The discrete approach for modelling fracture is followed, introducing discontinuous terms in the 144 K. De Proft et al. / Journal of Computational and Applied Mathematics 168 (2004) 135 – 144 Numerical Experimental 2000 P(N) 1500 1000 500 0 0 0.005 0.01 0.015 0.02 d (mm) Fig. 6. Load–deformation curve for experimental results and numerical simulation. displacement eld. The partition of unity concept is used for the numerical implementation. It was shown that this modelling technique yields results that are independent of the discretization. The behaviour of the discontinuities is modelled with a mode I model based on the Rankine yield surface. Since discontinuities are only embedded when the tensile strength of the material is reached, the deformations in the discontinuity are completely inelastic. A rigid plastic or rigid damage model is necessary to describe the behaviour. The model for cyclic loading used in this paper simply splits the inelastic behaviour into a recoverable and an unrecoverable part. Simulations showed that this simple model is capable to simulate the behaviour of limestone under cyclic loading. In the future, more complicated material models will be used and implemented so that the modelling is extended to mixed mode cracking. References [1] K. De Proft, L.J. Sluys, W.P. De Wilde, Combined experimental–numerical study to monotonic and cyclic behaviour of limestone, in: C.A. Brebbia, S.I. Nishida (Eds.), Damage and Fracture Mechanics, Wessex Institute of Technology, UK, 2003. [2] C.A. Duarte, J.T. Oden, H–p clouds and h–p meshless method, Numer. Methods Partial Di erential Equations 12 (1996) 673–705. [3] J. Oliver, On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations, Internat. J. Solids Struct. 37 (2000) 7207–7229. [4] R.H.J. Peerlings, R. De Borst, W.A.M. Brekelmans, M.G.D. Geers, Gradient-enhanced damage modelling of concrete fracture, Mech. Cohesive-Frictional Materials 3 (1998) 323–342. [5] G. Pijaudier-Cabot, Z. Bazant, Nonlocal damage theory, ASCE J. Eng. Mech. 113 (1987) 1512–1533. [6] J.C. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids, Comput. Mech. 12 (1993) 277–296. [7] G.N. Wells, L.J. Sluys, A new method for modelling cohesive cracks using nite elements, Internat. J. Numer. Methods Eng. 50 (2001) 2667–2682. [8] X.P. Xu, A. Needleman, Numerical simulations of fast crack growth in brittle solids, A new method for modelling cohesive cracks using nite elements, J. Mech. Phys. Solids 42 (1994) 1397–1434.