A DIFFUSE COHESIVE ENERGY APPROACH TO FRACTURE AND PLASTICITY: THE ONE-DIMENSIONAL CASE GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH Abstract. In the fracture model presented in this paper, the basic assump- tion is that the energy is the sum of two terms, elastic and cohesive, depending on the elastic and inelastic part of the deformation, respectively. Two variants are examined, a local model, and a non-local model obtained by adding a gradient term to the cohesive energy. While the local model only applies to materials which obey Drucker’s postulate and only predicts catastrophic failure, the non-local model describes the softening regime, and predicts two collapse mechanisms, one for brittle and one for ductile fracture. In its non-local version, the model has two main advantages over the mo- dels existing in the literature. The first is that the basic elements of the theory (yield function, hardening rule, evolution laws) are not assumed, but are de- termined as necessary conditions for the existence of solutions in incremental energy minimization. This reduces to a minimum the number of the indepen- dent assumptions required to construct the model. The second advantage is that, with appropriate choices of the analytical shape of the cohesive energy, it becomes possible to reproduce, with surprising accuracy, a big variety of ob- served experimental responses. In all cases, the model provides a description of the entire evolution, from the initial elastic regime to final rupture. 1. Introduction Since the pioneering paper of Francfort and Marigo [24] in 1998, the variational ap- proach to fracture has been the object of intense research. A reason for its fortune was the availability of numerical instruments supporting the theory, and capable of solving practical problems in fracture mechanics. One of such instruments is the regularization technique of Ambrosio and Tortorelli [3], initially conceived for the image segmentation problem [41], and applied to fracture problems by Bour- din, Francfort and Marigo [11]. The regularization consists in approximating the energy of a fractured body, which by its own nature is defined on discontinuous de- formations with bounded variation, by a family of functionals defined in a Sobolev space. This makes possible to solve fracture problems using standard finite element techniques. The original paper [24] and the subsequent developments were based on the brittle fracture model of Griffith [28], in which the presence of singular discontinu- ity surfaces for the displacement is allowed at the interior of an elastic body. In fact, before the appearing of [24], more sophisticated models, based on Barenblatt’s cohesive fracture hypothesis [5], were used in fracture mechanics. They were par- ticularly efficient in situations in which fracture is preceded by a regime of large Date: January 28, 2013. Key words and phrases. Non-local plasticity, Strain localization, Variational fracture, Ductile fracture, Incremental energy minimization. 1 2 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH inelastic deformation, like in the rupture of concrete and other non-metallic mate- rials [14, 33]. In this respect, a fundamental role was played by the cohesive crack model of Hillerborg et al. [34]. The development of the theories of rupture for concrete over the following decades is documented in the review papers [6, 7, 8, 13]. With the cohesive energy model, two basic aspects of energy minimization came to evidence. The first is the relevance of local, versus global, minimization. This point received large attention in the recent literature [16, 39], and is now generally accepted. The second aspect is the relevance of the role played by energy dissi- pation. Indeed, with a minimization constrained by a dissipation inequality it is possible to reproduce some sophisticated aspects of material response such as, for instance, elastic unloading in plasticity [19]. The present paper originates from a reflection on similarities and differences between the cohesive energy scheme and the Ambrosio–Tortorelli approximation. For both, the energy functional depends on the elastic deformation plus a second variable which provides a regularized representation of fracture. But, while the cohesive energy is usually assumed to be a surface energy, the regularizing term in [3] is a volume integral. The question naturally arises, of whether or not it is convenient to assume a cohesive energy concentrated on singular surfaces. In ductile materials, before rupture, the appearance of a process zone is observed, in which the material becomes weaker and more deformable. This phenomenon can be described by damage [4] or plasticity [17] theories. Cohesive energies diffused over the volume were recently considered in the modeling of brittle fracture [48] and damage [10, 44], and in the variational theory of fracture [26]. In particular, the paper [42] by Pham and Marigo deals with a model of energy minimization, in which rupture is preceded by progressive damage. In it, the one-dimensional energy functional has the form (1.1) ∫ l 0 ( w(u′(x), γ(x)) + θ(γ(x)) ) dx , where w is the bulk energy density, u′ is the axial deformation, and γ is an internal variable, representing the intensity of damage. In the present paper we consider the special case (1.2) ∫ l 0 ( w(ε(x)) + θ(γ(x)) ) dx , in which ε(x) = u′(x) − γ(x) is the elastic part of the deformation. Of the two energy densities w and θ, the first is assumed to be reversible, and the second is assumed to be totally dissipated. By assuming w independent of γ, we purposedly renounce any description of damage. Though our initial purpose was to construct something more general than a plasticity theory, to our surprise we realized that, as we shall see below, the energy (1.2) provides exactly the classical incremental theory of plasticity exposed, for example, in [32]. A weak point of the model is that it fails to describe the strain-softening response which, as well known from the literature, is associated with localization of the plastic deformation. These effects can be captured by adding a non-local term of the gradient type (1.3) ∫ l 0 ( w(ε(x)) + θ(γ(x)) ) dx+ 1 2 α ∫ l 0 γ′2(x) dx , DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 3 which successfully captures some specific effects of material response in a large variety of models, like the Van der Waals model for capillarity [15, 47], the Cahn– Hilliard model for phase transition [12], the Mumford-Shah image segmentation theory [41], some theories of liquid crystals [25, 22], and some numerical models for the rupture of concrete [18]. A non-local term was also used by Truskinovsky, [46], to construct a cohesive zone model starting from the energy of interatomic attraction (Lennard-Jones potential). The addition of a gradient term renders the expression (1.3) of the energy similar to the Ambrosio–Tortorelli functional (1.4) ∫ l 0 ( s2(x) + κα)w(u ′(x) ) dx+ 1 2 α ∫ l 0 ( s′2(x) + c α2 (1−s(x))2) dx. Indeed, in both cases the total energy is the sum of a basic energy and a pertur- bing gradient term, and in both cases the deformation is accompanied by a second independent variable. The main difference is that in (1.3) the second variable is the inelastic part of the deformation, while in (1.4) the scalar field s ∈ [0, 1] measures the fracture intensity in the regularized model. Also different is the nature of the smallness parameter α. In (1.4), α is the coefficient of a singular perturbation, and what is important is the asymptotic behavior of the solution when α → 0. In fracture theory, α is a fixed material constant, usually interpreted as an internal length determined by the microscopic structure of the material [7, 9, 44]. The value of α strongly influences the material’s response. Indeed, the ratio between internal length and the length l of the bar is an indicator of the brittleness of the fracture. A reduction of α has the same effect as an increase of l, that is, a transition from rupture preceded by a long regime of gradual weakening to sudden catastrophic failure. For fixed α and varying l, this is the well known size effect of fracture mechanics. The requirement of non-negativeness of the first variation leads to Kuhn-Tucker conditions very similar to, and sometimes identical with, the governing equations of the gradient theories of plasticity [2, 23, 31] and damage [10, 43, 44], which, starting from the pioneering paper [1], have been developed by several authors in the last two decades. In particular, the expression of the first variation of the energy given in our Section 3.2 is the one-dimensional counterpart of the expression of the virtual power currently assumed in the theories of strain-gradient plasticity [29, 30, 31]. An accurate review, analysis, and comparison of such theories can be found in the papers [37, 38]. Most theories include a number of assumptions which fix the basic physical aspects of plastic response: - the stress-strain elastic law, - the hardening rule, - the evolution laws, - the loading-unloading conditions, - the yield function. The underlying continuum can be either a classical or a generalized continuum. In a classical continuum the plastic deformation is regarded as a state variable, the Cauchy stress is the unique stress measure, and supplementary generalized forces are defined as the conjugates of the gradients of the plastic deformation. In a generalized continuum, higher-order stresses appear as conjugates of the higher- order deformation gradients. 4 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH The two cases are analyzed in [38] and [37], respectively. The number and variety of theories reported there shows by itself that there is yet no agreement about the choice of a fully satisfactory model. The present paper has the ambition of presenting a mathematically sound model, based on very simple assumptions, and providing a very flexible description of the entire evolution from the initial elastic regime to the final rupture. The basic assumptions are the forms of the elastic energy and of the dissipa- tive cohesive energy. Proceeding by incremental energy minimization, all classical hypotheses listed above are obtained as necessary conditions for a minimum. The analytical forms of the functions defined in the above hypotheses are already con- tained in the assumed form of the cohesive energy. In particular, the dissipative character assumed for the cohesive energy provides the basis for the description of elastic unloading, and the shape of the cohesive energy can be chosen to reproduce a big variety of responses, hardening or softening, and with the plastic strain both diffusing over the body or localizing on small regions. Surface cohesive energies have proved to be an efficient tool to formulate and solve fracture problems. Here we show that volume cohesive energies can be a very natural way for describing fracture as a terminal event in an evolution involving plasticity or damage. Of course, there is no claim of generality for the model pre- sented here. It can be generalized in several directions. We just like to compare the simplicity of the starting assumptions and of the incremental minimizing procedure, with the richness of the responses that can be obtained by acting on the analytical shape of the cohesive energy. Our analysis is restricted to the one-dimensional case. This allows us to get the basic results avoiding a number of technical complications, at the cost of neglect- ing some important aspects congenital to higher dimension. The paper is divided into four sections. After the present Introduction, Sections 2 and 3 deal with the local and the non-local model, respectively, and Section 4 reports some numerical simulations on the non-local model and their comparison with experiments. The local model is introduced in Subsection 2.1, and in Subsection 2.2 the equilibrium configurations are identified with the stationary points of the energy (1.2). Due to the unilateral character of the dissipation constraint, the Euler equation is re- placed by an inequality. This inequality provides an upper bound for the axial force, depending on the current value of the plastic deformation. In this way, the yield condition is obtained. The minimum properties of the stationary points are investigated in Subsection 2.3. The next Subsections deal with the quasi-static evolution of the inelastic deformation under varying load, along a path made of energy minimizers. For simplicity, we only consider a bar under a hard device, for which the load consists of displacements applied at the endpoints. Subsection 2.4 is devoted to incremental energy minimization. In a first-order approximation, the minimization consists in determining the direction of steepest descent under the given load increment. The result is that along this direction the inelastic deformation does not change in the elastic regime, that is, when the yield condition is satisfied as a strict inequality. On the contrary, in the inelastic regime the inelastic strain rate remains unspecified. For its determination, a second-order minimization is required. In it, some basic properties of classical plasticity, such as the hardening rule, the consistency condi- tion for the inelastic deformation at loading, and the property of elastic unloading, DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 5 are obtained as necessary conditions for a minimum. The minimization leads to a first-order differential problem, involving the current values of the inelastic de- formation and of the inelastic strain rate. This problem is studied in Subsection 2.5. In the inelastic regime, at each time t the response is determined by the sign of the second derivative of θ at the current deformation γt: the response is work- hardening for θ′′(γt) > 0, and perfectly plastic for θ′′(γt) = 0. For θ′′(γt) negative, the deformation concentrates at a single point, determining catastrophic rupture. The local model provides a classical incremental plasticity theory [32], obeying Drucker’s material stability postulate [21] (1.5) σ̇tγ̇t ≥ 0 , which does not allow for a strain-softening response. Within this restricted context, the model correctly describes a plastic response followed by brittle fracture in the sense of Griffith’s theory. The non-local model based on the energy (1.3) is discussed in Section 3. Due to the presence of a gradient term, more regularity and supplementary boundary conditions are required. The additional regularity assumptions are made in Sub- section 3.1, and two different options for the supplementary boundary conditions are considered in Subsection 3.2. The choice of the supplementary boundary con- ditions is a delicate problem, see, e.g., the discussion in [38], Sect. 3.4. Our choice of imposing a null inelastic deformation γ at the boundary is motivated by a closer adherence to the constraint devices used in laboratory tests. The minimum properties of the equilibrium configurations, still identified with the stationary points of the energy, are investigated in Subsection 3.3. Due to the stabilizing effect of the non-local term, the necessary conditions for a minimum are weaker than in the local model. In particular, moderate negative values of θ′′(γ(x)) become acceptable, and this renders admissible a strain-softening response. In the language of [37, 38], the non-local term acts as a localization limiter, since it opposes the localization of the inelastic deformation. The quasi-static evolution under a given load process is studied in Subsection 3.4. A first-order minimization still provides the direction of steepest descent for the energy functional. Along this direction, the inelastic strain rate γ̇t(x) is zero in the elastic zone, that is, at points x at which the yield condition is satisfied as a strict inequality. Surprisingly enough, this conclusion is partially contradicted in the second-order minimization, in which non-null inelastic strain rates are allowed inside the elastic zone, at points sufficiently close to the interface with the inelastic zone. The second-order minimization provides a set of Kuhn-Tucker conditions, which governs the quasi-static evolution of the inelastic deformation. Due to the presence of the non-local terms, closed-form solutions can be found only in some special cases. One of them is the onset of the inelastic regime, studied in detail in Subsection 3.5. In the determination of the incremental response at the onset, a fundamental role is played by the non-dimensional ratio l/li, where (1.6) li . = 2π √ α −θ′′(0) is the internal length of the material at the onset. Indeed, the response is work- hardening if l < li/2 and strain-softening if l > li/2. Moreover, for l < li the response is full-size, that is, the inelastic strain rate is diffused over the whole bar, 6 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH while for l > li the inelastic strain rate localizes on a region of length li. It is possible that catastrophic failure occurs just at the onset of the inelastic regime. This circumstance, which we call totally brittle fracture, depends on a second non- dimensional ratio l/lc, where (1.7) lc . = 2π √ α w′′(βc) is the characteristic length at the onset, and βc is the value of the load β at the onset. In fact, totally brittle fracture occurs when l/lc is less than a given function of l/li. A qualitative analysis of the response in the inelastic regime is made in Subsection 3.6. In it, internal length and characteristic length are considered as functions of the current inelastic and elastic deformation (1.8) lit = 2π √ α −θ′′t , lct = 2π √ α w′′(εt) , respectively, where θ′′t is an average of θ ′′(γt(x)) over the bar in the current inelastic deformation γt. In this way, by extrapolation of the exact results obtained at the onset, it is possible to obtain an approximate picture of the bar’s evolution up to final rupture. In particular, from this qualitative analysis two fundamental types of evolution emerge. They depend on the convexity-concavity properties of the derivative θ′ of the cohesive energy, and determine two very different collapse mechanisms. If θ′ is concave, an increasing load produces a progressive localization of the inelastic deformation, and a progressive increase of the negative slope of the force-elongation response curve. In the limit, an extreme localization and an infinite slope of the response curve produce catastrophic failure. On the contrary, if θ′ is convex, under increasing load the inelastic deformation diffuses over larger and larger zones, and the negative slope of the curve decreases. In the limit the slope tends to zero, and the bar elongates indefinitely at constant force. The two collapse mechanisms are typical of materials with crystalline struc- ture and with random structure, respectively. In Section 4 both are investigated by means of numerical simulations, with the purpose of reproducing the response curves of two real materials, concrete and steel. The correspondence of the simu- lations with the experimental curves is obtained by fitting a very small number of parameters, that is, the material constants necessary to represent the cohesive energy θ as a piecewise polynomial function. The result is impressive. It shows that an appropriate choice of the expression of the cohesive energy can indeed capture not only the overall behavior, but also many detailed features of the response curve, in the whole evolution from the initial natural state up to rupture. Roughly, the main results of the present study can be summarized as follows. - In its non-local version, the diffuse cohesive energy model captures both brit- tle and ductile fracture modes. In brittle fracture the inelastic deformation tends to localize on singular surfaces, and in ductile fracture it tends to spread over the whole bar. Brittle fracture may also take place without any preliminary regime of inelastic deformation. DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 7 - A convex cohesive energy produces a full-size work-hardening response. A concave energy may produce full-size or localized, work-hardening or strain- softening response, depending on the ratio between the current internal length and the length of the bar. - For θ concave, a convex θ′ produces brittle fracture, and a concave θ′ pro- duces ductile fracture. - Brittle fracture occurs when the current internal length reaches a critical value depending on the current characteristic length. Some of these results, reflecting an earlier stage of the present research, were an- nounced in the paper [20]. But even the present results are not complete, since many problems are left open. Specifically, it seems very difficult to reach the same level of detail in a two- or three-dimensional context. Also, a big challenge is how to find general, physically motivated correlations between the convexity-concavity properties of θ and θ′ and the microstructural properties of real materials. This would require a truly interdisciplinary effort, far beyond the domain of classical continuum mechanics. 2. The local model 2.1. Basic assumptions. Consider a straight bar of length l, homogeneous, with a constant cross section, subject to prescribed axial displacements at the endpoints (2.1) u(0) = 0 , u(l) = βl , and free of applied loads. In what follows, the term load will be used to denote the relative elongation β. A deformed configuration is described by the axial displace- ments u(x) of the points x of the bar’s axis, and the derivative u′(x) is a measure of the axial deformation at x. We make four basic assumptions. The first is that at every x the deformation can be split into the sum of an elastic part ε(x) and an inelastic part γ(x) (2.2) u′(x) = ε(x) + γ(x) , x ∈ (0, l) . A pair (ε, γ) of functions continuous in [0, l] is a configuration of the bar. By integration, from the boundary conditions (2.1) it follows that (2.3) β = ε̄+ γ̄ , with (2.4) ε̄ = 1 l ∫ l 0 ε(x) dx , γ̄ = 1 l ∫ l 0 γ(x) dx . Our second assumption is that the strain energy of the bar has the form (2.5) E(ε, γ) = ∫ l 0 ( w(ε(x)) + θ(γ(x)) ) dx , where w and θ are the volume densities of the elastic strain energy and of the cohesive energy, respectively. We assume that w is C2 and strictly convex, and that (2.6) w(0) = w′(0) = 0 . 8 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH Then w′ is strictly increasing, and w(ε) > 0 ∀ε = 0 , w′(ε) { > 0 ∀ε > 0 , < 0 ∀ε < 0 , w ′′(ε) > 0 ∀ε ∈ R . We also assume that θ is C2 and monotonic increasing, and that (2.7) θ(0) = 0 , θ′(0) > 0 , lim γ→+∞ θ(γ) < +∞ . The third assumption is that the elastic part of the deformation is related to the axial force σ by the constitutive equation (2.8) σ = w′(ε) , and the fourth assumption is that w can be stored, while θ is totally dissipated. That is, in every deformation process t → (εt, γt), at every x the cohesive power is non-negative: (2.9) θ′(γt(x)) γ̇t(x) ≥ 0 . Here and in the following, a superimposed dot denotes the right derivative with respect to the parameter t. Due to the rate-independent character of the model, t may be any monotonic non-decreasing function of the physical time. In what follows we consider the case of θ strictly monotonic, for which the dissipation inequality (2.10) γ̇t(x) ≥ 0 holds at all x and for all t. 2.2. Equilibrium. For a given configuration (ε, γ), a perturbation is a pair (δε, δγ) such that δε(x) = ε̇t(x) , δγ(x) = γ̇t(x) , in some deformation process with (εt, γt) = (ε, γ). While δε is unrestricted, δγ is subject to the condition (2.11) δγ(x) ≥ 0 ∀x ∈ (0, l) , due to the dissipation inequality (2.10). We say that (ε, γ) is an equilibrium confi- guration if the first variation of the energy (2.12) δE(ε, γ, δε, δγ) = lim λ→0+ 1 λ ( E(ε+λ δε, γ+λ δγ)− E(ε, γ)) = ∫ l 0 ( w′(ε(x)) δε(x) + θ′(γ(x)) δγ(x) ) dx is non-negative for all perturbations which preserve the length of the bar (2.13) δε̄+ δγ̄ = 0 , δε̄ = 1 l ∫ l 0 δε(x) dx , δγ̄ = 1 l ∫ l 0 δγ(x) dx . In particular, for perturbations with δγ=0 the first variation is δE(ε, γ, δε, 0) = ∫ l 0 w′(ε(x)) δε(x) dx . By (2.13), δε̄ = 0 for such perturbations. Then, a standard argument of the Calculus of Variations leads to the conclusion that the first variation is non-negative if and only if w′(ε(x)) is a constant. DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 9 By the constitutive equation (2.8), this implies that the axial force σ is constant. Moreover, by the strict convexity of w, the derivative w′ is strictly increasing and, therefore, invertible. Then, ε is constant as well (2.14) ε = (w′)−1(σ) . By (2.13), for constant ε the first variation reduces to δE(ε, γ, δε, δγ) = ∫ l 0 ( θ′(γ(x)) − w′(ε)) δγ(x) dx , and because δγ(x) is arbitrary non-negative, the condition (2.15) w′(ε) ≤ θ′(γ(x)) ∀x ∈ (0, l) follows. Conversely, for every configuration (ε, γ) obeying this condition the first variation is non-negative. Therefore, an equilibrium configuration is a configuration (ε, γ) with constant ε, which satisfies condition (2.15) for all x. An equilibrium configuration is equilibrated with the load β given by (2.3). Then inequality (2.15), rewritten in the form (2.16) w′(β − γ̄) ≤ θ′(γ(x)) ∀x ∈ (0, l) , shows that an equilibrium configuration can be defined alternatively as a pair (β, γ) which satisfies (2.11) and (2.16). There are many γ for which (β, γ) is an equilibrium configuration. Indeed, by the assumptions made on w and θ, all configurations (β, γ) with γ̄ > β satisfy inequality (2.16), because the right side of the inequality is non-negative, while for β < γ̄ the left side is negative. With every equilibrium configuration (β, γ) are associated the axial force (2.17) σ = w′(β − γ̄) , and the energy (2.18) Eeq(β, γ) . = lw(β − γ̄) + ∫ l 0 θ(γ(x)) dx . The set of all x at which inequality (2.16) is strict is the elastic zone (2.19) E(β, γ) .= {x ∈ (0, l) | w′(β − γ̄) < θ′(γ(x))} , and the complementary set J (β, γ) .= (0, l) \ E(β, γ) at which (2.16) is verified as an equality is the inelastic zone. We say that the bar is in the elastic regime if E(β, γ) = (0, l), and that it is in the inelastic regime otherwise. The equilibrium condition (2.16), rewritten in the form (2.20) σ ≤ θ′(γ(x)) ∀x ∈ (0, l) , shows that there is an upper limit for σ. In the language of Plasticity, this inequality is a yield condition, which imposes to the current stress σ to be not greater than the yield limit θ′. The difference θ′ − σ is the yield function. 10 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH 2.3. Energy minimizers. Among all inelastic deformations γ equilibrated with a given β, of interest are those which minimize the energy Eeq(β, · ). A necessary condition for a minimum at γ is that, in the expansion Eeq(β, γ + λδγ)− Eeq(β, γ) = λ ∫ l 0 ( θ′(γ(x)) − w′(β − γ̄)) δγ(x) dx + 1 2 λ2 ∫ l 0 ( θ′′(γ(x)) δγ2 + w′′(β − γ̄) δγ̄2) dx+ o(λ2) , the first variation be non-negative for all perturbations δγ and for all positive values of the smallness parameter λ. By (2.11) and (2.16), this is true if and only if (β, γ) is an equilibrium configuration. Therefore, all minimizers for Eeq(β, · ) are equilibrium configurations. A second necessary condition is that the second variation be non-negative for all δγ for which the first variation is zero. In particular, it must be non-negative for all δγ with support in the inelastic zone J (β, γ). For this, it is necessary that (2.21) θ′′(γ(x)) ≥ 0 a.e. x ∈ J (β, γ) . Indeed, if θ′′(γ(x)) is negative on a subset of J (β, γ) with non-null measure, a negative second variation is obtained by concentrating δγ on that set. The exis- tence of perturbed configurations (β, γ + δγ) with arbitrarily large negative energy characterizes a fractured configuration. Then, (2.21) is in fact a safety condition against fracture. The non-negativeness of the first variation joined with the stronger condition (2.22) θ′′(γ(x)) > 0 ∀x ∈ (0, l) , is sufficient for a local minimum at γ with respect to the L∞ norm (2.23) ‖ δγ ‖ = sup x∈(0,l) |δγ(x)| . Indeed, since θ is C2 and γ is continuous in the closed set [0, l], the map θ′′(γ( · )) is uniformly continuous. Then for all ε > 0 there is a δ > 0 such that |a| < δ =⇒ | θ′′(γ(x) + a)− θ′′(γ(x)) | < ε ∀x ∈ (0, l) . Consider the finite expansion of Eeq(β, · ) with Lagrange remainder Eeq(β, γ+ δγ) = Eeq(β, γ) + δEeq(β, γ, δγ) + 1 2 ∫ l 0 ( θ′′ ( γ(x) + ξ(x) δγ(x) ) δγ2(x) + w′′ ( β − γ̄ − ξ̄ δγ̄) δγ̄2) dx , with ξ̄ and all ξ(x) in (0, 1). For every perturbation δγ with ‖δγ‖ < δ, from the non-negativeness of the first variation and the positiveness of w′′ we have Eeq(β, γ+ δγ)− Eeq(β, γ) ≥ 1 2 ∫ l 0 ( θ′′(γ(x))− ε) δγ2(x) dx . If (2.22) holds, the map θ′′(γ( · )) is positive and continuous, and therefore it has a positive minimum c in [0, l]. It is then sufficient to take ε 0 there is a τo > 0 such that ‖γt+τ − γt‖ < δ ∀τ ∈ (0, τo) , where ‖ · ‖ is the L∞ norm (2.23). By consequence, in the time interval (0, τo) the evolution τ → γt+τ takes place inside the ball B(γt, δ) of L∞(0, l) of radius δ centered at γt. Assume that γt(·) is a deformation with constant value γ̄t, and that θ′′(γ̄t) is positive. Then it is possible to choose δ such that θ′′(γ(x)) > 0 for all γ in B(γ̄t, δ), so that the restriction of the energy to B(γ̄t, δ) is strictly convex for all τ < τo. Therefore, by Jensen’s inequality, Eeq(βt+τ , γt+τ ) ≥ Eeq(βt+τ , γ̄t+τ ) ∀τ ∈ (0, τo) , where γ̄t+τ is the homogeneous deformation defined as in (2.4), and equality holds only if γt+τ = γ̄t+τ . That is, the unique minimizer of the energy at time t+τ is a homogeneous deformation. Consequently, for every quasi-static evolution from a homogeneous deformation γ̄t with θ ′′(γ̄t) > 0 there is a time interval (0, τo) in which γt+τ is a homogeneous deformation. From here onwards we consider homogeneous deformation processes t → γt, and we identify the functions γt(·), γ̇t(·) with their constant values, which we denote by γt, γ̇t. For θ ′′(γt) > 0 and for δ and τo as above, the determination of a quasi-static evolution from (βt, γt) reduces to the following incremental minimization problem: for a given τ ≤ τo, find the homogeneous deformation γt+τ which minimizes the energy Eeq(βt+τ , · ) in B(γt, δ). An approximate solution is obtained by replacing βt+τ and γt+τ by the piecewise linear approximations (2.24) βt+τ ≈ βt + τβ̇t , γt+τ ≈ γt + τ γ̇t , and then determining γ̇t as the minimizer for the first-order approximation (2.25) Et+τ (γ̇) ≈ Et + τĖt(γ̇) of the energy Eeq(βt+τ , γt+τ ), where (2.26) Et = E eq(βt, γt) = l ( w(εt) + θ(γt) ) , Ėt(γ̇) = l ( σt (β̇t − γ̇) + θ′(γt) γ̇ ) , with σt = w ′(εt) and εt = βt− γt. This minimization determines the direction of steepest descent for the energy. Because Et and β̇t are known, the problem reduces to the minimization of the linear function (2.27) I(γ̇) = (θ′(γt)− σt) γ̇ , under the condition γ̇ ≥ 0. In the elastic regime, in which (θ′(γt)− σt) is positive, the minimum is zero and is attained at γ̇ = 0. In the inelastic regime, I(γ̇) is 12 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH identically zero because (θ′(γt) − σt) is zero. That is, the direction of steepest descent is indetermined. The determination of a quasi-static evolution in the inelastic regime then requires a second-order approximation of the energy (2.28) Et+τ (γ̇) ≈ Et + τĖt(γ̇) + 12 τ2Ët(γ̇) , where (2.29) Ët(γ̇) = l ( w′′(εt) (β̇t − γ̇)2+ θ′′(γt) γ̇2 ) . After elimination of the terms independent of γ̇, the problem reduces to the mini- mization of the quadratic function (2.30) J(γ̇) = ( θ′(γt)−σt−τ w′′(εt) β̇t ) γ̇ + 12 τ ( θ′′(γt) + w′′(εt) ) γ̇2 , under the condition γ̇ ≥ 0. The sum (θ′′(γt) + w′′(εt)) is positive by the strict convexity of w and the necessary condition (2.21) on θ′′(γ). Then, J is strictly convex. If the bar is in the elastic regime, θ′(γt)−σt > 0, the coefficient of γ̇ is positive for sufficiently small τ . Then the global minimum is zero, and is achieved at γ̇ = 0. In this case, the elastic incremental response (2.31) σ̇t = w ′′(εt) β̇t follows from the incremental force-elongation relation (2.32) σ̇t = w ′′(εt) (β̇t − ¯̇γt) , obtained by time differentiation of (2.17). This confirms the result of the first- order minimization, that there is no increase of inelastic deformation in the elastic regime. In the inelastic regime, θ′(γt) = σt, the minimum is still achieved at γ̇ = 0 at unloading, β̇t ≤ 0. Then, there is no increase of inelastic deformation at unloading. At loading, β̇t > 0, the minimum of J is achieved at (2.33) γ̇t = w′′(εt) θ′′(γt) + w′′(εt) β̇t , and, by (2.32), the inelastic response is (2.34) σ̇t = θ′′(γt)w′′(εt) θ′′(γt) + w′′(εt) β̇t . The different response at loading and unloading in the inelastic regime is a peculiar property of plastic behavior. At loading, from (2.33) and (2.34) it follows that (2.35) σ̇t = θ ′′(γt) γ̇t , whence (2.36) σt+τ ≈ σt + τσ̇t = θ′(γt) + τθ′′(γt)γ̇t ≈ θ′(γt+τ ) . Equation (2.35) is the hardening rule, which in classical plasticity specifies σ̇ as a function of the current value of γ. Equation (2.36) states that, within the approxi- mation (2.24), at time t + τ the yield condition holds as an equality. That is, the inelastic regime is preserved. This is the consistency condition of classical plastic- ity. At unloading, the response (2.31) expresses the property of elastic unloading. Classically, hardening rule, consistency, and elastic unloading are assumptions orig- inating from experimental observation. Here, they come as necessary conditions for a minimum in the second-order energy minimization. DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 13 Equation (2.34) provides the slope dσt/ dβt = σ̇t/β̇t of the force-elongation re- sponse curve at loading. For θ′′(γt) > 0 this slope is positive but smaller than the slope w′′(εt) at unloading. A positive slope at loading characterizes a work- hardening response. For θ′′(γt) = 0, the restriction to homogeneous deformations is not anymore possible. In the inelastic regime, σt = θ ′(γt), the function J(γ̇) reduces to J(γ̇) = −τw′′(εt) β̇t ¯̇γ + 12 τw′′(εt) ¯̇γ2 . The minimum is achieved at ¯̇γt = β̇t, which corresponds to σ̇t = 0, that is, to perfectly plastic response. In this case, the incremental minimization specifies the average ¯̇γt, but the punctual values γ̇t(x) are not determined. Thus, in this case the evolution from a homogeneous configuration need not be homogeneous. The case of θ′′(γt) negative, corresponding to a negative slope, that is, to a strain-softening response, is excluded by condition (2.21). Multiplication of (2.35) by γ̇t shows that Drucker’s postulate (1.5) is obeyed only if θ ′′(γt) is non-negative. This shows that the local model provides a valid response only for materials which obey Drucker’s postulate. The indeterminacy of the minimizers in a perfectly plastic response and the impossibility of reproducing the softening response are two major drawbacks of the local model. 2.5. Evolution from the natural configuration. Let us determine the quasi- static evolution for a load process t → βt, with β̇t > 0 for all t, from the natural configuration (β0, γ0) = (0, 0). At t=0 we have σ0 = 0 , θ ′(γ0) = θ′(0) > 0 , so that inequality (2.20) is strict, and γ̇0 is zero. By continuity, inequality (2.20) remains strict, and therefore γ̇t remains equal to zero, over a time interval (0, tc). Because γ0 = 0, γt is zero in this interval. Then an elastic regime takes place, with εt = βt , σt = w ′(βt) . This regime ends when, at some tc > 0, βt reaches the value (2.37) βc = (w ′)−1(θ′(0)) , at which, since σtc= w ′(βc)= θ′(0), inequality (2.20) becomes an equality. At t= tc, all points of the bar switch from the elastic to the inelastic zone. The response at t = tc is work-hardening if θ ′′(0) > 0, and perfectly plastic if θ′′(0) = 0. In the work-hardening regime, by (2.33), the evolution t → γt is determined by the solution of the differential equation (2.38) γ̇t = w′′(βt−γt) θ′′(γt) + w′′(βt−γt) β̇t , t ≥ tc , under the initial condition γtc = 0. A work-hardening regime persists as long as β̇t and θ ′′(γt) remain positive. For β̇t ≤ 0 elastic unloading takes place, and for θ′′(γt) = 0 the response becomes perfectly plastic. Fracture occurs as soon as γt becomes larger than (2.39) γr . = inf { γ ≥ 0 | θ′′(γ) < 0} . 14 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH Indeed, for γ > γr there are perturbed configurations with arbitrarily large negative energy. Because the energy suddenly drops from a finite value to −∞, the frac- ture has a catastrophic character. This fracture mode is called brittle fracture. If θ′′(0) < 0, fracture occurs without any intermediate inelastic regime. This is the case of totally brittle fracture. In conclusion, in the inelastic regime at loading the response at (βt, γt) is work- hardening if θ′′(γt) is positive, and perfectly plastic if θ′′(γt) = 0. For θ′′(γt) > 0, a quasi-static evolution from a homogeneous deformation is made of homogeneous deformations. The fracture is always catastrophic, and occurs as soon as γt becomes larger than γr. 3. The non-local model 3.1. Basic assumptions. The local model provides only a rough description of the bar’s behavior. In particular, it does not describe ductile fracture and strain- softening response. These phenomena are captured by adding to the energy (2.5) a non-local term proportional to the square of the derivative of γ E(ε, γ) = ∫ l 0 ( w(ε(x)) + θ(γ(x)) + 12 αγ ′2(x) ) dx , with α a positive constant. The new term requires a stronger regularity of γ. We assume that γ is C1 on [0, l] and piecewise C2, that is, that γ′′ is continuous except at a finite number of jump points xi, at which the left limit γ ′′(xi−) and the right limit γ′′(xi+) exist. The difference (3.1) [[γ′′]](xi) . = γ′′(xi+)− γ′′(xi−) is the jump of γ′′ at x. Like in the local model, the part of the energy depending on the inelastic defor- mation is supposed to be dissipative. That is, the cohesive power (3.2) θ′(γt(x)) γ̇t(x) + αγ′t(x) γ̇ ′ t(x) is assumed to be non-negative in every deformation process t → γt. For this, it is necessary that the product θ′(γt(x)) γ̇t(x) be non-negative at all x. Indeed, if it is negative at some xo, since γ̇ ′ t(xo) may have opposite signs in different deformation processes from γt(xo), for at least one choice the power (3.2) is negative. By the assumed strict monotonicity of θ′, the dissipation inequality γ̇t(x) ≥ 0 is necessary for dissipativity. In the following Subsection we show that, if the yield limit is non-negative, this condition is sufficient to guarantee the global dissipativity, that is, the non-negativeness of the integral of (3.2) over (0, l), in every equilibrium process. 3.2. Equilibrium. In the non-local model, an equilibrium configuration is a con- figuration (ε, γ) at which the first variation δE(ε, γ, δε, δγ) = ∫ l 0 ( w′(ε(x)) δε(x) + θ′(γ(x)) δγ(x) + αγ′(x) δγ′(x) ) dx is non-negative for all C1([0, l]) and piecewise C2 perturbations (δε, δγ), which satis- fy the dissipation inequality (2.11) and the boundary condition (2.13). Proceeding as in Subsection 2.2, using perturbations with δγ = 0 we find that σ and ε must be DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 15 constant over the bar. Then setting w′(ε) = σ, δε = −δγ̄, and integrating by parts, we get δE(ε, γ,−δγ̄, δγ) = ∫ l 0 ( θ′(γ(x))−σ−αγ′′(x)) δγ(x) dx + α[γ′(x) δγ(x)]l 0 . The non-negativeness of the first variation for all non-negative δγ requires (3.3) σ ≤ θ′(γ(x))− αγ′′(x) , in (0, l), and (3.4) γ′(l) δγ(l) ≥ 0 , γ′(0) δγ(0) ≤ 0 , at the boundary. Inequality (3.3) is the non-local version of the yield condition (2.20). The difference with the local model is that, while the yield limit θ′(γ(x)) of the local model depends only on the value of γ at x, the additional term αγ′′(x) introduces a dependence on the values taken by γ at the neighboring points. There are two ways for satisfying the boundary conditions (3.4). If δγ is allowed to take arbitrary positive values at the boundary, these conditions reduce to γ′(l)≥0 and γ′(0) ≤ 0, respectively. The alternative is to require that (3.5) γ(l) = γ(0) = 0 , and, therefore, to take perturbations δγ such that (3.6) δγ(l) = δγ(0) = 0 . We choose this second possibility, which better describes the standard experimental conditions. Indeed, while this choice keeps the inelastic deformations away from the boundary, the other choice favors the concentration of the inelastic deformation at the boundary. For a similar effect in different models, see [49], Fig. 9 and 10, [10], Fig. 3a, and the discussion in [38]. In laboratory tests this effect is carefully avoided, either by reinforcing the specimen’s end sections, or by weakening the central part of the bar, for example, with the creation of a notch. Thus, our choice allows for a comparison, at least qualitative, with the standard experimental curves. The alternative possibility of imposing a null derivative γ′ was chosen in [10, 43, 44]. This choice has the advantage of allowing for homogeneous inelastic equilibrium configurations, and this greatly facilitates the study of the inelastic regime. Unfortunately, due to the devices used to fix the bar’s ends in real experiments, this circumstance is hardly reproducible in practice. In conclusion, for the non-local model an equilibrium configuration is a pair (ε, γ), with ε a constant and γ a C1([0, l]) and piecewise C2 function, which satis- fies inequality (3.3) at the interior points and conditions (3.5) at the boundary. Assume that the dissipation inequality (2.10) holds, and that the yield limit θ′(γ) − αγ′′ is non-negative. After integrating over (0, l), an integration by parts yields 0 ≤ ∫ l 0 ( θ′(γ(x)) − αγ′′(x)) γ̇(x) dx = ∫ l 0 ( θ′(γ(x)) γ̇(x) + αγ′(x) γ̇′(x) ) dx− α[γ′(x)γ̇(x)]l 0 . On the right side the integral is the total cohesive power. Moreover, by (3.4), the remaining term is non-positive in any equilibrium process. This proves that, if the 16 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH yield limit is non-negative, the dissipation inequality (2.10) guarantees the global dissipativity in every equilibrium process. An equilibrium configuration (ε, γ) is equilibrated with the load β given by (2.3). Alternatively, an equilibrium configuration can be defined as a pair (β, γ), and the energy can be re-defined by (3.7) Eeq(β, γ) = lw(β − γ̄) + ∫ l 0 θ(γ(x)) dx + 1 2 α ∫ l 0 γ′2(x) dx . For an equilibrium configuration (β, γ) defined in this way, the elastic and inelastic zones are (3.8) E(β, γ) = {x ∈ (0, l) | w′(β − γ̄) < θ′(γ(x)) − αγ′′(x)} , J (β, γ) = (0, l) \ E(β, γ) , respectively. As in the local model, we say that the bar is in the elastic regime if E(β, γ) = (0, l), and that it is in the inelastic regime otherwise. 3.3. Energy minimizers. To determine the inelastic deformations γ which mini- mize the energy Eeq(β, · ), consider the expansion Eeq(β, γ + δγ)− Eeq(β, γ) = ∫ l 0 ( θ′(γ(x)) δγ(x) − w′(ε) δγ̄ + αγ′(x) δγ′(x)) dx + 1 2 ∫ l 0 ( θ′′(γ(x)) δγ2(x) + w′′(ε) δγ̄2 + α δγ′2(x) ) dx+ o(‖δγ‖2) , with ε = β−γ̄, and with perturbations δγ non-negative and satisfying the boundary conditions (3.6). A necessary condition for a minimum is the non-negativeness of the first variation. Because this requirement characterizes an equilibrium configuration, all minimizers are equilibrium configurations. Another necessary condition is that the second variation be non-negative for all δγ for which the first variation is zero. In particular, for all δγ with support in the inelastic zone J (β, γ). If (a, a + lj) is an interval in the inelastic zone, the second variation is non-negative only if the smallest eigenvalue of the eigenvalue problem (3.9) ∫ a+lj a ( θ′′(γ(x)) δγ2(x) + l lj w′′(ε) δγ̄2 + α δγ′2(x) ) dx = αρ ∫ a+lj a δγ2(x) dx , is non-negative for all perturbations δγ with support in this interval. This problem can be identified with problem (A.15) in the Appendix, with H as in (A.9) due to the presence of the constraint δγ(x) ≥ 0, and with a=0 , lj = L , θ ′′(γ(x))=αh(x) , δγ(x)=y(x) , δγ̄ = lj l ȳ . The last equality follows from the definitions ȳ = 1 lj ∫ a+lj a y(x) dx , δγ̄ = 1 l ∫ l 0 δγ(x) dx = 1 l ∫ lj 0 δγ(x) dx , which hold for δγ(x) = y(x) with support in (0, lj). Then, from the identification l lj w′′(ε) δγ̄2 = αω ȳ2 , DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 17 comes the relation ω = lj l w′′(ε) α . The smallest eigenvalue αρo has the upper bound (A.17). It provides the necessary condition for non-negativeness (3.10) αλ2oj + θ ′′ oj ≥ 0 , where λoj is a solution of (A.14), now rewritten in the form (3.11) λ2oj = lj l w′′(ε) α ψo(λoj lj) with ψo as in (A.13). Moreover, θ ′′ oj is defined as in (A.16), (3.12) θ′′oj = ∫ a+lj a θ′′(γ(x)) y2oj(x) dx∫ a+lj a y2oj(x) dx , with yoj given by (A.8) if ωL 2 = l3j w ′′(ε)/αl ≤ 4π2, and given by (A.12) if l3j w ′′(ε)/αl > 4π2. Because λoj is a decreasing function of lj , inequality (3.10) must be satisfied for the largest interval contained in the inelastic zone. To find a sufficient condition for a local minimum of Eeq(β, γ), write the finite expansion with Lagrange remainder Eeq ( β, γ+δγ ) = Eeq ( β, γ ) + ∫ l 0 ( θ′(γ(x)) δγ(x) − w′(ε) δγ̄ + αγ′(x) δγ′(x)) dx + 1 2 ∫ l 0 ( θ′′ ( γ(x) + ξ(x) δγ(x) ) δγ2(x) + w′′ ( ε − ξ̄ δγ̄) δγ̄2 + α δγ′2(x)) dx , with 0 ≤ ξ(x) ≤ 1 for all x and 0 ≤ ξ̄ ≤ 1. The first integral, which is the first variation, is non-negative. We wish to find conditions for which the second integral is non-negative as well. By the continuity of θ′′ and w′′, for every ε > 0 there is a positive δ such that sup x∈(0,l) ∣∣ θ′′(γ(x)+ξ(x) δγ(x)) − θ′′(γ(x)) ∣∣ < ε , ∣∣w′′(ε− ξ̄ δγ̄)− w′′(ε) ∣∣ < ε , for all δγ with ‖ δγ ‖ < δ. For all such δγ, the second integral is greater or equal than∫ l 0 ( θ′′(γ(x)) δγ2(x) + w′′(ε) δγ̄2 + α δγ′2(x) ) dx− ε ∫ l 0 ( δγ2(x) + δγ̄2 ) dx . The first integral is the same as in (A.15), with l = L , δγ = y , w′′(ε) = αω , θ′′(γ(x)) = αh(x) . Therefore, Eeq ( β, γ+δγ ) ≥ Eeq(β, γ)+ (αρo − 2 ε) ∫ l 0 δγ2(x) dx . By the arbitrariness of ε, γ is a minimizer if the smallest eigenvalue ρo is strictly positive. By (A.19), a sufficient condition for a minimum is (3.13) αλ2o + θ ′′ min > 0 , 18 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH with λo the solution of (A.14) for the given ω, and with θ ′′ min the smallest value of θ′′(γ(x)) in (0, l). Conditions (3.10), (3.13) improve the corresponding conditions (2.21), (2.22) of the local model. In particular, they show that moderate negative values of θ′′(γ(x)) are allowed in an energy minimizer. As remarked in the Appendix, the two conditions can be far away from each other. In the special case of lj = l and constant θ′′(γ(x)) they coincide, save for the fact that the second inequality is strict. 3.4. Quasi-static evolutions. For a given load process t → βt, a quasi-static evolution from a given deformation γ0 at t = t0 is a family t → γt of inelastic deformations such that (i) γt0 is equal to γ0, (ii) for all t ≥ t0, the function γt(·) is C1 and piecewise C2 in [0, l], (iii) for all t ≥ t0, γt is a local minimizer for Eeq(βt, · ). The regularity of the time derivative γ̇t is determined by Hadamard’s kinematic compatibility condition [[γ̇′t]](x) = −[[γ′′t ]](x) dx dt , where dx/dt is the velocity of a moving jump of γ′t. This condition tells us that to a moving jump of γ′′t corresponds a moving jump of γ̇′t. Therefore, if γt is C1 and piecewise C2, γ̇t is only continuous, piecewise C 1, and piecewise C2. Like in the local model, item (iii) implies that all configurations (βt, γt) are equilibrium configurations. Note that, because of the boundary conditions (3.5), γt is homogeneous only if γt = 0. Therefore, while in the local model it is sometimes possible to restrict the analysis to quasi-static evolutions made of homogeneous deformations, here there are no homogeneous deformations besides γt = 0. Repeating the procedure followed for the local model, fix a time step τ and consider the piecewise linear approximations (2.24) of t → βt and t → γt. For them, we first determine the direction γ̇t of steepest descent, by minimizing the first-order approximation (2.25) of the energy. The problem reduces to the minimization of the linear functional I(γ̇) = ∫ l 0 ( θ′(γt(x)) γ̇(x)− σt ¯̇γ + αγ′t(x) γ̇′(x) ) dx = ∫ l 0 ft(x) γ̇(x) dx , where σt=w ′(βt−γ̄t), and (3.14) ft(x) . = θ′(γt(x))−σt − αγ′′t (x) is the yield function. The domain of I( · ) is the set of all γ̇ which satisfy the dissipation inequality (2.10) and the boundary conditions (3.15) γ̇(0) = γ̇(l) = 0 , which follow from (3.6). In the first variation δI(γ̇t, δγ̇) = ∫ l 0 ft(x)) δγ̇(x) dx , the perturbations δγ̇ satisfy the boundary conditions (3.16) δγ̇(0) = δγ̇(l) = 0 , DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 19 and the inequality (3.17) γ̇t(x) + δγ̇(x) ≥ 0 ∀x ∈ (0, l) . This is because, in a quasi-static evolution, (γ̇t + δγ̇) is the derivative of γt in some deformation process starting from (βt, γt), and therefore it obeys the boundary conditions (3.15) and the dissipation inequality (2.10). The non-negativeness of the first variation requires that ft(x) δγ̇(x) ≥ 0 a.e. x ∈ (0, l) . By (3.17), at points at which γ̇t(x) = 0 the perturbation δγ̇(x) is arbitrary non- negative, and therefore ft(x) must be non-negative. At points at which γ̇t(x) > 0 the perturbation may have any sign, and therefore ft(x) must be zero. Together with the dissipation inequality (2.10) and the yield condition (3.3), this leads to the Kuhn-Tucker conditions (3.18) γ̇t(x) ≥ 0 , ft(x) ≥ 0 , ft(x) γ̇t(x) = 0 , a.e. x ∈ (0, l) , as necessary conditions for a minimum of I at γ̇t. The last equality, called the complementarity condition, requires that γ̇t(x) = 0 if ft(x) > 0, that is, at all points in the elastic zone E(βt, γt). For all γ̇t which satisfy conditions (3.18), I(γ̇t) is zero and Ėt(γ̇t) is equal to l σt β̇t. Thus, like in the local model, the direction γ̇t of steepest descent is not determined in the inelastic zone. For this, it is necessary to minimize the second- order approximation (2.28) of the energy. This reduces to the minimization of the functional (3.19) J(γ̇) = ∫ l 0 ft(x) γ̇(x) dx + 1 2 τ (∫ l 0 ( θ′′(γt(x)) γ̇2(x)+w′′(εt) ¯̇γ2+αγ̇′2(x) ) dx− 2l w′′(εt) β̇t ¯̇γ ) , with respect to all perturbations δγ̇ which satisfy conditions (3.16) and (3.17). For the minimizers of J we keep the same symbol γ̇t used for the minimizers of I. While a minimizer for I is a direction of steepest descent, a minimizer for J will be called a continuation of the quasi-static evolution at t. The first variation of J is δJ(γ̇t, δγ̇) = ∫ l 0 ft(x) δγ̇(x) dx + τ ∫ l 0 ( θ′′(γt(x)) γ̇t(x) δγ̇(x)− σ̇t δγ̇(x) + α γ̇′t(x) δγ̇′(x) ) dx = ∫ l 0 ft(x) δγ̇(x) dx + τ (∫ l 0 ( θ′′(γt(x)) γ̇t(x)−σ̇t−αγ̇′′t (x) ) δγ̇(x) dx −α ∑ i [[γ̇′t]](xi) δγ̇(xi) ) = ∫ l 0 (ft(x) + τ ḟt(x)) δγ̇(x) dx − τα ∑ i [[γ̇′t]](xi) δγ̇(xi) . The derivative γ̇′′t (x) which appears in ḟt(x) is the regular part of the distributional derivative of γ̇′t, while the singular part is given by the jumps [[γ̇′t]](xi) at the jump points xi. Necessary conditions for a minimum at γ̇ = γ̇t are (3.20) ( ft(x) + τ ḟt(x) ) δγ̇(x) ≥ 0 a.e. x ∈ (0, l) [[γ̇′t]](x) δγ̇(x) ≤ 0 ∀x ∈ (0, l) . 20 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH The jump condition have been imposed to all x and not only to the jump points, since [[γ̇′t]](x) is zero at all x which are not jump points. By (3.17), at points at which γ̇t(x) = 0 a perturbation δγ̇(x) is arbitrary po- sitive, and therefore ( ft(x)+ τ ḟt(x) ) must be non-negative and [[γ̇′t]](x) must be non-positive. At points at which γ̇t(x) > 0 a perturbation may have any sign, and therefore the same terms must be zero. Then, the sets of Kuhn-Tucker conditions for ft (3.21) γ̇t(x) ≥ 0 , ft(x) + τ ḟt(x) ≥ 0 ,( ft(x) + τ ḟt(x) ) γ̇t(x) = 0 , a.e. x ∈ (0, l) , and for the jumps of γ̇′t (3.22) γ̇t(x) ≥ 0, [[γ̇′t]](x) ≤ 0 , [[γ̇′t]](x) γ̇t(x) = 0 , ∀x ∈ (0, l) , follow. At every jump point xo for γ̇ ′ t, [[γ̇′t]](xo) = γ̇ ′ t(xo+)−γ̇′t(xo−) = lim ε→0+ γ̇t(xo+ ε)−γ̇t(xo) ε − lim ε→0− γ̇t(xo+ ε)−γ̇t(xo) ε . For γ̇t(xo) = 0 the right-hand side is non-negative, and therefore [[γ̇ ′ t]](xo) is non- negative. But it is also non-positive by (3.22)2. Hence, a point at which γ̇t(xo) = 0 cannot be a jump point for γ̇′t. On the other hand, when γ̇t(xo) is positive, [[γ̇ ′ t]](xo) is zero by the complementarity condition (3.22)3. Then [[γ̇ ′ t]](xo) is zero in all cases, that is, there is no creation of new jump points for γ̇′t in the second-order minimization. In particular, if γ′t is continuous at the initial time t0, it remains continuous at all subsequent t. In what follows, we only consider processes in which γ̇′t is continuous instead of piecewise continuous. The assumed continuity of γ̇′t requires a couple of comments. The first is that it does not imply the continuity of γ′′t . Indeed, by Hadamard’s condition, this is true only if dx/dt = 0. That is, a discontinuity of γ′′t is possible, provided that it does not change its position with t. The second comment is that the continuity of γ′′t has been deduced only for the discretized problem. It need not be preserved in the limit when the time step τ tends to zero. In fact, in the numerical simulations in Section 4 we will find a situation in which the continuity is not preserved in the limit. The solutions of (3.21) and (3.22) depend on the time step τ , and do not imply the conditions obtained from the first-order minimization. Specifically, the com- plementarity conditions (3.18)3 does not hold in the second-order minimization. Therefore, the solutions for τ > 0 are not, in general, directions of steepest descent. Another necessary condition for a minimum of J at γ̇t is that the second variation (3.23) δ2J(δγ̇) = ∫ l 0 ( θ′′(γt(x)) δγ̇2(x) + w′′(εt) δ¯̇γ2+ α δγ̇′ 2 (x) ) dx be non-negative for all δγ̇ for which the first variation δJ(γ̇t, δγ̇) = ∫ l 0 ( ft(x) + τ ḟt(x) ) δγ̇(x) dx is zero. For x in the elastic zone, ft(x) is positive, and (ft(x) + τ ḟt(x)) is positive for sufficiently small τ . Then, by (3.20)1, for sufficiently small τ the first variation is zero only for perturbations with δγ̇(x) = 0 almost everywhere in the elastic zone. DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 21 In the inelastic zone, where ft(x)=0 by definition, conditions (3.21) reduce to (3.24) γ̇t(x) ≥ 0 , ḟt(x) ≥ 0 , ḟt(x) γ̇t(x) = 0 , a.e. x ∈ J (βt, γt) , and, by (3.20)1, the first variation is zero only if (3.25) ḟt(x) δγ̇(x) = 0 almost everywhere in the inelastic zone. In any interval at which γ̇t(x) = 0, we have ḟt(x) = θ ′′(γt(x)) γ̇t(x)− σ̇t − αγ̇′′t (x) = −σ̇t . Because ḟt(x) ≥ 0 by (3.24)2, such intervals do not exist in a hardening response, σ̇t > 0, and may exist only with ḟt(x) = 0 in a perfect plastic response, and only with ḟt(x) > 0 in a softening response. Then in a softening response condition (3.25) requires δγ̇(x) = 0. Leaving aside the exceptional case of a perfectly plastic response, we conclude that the first variation is zero only if δγ̇(x) is zero in the elastic zone and at almost all points of the inelastic zone at which γ̇t(x) = 0. Therefore, it is sufficient to consider perturbations δγ̇ with support in the part of the inelastic zone at which γ̇t > 0. Let (a, a+ lj) be an interval in this region. A necessary condition for a minimum is that the integral (3.26) ∫ a+lj a ( θ′′(γt(x)) δγ̇2(x) + l lj w′′(εt) δ¯̇γ2 + α δγ̇′ 2 (x) ) dx be non-negative for all perturbations δγ̇ with δγ̇(a) = δγ̇(a + lj) = 0. This corre- sponds to the non-negativeness of the smallest eigenvalue ρ1 of problem (A.15) in the Appendix, with H as in (A.5), and with a=0 , lj = L , θ ′′(γt(x))=αh(x) , lj l w′′(εt)=αω , δγ̇(x)=y(x) , δ¯̇γ = lj l ȳ . The last equality follows from the definitions ȳ = 1 lj ∫ lj 0 y(x) dx , δ¯̇γ = 1 l ∫ l 0 δγ̇(x) dx = 1 l ∫ lj 0 δγ̇(x) dx , which hold for δγ̇(x) = y(x) with support in (0, lj). Then the relation between w′′(εt) and ω comes from the identification l lj w′′(εt) δ¯̇γ2 = αω ȳ2 . The upper bound (A.17) for ρ1 provides the necessary condition (3.27) αλ2oj + θ ′′ oj ≥ 0 , where λoj is the smallest eigenvalue of problem (A.1). It is the solution of (A.7), now rewritten in the form lj l w′′(εt) α l2j = λ3oj l 3 j λoj lj − 2 tanλoj lj/2 , and θ′′oj is as in (3.12), with yoj given by (A.8). Because the first variation δJ(γ̇t, yoj) is zero and J is quadratic, one has J(γ̇t + c yoj) = J(γ̇t) + 1 2 τc2 δ2J(γ̇t, yoj) = J(γ̇t) + 1 2 τc2ρ1 ∫ l 0 y2oj(x) dx , 22 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH for every positive constant c. Therefore, ρ1 ≥ 0 is a necessary condition for a minimum of J . Moreover, for any other perturbation with support in (0, lj), since the first variation is non-negative and ρ1 is the smallest eigenvalue, J(γ̇t+δγ̇) ≥ J(γ̇t) + 1 2 τδ2J(γ̇t, δγ̇) ≥ J(γ̇t) + 1 2 τρ1 ∫ l 0 δγ̇2(x) dx , Therefore, ρ1 ≥ 0 is also a sufficient condition for a minimum. The lower bound (A.19) for ρ1 provides the explicit sufficient condition (3.28) αλ2oj + θ ′′ min ≥ 0 , with θ′′min = inf x∈(0,lj) θ′′(γt(x)) . For ρ1 negative, J(γ̇t + c yoj) can take unlimited negative values. Like in the local model, this event corresponds to brittle fracture. Therefore, ρ1 ≥ 0 is a necessary and sufficient condition for a minimum of J , and a negative ρ1 corresponds to brittle fracture. It is of interest to see in which cases the continuations are elastic, γ̇t = 0. From (3.19), J(γ̇) = ∫ l 0 ( ft(x)−τ w′′(εt) β̇t ) γ̇(x) dx + 1 2 τ ∫ l 0 ( θ′′(γt(x)) γ̇2(x) + w′′(εt) ¯̇γ2 + αγ̇′2(x) ) dx ≥ ∫ l 0 ( ft(x)−τ w′′(εt) β̇t ) γ̇(x) dx + 1 2 τ αρ1 ∫ l 0 γ̇2(x) dx ≥ ∫ l 0 ( c γ̇(x) + 1 2 ταρ1 γ̇ 2(x) ) dx , where c . = fmin − τw′′(εt) β̇t , fmin .= inf x∈(0,l) ft(x) . Because J(0) = 0, γ̇ = 0 is a minimizer if the right-hand side of the inequality is non-negative. In particular, γ̇ = 0 is a minimizer if both ρ1 and c are non-negative. There are two remarkable cases in which c is positive. The first is the case of β̇t < 0. That γ̇t = 0 for β̇t < 0 is the non-local version of the property (2.31) of elastic unloading for the local model. The second is the case of a strictly elastic regime, that is, of configurations (βt, γt) for which fmin > 0. Indeed, in this case c is positive for sufficiently small τ . This is a relaxed version of the condition found in the local model, that there is no increase of inelastic deformation in the elastic regime. Here, the same is true only in a strictly elastic regime. The foregoing discussion can be summarized as follows: (3.29) ρ1 ≥ 0 ⇐⇒ there are minimizers , ρ1 < 0 ⇐⇒ brittle fracture . ρ1 ≥ 0 and β̇t < 0 =⇒ γ̇t = 0 is a minimizer , ρ1 ≥ 0 and fmin > 0 =⇒ γ̇t = 0 is a minimizer . The last two implications tell us that, in a quasi-static evolution, fracture may occur only at loading, β̇t ≥ 0, and in a regime which is not strictly elastic. DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 23 To determine the explicit form of a continuation γ̇t when c < 0 is not as easy as it was in the local model. Indeed, in the region at which γ̇t(x) = 0, the in- troduction of the non-local term transforms the incremental equilibrium equation σ̇t = θ ′′(γt(x)) γ̇t(x) into the differential equation σ̇t = θ ′′(γt(x)) γ̇t(x)− αγ̇′′t (x) . This region is not known a priori, since its determination is a part of the solution of the problem. Moreover, only in some special cases a closed-form solution is available. One of them is the onset of the inelastic regime, studied in the next Subsection. 3.5. The onset of the inelastic regime. Consider a load process t → βt from the natural configuration (βt0 , γt0) = (0, 0), with β̇t > 0 for all t. At t = t0 we have ft0(x) = θ ′(γt0(x))−σt0−αγ′′t0(x) = θ′(0) > 0 , so that the whole bar is in a strictly elastic regime. By continuity, this regime persists over a finite time interval (t0, tc). In this interval γ̇t = 0 by (3.29)4, and γt = 0 because of the initial condition γt0 = 0. The deformation of the bar is homogeneous, with (3.30) u′(x, t) = εt = βt , σt = w′(βt) , ft = θ′(0)− w′(βt) . The elastic regime ends at the time tc at which β reaches the critical value βc given by (2.37). At this time ft becomes zero, and all points of the bar switch from the elastic to the inelastic zone. This is the onset of the inelastic regime. To determine the continuation γ̇t at t = tc, we start from the complementarity condition (3.24)3 (3.31) 0 = ḟ(x) γ̇(x) = θ′′(0) γ̇2(x)− σ̇ γ̇(x)− αγ̇′′(x) γ̇(x) . Here and in the following, for simplicity, we omit all subscripts tc. Integrating over (0, l) and recalling that γ̇′ has no jumps if γ̇ is a minimizer for J , we get σ̇ ∫ l 0 γ̇(x) dx = ∫ l 0 ( θ′′(0) γ̇2(x) + α γ̇′2(x) ) dx . If θ′′(0)≥ 0, the right-hand side is positive and therefore σ̇ is positive. Moreover, γ̇(x) is strictly positive almost everywhere in (0, l). Indeed, if γ̇(x) = 0 in some interval (a, b), from inequality (3.24)2 on that interval we have 0 ≤ ḟ(x) = −σ̇ , in contradiction with the positiveness of σ̇. Therefore, γ̇ is a solution of the diffe- rential problem (3.32) θ′′(0) γ̇(x) − σ̇ − α γ̇′′(x) = 0 ∀x ∈ (0, l) , γ̇(0) = γ̇(l) = 0 , subject to the dissipation condition γ̇(x) ≥ 0. For θ′′(0) > 0 the solution is (3.33) γ̇(x) = σ̇ θ′′(0) ( 1− coshκ (l/2− x) coshκl/2 ) , with κ = (θ′′(0)/α)1/2. By integration over (0, l), ¯̇γ = σ̇ θ′′(0) ϕ(κl) , ϕ(κl) . = 1− tanhκl/2 κl/2 , 24 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH and, because σ̇ = w′′(βc)(β̇ − ¯̇γ), (3.34) σ̇ = θ′′(0)w′′(βc) θ′′(0) + ϕ(κl)w′′(βc) β̇ . Comparing with the solution (2.34) of the local model, we see that the non-local effect is concentrated in the factor ϕ(κl). Because this factor is positive and less than one, the slope σ̇/β̇ of the response curve (σ, β) is positive and greater than the slope predicted by the local model, but smaller than the slope w′′(βc) at unloading. For θ′′(0) = 0, the solution is (3.35) γ̇(x) = σ̇ 2α x (l − x) . By integration over (0, l), (3.36) ¯̇γ = σ̇ l2 12α , and from σ̇ = w′′(βc)(β̇ − ¯̇γ), (3.37) σ̇ = 12αw′′(βc) 12α+ l2w′′(βc) β̇ . Thus, for θ′′(0) = 0 the perfectly plastic response σ̇ = 0 of the local model is replaced by a work-hardening response. For θ′′(0) < 0, the solution is (3.38) γ̇(x) = w′′(βc) (β̇ − ¯̇γ) θ′′(0) ( 1− cos k(l/2− x) cos kl/2 ) , with k = (−θ′′(0)/α)1/2. Integrating over (0, l) we find (3.39) ¯̇γ = w′′(βc)(β̇ − ¯̇γ) θ′′(0) ψ(kl) , ψ(kl) = 1− tan kl/2 kl/2 , and, therefore, (3.40) ¯̇γ = ψ(kl) w′′(βc) θ′′(0) + ψ(kl)w′′(βc) β̇ . Again, the difference with (2.33) is due to a single factor, which now is ψ(kl). This factor is negative for kl < π and positive for π < kl < 2.861π, with a jump from −∞ to +∞ at kl = π. The solution (3.38) satisfies the dissipation inequality (2.10) for all x only if kl ≤ 2π. Indeed, for kl > 2π, γ̇(x) takes negative values near the boundary. For such kl, let us consider the possibility of localized solutions, that is, of solutions of equation (3.31) which satisfy the differential equation (3.32) in a subinterval (a, a + li) of (0, l) and are zero outside. Without loss of generality, we take the interval (0, li), with li < l. The boundary conditions are (3.41) γ̇(0) = γ̇(li) = 0 , γ̇ ′(li) = 0 , where γ̇(0) = 0 is condition (3.15)1, and the two remaining conditions are due to the continuity of γ̇ and γ̇′, which are zero on (li, l). The first two conditions determine a solution of the form (3.38) with l replaced by li, and the third determines the length (3.42) li = 2π/ k . DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 25 Therefore, the solution is (3.43) γ̇(x) = w′′(βc) (β̇ − ¯̇γ) θ′′(0) (1− cos kx) , x ∈ (0, li) . It holds for li < l, that is, for kl > 2π. For kl = 2π, it coincides with (3.38). By integration over (0, li) it follows that l¯̇γ = li w′′(βc) (β̇ − ¯̇γ) θ′′(0) . and, therefore, (3.44) ¯̇γ = w′′(βc) w′′(βc) + θ′′(0) l/ li β̇ . It is convenient to introduce the positive constant (3.45) ψf . = − θ ′′(0) w′′(βc) . With this constant, and with ψo as in (A.13), equations (3.40) and (3.44) take the common form (3.46) ¯̇γ = ψo(kl) ψo(kl)− ψf β̇ . For kl < π, ψo(kl) is negative. Therefore, the condition ¯̇γ ≥ 0 is satisfied. For kl > π the numerator becomes positive, and ¯̇γ is positive only if (3.47) ψo(kl) > ψf . From (3.46), using (2.32), we get the incremental response law (3.48) σ̇ = θ′′(0) ψo(kl)− ψf β̇ . The slope σ̇/β̇ of the response curve is positive if kl < π, and negative if kl > π and condition (3.47) holds. Therefore, at the onset of the inelastic regime, for θ′′(0) negative the continuation is work-hardening, if kl < π , strain-softening, if kl > π and ψo(kl) > ψf . In the separating case kl = π the continuation is perfectly plastic. For kl > π and ψo(kl) ≤ ψf , by (3.47), there is no continuation obeying the dissipation inequality (2.10). We already know that in the solution (3.40) γ̇(x) is zero only at the boundary, while in the solution (3.43) γ̇(x) is zero on a portion of (0, l) of positive length (l − li). We say that the first is a full-size solution, and that the second is a localized solution. Therefore, for θ′′(0) negative, at the onset of the inelastic regime the continuations are full-size, if kl ≤ 2π , localized, if kl > 2π . From (3.34) and (3.37) it is clear that for θ′′(0) ≥ 0 all continuations are work- hardening and full-size. The continuations (3.38), (3.43) have been determined using the Kuhn-Tucker conditions (3.21), (3.22), which are necessary for a minimum for the functional J defined in (3.19). It remains to check whether or not these continuations are indeed 26 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH minimizers. At the onset, ft is zero. Then, neglecting the factor τ/2, J has the form J(γ̇) = ∫ l 0 ( θ′′(0) γ̇2(x)+w′′(βc) ¯̇γ2+αγ̇′2(x)− 2w′′(βc) β̇tγ̇(x) ) dx , and direct computation yields J(γ̇+δγ̇) = J(γ̇)+2 ∫ l 0 ḟt(x) δγ̇(x) dx+ ∫ l 0 ( θ′′(0) δγ̇2(x)+w′′(βc) δ¯̇γ2+α δγ̇′2(x) ) dx . On the right side, the first integral is non-negative by (3.20)1, and the second integral has the form (3.26). Then we can use the results of Section 3.4. In the present case, both constants θ′′oj , θ′′min which appear in the bounds (3.27), (3.28) are equal to θ′′(0), and the two bounds coincide. They provide the necessary and sufficient condition αλ2o ≥ −θ′′(0), that is, (3.49) λ2o ≥ k2 , where λ2o is the smallest eigenvalue of problem (A.1), with H as in (A.5). Consider first the case π < kl ≤ 2π, for which the continuation (3.38) is positive over (0, l). In this case, λo is the solution of equation (A.7) with L = l , ω = w′′(βc)/α , y(x) = δγ̇(x) , that is, (3.50) w′′(βc) α l2 = λ2ol 2 ψ(λol) . Recalling the definition (3.45) of ψf , we have (3.51) w′′(βc) α = −θ′′(0) αψf = k2 ψf , and, therefore, ψf = k2 λ2o ψ(λol) . If k ≤ λo, then ψ(λol) ≤ ψ(kl) = ψo(kl) , because ψ is a decreasing function and, by (A.13), coincides with ψo for π < kl ≤ 2π. Then, k ≤ λo implies inequality (3.47). Conversely, if k > λo the above inequality is reversed. That is, for π < kl ≤ 2π the condition (3.49) for a minimum coincides with the condition (3.47) for a non-negative ¯̇γ. Now consider the case kl > 2π, for which the continuation (3.43) is positive only on (0, li), with li = 2π/k < l. In this case, λo is the solution of equation (A.7) with L = li , ω = li l w′′(βc) α , y(x) = δγ̇(x) , ȳ = l li δ¯̇γ , as shown in Subsection 3.4. Using (3.51), (A.7) takes the form li l k2 ψf l2i = λ2o l 2 i ψ(λoli) , and recalling that for kl > 2π, by (A.13), ψo(kl) = 2π kl = li l , DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 27 hardening full-size softening full-size softening localized brittle fracture (1) (2) c ll i ll 0 ½ 1 0 1 Figure 1. The response at the onset of the inelastic deforma- tion. The three regions above the curve correspond to hardening full-size, softening full-size, and softening localized response. The region below the curve corresponds to totally brittle fracture. The lines (1), (2) represent families of bars made of the same material and with varying length we finally get ψf = k2 λ2o ψo(kl)ψ(λoli) . If k ≤ λo, then λoli ≥ kli = 2π. Therefore, ψ(λoli) ≤ ψ(2π) = 1 , and inequality (3.47) follows. If k > λo, then 2π > λoli, and the reversed inequality holds. Thus, the coincidence of (3.47) with (3.49) is proved for all kl > π. The incremental response law (3.48) shows that, when these inequalities are satisfied as equalites, the slope σ̇/β̇ of the response curve becomes infinite. This is the totally brittle fracture, that is, the catastrophic fracture at the onset of the inelastic deformation already described by the local model. The difference is that, while in the local model brittle fracture occurs as soon as θ′′ ceases to be positive, here it occurs at negative values of θ′′. Therefore, the sharp inequality (3.47) is a safety condition against totally brittle fracture. It is commonly asserted that the non-local energy term introduces an internal length of the material. Here, this length can be identified with the length li defined by (3.42), which is the length of the localization zone in a localized continuation. As 28 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH seen before, the ratio l/ li = kl/2π determines the quality of the response, hardening or softening, full-size or localized. The present analysis suggests the definition of a second characteristic length (3.52) lc = 2π √ α w′′(βc) , with which the no-fracture condition (3.47) takes the form ψo ( 2π l/ li ) > l2c/ l 2 i , that is, (3.53) l/ lc > l/ li√ ψo(2π l/ li) . = Ψ ( l/ li ) . The parameters which determine the response at the onset are the bar length l and the material constants w′′(βc), θ′′(0), and α. With the above definitions, they reduce to the ratios l/ lc and l/ li. A complete representation of the response is given in Fig. 1. In it, each point of the plane represents a bar. The points located below the curve l/ lc = Ψ(l/ li) correspond to bars which undergo totally brittle fracture at the onset of the inelastic deformation. The region above the curve is divided into three subregions, corresponding to full-size and hardening (l/ li < 1/2), full-size and softening (1/2 < l/ li < 1), and localized and softening (l/ li > 1) continuations, respectively. Straight lines from the origin represent families of bars made of the same material and with l growing with the distance from the origin. All families have full-size hardening continuations for small l, softening continuations for intermediate l, and catastrophic failure for large l. Moreover, for families with slope li/ lc > 1 the softening continuations may be both full-size and localized, while for those with li/ lc < 1, corresponding to straight lines which do not cross the third region above the curve, the softening continuations can only be full-size. 3.6. Collapse mechanisms. In the non-local model there are two possible collapse mechanisms, brittle fracture and ductile fracture. Both are characterized by the slope σ̇/ β̇ of the response curve, which is −∞ for brittle fracture and zero for ductile fracture. Brittle fracture coincides with the unique collapse mode provided by the local model. As shown in the next Section, the two collapse modes are observed in steel and in concrete, respectively. In a quasi-static evolution t → (βt, γt), assume that θ′′(γt(x)) is negative at all x, and that it can be approximated by a negative constant θ′′t . Define the current values of k and ψf (3.54) kt . = √ −θ′′t α , ψft . = −θ′′t w′′(εt) . With them, the relations (3.46), (3.48), obtained for the onset of the inelastic regime, can be extrapolated to the subsequent evolution (3.55) ¯̇γt = ψo(ktl) ψo(ktl)− ψft β̇t , σ̇t = θ′′t ψo(ktl)− ψft β̇t , and condition (3.47) against brittle fracture can be replaced by (3.56) ψo(ktl) > ψft . DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 29 This extrapolation helps us to understand some general features of the bar’s re- sponse in the softening regime. Of course, the conclusions are purely qualitative, since the exact values of θ′′t and kt are not known. More precise results can only be obtained from numerical simulation. Assume that −θ′′ is monotonic increasing. For all x, the function t → −θ′′(γt(x)) is increasing as well, because t → γt(x) is increasing by the dissipation inequality (2.10). Then the average −θ′′t is increasing, kt is increasing, and ψo(ktl) is decrea- sing. Supposing w′′(εt) substantially independent of t, ψft is increasing. Therefore, the positive denominator in (3.55) decreases with t. When it becomes zero, the negative slope of the response curve becomes infinite, and catastrophic failure takes place. On the contrary, for −θ′′ monotonic decreasing, −θ′′t decreases with t, and the slope |σ̇t/ β̇t| in (3.55) decreases as well. If −θ′′t tends to zero, the slope tends to zero, and this corresponds to ductile fracture. Thus, brittle fracture takes place for −θ′′ increasing, that is, for θ′ concave, and ductile fracture takes place for θ′ convex. In spite of some early statements and of more recent observations, see [33] and ([38], Fig. 16), respectively, this correlation between fracture modes and the concavity-convexity properties of θ′ did not receive adequate attention in the literature. This correlation is illustrated by the picture in Fig. 1. If θ′′t increases with t, the internal length lit = 2π/ kt decreases and the ratio l/ lit increases, while l/ lc stays approximately constant. The point representing the current status of the bar moves horizontally to the right, tending to brittle fracture which, at the onset, occurs at the curve Ψ. On the contrary, if θ′′t decreases with t the point moves to the left, tending to the perfectly plastic response which, at the onset, takes place at the line l/li = 1/2. The function θ′ cannot be concave for all γ > 0. Indeed, concavity implies θ′(γ) ≤ θ′(0) + γ θ′′(0) , and for θ′′(0) < 0 the assumption θ′(γ) > 0 is violated for sufficiently large γ. On the other hand, the assumption (2.7)3 that θ has a finite limit at +∞ requires that lim γ→+∞ θ ′′(γ) = 0 . Then−θ′′ has a maximum at a finite γ. Brittle fracture occurs only if this maximum produces a kt sufficiently large to violate inequality (3.56). 4. Numerical simulations For the non-local model, the quasi-static evolution in the inelastic regime has been investigated by performing a number of numerical simulations. An iterative pro- cedure was used: after fixing a time step τ , for each t = nτ the continuation γ̇t at (βt, γt) was determined approximately by solving a sequence of discretized minimum problems, with the data taken from the solution at t = (n− 1) τ . In the case of loading, β̇ > 0, the parameter t is identified with the load β, so that β̇ = 1, and the time step τ is in fact a load step. The discretized problem consists in minimizating a quadratic approximation of the energy Eeq(βt+τβ̇, γt+τ γ̇). For 30 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH convenience, we refer to the more general functional (4.1) F (β, γ, τ, δγ) = ∫ l 0 ( θ(γ(x)) + w(β − γ̄) + 1 2 αγ′2(x) ) dx τ ∫ l 0 ( θ′(γt(x)) δγ(x)+w′(β−γ̄)(1−δγ(x)) + αγ′(x) δγ′(x) ) dx + 1 2 τ2 (∫ l 0 ( θ′′(γ(x)) δγ2(x) + α δγ′2(x) ) dx+ lw′′(β − γ̄)(1− δγ̄)2 ) , such that F (βt, γt, 0, 0) coincides with E eq(βt, γt), and F (βt, γt, τ, γ̇) coincides with the second-order approximation (2.28) of Eeq(βt+τβ̇, γt+τ γ̇) for the given τ . 4.1. The algorithm. In the solution procedure an iterative algorithm, imple- mented in an ad hoc finite element code, was used. At each iteration step, for a given initial configuration (β0, γ0), with γ0 a local minimizer for F (β0, · , 0 , 0), the confi- guration (β0+ τ, γ0+ τ γ̇0) is determined by minimizing the function F (β0, γ0, τ, · ), constrained by the dissipation inequality (2.10) and by the boundary conditions (3.6). At each load step, the solution accuracy is refined through an iterative pro- cedure acting on the quadratic approximation F2(β0, γ, δγ) . = F (β0, γ, 0 , 0) + ∂ ∂(δγ) F (β0, γ, 0 , 0) δγ + 1 2 ∂2 ∂(δγ)2 F (β0, γ, 0 , 0) δγ 2 of F (β0, γ, 0 , δγ). The algorithm runs as follows. 1. Incremental step (i) Compute γ̇0 = argmin { F (β0, γ0, τ, · ) | δγ(x) ≥ 0, δγ(0) = δγ(l) = 0 } . (ii) Set β1 = β0 + τ , γ 0 1 = γ0 + τ γ̇0 . 2. Iterative refinement (i) Compute δγ1 = argmin { F2(β1, γ 0 1 , · ) | δγ(x) ≥ γ0(x)−γ01(x), δγ(0) = δγ(l) = 0 } , and set γ11 = γ 0 1 + δγ 1. (ii) Compute δγi = argmin { F2(β1, γ i−1 1 , · ) | δγ(x) ≥ γ0(x)−γi−11 (x), δγ(0) = δγ(l) = 0 } , and set γi1 = γ i−1 1 + δγ i. (iii) Stop when the L2 norm of (γi1 − γi−11 ) is less than a given tolerance γ̂. 3. End (i) Take as γ1 the last γ i 1 in the preceding iteration. (ii) If the L2 norm of (γ1 − γ0) is less than a given tolerance γ̃, perform a new incremental step from (β1, γ1). Otherwise, repeat the computation with τ replaced by τ/2. DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 31 The last control avoids the overcoming of energy barriers due to exceedingly large incremental steps. In the numerical code, the spatial variable x is discretized using linear finite ele- ments. The quadratic programming problems involving the minimization of F2 are solved using the projection method [27, 45] implemented in the quadprog.m function of Matlab. The code generates a mesh refinement when the number of elements in the inelastic zone is smaller than a certain number (100 in the simulations presented below). In this case, each element is split into two sub-elements. 4.2. Choice of the energy densities. For the elastic strain energy density we take the quadratic expression (4.2) w(ε) = 12 EAε 2, where the axial stiffness EA is the product of the Young modulus E of the material by the area A of the cross section. For the cohesive energy density we take the piecewise cubic function (4.3) θ(γ) = Ai +Biγ + 1 2 Ciγ 2 + 1 6 Diγ 3 , γ ∈ [γi−1, γi ] , where { [γi−1, γi ] , i = 1, . . . , n, } is a subdivision of the interval [γ0, γn] into n adjacent subintervals, with γ0 = 0 and γn ≤ +∞. The coefficients Ai, Bi, Ci and Di satisfy the 3(n− 1) conditions Bi+1 = Bi − 3 (Ai+1 −Ai) γ−1i , Ci+1 = Ci + 6 (Ai+1 −Ai) γ−2i , Di+1 = Di − 6 (Ai+1 −Ai) γ−3i , i = 1, . . . (n− 1) , which guarantee the continuity of θ, θ′ and θ′′ at γi. We fix A1 = θ(0) = 0, while B1 = θ ′(0), the value of the axial force at the onset of the inelastic regime, is identified on the experimental curve. The remaining constants are determined by fixing n + 1 coefficients, for example, C1, D1, and Aj for j = 2, ..., n. In the following simulations these constants are selected with the purpose of reproducing the response curves of two specific experimental tests, one on a steel bar and one on a concrete specimen. 4.3. Simulation 1: the tensile response of a steel bar. For the first simula- tion, the data were obtained from a test on a steel bar, stretched to rupture under controlled end displacements. The test was made in Ancona, at the Laboratorio Prove Materiali e Strutture of the Università Politecnica delle Marche. The speci- men was a ribbed bar, made of B450C (FeB44k) steel, with diameter φ = 16 mm and length l = 200 mm. The axial stiffness EA = 42× 103 kN is the product of the Young modulus of steel, E = 210 kN/mm2, by the cross- sectional area A = πφ2/4. The axial force at the onset of the inelastic deformation, measured on the response diagram, is B1 = 109.5 kN . 32 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH With the purpose of investigating the size effect, simulations with three different lengths l = 100, 200, 300 mm , were made. The non-locality parameter α was taken equal to α = 10 kNmm2. This is a tentative value. The criteria for choosing the value of α are still under investigation. For the time step τ , the mesh size h, and the tolerances γ̂ and γ̃ introduced in Subsection 4.1, the values τ = 10−4, h = 0.5mm , γ̂ = 10−6 , γ̃ = 10−3 , were chosen. These values provide a good accuracy and, at the same time, keep a moderate size for the discretized problem. For the cohesive energy density a function θ1 of the form (4.3) was taken, with n = 4 , γi = (0.10, 1.62, 2.01, 10) , C1 = 380 kN , D1 = −3800 kN , Ai = (0, −0.61, 499.69, 223.74) kN . This energy, plotted in Fig. 2a, is convex for γ < 0.1 and concave for γ > 0.1. 0 0.5 1 1.5 2 2.5 0 50 100 150 200 250 0 0.5 1 1.5 2 2.5 0 20 40 60 80 100 120 140 0.10 2.101.62 (a) (b) 1 ' 1 Figure 2. Simulation 1. The cohesive energy θ1 (a), and its de- rivative θ′1 (b) In Fig. 3, the experimental response curve (dotted line) is compared with the response curves of simulations with lengths l = 100, 200, 300mm. In the experimen- tal curve, the initial hardening regime is followed by a softening regime, in which the negative slope of the curve increases with β, and becomes infinite at rupture. DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 33 In the model, this behavior is reproduced by taking θ1 initially convex, and then concave with growing concavity. The response curve of the simulation with l = 200 mm, which is the specimen’s length, is very close to the experimental curve. Of course, this is due to an ap- propriate choice of the material constants. Anyway, it is remarkable that a small number of constants is sufficient to reproduce the evolution of the bar’s deforma- tion from the natural configuration to final rupture. Also remarkable is the perfect reproduction of the elastic unloading, shown by the dashed line in the figure, and obtained by reversing the sign of β̇ up to complete unloading, σ = 0, and then reloading. As shown in the figure, the reloading curve traces back the unloading line and then continues along the loading curve of monotonic loading. In the three simulations shown in the figure, the hardening parts of the response curves are identical. The hardening regime ends at β ≈ 0.102, where the force σ reaches its maximum and a softening regime starts. Shortly after, strain localization begins, and the three curves separate. The figure shows that for longer bars the softening regime produces larger negative slopes, and rupture occurs at smaller β. This size effect agrees with the theoretical preview (3.55)2. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 20 40 60 80 100 120 140 β [kN] l=300 l=200 l=100 Figure 3. Simulation 1. The experimental response curve (dotted line), and the response curves of simulations with different values of l. The dashed line shows the model’s response to unloading and reloading The evolutions of γ and γ̇ in the simulation with l = 200 mm are shown in Fig. 4. In the initial hardening regime, represented in the upper part of the picture, both γ and γ̇ are constant, except near the boundary, where γ(x) = γ̇(x) = 0 as required by the boundary conditions. The regime of full-size softening response is very short. Indeed, it starts at β ≈ 0.102 and ends at β ≈ 0.103 when, as shown in Fig. 4b and 4d, localization takes place. For larger β, the deformation concentrates at shorter and shorter central zones of the bar. Note the strict similarity of the diagrams of Fig. 4c, reproduced in Fig. 5a, with the results of experiments made long ago by Miklowitz 34 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH [40], shown in Fig. 5b. Also, the evolution curves for γ are qualitatively the same as the curves in Fig. 17.1 of [38]. 0 20 40 60 80 100 120 140 160 180 200 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 20 40 60 80 100 120 140 160 180 200 0 0.5 1 1.5 2 75 80 85 90 95 100 105 110 115 120 125 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 75 80 85 90 95 100 105 110 115 120 125 0 100 200 300 400 500 600 700 x =0.103 =0.040 =0.080 =0.100 =0.103 =0.040 =0.080 =0.100 x xx ( b )( a ) ( d )( c ) =0.104 =0.110 =0.115 =0.120 =0.125 =0.126 =0.104 =0.110 =0.115 =0.120 =0.125 =0.126 Figure 4. Simulation 1 with l = 200 mm. Evolutions of γ (a), (c), and γ̇ (b), (d) for increasing β 80 85 90 95 100 105 110 115 120 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) (b) u =11.40meas u =12.40meas u =12.84meas u =13.15meas u =10.40meas x MINIMUM SECTION OF NECK LONGITUDINAL AXIS OF BAR U N IT C O N V E N T IO N A L S T R A IN PLOT ALONG ORIGINAL POSITIONS IN GAGE LENGTH Figure 5. The diagrams of Fig. 4c, compared with the experi- ments on a cylindrical bar of medium-carbon steel (from [40]) 4.4. Simulation 2: the tensile response of a concrete specimen. The aim of our second simulation was to reproduce the experimental force-elongation curve of a prismatic specimen made of lightweight concrete, taken from Fig. 38 of [35]. The specimen is prismatic, with dimensions 50 × 50 × 150 mm, Young modulus E = 18 kN/mm2, and the maximum aggregate size is 8 mm. In the simulation we take EA = 45× 103 kN , B1 = 6.3 kN, DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 35 with B1 the axial force at the onset of the inelastic regime, read on the experimental curve. We consider the three lengths l = 100, 150, 200 mm . For the non-locality parameter α we take α = 3497 kNmm2 . This value corresponds to an internal length li = 2π/k = 8 mm, equal to the maximum aggregate size. For the time step τ , the mesh size h, and the tolerances γ̂ and γ̃ we choose τ = 10−7, h = 0.25mm , γ̂ = 10−6, γ̃ = 10−3. For the cohesive energy we take the function θ2 represented in Fig. 6. This is an energy of the form (4.3), with n = 4 , γi = (0.7, 1.2, 65.0, 100.0)× 10−4, and with C1 = 18× 103 kN , D1 = −25.7143× 107 kN , Ai = (0, −0.12, −0.25, 154.27)× 10−4 kN . This energy is convex for γ < 0.7× 10−4 and concave for γ > 0.7× 10−4. 0 1 2 3 4 5 6 7 8 x 10 −3 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0 1 2 3 4 x 10 −4 0 0.5 1 1.5 2 2.5 x 10 −3 0 1 2 3 4 x 10 −4 6 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7 0.7 concaveconvex (a) (b) 2 ' 2 ' 2 2 0 1 2 3 4 5 6 7 8 x 10 −3 0 1 2 3 4 5 6 7 Figure 6. Simulation 2. The cohesive energy θ2 (a), and its de- rivative θ′2 (b) The force-elongation response curves for the three considered values of l are shown in Fig. 7. The correspondence with the experimental curve is good, but not as good as in the preceding simulation. This may be due to the fact that in the 36 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH experimental curve in [35] the elongation βmlm was measured on a central zone of the bar, of length lm = 35 mm, while here the total elongation was evaluated with the approximate formula βl ≈ βmlm + (l − lm)(β − γ̄) , obtained by neglecting the inelastic deformation outside the central zone. Another cause of discrepancy is that in the experiment the specimen had a central notch, not considered in the simulation. Generally, the data fitting was more difficult than in in Simulation 1, where all data came from our own experiment and not from one taken from the literature. As in the preceding simulation, the three numerical curves coincide in the elastic regime and in the subsequent hardening and full-size softening regimes. They start to differ at the onset of localization, at β ≈ 2.23× 10−4. In the localized softening regime, in accordance with the size effect, the negative slope of the response curve is larger for longer bars. The curve for l = 200 mm exhibits brittle fracture at the onset of the localization. This agrees with the safety condition (3.56) against brittle fracture, which now takes the form 2π ktl > αk2t EA . Indeed, in a localized solution ψo(ktl) is equal to 2π/ktl by (A.13), and ψft is is equal to αk2t/EA by (3.54) and (4.2). Setting kt = 2π/lt, with lt = 8 mm equal to the maximum aggregate size, condition (3.56) becomes l < 2πEA k3t α = EA l3t 4π2α ≈ 167 mm . For the bars with l = 100 mm and 150 mm, the convex shape of θ′2 determines a ductile fracture mechanism. For l = 150 mm, as shown in Fig. 8, at β ≈ 2.22×10−4 the inelastic strain rate γ̇ begins to localize and the force starts to decrease. The length of the localized zone attains a minimum of about 8 mm at β ≈ 2.25 ×10−4, and then it steadily increases. Numerical investigations not reported here show that the localization zone eventually diffuses over the whole bar. For the bar with l = 200 mm, catastrophic failure occurs at β ≈ 2.24× 10−4. The evolutions of γ and γ̇ are shown in Fig. 8. The dissymmetries of the diagrams of γ̇ in Fig. 8d are purely numerical, and their nature is random. This is confirmed by the symmetry of their time integral γ in Fig. 8c. For β ≈ 10−3 the inelastic deformation γ is concentrated on a portion of the bar of length of about 20 mm. This agrees with an estimate given in [9], according to which the characteristic length of concrete is about 2.7 times the maximum aggregate size. By effect of the expansion of the inelastic zone the force decreases, and tends to a horizontal asymptote. In [36, 37], the expansion of the inelastic zone is considered as unphysical, and the stress locking, that is, the convergence of the force to a positive limit, is considered as a proof of this. While the horizontal asymptotes of the response curve shown in [36], Fig. 3, are about 80% of the maximum σ, a numerical investigation, not reported here, shows that in the present simulation the horizontal asymptote is only about 2% of the maximum σ. This modest residual force cannot be eliminated. Indeed, looking at the diagrams of γ in Fig. 8c we see that, for β > 4.5 × 10−4, the largest value γ(l/2) of γ is greater than 6× 10−3. Therefore, as shown by the constitutive curve DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 37 0 0.2 0.4 0.6 0.8 x 10 −3 0 1 2 3 4 5 6 7 8 β [kN] l=150 l=100 l=200 2.22 2.24 2.26 2.28 2.3 x 10 −4 6.6 6.7 6.8 6.9 l=150 l=100 l=200 Figure 7. Simulation 2. The experimental (dotted line) and the numerical response curve (solid line) of Fig. 6b, the corresponding value of θ′ is zero. Because the point x = l/2 is in the inelastic zone, from the yield condition (3.2) with θ′ = 0 it follows that σ = −αγ′′(l/2) , with γ′′(l/2) strictly negative because γ has a maximum at x = l/2. 0 50 100 150 0 1 2 3 4 5 6 7 8 x 10 −5 0 50 100 150 0 5 10 15 20 x x ( b ) ( a ) =1.90*10 -4 =2.22*10 -4 =2.23*10 -4 =2.24*10 -4 =2.25*10 -4 =1.90*10 -4 =2.22*10 -4 =2.23*10 -4 =2.24*10 -4=2.25*10 -4 0 50 100 150 0 0.002 0.004 0.006 0.008 0.01 0.012 0 50 100 150 0 5 10 15 20 25 30 35 40 45 50 x x ( d ) ( c ) =11.00*10 -4 =8.00*10 -4 =5.00*10 -4 =3.00*10 -4 =11.00*10 -4 =8.00*10 -4 =5.00*10 -4 =3.00*10 -4 Figure 8. Simulation 2. Evolution of the inelastic deformation for l = 150 mm 38 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH Appendix A. Solution of some eigenvalue problems Consider the eigenvalue problem (A.1) ∫ L 0 ( y′2(x) + ω ȳ2 ) dx = λ2 ∫ L 0 y2(x) dx , y ∈ H, with ω a given positive constant, and (A.2) ȳ = 1 L ∫ L 0 y(x) dx . The eigenfunctions are the solutions of the Euler equation (A.3) y′′(x) − ω ȳ + λ2y(x) = 0 , and have the form (A.4) y(x) = A sinλx +B cosλx+ ω λ−2ȳ . We wish to determine the smallest eigenvalue λ2o and the corresponding eigenfun- ctions yo, for two different choices of the domain H . (i) Let H be the set (A.5) H = { y ∈ H1(0, L) | y(0) = y(L) = 0} . In this case, the constants A,B are (A.6) A = −ω λ−2ȳ tanλL/2 , B = −ω λ−2ȳ . Moreover, integrating (A.4) over (0, L) we get L ȳ = ω λ−3ȳ ( (cosλL− 1) tanλL/2− sinλL+ λL) = ω λ−3ȳ (λL− 2 tanλL/2) , that is, (A.7) ω L2 = λ3L3 λL− 2 tanλL/2 = λ2L2 ψ(λL) , with ψ as in (3.39)2. The right-hand side is an increasing function of λL, and is equal to zero at λL = π and to +∞ at λL ≈ 2.861 π. Therefore, for every ωL2 > 0 there exists a unique λL in (π, 2.861π) which satisfies (A.7). It determines the smallest eigenvalue λ2o for given ω and L. The corresponding eigenfunctions are the scalar multiples of (A.8) yo(x) = cosλo(L/2− x)− cosλoL/2 . The problem has been formulated in H1(0, L), which is the largest space in which the integrals which appear in (A.1) exist. However, the eigenfunction yo belongs to C∞[0, L]. Then the solution does not change if H is taken to be any space included between C∞[0, L] and H1(0, L). (ii) Now let H be the set (A.9) H = { y ∈ H2(0, L) | y(0) = y(L) = 0 , y(x) ≥ 0 ∀x ∈ (0, L)} . In (A.8), yo satisfies condition y(x) ≥ 0 only for λoL ≤ 2π. By (A.7), this corre- sponds to ωL2 ≤ 4π2. For larger values of ωL2, yo(x) becomes negative near the boundary. DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 39 For ωL2 > 4π2 we look for solutions y which are of the form (A.4) in a subinterval (0, Ly) of (0, L) and zero outside. To belong to H 2(0, L), they must satisfy the conditions y(0) = y(Ly) = 0 , y ′(Ly) = 0 . By (A.4), the first two conditions are satisfied by A = −ω λ−2ȳ tanλLy/2 , B = −ω λ−2ȳ , and the third condition is satisfied if A = 0. This implies tanλLy/2 = 0, that is, (A.10) λLy = 2π . Because Ly is smaller than L, this solution is possible only for λL > 2π, that is, for ωL2 > 4π2. But this is exactly the range not covered by the solution (A.7), (A.8). For A = 0 and B = −ω λ−2ȳ, the integration of (A.4) yields (A.11) ωL2 = λ2L3 Ly = λ3L3 2π . Therefore, for all ωL2 > 4π2 the smallest eigenvalue is the solution λ2o of this equation. The corresponding eigenfunction (A.12) yo(x) = { 1− cosλox if x ≤ 2π/λo , 0 if x > 2π/λo , satisfies the condition yo(x) ≥ 0 for all x. Therefore, λ2o is the smallest eigenvalue, and yo is a corresponding eigenfunction. In conclusion, the smallest eigenvalue of problem (A.1) with H as in (A.9) is the solution of (A.7) if ωL2 ≤ 4π2, and of (A.11) if ωL2 > 4π2. By defining (A.13) ψo(λL) . = ⎧⎨ ⎩ ψ(λL) if 0 ≤ λL ≤ 2π , 2π λL if 2π < λL , both (A.7) and (A.11) take the form (A.14) λ2 = ω ψo(λL) , and the smallest eigenvalue λ2o is the unique solution of this equation. The function ψo is C ∞ in (π,+∞), except at the point λL = 2π, where it is continuous and has a continuous derivative. Now consider the eigenvalue problem (A.15) ∫ L 0 ( y′2(x) + ω ȳ2 + h(x) y2(x) ) dx = ρ ∫ L 0 y2(x) dx , y ∈ H , with ω a given positive constant, and h an integrable function. We wish to deter- mine upper and lower bounds for the smallest eigenvalue ρo, both for H as in (A.5) and as in (A.9). For H as in (A.5), let λ2o be the smallest eigenvalue of problem (A.1) and let yo be a corresponding eigenfunction. Define (A.16) ho . = ∫ L 0 h(x) y2o(x) dx∫ L 0 y2o(x) dx . 40 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH Then, by (A.1),∫ L 0 ( y′o 2 (x) + ω ȳ2o + ho y 2 o(x) ) dx = (λ2o+ ho) ∫ L 0 y2o(x) dx . That is, (λ2o+ho) is an eigenvalue of (A.15), and yo is a corresponding eigenfunction. Therefore, an upper bound for ρo is (A.17) λ2o+ ho ≥ ρo . To get a lower bound, set (A.18) hmin = inf x∈(0,L) h(x) . Then the left-hand side of (A.15) is greater or equal than∫ L 0 ( y′2(x) + ω ȳ2 + hmin y2(x) ) dx ≥ (λ2o+ hmin) ∫ L 0 y2(x) dx ∀y ∈ H , and the lower bound (A.19) ρo ≥ λ2o+ hmin follows. Depending on the given function h, the two bounds can be far away from each other. On the contrary, they coincide if h is a constant, h(x) = h̄. Indeed, in this case ho = hmin = h̄ and, therefore, (A.20) ρo = λ 2 o + h̄ . For H as in (A.9) the same bounds hold, with λ2o and yo replaced by the smallest eigenvalue of problem (A.1) and by a corresponding eigenfunction, respectively. References [1] E.C. Aifantis, On the microstructural origin of certain inelastic models, J. of Engng. Materials and Technology, ASME 106: 326-330, 1984 [2] E.C. Aifantis, On the role of gradients in the localization of deformation and fracture, Int. J. Engng. Sci. 30: 1279-1299, 1992 [3] L. Ambrosio, V.M. Tortorelli, Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence, Comm. Pure Appl. Math. 43: 999-1036, 1990 [4] J.-F. Babadjian, A quasistatic evolution model for the interaction between fracture and da- mage, Arch. Rat. Mech. Analysis 200: 945-1002, 2011 [5] G.I. Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, Adv. Appl. Mech. 7: 55-129, 1962 [6] Z.P. Bažant, E.P. Chen, Scaling of structural failure, Appl. Mech. Rev. 50: 593-627, 1997 [7] Z.P. Bažant, M. Jirásek, Nonlocal integral formulations of plasticity and damage: survey of progress, ASCE J. of Engrg. Mechanics 128: 1119-1149, 2002 [8] Z.P. Bažant, J.L. Le, Size effect on strength and lifetime distributions of quasibrittle struc- tures, Proc. ASME Int. Mech. Engrg. Congress IMECE2009, pp 1-9, 2009 [9] Z.P. Bažant, G. Pijaudier-Cabot, Measurement of characteristic length of nonlocal continu- um, ASCE J. of Engrg. Mechanics 115: 755-767, 1989 [10] A. Benallal, J.-J. Marigo, Bifurcation and stability issues in gradient theories with softening, Modelling and Simulation in Materials Scienc and Engineering 15: S283-S295, 2007 [11] B. Bourdin, G.A. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture, J. Mech. Phys. Solids 48: 797-826, 2000 [12] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy, J. Chem. Physics 28: 258-267, 1958 [13] A. Carpinteri, Application of fracture mechanics to concrete structures, ASCE J. Struct. Division 108: 833-848, 1982 DIFFUSE COHESIVE ENERGY FOR FRACTURE AND PLASTICITY 41 [14] A. Carpinteri, R. Massabò, Continuous vs discontinuous bridged-crack model for fiber- reinforced materials in flexure, Int. J. Solids Structures 34: 2321-2338, 1997 [15] P. Casal, La théorie du second gradient et la capillarité, C. R. Acad. Sc. Paris 274: A1571- 1574, 1972 [16] G. Dal Maso, R. Toader, A model for the quasi-static growth of brittle fractures based on local minimization, Arch. Rat. Mech. Analysis 162: 101-135, 2002 [17] G. Dal Maso, R. Toader, Quasi-static crack growth in elasto-plastic materials: the two- dimensional case, Arch. Rat. Mech. Analysis 196: 867-906, 2010 [18] R. De Borst, J. Pamin, Gradient plasticity in numerical simulation of concrete cracking, Eur. J. Mech. A/15: 295-320, 1996 [19] G. Del Piero, L. Truskinovsky, Elastic bars with cohesive energy, Cont. Mech. Thermody- namics 21: 141-171, 2009 [20] G. Del Piero, G. Lancioni, R. March, Diffuse cohesive energy in plasticity and fracture, Technische Mechanik 32: 174-188, 2012 [21] D.C. Drucker, A more fundamental approach to plastic stress strain relations, Proc. 1st U.S. Nat. Congr. Appl. Mech., ASME, pp 487-491, 1951 [22] J.L. Ericksen, Liquid crystals with variable degree of orientation, Arch. Rational Mech. Anal. 113: 97-120, 1990 [23] N.A. Fleck, J.W. Hutchinson, A reformulation of strain gradient plasticity, J. Mech. Phys. Solids 49: 2245-2271, 2001 [24] G.A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46: 1319-1342, 1998 [25] F.C. Frank, On the theory of liquid crystals, Discuss. Faraday Soc. 25: 19-28, 1958 [26] F. Freddi, G. Royer-Carfagni, Regularized variational theories of fracture: a unified approach, J. Mech. Phys. Solids 58: 1154-1174, 2010 [27] P.E. Gill, W. Murray, M.H. Wright, Practical optimization, Academic Press, 1981 [28] A.A. Griffith, The phenomena of rupture and flow in solids, Phil. Trans. Roy. Soc. A221: 163-198, 1920 [29] P. Gudmundson, A unified treatment of strain gradient plasticity, J. Mech. Phys. Solids 52: 1379-1406, 2004 [30] M.E. Gurtin, On a framework for small-deformation viscoplasticity: free energy, microforces, strain gradients, Int. J. of Plasticity 19: 47-90, 2003 [31] M.E. Gurtin, L. Anand, A theory of strain-gradient plasticity for isotropic, plastically irro- tational materials. Part I: Small deformations, J. Mech. Phys. Solids 53: 1624-1649, 2005 [32] R. Hill, The Mathematical Theory of Plasticity, Oxford University Press, Oxford, 1950 [33] A. Hillerborg, Application of the fictitious crack model to different types of materials, Int. J. of Fracture 51: 95-102, 1991 [34] A. Hillerborg, M. Modéer, P.E. Peterson, Analysis of crack formation and crac growth in concrete by means of fracture mechanics and finite elements, Cement and Concrete Research 6: 773-782, 1976 [35] D.A. Hordijik, Tensile and tensile fatigue behaviour of concrete; experiments, modelling and analyses, Heron 37: 1-79, 1992 [36] M. Jirásek, Nonlocal models for damage and fracture: comparison of approaches, Int. J. Solids Structures 35: 4133-4145, 1998 [37] M. Jirásek, S. Rolshoven, Localization properties of strain-softening gradient plasticity mod- els. Part I: strain-gradient theories, Int. J. Solids Structures 46: 2225-2238, 2009 [38] M. Jirásek, S. Rolshoven, Localization properties of strain-softening gradient plasticity mod- els. Part II: theories with gradients of internal variables, Int. J. Solids Structures 46: 2239- 2254, 2009 [39] J.-J. Marigo, L. Truskinovsky, Initiation and propagation of fracture in the models of Griffith and Barenblatt, Cont. Mech. Thermodynamics 16: 391-409, 2004 [40] J. Miklowitz, The influence of the dimensional factors on the mode of yielding and fracture in medium-carbon steel. II. The size of the round tensile bar, J. Appl. Mech 17: 159-168, 1950 [41] D. Mumford, J. Shah, Optimal approximations by piecewise smooth functions and associated variational problems, Comm. Pure Appl. Math. 42: 577-685, 1989 [42] K. Pham, J.-J. Marigo, Approche variationnelle de l’endommagement. I. Les concepts fon- damentaux, C. R. Mécanique 338: 191-198, 2010 42 GIANPIETRO DEL PIERO, GIOVANNI LANCIONI, AND RICCARDO MARCH [43] K. Pham, J.-J. Marigo, Approche variationnelle de l’endommagement. II. Les modèles à gradient, C. R. Mécanique 338: 199-206, 2010 [44] K. Pham, H. Amor, J.-J. Marigo, C. Maurini, Gradient damage models and their use to approximate brittle fracture, Int. J. Damage Mechanics 20: 618-652, 2011 [45] E. Polak, Computational methods in optimization, Academic Press, 1971 [46] L. Truskinovsky, Fracture as a phase transition, In: R.C. Batra, M.F. Beatty eds., Con- temporary Research in the Mechanics and Mathematics of Materials, CIMNE, Barcelona 1996 [47] J.D. Van der Waals, The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density (in Dutch), Verh. Kon. Akad. Veten. 33, N. 8, 1893 [48] K.Y. Volokh, Nonlinear elasticity for modeling fracture of isotropic brittle solids, J. Appl. Mech. 71: 141-143, 2004 [49] T. Yalcinkaya, W.A.M. Brekelmans, M.G.D. Geers, Deformation patterning driven by rate dependent non-convex strain gradient plasticity, J. Mech. Phys. Solids 59: 1-17, 2011 Dipartimento di Ingegneria, Università di Ferrara, Via Saragat 1, 44100 Ferrara, Italy E-mail address: dlpgpt@unife.it Dipartimento di Ingegneria Civile, Edile e Architettura, Università Politecnica delle Marche, Via Brecce Bianche 1, 60131 Ancona, Italy E-mail address: g.lancioni@univpm.it Istituto per le Applicazioni del Calcolo, CNR, Via dei Taurini 19, 00185 Roma, Italy E-mail address: r.march@iac.cnr.it

5

Abstract

Extracted data

5

Abstract

Extracted data

Related items

Related items