Bifurcation phenomena in cohesive crack propagation

Compr ~ er s Printed & Slr uct ur es in Great Vol. 44, No. l/ 2, pi. 55-62, 1992 C04 5-194 9192 Britain. 8 1992 $5.00 Per gamon + 0.00 Press Ltd BIFURCATION PHENOMENA IN COHESIVE CRACK PROPAGATION s. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIH V ALEN T E Dipartimento di Ingegneria Strutturale, Politecnico di Torino, 10129 Turin, Italy Abstract-The ultimate behaviour of several materials, such as concrete, rocks, bricks, fibre-reinforced composites, is characterized by the localization of a non-linear zone within a very narrow band whilst the rest of the material retains its linear behaviour. The cohesive model describes this band by means of a fictitious crack carrying softening stresses which are decreasing functions of the displacement discontinuity. This oaoer describes an extension of such a model which enables the possible bifurcations of the e&ilibr&n-path to be detected and analysed. et al. [17l in the flexural behaviour of an elastic-perfectly plastic beam, by Nemat-Nasser et al. [18] in the interaction of two linear elastic cracks, and, more recently, by Rots et al. [19] in the analysis of the direct tension test by means of the smeared crack model. This paper describes an extension of the cohesive crack model which enables the possible bifurcations of the equilibrium path to be detected and analysed. 1 . I N T RODU CT I ON The ultimate behaviour of several materials, such as concrete, rocks, bricks, fibre-reinforced composites, is characterized by the localization of a non-linear zone within a very narrow band, whilst the rest of the material retains its linear behaviour. The cohesive model describes this band as a fictitious crack (or process zone) where the material, albeit damaged, can still transfer stresses. The stresses involved are described as decreasing (strain-softening) functions of the relative displacements of the opposing surfaces of the fictitious crack (displacement discontinuity). The process zone originates perpendicularly to the principal tensile stress when the latter reaches its ultimate value u,. The point at which this condition occurs is called the fictitious crack tip: it represents the boundary between the integral and the damaged material. The point at which the displacement discontinuity reaches the critical value w,, beyond which no stresses are transferred, is called the real crack tip; this point divides the stress-free crack from the process zone (Fig. 1). According to the foregoing assumptions, there are no stress singularities in this model. The cohesive model was initially proposed by Barenblatt [l] and, independently, by Dugdale [2]. Later on, this model was reconsidered by Bilby et al. [3], Willis [4] and Rice [S]. Subsequently, it was reproposed by Wnuk [6] under the name of ‘final stretch model’ and by Hillerborg et al. [7j, under the name of ‘fictitious crack model’. Recently, using the length of the fictitious crack as a control variable in numerical simulations (‘fictitious crack length control scheme’), Carpinteri et al. have provided an explanation of phenomena such as catastrophic collapse (or snap-back) and the ductile to brittle failure transition governed by the brittleness number s,, which depends on the material properties and on the structural size. In particular, [8,9] deal with mode I deflection problems, while in [lo-161 the method has been extended to mixed mode problems. On the other hand, bifurcation phenomena due to strain-softening were originally described by Maier 2 . T H E COH ESI V E CRACK M ODEL I N M lX ED M ODE CON DI T I ON S AN D CRACK LEN GT H CON T ROL By neglecting tangential cohesive stresses and assuming, for normal cohesive stresses, a linear softening law, we can write (Fig. 2) 6,+ = -a; 6,+=-o- = 0, = a,(1 - w/w,), forO<w<w,andi>O (la) for w > w,, (lb) c =a,=o, where u, stands for the ultimate tensile strength of the material, w is the normal component of the displacement discontinuity, w, is the critical value of w, the superscript ‘+’ denotes the positive side of the crack (Fig. 1) and the superscript ‘-’ the negative one, the dot denotes the derivation with respect to time (evolutionary problem). The energy dissipated in the process which gives rise to a unit of stress-free surface is called fracture energy, G,, and turns out to be a property of the material G,= S’< u=dw. s0 (2) Having assumed the constitutive law (la), we get G,= 6, wJ2. The crack irreversibility condition can be written as follows: i 20. 55 (3) S. VALENTE 56 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA Fictitious crack or process zone Fig. 1. The cohesive model represents the process zone in the form of a fictitious crack. the solution obtained can be regarded as acceptable, it is necessary to calculate the values of w in the process zone. If at any point w > w,, then, based on the constitutive law (1b), the cohesive stresses vanish and the process zone decreases. Similarly, if at any point i c 0, then, based on the constitutive law (3), the process zone is modified locally so as to enable (4) rigid unloading (a = 0) to take place. If a modififJ/J,< QU, cation occurs in the process zone, then it becomes necessary to reassemble the matrix (K - C) and to where aP, stands for the principal tensile stress. repeat the procedure starting from eqns (5). The By subdividing the domain into a finite number of iteration chain ends when the solution obtained does elements and employing the principle of virtual work, the following system of linear equations can be not bring about any variation in the process zone. At this point it is possible to have the fictitious obtained [lo-12, 161 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA crack propagate, by a pre-determined length, perpendicularly to the direction of the principal (K-C)u=f,+lf,, (5) Equation (3) entails that, if the crack tends to close, rigid unloading takes place which preserves the opening displacement (G = 0), whilst eqn (la) is no longer applicable. In order to prevent new cracks from forming, at each point in the domain, the following condition has to be fulfilled where K = stiffness matrix (n x n), C = symmetrical strain-softening matrix, K - C = effective stiffness matrix (n x n), u = vector of the n unknown nodal displacements, f, = load vector depending on the crack length, f, = external load vector and L = external load multiplier. The (n + 1)th unknown 1 is determined by setting that at the fictitious crack tip (the centre of gravity of the dashed element in Fig. 3) the principal tensile stress reaches the value Q, ~~-(6,+~~)uu+~xu~-~*=0. I t (6) Since ox, au, T can be expressed as a function of u, eqns (5) and (6) make up a system of (n + 1) equations with (a + 1) unknowns, which must be solved at each step in the growth of the crack. Before Opening (w) Fig. 2. Stress vs crack opening displacement constitutive I_._ raw. Bifurcation phenomena in cohesive crack propagation 57 This initial behaviour is common to all problems, i.e. it does not depend on the geometrical and mechanical characteristics of the structure. As the crack grows, two different situations may occur a, > 0, wmx= WC3 for zyxwvutsrqponmlkjihgfedcbaZYXWV 63) %*x= WC9 for Fig. 3. Finite element rosette at the fictitious crack tip. tensile stress, that is to say, in the direction given by 9 = arctan [h/(a, - a,)]/2. (7) In order to follow the new position of the fictitious crack tip, the finite element rosette (Fig. 3) translates and rotates and a new mesh is automatically generated by the computer. Therefore, the solution process for the subsequent crack growth step can start. In other words, the numerical simulations are controlled by means of the length of the fictitious crack, whilst all the other parameters (u, 2, process zone) follow according to eqns (I), (3) and (5)-(7). It should be noted that when rr, = a,,, Mohr’s circle has degenerated to a point; hence it is not possible to determine from eqn (7) the subsequent direction of growth of the crack. This technique has been called ‘fictitious crack length control scheme’ [IO-12, 161. 3. SOLUTION OF A SINGLE CRACK GROWTH STEP During the evolution of the crack, the linear system (5) may go through different stages. Initially, the process zone is small and therefore the matrix C (negative definite) is small compared with matrix K (positive definite). As a result, matrix (K-C) is positive definite and all its eigenvalues ai are real and positive (0 < a, < a2. . ’ < an). 80 I- = 0 and a2 > 0, (9) where w,,,,, is the maximum value of w in the process zone. By way of exemplification, numerical simulations concerning the three-point bending test (Fig. 4), performed by taking into account a wide range of geometrical and mechanical ratios, have shown that in this case situation (8) invariably occurs. It involves a reduction in the process zone, with the ensuing increase of a,, so that a, never becomes negative. Thus it is possible to simulate the test through the end, until structural collapse, by always using (K - C) positive definite matrices. System (5) always turns out to be well conditioned, even in the event of catastrophic collapse (or snap-back) taking place. If we now take into consideration mixed mode problems, for the geometrical and mechanical ratios analysed in [IO-16], we find that situation (9) invariably occurs. Between situations (8) and (9) there is a significant difference: the former, once it has been reached, tends to persist during the following steps. The latter, instead, represents a critical point which can be easily overcome in the course of the numerica simulation. Considering the importance of this situation, which may represent the start of bifurcation paths, it is advisable, as soon as a, becomes negative, to set back the fictitious crack tip, in order to approximate the condition t(, = 0 as closely as possible. Once this point has been analysed, going on with the growth of the crack, we get ai c 0. Matrix (K-C) is now invertible again, although it is no longer positive definite. Thus, it can be concluded that during crack growth the matrix (K-C) may, in some cases, become singular. b 9 a1 t Fig. 4. Layout of the three-point bending test. S. VALENTE 58 3.1. The (K - C) matrix is non-singular In this case it is possible to solve the system (5) for two right-hand sides (K - C)u, = f,, (W (K-C)u,=f,. (1Ob) By denoting with 8, the stress at the fictitious crack tip (centre of gravity of the element dashed in Fig. 3), corresponding to displacements u,, and with ur the stress corresponding to displacements ur, we can write a=[a,,u,,r]r=b,+au*. (11) By substituting eqn (11) into (6) we get a secondorder equation in 1, which, if it has any solution, makes it possible to determine 1. Thus, the solution, in terms of displacements, can be written u=u,+au,. Equation (18) represents a system of n independent equations. Since we have assumed that all the eigenvalues are distinct, when K - C becomes singular, just one eigenvalue vanishes, i.e. tlk = 0. Writing eqn (18) for i = k, noting that the first term vanishes and solving with respect to 1, we get 1 = -(Gf,)/(qkTfr). (19) By solving eqn (18) in respect of vi, we get vi = q’(f, + J.f,)/a,, for i = 1,2, . . . , n, i #k. (20) The state of stress u at the fictitious crack tip can be written as a function of the contributions ui, pertaining to the eigenvector qi, in the following form u =[a,,u,,7]T= &Ii. i=l (21) (12) By substituting eqn (20) into (21) it is possible to make the state of stress u depend solely on the 3.2. The (K - C) matrix is singular unknown ok. In order to impose the fictitious crack Assuming that the eigenvalues xi of K-C are propagation conditions it is now sufficient to substidistinct, the corresponding eigenvectors q, are unique tute eqn (21) into (6) so as to obtain a second-order (within scalar multipliers). Then we have equation in vk, which, if it has any real solution, makes it possible to determine v,. By substituting the generalized displacements v (K - C)qi = aiqi. (13) into eqn (16), we obtain the solution expressed in terms of nodal displacements u. Since K - C is symmetrical, the matrix Q, defined as If $f2 = 0, it is necessary to check zyxwvutsrqponmlkjih qkTf, . If the latter follows: zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA is not zero, it becomes necessary to move back the fictitious crack tip and then make it grow at a slower (14) Q=h,,qz,...,9nl, rate. An infinite value of 1, in fact, is not compatible with eqn (4). If, instead, we find turns out to be orthonormal. Hence it is possible to normalize the eigenvectors so as to obtain q;r, = df* = 0, q$=O for i #j and q,qT= 1. (15) Let us now consider new displacements v related to the nodal displacements by the following linear transformation (singular value decomposition) u=Qv. (16) v is called the generalized displacement vector. By substituting eqn (16) into (5) and pre-multiplying the latter by Qr we get Qr(K - C)Qv = QT(f, + lf2). (17) By pre-multiplying eqn (13) by $, using eqns (15), and substituting it into eqn (17), we can write cl++= q’(f, + If,), for i = 1,2, . . . , n. (18) (22) it is useful to re-write eqn (18) for i = k ov, = 0. (23) Regardless of the v, value selected, by substituting vi, given by eqns (20), and uk into eqn (21), and then eqn (21) into eqn (6), we obtain a second-order equation in 1, which, if it has any real solution, makes it possible to determine 1. Equation (23) represents a bifurcation point, as two solutions are possible, i.e. vk = 0 and vk # 0. In order to reduce the computational effort in the calculation of the eigenvalues and eigenvectors, it is useful to reduce the order of the linear system by means of static condensation. This problem will be taken up again below, by providing a numerical example. The case of multiple eigenvalues is not discussed in this paper. Bifurcation phenomena in cohesive crack propagation 59 4b 2b , ,0.4b substructure B b L Fig. 5. Layout of the four-point shear test. zyxwvutsrqponmlkjihgfedcbaZYXWVU 4 FOUR-POINT SHEAR TEST where t stands for the thickness of the specimen, while the following dimensionless mechanical parameters have been taken into account As an example, let us consider a numerical simulation of the four-point shear test on a specimen with a double notch (Figs 5 and 6). This is a structure with sE= Gf/a,b = 0.00025, v = 0.1, two axes of symmetry, in which the external loads 6, = a,,‘E = 0.741 x IO-*. (25) and the support reactions make up a system of balanced external forces having a polar symmetry In the numerical simulation of this test, performed as around the centre of the specimen. described in the previous paragraphs, it is possible to According to the polar s~rnet~, the appli~tion identify the following stages. points of loads i;l and zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA Fz undergo the same displacement 6, = 6, (Fig. 5). Numerical and experimental 4.1. Initial stage @ ositivedejkite K - C) results, relating to different values of the geometrical This first stage involves the propagation of two and mechanical ratios of the specimen, are given cracks growing symmetrically from the two notches. in [lo-13, 161. Because of the symmetry, one equation only, and In the example described below, the following namely eqn (6), can impose the propagation congeometrical ratios have been considered (Fig. 5) ditions of both fictitious cracks. The point, which describes the state of structure, follows the portion I/b -4, c/b = 0.8, a,/b = 0.2, t/b = 1, (24) O-A of the curves illustrated in Figs 8-l 1. 4.2. The K -C matrix is singular In this case (CQ=O), the eigenvector qk, by neglecting a rigid motion, is seen to have a polar anti-symmetry (Fig. 7). This means that the eigen____I_----_-,----___‘~‘~~ I I I I I I zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONM I zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPO I I I I I I I I I I I Fig. 6. Example of a finite element mesh. Fig. 7. Eigenvector related to the vanishing eigenvaiue. S. VALENTE .g 0 02 Dimensionless deflection, 104 (5 6, + 6,) / 6 b Fig. 8. Total load vs the displacement of its application point. C I I II,,1 I I I -10 8 -6 -7 -6 -5 -4 -3 -2 -1 Dimensionless 123456 deflections, zyxwvutsrqponmlkjihgfedc 104 6, / b , 104 62 / b Fig. 9. F,, 6, and F2. 6, diagrams. The point describing the state of the structure follows displacements, in two polar-symmetrical points, are the portion A-C of the curves plotted in Figs 8-l 1. the same, in terms of both module and sign. Since the It should be noted that in the symmetrical case 6, is external loads, applied on two polar-symmetrical always positive (downward); following the loss of points, share the same module but have opposite symmetry, it changes sign. The differences in the signs, we find s,‘f, = 0. Since this is the first eigenvalue to vanish, the cracking path up to this point is cracking trajectories observed in the two cases, on the unique and polar-symmetrical, and hence df, = 0. other hand, are negligible; the maximum loads are nearly the same. Equation (22) being satisfied, we have met a bifurAs shown in Fig. 8, loss of symmetry in the crack cation point where two solutions are possible, i.e. propagation can cause snap-back instability. In order vk = 0 and vk # 0. 4.2.1. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA Symmetrical propagation. Solution vk = 0 to obtain a stable crack, the crack mouth opening displacement (CMOD) has to be chosen as the drivrules out any contribution on the part of the eigenvalue ok to the displacement field. For all subsequent ing parameter. As shown in Fig. 10, in fact, CMOD crack growth steps, the polar-symmetry condition is is a monotonic increasing function of time, during the irreversible fracture propagation process. In the eximposed and therefore the (K - C) matrix immediately returns non-singular. After a small growth of perimental test, a small non-symmetric imperfection has to be introduced in order to determine the main the crack, the (K - C) matrix becomes singular again crack (the lower one, in this numerical example) (CQ= 0). The product df2 is different from zero, there before applying the load. It is therefore possible to are no further bifurcations and the analysis is carried apply the displacement tranducer on the main crack on with (K -C) no longer positive definite without mouth, in order to close the control loop of the encountering any additional singularity. The point testing machine on such a displacement. that describes the state of the structure follows the In this way, both crack patterns, symmetric and portion A-B of the curves plotted in Figs 8-11. 4.2.2. Non-sy mmetrical propagation. Without loss non-symmetric, are stable. Therefore, a question of generality, we can normalize the eigenvector Q as arises: which is more stable? According to the criillustrated in Fig. 7, i.e. the eigen-displacement field terion proposed by Nemat-Nasser [20], which is deshows an opening discontinuity along the lower crack rived from Gibbs’ statement of the second law of and a closing discontinuity, with the same absolute thermodynamics, the load being the same, the state value, along the upper crack. involving a smaller quantity of stored elastic energy In these hypotheses, regardless of the absolute is more stable. In order to apply this criterion, the value of v,, the assumption of v, > 0 results in It < 0 elastic strain energy stored in the structure has been along the upper crack. In order to follow the trend of plotted in Fig. 11, as a function of thm; >tal load. This the evolutionary process towards a loss of symmetry, diagram shows that the non-symmetric solution turns it is therefore necessary to modify the constitutive law out to be slightly more stable. along the upper crack (@ = 0). This modification The amount of elastic stored is given by the work makes the (K-C) matrix invertible. But then the done by the external loads on the structure minus the matrix becomes singular again (CQ= 0) following a energy dissipated along the cracks. If we examine the small growth of the lower crack. The product q:f, is evolution of the total load vs the displacement of its different from zero, no further bifurcation occurs and application point (Fig. 8), we find that the area lying the analysis is carried on with (K-C) no longer under the O-A-B curve (symmetrical case) is greater positive definite without meeting any further singuthan the one under the &A-C curve (non-symmetrilarity. During this stage, the tensile stress at the tip cal case). This area is proportional to the work of the upper fictitious crack decreases monotonically. performed by the external loads, which therefore is Bifurcation phenomena in cohesive crack propagation 61 w n tf + 1 2 + co.75 $ -0 ln 0.5 P 6 ‘g 0.25 E” a 0 123456789 Dimensionless crack mouth opening disp., 104 CMOD / b Fig. 10. Total load vs crack mouth opening displacement (CMOD). author wishes to express his gratitude to Professor Albert0 Carpinteri for the stimulating discussions on the topics dealt with in this paper. This study was carried out with the financial support of the Ministry of Education (M.P.I.) and the National Research Council (C.N.R.). Acknowledgements-The REFERENCES I. G. I. Barenblatt, The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses. Axially-symmetric cracks. J. Appl. Math. Mech. 23, Dimensionless load, (F, + F,) / uub t 622-636 (1959). 2. D. S. Dugdale, Yielding of steel sheets containing slits. J. Mech. Phys. Solids 8, 100-104 (1960). Fig. I I. Elastic strain energy stored in the structure vs total 3. B. A. Bilby, A. H. Cottrell and K. H. Swinden, The load. spread of plastic yield from a notch. Proc. R. Sot. A272, 304-314 (1963). 4. J. R. Willis, A comparison of the fracture criteria of Griffith and Barenblatt. J. Mech. Phys. Solids 15, greater in the symmetrical case. But the amounts of 151-162 (1967). energy dissipated along the two cracks are greater in 5. J. R. Rice, A path independent integral and the approxithe symmetrical case too. This is the reason why the mate analysis of strain concentration by notches and elastic energy stored is comparable in the two cases, cracks. J. Appl. Mech. 35, 379.-386 (1968). even if the two terms of the balance are much 6. M. P. Wnuk, Quasi-static extension of a tensile crack contained in a viscoelasticplastic solid. J. Appl. Mech. different. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA 41, 234-242 5. CONCLUSIONS 1. An extension of the cohesive crack model, based on the singular value decomposition of the matrix (K - C) has been presented. It makes it possible to detect and analyse the possible bifurcations of the equilibrium path. 2. Bifurcation can cause snap-back instability. 3. If two different states of cracking are both stable, the criterion proposed by Nemat-Nasser [20] can be used in order to determine which state is more stable. 4. Applying such a criterion in a numerical example (four-point shear test), it appears that the non-symmetrical cracking state is somewhat more stable (just slightly less elastic energy) than the symmetrical state. The maximum loads, in the two cases, are nearly the same. (1974). 7. A. Hillerborg, M. Modeer and P. E. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement Concrere Rex 6, 773-782 (1976). 8. A. Carpinteri, Interpretation of the Griffith instability as a bifurcation of the global equilibrium. NATO Advanced Research Workshop on Application of Fracture Mechanics lo Cementitio& Cdhposites, -Evanston (Illinois). Seotember 4-7. 1984: (Edited bv S. P. Shah). pp. 284y318. Martinus Nijhoff (1985). 9. A. Carpinteri, Softening and snap-back instability in cohesive solids. Int. J. Numer. Meth. Engng 28, 1521-1537 (1989). 10. A. Carpinteri and S. Valente, Size-scale transition from ductile to brittle failure: a dimensional analysis approach. In Cracking and Damage, (Mited by J. Mazars and Z. P. Bazant), pp. 477-490. Elsevier (1988). 11. A. Carpinteri and S. Valente, Mixed mode crack propagation and ductile-brittle transition by varying structural size. Atti del IX Congress0 Nazionale AIMETA, pp. 107-I 10, Bari, Italy (1988). II 62 S. VALENTE 12. A. Carpinteri, S. Valente and P. Bocca, Mixed mode 16. P. Bocca, A. Carpinteri and S. Valente, Mixed mode cohesive crack propagation. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA Seventh in!ernational fracture of concrete. Int. J. So&& Structures 21, Conference on Fracture (ICF7), Houston (Texas), 1139-1153 (1991). pp. 2243-2257 (1989). 17. G. Maier, A. Zevelani and J. C. Dotreppe, Equilibrium 13. A. Carpinteri and S. Valente, Design and analysis branching due to flexural softening. J. Engng M ech. of orthotropic composite materials through a mixed Div., ASCE 99, 897-901 (1973). mode cohesive crack simulation. Proc. 3rd Europ. 18. S. Nemat-Nasser, Y. Sumi and L. M. Keer, Unstable Conf. on Composite M aterials, pp. 309-314. Elsevier growth of tension crack in brittle solids: stable and un(1989). stable bifurcations, snap-through, and imperfection sen14. P. Bocca, A. Carpinteri and S. Valente, Evaluation of sitivity. Inr. J. Solids Structures 16, 1017-1035 (1980). concrete fracture energy through a pull-out testing 19. J. G. Rots, D. A. Hordijk and R. de Borst, Numerical procedure. Proc. Int. Conf. on Recent Developments on simulation of concrete fracture in direct tension. Proc. the Fracture of Concrete and Rock, pp. 347-356, Cardiff 4th Int. Conf. on Numerical M ethods in Fracture M ech(1989). anics, pp. 457-471. Pineridge Press, Swansea (1987). 15. P. Bocca, A. Carpinteri and S. Valente, Size effects in 20. S. Nemat-Nasser, The second law of thermodynamics the mxed mode crack propagation: softening and snapon non-collinear crack growth. Proc. of the Third back analysis. Engng Fracture M ech. 35, 159- 170 ASCEIEM D Specialty Conference, pp. 4499452. (1990). Austin, Texas (1979).