One-Dimensional Failure Modes for Bodies with Non-convex Plastic Energies

In this paper, a complete picture of the different plastic failure modes that can be predicted by the strain gradient plasticity model proposed in Del Piero et al. (J. Mech. Mater. Struct. 8:109–151, 2013) is drawn. The evolution problem of the elasto-plastic strain is formulated in Del Piero et al. (J. Mech. Mater. Struct. 8:109–151, 2013) as an incremental minimization problem acting on an energy functional which includes a local plastic term and a non-local gradient contribution. Here, an approximate analytical solution of the evolution problem is determined in the one-dimensional case of a tensile bar. Different solutions are found describing specific plastic strain processes, and correlations between the different evolution modes and the convexity/concavity properties of the plastic energy density are established. The variety of solutions demonstrates the large versatility of the model in describing many failure mechanisms, ranging from brittle to ductile. Indeed, for a convex plastic energy, the plastic strain diffuses in the body, while, for a concave plastic energy, it localizes in regions whose amplitude depends on the internal length parameter included into the non-local energy term, and, depending on the convexity properties of the first derivative of the plastic energy, the localization band expands or contracts. Complex failure processes combining different modes can be reproduced by assuming plastic energy functionals with specific convex and concave branches. The quasi-brittle failure of geomaterials in simple tension tests was reproduced by assuming a convex-concave plastic energy, and the accuracy of the analytical predictions was checked by comparing them with the numerical results of finite element simulations.


Introduction
In our previous paper [1] we presented two models, local and nonlocal, to describe fracture in the plastic regime. The local model is unable to describe softening, since an instability arises resulting in the extreme localization of the plastic deformation. This fact was noticed since the 1960's and since then it has been the object of a large literature because of some peculiar aspects such as mesh-dependence of the solutions in the finite-element computations and size effect. An accurate description and ample references are given by Jirásek and Bažant in their book [2]. Among the many localization limiters devised to contain the localization within finite regions of the body, very efficient came out to be the ones making use of the idea of a nonlocal continuum. In [2] three models of a nonlocal continuum are considered: integral models, and gradient models involving either the strain gradient or the gradient of an internal variable. Strain gradient theories, strictly related with the generalized continua [3], are compared in [4], paying attention on the models ability to localize strains. Analogous comparative studies are proposed in [5,6] where different plastic strain gradient theories are considered. In [7], different nonlocal models involving a damage internal variable are analysed, focusing on their ability to reproduce softening processes.
The model presented in [1] and resumed here is a model involving the gradient of a state variable, the plastic strain. It is developed in the simple one-dimensional setting where the key aspects of fracture phenomena can be easily extracted. The peculiarity of the formulation is that of being based on an energetic approach, so that the localization limiter has the form of a supplementary energy term. Our starting point is an energy of the form where and γ are the terms of the additive decomposition of the one-dimensional deformation displacement gradient u into an elastic and a plastic part, u being the axial displacement of the bar (0, l). The functions w and θ are energy densities depending on the two parts of the deformation, and the last term in the integral is the nonlocal term, the nonlocality being due to the fact that the derivative γ (x) somehow involves the values taken by γ in a neighborhood of x. This model has been shown to be insufficient to properly capture the final stage of material failure leading to rupture [8].
Nevertheless we keep the form (1) of the energy for simplicity, not after observing that more refined effects can be captured taking nonlocal energy terms of a more complicated form. Extension to multidimensional setting is also possible, and different gradient plasticity theories can be recovered. The theory of Aifantis [9,10] is obtained when local and nonlocal plastic energies are taken as function of the cumulated plastic strain, and the theory of Gurtin and Anand [11] turns out when a dependence on the plastic strain tensor is assumed (see [12] for a comparative study and [13] for a detailed description of the gradient plasticity theories). The method followed in [1] to determine the response to a given load is the incremental energy minimization. The load is taken to be an imposed stretch β, such that βl is the difference between the axial displacements at the bar's ends, and the model follows the evolution of the bar's deformation from the initial unloaded state up to rupture. The model succeeds in capturing the many phenomenological aspects of material failure reflected in the very different fracture modes of materials like, for example, steel and concrete, describing them with the aid of a small number of material parameters. The purpose of the present paper is to simplify our previous model, in order to provide a simple but unitary view of the very complicated aspects observed in fracture in the plastic regime. For this we introduce a basic simplification, consisting in replacing the elastic energy density w by its quadratic approximation at = 0 and the plastic energy density θ by its quadratic approximation at a representative valueγ of γ (x). Of course the approximation is exact when both w and θ are quadratic, and is accurate when either w and θ are close to quadratic and when both (x) and the difference (γ − γ (x)) are small.
Just as in our previous paper [1], the equilibrium configurations are identified with the pairs ( , γ ) for which the first variation of the energy is non-negative, and the stable equilibrium configurations for which the second variation is non-negative (when the first variation is zero). As a result of our simplifying assumption, the complete landscape of the stable equilibrium configurations for each value of β is provided. When β varies, the body evolves within these configurations, according to the usual scheme of plastic loading and elastic unloading.
An analytical solution of the evolution problem is found, under the hypothesis of homogeneous θ , in the case of a convex-concave plastic energy density. The solution allows to draw a complete picture of the different modes of plastic strain evolution that can be predicted, and to establish correlations between the shape of the plastic energy and the failure mode. Indeed, different phases are encountered in the strain process, characterized by specific solutions, ranging from the initial homogeneous strain growth in regime of strain hardening, to the final strain localization under stress softening. The progressive plastic strain localization is induced by the concave branch of the plastic energy, and it interprets the growing of a cohesive fracture. 1 Non-convex plastic energies were also used in [8,14,15] to describe different localization phenomena: McReynolds' plastic waves propagation and necking in tensile metallic bars [8], and plastic slip patterns in single crystal plasticity [14,15]. Non-convexity of surface energy is a crucial requirement in fracture mechanics as well. In [16], theories of cohesive fracture were deduced from non-convex elasticity by considering a double-well elastic potential typical of phase transition. Moving on the same research line, in [17,18], micro-and macro-cracking were reproduced by assigning specific non-convex shapes (convex-concave, bi-modal) to the surface energy included into the fracture model.
The analytical solution was specialized to the case of tensile geo-materials to reproduce the peculiar features observed in the rupture of quasi-brittle materials subjected to tensile loadings [19][20][21]. While quasi-brittle fracture is usually described by damage theories [22][23][24], here it is captured by properly tailoring the functional of the plastic energy density. The accuracy of the analytical solution is checked by comparing it with experimental results from [19] and finite element numerical simulations.
The paper is organized as follows. In Sect. 2, the one-dimensional model is formulated, and the equilibrium and the evolution equations are deduced as necessary conditions for certain energy minimum problems to be satisfied. In Sect. 3, the evolution problem is solved for convex-concave plastic energies. Different solutions are determined associated to the different stages of the plastic strain evolution. In Sect. 4, the accuracy of the analytical solution is tested by comparing it with the experimental data of tensile tests on concrete specimens. The analytical solution is also compared to the numerical solution resulting from a finite element simulation. Finally, concluding remarks are drawn in Sect. 5.

General Assumptions
For the strain energy of the bar we assume the form (1) with and γ the elastic and plastic parts of the deformation, w and θ functions of class C 2 , and α a given positive material constant. The pair ( , γ ) is a configuration of the bar. A deformation process is a one-parameter family t → ( t , γ t ) of configurations, and the parameter t is the time. We assume that w is strictly convex and with and that θ is strictly increasing and with We also assume that the function θ is dissipative, that is, that in every deformation process and at every point x of the bar the power θ (γ t (x))γ t (x) is non-negative. In view of the strict monotonicity of θ , this is equivalent to assume thaṫ in all deformation processes and at all x. Then in any deformation process γ t is a nondecreasing function of time.
No external forces are supposed to act on the bar, so that the strain energy is in fact the total energy. We also prescribe longitudinal displacements u t (0), u t (l) at the bar's ends, to which corresponds the total elongation We say that β t is a load, and that the map t → β t is a load process. We only consider the case of positive loads.

Equilibrium
A perturbation of a configuration ( , γ ) is a pair (δ , δγ ) in which δ and δγ are the initial time derivatives of some deformation process starting from ( , γ ). While δ is unrestricted, δγ is subject to the condition due to the dissipation inequality (4). We say that ( , γ ) is an equilibrium configuration if the first variation of the energy 2 is non-negative for all perturbations (δ , δγ ) which satisfy (6) and leave unchanged the length of the bar. In view of the decomposition (2) of u , from the boundary condition (5) rewritten in the form follows that the configurations ( , γ ) and ( + δ , γ + δγ ) correspond to the same length of the bar if In particular, for perturbations with δγ = 0 the first variation reduces to and δ is zero. In this case, a standard argument of the Calculus of Variations leads to the conclusion that the first variation is non-negative if and only if the map x → w ( (x)) is constant. Since by the assumed strict convexity of w the derivative w is strictly increasing, is the same at all x. By (8), the first term on the right-hand side of (7) becomes This eliminates the dependence of the first variation (7) on δ . Moreover, by the expression (9) of σ , the dependence on reduces to a dependence on β. Then for an equilibrium configuration, for every given β the first variation becomes a function of γ and δγ which, after an integration by parts, takes the form From the non-negativeness of δE and from the arbitrariness of δγ ≥ 0, it follows that These are necessary conditions for the non-negativeness of the first variation, to be satisfied at all interior points of the bar and at the endpoints, respectively. In turn, these conditions imply the non-negativeness of the first variation. Therefore, inequalities (10) are both necessary and sufficient for equilibrium. Moreover, in view of the dependence (9) of σ on β, inequalities (10) characterize the configurations equilibrated with the load β. Inequality (10) 1 is the yield condition which imposes an upper bound to the axial force σ . This condition is nonlocal, since the bound imposed at x involves the values of γ at the neighboring points through the derivative γ (x). The set of all points of the bar at which the yield condition is verified as an equality is the plastic zone. Then, after setting we have f (x) = 0 in the plastic zone and f (x) > 0 outside.

Incremental Minimization
For a given load process t → β t , consider the problem of finding a deformation process t →( t , γ t ) such that at each t the configuration ( t , γ t ) is a local minimizer of the energy in the class of all configurations equilibrated with the load β t . Since in an equilibrium configuration t is determined by β t and γ t and β t is given, the problem reduces to the determination of a minimizing process t →γ t . This problem can be solved by incremental energy minimization. Assuming that a local minimizer γ t at the time t is known, for a sufficiently small τ > 0 one looks for a local minimizer γ t+τ , close to γ t , in the class of all configurations equilibrated with the load β t+τ . Using the expansions β t+τ = β t + τβ t + o(τ ) and γ t+τ = γ t + τγ t + o(τ ), setting and neglecting the terms o(τ ), the energy of (β t+τ , γ t+τ ) takes the form For an arbitrarily small perturbation δγ ofγ t , we have and subtracting the previous equation we get For the first-order terms in δγ , integrating by parts we find with f t and σ t as in (11) and (12) 1 and, by time differentiation, By the arbitrariness of δγ , local necessary conditions for a minimum atγ t+τ are for all perturbations δγ such that, by (6), At points at whichγ t (x) is positive, both positive and negative values of δγ (x) are allowed. Then, by condition (16) ) must be zero. This alternative is expressed by the complementarity condition At the endpoint x = 0, assume that γ t (0) =γ t (0) = 0 at some t . Then δγ (0) ≥ 0 by (17) andγ t (0) ≤ 0 by (16) 2 . But in the expansioṅ γ t (x) must be non-negative by the dissipation inequality. Thenγ t (0) = 0 impliesγ t (0) ≥ 0. Therefore,γ t (0) must be zero. If γ t (0) = 0 andγ t (0) > 0, both positive and negative values of δγ (0) are allowed by (17), and (16) 2 requires thatγ t (0) = 0. Thus, in any case γ t (0) = 0 at some t impliesγ t (0) = 0. Then γ t+τ (0) = γ t (0) + τγ t (0) = 0. Then repeating the above with t replaced by t + τ we getγ t+τ (0) = 0, and so on.
A similar conclusion holds for the endpoint x = l. Therefore, if in a minimizing process t →γ t the Neumann boundary conditions γ t (0) = 0 and γ t (l) = 0 hold at some t , then the conditionsγ hold at all subsequent t .
In what follows, we make some simplifying assumptions: -the initial configuration is the natural configuration ( 0 , γ 0 ) = (0, 0), -the load is increasing at a constant rateβ, -the elastic energy density is quadratic From them it follows that, at all t , -the Neumann conditions (19) hold, -β t is a positive constantβ, -the axial force σ t is Under these assumptions and with the definitions (15), at all x in (0, l) the complementarity condition (18) The minimizing process t → γ t is determined solving the initial-boundary value problem consisting of this alternative, of the boundary conditions (19) and of the initial condition γ 0 = 0. The solution strongly depends on the shape of the function θ . Here we consider a special class of functions, which we call convex-concave, which seem to reproduce quite well the behavior of a large class of materials.

Solutions for Convex-Concave Energies
We say that a function γ →θ(γ ) is convex-concave if there is a γ p > 0 such that θ is convex for γ < γ p and concave for γ > γ p . Then θ (γ ) is positive for γ < γ p and negative for γ > γ p (see graphs of Fig. 1). Moreover, by the assumptions made in Sect. 2.1, γ →θ (γ ) decreases up to a negative minimum θ m at some γ m > γ p and then increases again, tending to zero for γ → +∞. For energies θ of this form, the minimizing process t →γ t comes out to be made of a sequence of different regimes.

The Initial Elastic Regime
At the initial time t = 0 we have γ 0 = β 0 = 0 by assumption, σ 0 = 0 by (20), and f 0 (x) = θ (0) by (11). The support ofγ 0 consists of intervals (a, b), in each of which equation (21) 1 has the form Moreover, at the endpoints we haveγ and since both K and θ (0) are positive, the inequality follows. Since K, θ (0) andβ are given, this inequality can be violated choosing a sufficiently small τ . This renders the first alternative in (21) unavailable. By consequence,γ 0 (x) must be zero all over the bar. Then γ t (x) = γ 0 (x) + tγ 0 (x) = 0 at the subsequent t , at which inequality (23) Since plastic strains can develop from now on, the instant t o identifies the onset of the plastic regime. Subscript letter o stands for onset.

The Uniform Hardening Regime
At t = t o we have f o (x) = 0, and in each of the intervals (a i , b i ) which form the support oḟ γ o equation (21) 1 becomes Integrating over (a i , b i ) under the boundary conditions (22) and summing over all intervals (a i , b i ), Subtracting to equation (25), we get whereγ  interval (a, b), and it coincides with (0, l). Then l * o = l, and from (26) we haveγ At the instants which follow t o the bar keeps a uniform regime, with x →γ t (x) constant and l * t = l. Moreover, by the complementarity condition (18), That is, the plastic zone, which for t < t o was empty, now occupies the whole bar. Then the regime starting at t o is a uniform plastic regime. For t > t o , equation (29) holds with θ (0) replaced by θ (γ t ) Moreover, the relationσ shows that the slopeσ t /β is positive as long as θ (γ t ) is positive. That is, the regime starting at t o is a uniform hardening regime.

The Uniform Softening Regime
The hardening regime ends at the time t p at which θ (γ t ) becomes zero. At this instant, γ = γ p , f p is zero, and the axial force attains the peak value σ p = θ (γ p ). This explains the choice of the subscript letter p (peak) labelling this time instant. Equation (21) 1 becomes Due to the boundary conditions (19), the integration over (0, l) yieldsγ p =β. Thenγ p (x) = 0, and from (15) we getḟ p (x) = 0, so that f t (x) = 0 at the instants t following t p . Therefore, the regime starting at t = t p is again a uniform plastic regime. At the instants t following t p , equation (21) 1 holds with γ t (x) = γ t and f t (x) = 0 at all x in (0, l). Then integrating over (0, l) we geṫ andσ t is as in (31). But, since now θ (γ t ) is negative, the slopeσ t /β of the response curve becomes negative, that is, a uniform softening regime takes place. Moreover, forγ ( With the boundary conditions (19) this differential equation has the null solution, except at κ t l = π , where it admits the trigonometric solutioṅ with A an arbitrary constant. The uniform solutionγ t =γ t , which is the only possible one for κ t l < π, for κ t l > π becomes unstable. Indeed, by (15), (31) and (32), But f t (x) = 0 andḟ t = 0 implyḟ t (x) = 0. Then in (14) the zero-and first-order terms in δγ vanish, and the condition of non-negativeness of the right-hand side of (14) reduces to In particular, for the perturbation δγ (x) = ε cos πx/l, with ε a small constant, we have a condition violated for −θ (γ t ) > απ 2 /l 2 . Equation (32) satisfies the dissipation inequalityγ t ≥ 0 as long as the denominator is positive. In the uniform softening regime, θ (γ t ) decreases from θ (γ p ) = 0 to the minimum θ m , and then increases tending to zero for γ t → +∞. If K + θ m is negative, there is a θ r > θ m for which K + θ r = 0. Then K + θ (γ t ) is negative for all γ t such that θ (γ t ) < θ r , and for all such γ t the dissipation inequality is violated. Moreover, by (31), when θ (γ t ) tends to θ r the slope of the response curve tends to −∞. The impossibility of reaching the γ t larger than γ r and the very quick reduction of the axial force when γ t approaches γ r are interpreted as signs of the outcoming brittle rupture of the bar.
In conclusion, the uniform softening regime continues up to γ t = +∞, except in two cases. The first is the case θ m < −K, in which brittle fracture takes place. The second is the case θ m < −απ 2 /l 2 , in which the constant solution becomes unstable. This second case is studied in the next section.

The Localized Softening Regime
If απ 2 /l 2 is smaller than −θ r and −θ m , the uniform softening regime ends at the time t c at which At t c there is a bifurcation of equilibrium, and the unstable uniform solution is replaced by a localized solution, in which the support ofγ t is strictly contained in (0, l). Since a phase of contraction of the plastic region initiates, the subscript letter c is used to define this time instant. The present study is restricted to the case in which the support ofγ t is a single interval (a, b), with a = 0 and b = l * t < l to be determined for each t > t c . At t = t c the localized regime starts with the solution (33) This is a second-order differential equation with a variable coefficient θ (γ t (x)). Leaving apart any attempt for a closed-form solution, we prefer to focus on the general behavior of the solutions. For this we make a drastic simplification: at each t , we replace the function x →θ (γ t (x)) by a constant, which we denote by θ t . 4 With this simplification and using the boundary conditionsγ by integration of (35) over (0, l * t ) we get and, therefore,γ Substituting into (35) and settingγ the differential problem reduces to Its closed-form solution iṡ Here, again, the boundary conditionγ t (l *

The Expansion of the Plastic Zone
In the expansive regime, at each time t the plastic zone (0, l t ) is known from the preceding step, and the support (0, l * t ) ofγ t is unknown but larger than (0, l t ). This situation requires a full application of the incremental minimization technique developed in Sect. 2.3.
We then select a sufficiently small τ > 0, and minimize the energy (13) under the approximation θ (γ t (x)) = θ t . Then the alternative (21) becomes where f t is a non-negative function, known from the previous minimizing steps and null over (0, l t ). Then integrating over (0, l * t ) we get Dividing by l * t and subtracting to (38) 1 we obtain the differential equation withγ defined as in (28). Moreover, from (22) and the continuity ofγ t andγ t at x = l * t we haveγ It is convenient to splitγ t into the sum of two functionsγ 1 ,γ 2 with support on (0, l * t ) and (l t , l * t ), respectively, satisfying the equations and the conditionṡ Forγ 1 withγ 1 (0) = 0, the solution iṡ with A an arbitrary constant. Forγ 2 we look at a solution of the forṁ with g t a function to be determined. Sincė the differential equation (40) 2 is satisfied by The conditionγ 2 (0) = 0 is identically satisfied, and fromγ 2 (0) = 0 we have g t (0) = 0. Therefore, The conditions (41) at x = l * t yield and from (44) and (45) we finally get follow t m we have h t > 0 and, for coherence with the sign of the term on the left-hand side of (50), the parenthesis must be positive, so that we have tan κ t l t < 0. Since κ t l t < π, it then follows that κ t decreases, l t increases and the product κ t l t decreases.
To conclude, a typical stress-strain curve that characterizes the evolution process described in this section is drawn in Fig. 2. It is constituted by five different branches, which correspond to the different phases of the stretching process. Red dots indicate the points of transition between phases associated to the instants (t o , t p , t c , t m ).

Quasi-Brittle Failure in Geomaterials
In this section, the analytical results of Sect. 3 are implemented in a numerical algorithm to reproduce the failure process observed in tensile tests on geomaterials. The accuracy of the analytical predictions are checked by comparing them with the experimental stressstrain curves of simple tension tests on lightweight concretes in [19], and with the numerical results of finite element simulations.

Semi-Analytical Algorithm
According to the incremental minimization procedure of Sect. 2.3, we discretize time into intervals of length τ , and, within each time step t → t + τ , we estimate the plastic strain γ t+τ supposing that γ t is known.
Given γ t , we estimate f t (x) by (11). If γ t is not homogeneous, θ t (γ t (x)) is approximated by the constant θ t = θ (γ t ) whereγ t is a representative value of γ t (x). We consider the following three alternative expressions forγ t : with (0, l t ) the plastic zone. The plastic strain rateγ t is determined by following the scheme sketched below.
• Elseif θ t < 0, three different solutions are possible, depending on the values of l, l t , and -If l * t ≥ l,γ t is homogeneous, and it is given by (30). -Elseif l * t ≤ l t < l,γ t localizes, and the plastic zone contracts;γ t has the expression (37), withγ t as in (36) 1 .

-End • End
Plastic strain and axial force at t + τ are This computation is repeated at each time step. The geometrical and constitutive data that must be assigned to perform the computation are l,β, τ , K, θ = θ(γ ) and α, and they are assigned in the next section.

Parameters Calibration
In this section, the constitutive parameters of the model are calibrated to reproduce the experimental results of [19], where simple tension tests are performed on prismatic specimens made of lightweight concrete.
We assign K = 18 · 10 3 MPa, corresponding to the Young's modulus of lightweight concrete [19]. For the plastic energy density, we take the piece-wise cubic polynomial with θ i = θ(γ i ), and i = 1, . . . , n. We consider n = 4. The four branches are defined in intervals between the points (0, γ c , γ m , γ u , ∞), as shown by the graphs of Fig. 3. The strain γ c is that at which the localization regime initiates (see Sect. 3.4), and the strain γ m is that at which θ attains the minimum value.
Some experimental data are needed to calibrate the parameters in (52). They are reported in Fig. 4, where the envelope of the experimental curves 7 is plotted. They are the stress σ o = 2.50 MPa at the end of the linear elastic phase, and the pair u p = 8 · 10 −3 mm and σ p = 3.25 MPa at peak stress. In Fig. 4, the horizontal axis measures the relative displacement u meas Fig. 3 Graphs of the piecewise cubic function θ = θ(γ ), with its first and second derivative Fig. 4 Experimental stress-displacement curve from [19]. Displacement u meas refers to the measuring length l meas = 35 mm within a measure length of 35 mm across the crack, and thus the corresponding stretching is β = u meas /35. With these data, the parameters γ i , θ i , θ i and θ i are determined as follows.
At γ = 0, θ 0 = σ o from (24) and θ 0 = 0. The values of θ 0 , γ c and θ c are where γ p = u p /35 − σ p /K is the plastic strain corresponding to the peak stress. The above expressions are determined by solving the equations The third equation is (34), defining the plastic strain at the onset of localization. Parameters γ m and θ m are related to the instant t m of transition from contraction to expansion of the plastic zone, and it corresponds to the point of the stress-strain curve separating the concave and the convex softening branches. Accordingly, the pair (γ m , θ m ) controls the shape of the stress-strain softening tail. We assign the values γ m = 2 · 10 −4 and θ m = −5 · 10 3 . The strain γ u satisfies the equation θ (γ u ) = 0, and its value is γ u = γ m − 2θ m /(θ m + θ u ). Finally, with i = c, m, u, which guarantee the continuity of θ and θ at γ i . A further datum needed to calibrate the non-local coefficient α is the size of the band where fracture localizes, which is about 2.7 times the maximum aggregate size, according to the estimates proposed in [25]. Since the concrete used for experiments in [19] has the maximum aggregate size equal to 8 mm, we assume the fracture bandwidth l f ract = 22 mm, and we suppose that it represents the minimum size of the localization zone. Since the smallest zone where plastic strains localize is attained at t m , and, its half-bandwidth is l * m = π −α/θ m from (37), we take l * m = 11 mm. By inverting the above relation, we assign α = −l * 2 m θ m /π 2 .
Finally, we assume that the bar length is l = 17.5 mm, corresponding to one half of the length of the specimen zone where strains are measured in experiments.

Results
The stress-strain curves obtained by implementing the semi-analytical algorithm of Sect. 4.1 for the three possible choices (51) ofγ t are plotted in Fig. 5 (red and blue curves), and they are compared with the experimental envelope, and with the numerical curve of a finite elements simulation (black curve). The FEM code used for the simulation is based on a staggered scheme which, at each time increment, iterates the two-step minimization of the energy functional (3), first with respect to u, keeping γ constant, and then with respect to γ , for constant u. Since minimization with respect to γ is a constrained non-linear problem, it is solved by a sequential quadratic programming scheme (see [8] for details on the finite element implementation).
Profiles of γ given by the analytical solution and those of the numerical simulation are plotted in Fig. 6 at different instant of the evolution process. The analytical solutions reproduce all the phases of the evolution of γ , which are described in the following. The differences between analytical and numerical results are also commented. Initially, γ grows homogeneously. Strain localization initiates when the representative strainγ t reaches the valueγ t = γ t = γ c . The phase of localization on smaller and smaller zones is related to the concave decreasing branch of the stress-strain curves of Fig. 5, and it terminates whenγ t = γ m . Analytical and numerical profiles ofγ at the initial and final instants of the localization process are compared in Fig. 7(a). We notice that the minimum length of the localization zone is about 11 mm (profiles for β = 2.8) in accordance with the length l * m = 11 mm assigned in Sect. 4.2. Then, the evolution proceeds with a progressive enlargement of the plastic zone, and the related branch of the stress-strain curve is the final convex tail. Profiles at two different instants of the enlargement phase are plotted in Fig.  7(b) and (c). We observe that the analytical profiles obtained with the assumption (51) 1 overestimate the size of the localization zone, which, for β > 5.7, cover the whole bar with homogeneousγ . On the contrary, the size of the localization zone is underestimated when the assumption (51) 2 is made.
The differences between the analytical and numerical profiles of Figs. 6 and 7 translate into differences among the stress-strain curves of Fig. 5. The analytical stress-strain curves exhibit an evident dependence on the choice of the representative valueγ t used into the algorithm. By comparing the blue and red curves, we notice that the assumptions (51) 1 and (51) 2 plays opposite effects. For the assumption (51) 1 (blue curve), the slope of the softening branch is overestimated in the initial concave part and it is underestimated in the next convex tail, while the opposite trend is observed if the assumption (51) 2 is made (red curve). We explain these differences by noticing that, from (36),σ is proportional to θ t . Then, θ (max{γ t }) < θ (lγ t /l t ) < 0, when bothγ t = max{γ t } andγ t = lγ t /l t belong to the set (γ c , γ m ) (phase of localization related to the concave stress-strain branch), because max{γ t } > lγ t /l t and the curve θ (γ ) is decreasing (see Fig. 3). As a result, the softening concave branch associated toγ t = max{γ t } (blue curve) is more sloped than that associated toγ t = lγ t /l t (red curve). On the contrary, 0 > θ (max{γ t }) > θ (lγ t /l t ), if max{γ t } and lγ t /l t are larger than γ m (phase of enlargement, characterized by the convex stress-strain branch), since θ (γ ) is increasing, and, thus, the blue softening convex branch related tõ γ t = max{γ t } has a lower slope than the red one related toγ t = lγ t /l t .

Conclusions
The incremental minimization problem proposed in [1] for the evolution of strain in nonlocal elasto-plastic bodies was analytically solved in the simple case of a bar subjected to an increasing stretching. A convex-concave plastic energy density θ was assigned to the material, in addition to a non-local contribution depending on the gradient of the plastic strain γ , and different solutions were found, describing specific evolution modes, and characterizing the different phases of the deformation process. Correlations were established between solutions and convexity/concavity properties of θ and of its first derivative θ .
A convex or a concave function θ makes the evolution stress-hardening or strainsoftening, respectively. In hardening regime, γ evolves homogeneously, spreading in the whole bar. In softening regime, the evolution of γ is still homogeneous if θ is moderately concave, but, for a sufficiently large concavity (i.e., large values of |θ |), γ localizes in a bar portion whose width is proportional to |θ |, and the localization zone shrinks or expands depending on the concavity/convexity of θ . Brittle failure is also captured by assuming a function θ whose concavity is large enough in comparison to the elastic modulus of the material.
Since the analytical solutions are determined under the hypothesis of homogeneous θ at each time step of the incremental evolution, the effects of this assumption on the solution accuracy was tested by comparing the analytical predictions with those of finite-element simulations. Criteria for parameters calibration were proposed, then the failure process experienced by concrete specimens in simple tension tests was reproduced. The main stages of the softening process observed in experiments were captured by the analytical solutions, and the influence of the different choices of the constant θ at each time step was analysed. By comparing the response curves of Fig. 5, it was found that the analytical solution also gives quantitatively accurate estimates, provided that the representative strainγ is properly chosen.
To conclude, the proposed study has pointed out the versatility of the gradient plasticity model [1] in describing a large variety of different plastic evolution modes, and it has led the way to reproduce complex failure processes by properly assigning the shape of the plastic energy functional. A further step of the research will be to extend the formulation to the multi-dimensional context.