Microsoft Word - numero 18 articolo 1.doc G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 5 The variational theory of fracture: diffuse cohesive energy and elastic-plastic rupture Gianpietro Del Piero Dipartimento di Ingegneria, Università di Ferrara, Ferrara (Italy) dlpgpt@unife.it Giovanni Lancioni Dipartimento di Architettura, Costruzioni e Strutture, Università Politecnica delle Marche, Ancona (Italy) Riccardo March Istituto per le Applicazioni del Calcolo, CNR, Roma (Italy) ABSTRACT. This communication anticipates some results of a work in progress [1], addressed to explore the efficiency of the diffuse cohesive energy model for describing the phenomena of fracture and yielding. A first local model is partially successful, but fails to reproduce the strain softening regime. A more robust non-local model, obtained by adding an energy term depending on the deformation gradient, describes many typical features of the inelastic response observed in experiments, including strain localization and necking. Fracture occurs as the result of extreme strain localization. The model predicts different fracture modes, brittle and ductile, depending on the analytical form of the cohesive energy function. SOMMARIO. Nella presente comunicazione si anticipano alcuni risultati del lavoro [1], in preparazione, in cui si analizza l’efficienza del modello dell’energia coesiva diffusa nel descrivere i fenomeni di plasticizzazione e rottura. Un primo modello locale dà risultati parzialmente soddisfacenti, ma si rivela incapace di descrivere il regime di strain softening. Un più robusto modello non locale, ottenuto aggiungendo un termine energetico dipendente dal gradiente della deformazione, riesce a descrivere molte peculiarità della risposta anelastica osservate sperimentalmente, incluse la localizzazione della deformazione e il necking. La frattura avviene per estrema localizzazione. Il modello prevede la possibilità di diversi modi di frattura, duttile o fragile, a seconda della forma analitica ipotizzata per l’energia coesiva. KEYWORDS. Fracture mechanics; Energy methods; Elastic-plastic rupture; Variational theory of fracture. INTRODUCTION leading idea in Fracture Mechanics is that the growth of a fracture surface is decided by the competition between the strain energy lost in, and the amount of fracture energy required to, the creation of a new fracture surface [2]. However, it was soon clear that this idea is appropriate to the description of brittle fracture, typically exhibited in large-size bodies, but is unfit to describe the ductile fracture modes which prevail in small-size bodies. This size effect is captured by the cohesive energy models [3, 4]. A http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 6 Other basic problems, such as the formation of a fracture in an initially unfractured body, were neglected for long time, and ony recently came to the attention of the scientific community. The use of variational techniques for energy minimization was introduced first for Griffith’s fracture model [5], and extended later to cohesive energy models [6, 7]. From the engineering side, it became clear that an appropriate choice of the analytic form of the cohesive energy may lead to a unified description of the phenomena of yielding, damage, and fracture [8, 9]. A very recent research trend consists in questioning the surfacic character of the fracture energy, and in exploring the alternative possibility of an energy diffused over the volume. An example is the paper [10], in which, in the context of a damage model, a rather detailed description of the process zone preceding the final rupture is obtained. In this communication a diffuse cohesive energy is considered, and fracture is regarded as the extreme stage of the localization of inelastic deformation. In a one-dimensional context, we assume that at every point of the bar the energy density is the sum of an elastic and a cohesive part. An incremental energy minimization is performed, under the only assumption that the cohesive energy is totally dissipative. With this simple model, a rather accurate description of the inelastic response of the bar, from the onset of the inelastic regime up to rupture, is obtained. The strain-softening case is not adequately described, but this inconvenience is repaired with the introduction of a non-local energy term of the gradient type. The traditional elements of elastic-plastic response, such as the yield condition, the flow law, the hardening law, the elastic unloading, come out of the model as necessary conditions for an energy minimum. We also find that a convex shape of the cohesive energy favorizes a uniform distribution of inelastic strain and a work-hardening response, while a concave shape is responsible of the localization of the inelastic strain, which, in turn, produces the phenomena of strain softening and necking. Rupture takes place when the localization becomes extreme. If the cohesive energy is approximated by a piecewise polynomial function, its expression becomes a function of a small number of material parameters. The first results of a series of numerical simulations show that an appropriate choice of these parameters provides a good approximation of the experimental response curves. This seems to confirm the efficiency and flexibility of the proposed model. THE LOCAL MODEL onsider a bar of length l , with constant cross section, free of external loads, and subject to the axial displacements u(0) = 0 , u(l) = l , (1) at the endpoints x  0 and x  l. The bar’s deformation is measured by the derivative uof the axial displacement u. We assume that, at every point x of the bar’s axis, the deformation u(x) can be split into the sum of an elastic part (x) and of an inelastic part (x): u(x) = (x) + (x) , (2) and that the total energy of the bar has the form   l dxxxwE 0 )))(())(((),(  (3) where w and  are the volume densities of the elastic strain energy and of the cohesive energy, respectively. We also assume that w can be recovered, while  is totally dissipated. The last assumption requires that the cohesive power ((x,t)) ),( tx be non-negative for all x and at all instants t, and if we assume that  is an increasing function of  we get the dissipation inequality 0),( tx (4) to be satisfied for all x and at all t. In the spirit of the Calculus of Variations, the equilibrium configurations of the bar are identified with the stationary points of the energy 0))())((')())(('(),,,( 0   dxxxxxwE l  (5) C http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 7 with the inequality due to the presence of the dissipation inequality (4), which requires the restriction  (x)  0 (6) on the perturbations . The analysis of the first variation leads to the following characterization of the equilibrium configurations [DLM] :  The elastic deformation (x) and the axial force   = w((x)) (7) are constant all over the bar,  The axial force cannot exceed a limit depending on the current inelastic deformation:   ((x)) . (8)  The equilibrium condition (8) is in fact a yield condition. It is remarkable that it has not been postulated, but deduced from the variational procedure. The energy minimization is insufficient to determine the evolution of the deformation under a loading process  = (t). To do this, it is necessary to formulate an incremental equilibrium problem, in which at any time t the deformations (t) and (x,t) are supposed to be known, and the unknowns are the deformation increments ),(),( txt   , produced by a given load increment )(t . One considers the expansion )())(),(())(),(())(),(())(),(( 22 2 1  ottEttEttEttE   (9) and minimizes for sufficiently small to legitimate the truncation at the second-order term. The term of order zero being known, a first-order approximation involves the minimization of the function  dxxttxttlttE l ))())(),(('()()())(),(( 0     (10) In it, everything is known except )(x . Then E is an affine function of  , and a proper use of the dissipation inequality leads to the Kuhn-Tucker conditions 0)( x , 0))(('  x , 0)()))(('(  xx   (11) as necessary conditions for a minimum. The first two conditions are (4) and (8), respectively. The third condition, the complementarity condition, states that when  is positive the force  must lie on the border of the yield surface, and that when  lies at the interior of the yield surface then  must be zero. Again, it is remarkable that these distinctive aspects of plastic response have not been postulated, but deduced from incremental energy minimization. Conditions (11) are still insufficient to determine the evolution of the deformation. This is done by minimizing the second-order term of (9) dxxtxtxttlttE l ))()()()),((''()()())(),(( 0     (12) which is a quadratic function of  . The new minimization provides a second set of Kuhn-Tucker conditions 0)( x , 0)())((''   xx , 0)())())((''(  xxx   , x (13) to be satisfied at the portion  of the bar at which (8) is satisfied as an equality. Indeed, we already know that  = 0 out of . In (13), the first condition is again the dissipation inequality (4). The second condition provides a relation between the increments of the force and of the inelastic deformation. By the complementarity condition (13)3 , inequality (13)2 is satisfied as an equality when  > 0. This equality provides the flow rule for the inelastic strain rate in the regime of inelastic deformation. Condition (13)3 is the consistency condition, which says that the inelastic deformation cannot increase when a http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 8 point x tends to re-enter the elastic range (8). This condition determines the property of elastic unloading. Again, these typical properties of plastic response are not assumed, but come as results of the second-order incremental minimization. Together with the necessary condition  ((x))  0 (14) imposed by the non-negativeness of the second variation, conditions (11) and (13) determine the solution ( ,  ) of the incremental problem. Here we summarize some properties of the continuations. For a more detailed analysis we refer to the main paper [1].  If at the time t = 0 we start with (0)=0 from the natural configuration (x) = (x) = 0, from the constitutive assumption (0) > 0 we have that, initially, inequality (8) is strict. Then  is zero, and a purely elastic evolution occurs. This elastic regime ends when, with growing , the force  reaches the limit value  =(0), corresponding to the boundary of the elastic range. This event marks the onset of the inelastic regime.  The response in the inelastic regime strongly depends on the sign of (0). If (0) is positive, the deformation is homogeneous, that is, both  and  are constant all over the bar. If 0 , the inelastic deformation is given by the solution of the differential equation )( ))()((''))(('' ))()(('' )( t ttwt ttw t         (15) joined with an appropriate initial condition. The slope of the force-elongation response curve is given by )( ))()((''))(('' ))()((''))(('' )( t ttwt ttwt t         (16) The slope is positive, that is, the inelastic regime starts with a work-hardening response. If 0 , elastic unloading takes place.  If (0) is zero, the initial response is perfectly plastic. The force is constant in time, 0)( t , the total inelastic deformation rate is equal to the total elongation rate )(tl , while the punctual distribution of )(t remains undetermined.  Finally, if (0) is negative the necessary condition for equilibrium (14) is violated. The deformation tends to concentrate on a very small portion of the bar, and the slope )(/)( tt   tends to . This is the model’s representation of the catastrophic failure of the bar. The work-hardening regime ends if, with growing inelastic deformation, ((t)) becomes equal to zero. Of course, this may or may not be possible, according to the shape assumed for the function   Figure 1: Force-elongation response curves at the onset of the inelastic regime, for different signs of (0)  The main result of the local model is that the initial part of the inelastic response is determined by the sign of the second derivative  (0). That is, by the initial convexity or concavity of the cohesive energy. The predictions of the local model http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 9 at the onset of the inelastic regime are summarized in Fig.1. The subsequent part of the response curve depends on the analytic expression of the cohesive energy away from the origin. Anyway, the picture shows a major defect of the model, that is, the incapability of reproducing the regime of strain softening, in which the response curve exhibits a negative slope. Indeed, catastrophic failure occurs as soon as  (0) takes a negative value.  THE NON-LOCAL MODEL he reproduction of strain softening becomes possible in the non-local model obtained by adding to the energy (3) a term proportional to the square of the gradient of the inelastic deformation   l dxxxwE 0 )))(())(((),(  +  l dxx 0 2 2 1 )(' (17) where  is a small positive constant. The addition of the new term brings additional terms to the yield condition (8)    ((x)) (x) (18) and to the Kuhn-Tucker conditions (11) 0)( x , 0)(''))(('  xx  , 0)())(''))(('(  xxx   (19) and (13) 0)( x , 0)('')())((''  xxx   , 0)())('')())((''(  xxxx   , x (20) It also adds a positive term to the second variation. Then the necessary condition (14) is relaxed, and this is the reason why a description of strain-softening becomes possible. At the onset of the inelastic regime, due to the additional boundary conditions required by the supplementary gradient term, the inelastic deformation is not anymore constant all over the bar. The choice of the additional boundary conditions is a delicate point. Our choice of taking  (l ) = (0) = 0 (21) is discussed in detail in the paper [1]. With these conditions, for the inelastic strain rate we get the expression xk xl x sinh 2 tanh 2 tanh )0('' )(              , (22) with  =((0) /)1/2, if (0) > 0, and )( 2 )( xlxx       (23) if (0) = 0. In both cases, the added energy term has the effect of increasing the slope of the force-elongation response curve. In particular, for (0) = 0 a hardening response takes the place of the perfectly plastic response predicted by the local model. For (0) < 0, there are two types of solutions: the full-size solution xk lkxk x sin 2 tan 2 tan )0('' )(           , (24) with k =((0) /)1/2, if kl < 2, and the localized solution  ))(cos1( )0('' )( axkx     , (25) T http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 10 if kl > 2.The latter are concentrated on an interval (a, a + ly) of length ly = 2 /k < l. The initial slope of the inelastic response curve is positive for kl < , and negative for kl > . It is usual to consider the material constant li =  / k (26) as an internal length of the material. From the preceding analysis it follows that the inelastic regime starts with  a full-size solution and a work-hardening response if l < li ,  a full-size solution and a strain-softening response if li < l < 2 li ,  a localized solution and a strain-softening response if l > 2 li . Catastrophic failure again occurs when the slope of the response curve becomes . This never happens for l < li . For li < l < 2 li , this happens when 0)('' 2/ 2/tan 1)0(''        cw kl kl  (27) where c is the elongation at the onset of the inelastic regime. For l > 2 li , this happens when 0)('')0(''  c y w l l  (28) The occurrence of catastrophic rupture at the onset is the response predicted by Griffth’s theory of brittle fracture. For all remaining situations, the subsequent response depends on the behavior of  away from the origin. Depending on the form assumed for as a function of , it may happen that an initially full-size solution localizes for some  > c , and in some cases it may also happen that a localized solution becomes full-size. Catastrophic fracture may take place in both localized and full-size solutions. NUMERICAL SIMULATIONS ith an appropriate choice of the expression of , the non-local model gives the possibility of describing many experimental situations in which rupture is preceded by more or less extended regimes of inelastic deformation, as well as many intermediate situations between the extreme cases of a totally brittle and a totally ductile response. However, the non-homogeneous character of the solutions in the non-local model makes hopeless any attempt of describing the evolution of the inelastic deformation by a simple differential equation, like Eq. (15) in the case of the local model. For this reason, we made a series of numerical simulations. Here we present the current status of our study, which is still in progress. The purpose of our simulations was to determine the shape of the function  giving the best reproducion of the response curves of two experimental tests, one on a steel bar and one on a concrete specimen. For  we chose a piecewise polynomial C2 representation. This was obtained by subdividing the domain {  0 } into N intervals, in each of which  was taken to be a third-order polynomial, and by imposing the continuity of the function and of its first and second derivatives at the nodal points of the subdivision. The number and position of the nodal points determine the accuracy of the approximation. The advantage of a better approximation obtained with large N is paid with the disadvantage of working with a larger number of material constants. For the steel bar we made two series of simulations, one with N = 3 and one with N = 6. The results of the first series are shown in Fig. 2. In it, we see the response curves obtained numerically for a bar length l equal to 40, 80, 160 mm. They show that a short bar is more ductile than a long bar, since for short bars catastrophic failure takes place for larger . This is a manifestation of the size effect. In the same figure, the curve for l = 80 mm is compared with the experimental curve, represented by the dotted line, for a bar of the same length. We see a very good agreement, except for the initial horizontal plateau exhibited by the experimental curve and not reproduced in the simulation. The second series of simulations was done with the purpose of eliminating this discrepancy, by increasing the number N of parameters. The result is shown in Fig. 3, where we see an almost perfect agreement between experiment and simulation. W http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 11  Figure 2: Response curves for a steel bar in numerical simulations with N = 3 and different values of l, and experimental response curve (the dotted line) for l = 80 mm.   Figure 3: Force-elongation response curve for a steel bar in a numerical simulation with N = 6 and l = 80 mm, compared with the experimental response curve (the dotted line) (a). Detail of the plateau at the onset of the inelastic regime (b) For the concrete specimen, at present the simulations are still in progress. A first result is shown in Fig. 4. While the general trend of the response is well reproduced, just after the onset of the inelastic deformation a sharp change of the slope occurs in the simulation but not in the experiment. A larger N is probably necessary to eliminate this discrepancy. The distribution of the inelastic deformation (x) along the bar’s axis is shown in Fig. 5 for different values of . The resemblance with the measurements by Miklowitz made many decades ago [11], is impressive. Both diagrams show a progressive concentration of inelastic deformation at the central zone of the bar. If we imagine that, by a sort of Poisson effect, axial elongation is accompanied by a proportional transversal contraction, we have that a concentration of axial deformation is accompanied by the necking of the cross section. Therefore, the proposed model provides a good description of the phenomenon of necking observed in bars subjected to a tensile load. http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 12  Figure 4: Force-elongation response curve for a concrete specimen in a numerical simulation with N = 4 compared with the experimental curve (the dotted line)  Figure 5: The inelastic deformation  as a function of x, for different values of the total elongation . The first diagram shows the results of the numerical simulations, and the second shows the experimental measurements reported in [11]. REFERENCES [1] G. Del Piero, G. Lancioni, R. March, A diffuse cohesive energy approach to fracture and plasticity: the one- dimensional case. In preparation (2011). http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true G. Del Piero et alii, Frattura ed Integrità Strutturale, 18 (2011) 5-13; DOI: 10.3221/IGF-ESIS.18.01 13 [2] A.A. Griffith, Phil. Trans. Roy. Soc. A221, 163-198 (1920). [3] D. Dugdale, J. Mech. Phys. Solids, 8 (1960)100. [4] G.I. Barenblatt, Adv. Appl. Mech. , 7 (1962) 55. [5] G.A. Francfort, J.-J. Marigo, J. Mech Phys. Solids, 46 (1998) 1319. [6] J.-J. Marigo, L. Truskinovsky, Cont. Mech. Thermodyn, 16 (2004) 391. [7] M. Charlotte, J. Laverne, J.-J. Marigo, Eur. J. Mech., A/25 (2006) 649. [8] G. Del Piero, L. Truskinovsky, Cont. Mech. Thermodynamics, 21 (2009) 141. [9] L. Freddi, G. Royer-Carfagni, J. Mech. Phys. Solids, 58 (2010) 1154. [10] K. Pham, H. Amor, J.-J. Marigo, C. Maurini, Int. J. Damage Mechanics, 20 (2011) 618. [11] J. Miklowitz, In: Proc. Nat. Meeting of the Applied Mechanics Division, ASME, Ann Harbor, Mich, (1949) 159. http://www.gruppofrattura.it http://dx.medra.org/10.3221/IGF-ESIS.18.01&auth=true