Bounds on the elastic threshold for problems of dissipative strain-gradient plasticity
BBounds on the elastic threshold for problems ofdissipative strain-gradient plasticity
B. D. Reddy ∗ and S. Sysala † April 27, 2020
Abstract
This work is concerned with the purely dissipative version of a well-established model ofrate-independent strain-gradient plasticity. In the conventional theory of plasticity theapproach to determining plastic flow is local, and based on the stress distribution in thebody. For the dissipative problem of strain-gradient plasticity such an approach is notvalid as the yield function depends on microstresses that are not known in the elasticregion. Instead, yield and plastic flow must be considered at the global level. This workaddresses the problem of determining the elastic threshold by formulating primal anddual versions of the global problem and, motivated by techniques used in limit analysisfor perfect plasticity, establishing conditions for lower and upper bounds to the threshold.The general approach is applied to two examples: of a plate under plane stress, andsubjected to a prescribed displacement; and of a bar subjected to torsion.
Keywords: dissipative strain-gradient plasticity, elastic threshold, duality, lower andupper bounds, penalization method, finite element method ∗ Department of Mathematics and Applied Mathematics, University of Cape Town, 7701 Rondebosch,South Africa. Email: [email protected]. BDR acknowledges support for his work from the NationalResearch Foundation, through the South African Chair in Computational Mechanics, SARChI Grant47584. † Institute of Geonics of the Czech Academy of Sciences, Department of Applied Mathematics andComputer Science & Department IT4Innovations, 708 00 Ostrava-Poruba, Czech Republic. E-mail:[email protected]. SS acknowledges support for his work from the Czech Science Foundation(GA ˇCR) through project No. 19-11441S. a r X i v : . [ phy s i c s . c l a ss - ph ] A p r Introduction
There has been sustained interest over the last few decades in the development and assess-ment of models of elastoplastic behaviour that take account of size-dependent responsesat the microscale. A central and popular direction has been the development of models inwhich gradients of the plastic strain enter the formulation, as a means of accounting forthe relationship between size-dependence and the development of geometrically necessarydislocations in the context of non-homogeneous deformations. Some representative ex-amples of such strain-gradient theories include [9, 14, 15] (see also the references in theseworks).The works [27, 25, 26] have focused on analyses of well-posedness for particular formula-tions of strain-gradient plasticity. These analyses have given attention typically to whatare referred to as energetic as well as dissipative models. Energetic formulations refer tomodels of strain-gradient plasticity in which the gradient terms arise from a recoverablefree energy, and appear as a back-stress in the yield function. Dissipative formulations,on the other hand, are characterized by an extension of the conventional flow relationin which plastic strain rates as well as their gradients appear in the normality relationof associative plasticity. The two models lead respectively to size-dependent behaviourin the form of an increase in strain-hardening with length scale, and in the second casestrengthening, that is, an increase in initial yield with increase in length scale. The works[4, 5] provide further examples of theoretical investigations of both kinds of behaviour.Some features of rate-independent versions of dissipative problems have attracted par-ticular attention. First, the generalization of the classical flow relation with a normalityrelation leads to a yield function expressed in terms of generalized microstresses which,because they are not known in the elastic region, cannot predict yield in the conventionalway (see also [12]). On the other hand, it has been shown [27, 25] that the appropriateway to interpret yield and plastic flow lies through a global form of this relation, in whichit then becomes possible to formulate the flow relation in terms of the Cauchy stress.This is most readily carried out for the kinematic form, using the dissipation function.In [2] the matter of determining a dual form of such a global flow relation is investigatedfurther, where it is shown that a closed-form expression for such a relation, in terms of aglobal yield function, is elusive.A second feature of dissipative problems relates to the response to non-proportional load-ing in the plastic range. This leads to an elastic gap: that is, an elastic response followingthe application of such loading. The phenomenon was first presented in [11], and has2een further studied in various works, for example [1, 13, 21, 22, 23].The objective of this work is to address, in the context of the dissipative problem ofstrain-gradient plasticity, the problem of determining the elastic threshold: that is, inrelation to a scalar load factor, the stage at which incipient plastic behaviour takes place.The basis for the study is the dissipative version of the rate-independent theory presentedand analyzed in [2, 27, 25].The approach taken in this work is to draw on the methods of limit analysis, a well-established area of investigation in conventional perfect plasticity. The link derives fromthe formulation, in limit analysis, of optimization problems that yield lower and upperbounds to collapse states [29, 3, 18, 28]. This correspondence is noted also in [12, Section7], and is explored in detail in [24], for a model in which size-dependence is through thegradient of a scalar function of plastic strain. In the current work it will be seen that theelastic threshold for problems of dissipative gradient plasticity may likewise be formulatedas lower- and upper bound problems.The rest of this work is organized as follows. In Section 2 the rate-independent problem ofdissipative strain-gradient plasticity is presented, using the notions of a yield function andassociative flow relation. The global form of the flow relation is derived in its kinematicform, that is, using the dissipation function or support function associated with the yieldfunction. Section 3 makes precise the notion of plastically admissible stress fields in thecontext of strain-gradient plasticity. Inspired by classical limit load analysis for perfectplasticity, we obtain necessary and sufficient conditions for stress fields to be plasticallyadmissible. This paves the way, in Section 4, to derive expressions for lower and upperbounds to the elastic threshold, in terms of a parameter characterizing monotonic loading.Section 5 introduces approximations in the form of penalizations of the highly nonlinearproblems for the bounds. These theoretical results are applied to two example problems:in Section 6, that of a square domain in plane stress and subjected to a prescribed dis-placement; and in Section 7, the problem of torsion of a circular rod. For both problemsapproximations to the bounds are sought first by assuming simple forms for the minimiz-ers or maximizers; and secondly, the problems are solved numerically, using finite elementapproximations in combination with penalized forms of the problems.3
The problem setting
Vector- and 2nd-order tensor- or matrix-valued functions will be written in lower caseboldface form. The scalar product of two tensors or matrices σ and τ will be denoted by σ : τ , and in component form, relative to an orthonormal basis, by σ ij τ ij , the summationconvention on repeated indices being invoked. Furthermore, we use upper case boldfaceletters to denote third-order tensors. For two such quantities Π and T , the inner productis denoted by Π ◦ T , or in component form Π ijk T ijk .Consider a homogeneous elastic-plastic body that occupies a bounded domain Ω ∈ R with boundary ∂ Ω. For a prescribed body force b ∈ [ L (Ω)] the equilibrium equationreads − div σ = b , (2.1)where σ is the symmetric Cauchy stress tensor. The infinitesimal strain ε is defined as afunction of the displacement u by ε := ε ( u ) = ( ∇ u + ( ∇ u ) (cid:62) ) . (2.2)The strain ε is decomposed into elastic and plastic components ε e and ε p , respectively,according to ε = ε e + ε p . (2.3)The material is assumed to be isotropic, with the elastic relation given by σ = C ε e := λ (tr ε e ) I + 2 µ ε e . (2.4)Here C is the fourth-order elasticity tensor, and λ and µ are the two Lam´e coefficients.This investigation is based on a model of strain-gradient plasticity [25, 27] in which theclassical notions of a convex yield surface and normality flow relation are extended to thegradient case. This model takes as a point of departure the important work [15], whichassumes viscoplastic behaviour without a yield surface. The theory makes use of second-and third-order microstresses π and Π respectively. The quantity π is symmetric anddeviatoric, and Π is symmetric and deviatoric in its first two indices, in the sense thatΠ ijk = Π jik , Π ppk = 0. We will also require the deviatoric stress σ D := σ − (tr σ ) I . Inaddition to the equation of macroscopic equilibrium the Cauchy stress and microstressessatisfy the microforce balance equation σ D = π − div Π or, in index form σ Dij = π ij − Π ijk,k . (2.5)4he local form of the flow relation is posed in terms of a yield function f (cid:96) ( π , Π ), dependenton a length parameter (cid:96) . The flow relation is then an extension of the classical normalityrelation and is given by ( ˙ ε p , ∇ ˙ ε p ) = γ (cid:18) ∂f (cid:96) ∂ π , ∂f (cid:96) ∂ Π (cid:19) ,γ ≥ , f (cid:96) ( π , Π ) ≤ , γf (cid:96) ( π , Π ) = 0 . (2.6)We will work with f (cid:96) ( π , Π ) = (cid:112) | π | + (cid:96) − | Π | − Y, (2.7)where Y > E , defined by E = { ( π , Π ) | f (cid:96) ( π , Π ) ≤ } . (2.8)The flow relation may be recast with the microstresses as dependent variables by intro-ducing the support function associated with E : in the present context this is referred toas the local dissipation function, denoted by D (cid:96) , and given by D (cid:96) ( q , ∇ q ) = sup ( π , Π ) f (cid:96) ( π , Π ) ≤ [ π : q + Π ◦ ∇ q ] . (2.9)For f (cid:96) given by (2.7), the local dissipation function takes the form D (cid:96) ( q , ∇ q ) = Y (cid:112) | q | + (cid:96) |∇ q | in Ω . (2.10)Then from a basic result in convex analysis, the flow relation (2.6) is equivalent to theinclusion ( π , Π ) ∈ ∂D (cid:96) ( ˙ ε p , ∇ ε p ); (2.11)here ∂D (cid:96) is the subdifferential of D (cid:96) , so that (2.11) reads D (cid:96) ( q , ∇ q ) ≥ D (cid:96) ( ˙ ε p , ∇ ε p ) + π : ( q − ˙ ε p ) + Π ◦ ∇ ( q − ˙ ε p ) ∀ q . (2.12) A note on alternative forms for f (cid:96) and D (cid:96) . A more general family of local dissipationfunctions may be defined by D (cid:96),r ( q , ∇ q ) := Y [ | q | r + ( (cid:96) |∇ q | ) r ] /r (2.13)5or real r satisfying 1 ≤ r < ∞ . The case r = ∞ corresponds to the function D (cid:96), ∞ ( q , ∇ q ) := Y max {| q | , (cid:96) |∇ q | ] . (2.14)The yield function corresponding to D (cid:96),r , 1 < r ≤ ∞ reads f (cid:96),r (cid:48) ( π , Π ) := (cid:104) | π | r (cid:48) + ( (cid:96) − | Π | ) r (cid:48) (cid:105) /r (cid:48) − Y, r + 1 r (cid:48) = 1 . (2.15)The case r = 1 is also of practical significance [8]. The yield function f (cid:96), ∞ correspondingto D (cid:96), is shown in [25] (see also [16, Section 4.3.2] to be given by f (cid:96), ∞ ( π , Π ) = max {| π | , (cid:96) − | Π |} − Y . (2.16)We will however focus in this work on the case r = 2, for which all significant features arelikely to be present. Initial and boundary conditions.
The initial conditions of the problem are u = , ε p = in Ω × { } . (2.17)For the boundary conditions associated with the equilibrium equation the boundary ∂ Ωis partitioned into complementary parts ∂ Ω u and ∂ Ω t on which macroscopic Dirichlet andNeumann boundary conditions respectively are prescribed. These are u = ¯ u on ∂ Ω u × (0 , T ) , σn = ¯ t on ∂ Ω t × (0 , T ) , (2.18)where ¯ u and ¯ t are respectively a prescribed displacement and traction, and n denotes theoutward unit normal to the boundary ∂ Ω.Boundary conditions associated with the microforce balance equation are imposed oncomplementary parts ∂ Ω H and ∂ Ω F . The conditions on these two parts are referred torespectively as microhard and microfree boundary conditions, and are given by ε p = on ∂ Ω H × (0 , T ) , Π n = on ∂ Ω F × (0 , T ) . (2.19)The initial-boundary value problem for gradient plasticity is given by equations (2.1),(2.3), (2.2), (2.4), (2.5), (2.6) together with the initial and boundary conditions (2.17),(2.18), (2.19). 6 emark. The model presented here is a dissipative one, in the sense that recoverableenergetic contributions to terms involving the plastic strain gradient are omitted. Atreatment of the full problem may be found in [25, 27, 16]. The structure of the dissipativemodel presents difficulties in its local form, and is approached as a global problem.
The global form of the flow relation.
It is clear from (2.6) or (2.12) that, in contrastto the case of conventional plasticity, it is in fact not possible to use the local yieldcondition to determine pointwise whether the elastic limit has been reached. This is thecase as π and Π are not known a priori (cf. the case of a rigid-plastic body, for which itis not possible to determine the stresses in the rigid region, and hence to determine whenyield occurs locally).This dilemma may be resolved by approaching plastic flow as a global notion. For thispurpose we construct a weak formulation of the microforce balance equation. We definethe following space of tensor fields: W = { q : Ω → R × sym | q ij = q ji , q ii = 0 , q ij ∈ L (Ω) , q ij,k ∈ L (Ω) , q = on ∂ Ω H } . The L -space of all pairs ( π , Π ) will be denoted by X . We take the inner product of(2.5) with q − ˙ ε p ∈ W , integrate, integrate the term involving Π by parts, and use theboundary condition (2.19) , to obtain (cid:90) Ω [ π : ( q − ˙ ε p ) + Π ◦ ∇ ( q − ˙ ε p )] dx = (cid:90) Ω σ D : ( q − ˙ ε p ) dx ∀ q ∈ W. (2.20)Next, integrate (2.12) and use (2.20) to eliminate the terms involving π and Π from theresulting global inequality. The result is the global flow relation (cid:90) Ω D (cid:96) ( q , ∇ q ) dx ≥ (cid:90) Ω D (cid:96) ( ˙ ε p , ∇ ε p ) dx + (cid:90) Ω σ : ( q − ˙ ε p ) dx ∀ q ∈ W . (2.21)
Remark.
The variational inequality (2.21) together with the weak form of the equilibriumequation has a unique solution in appropriate Sobolev spaces provided that there is somehardening present [25]. For example, kinematic hardening would result in the replacementof the second term on the righthand side of (2.21) by σ − κ ε p , where κ ∈ L ∞ (Ω) is requiredto satisfy κ ( x ) ≥ ¯ κ >
0. We will omit the hardening term in what follows, as it will notbe central to the focus in what follows, viz. a global form of the flow relation.7
Plastically admissible stresses
Within this section, we fix t ∈ (0 , T ). From the classical problem (2.1)–(2.6), (2.17)–(2.19) we say that the stress field σ is plastically admissible if there exists a pair ( π , Π )satisfying π − div Π = σ D in Ω , Π n = on ∂ Ω F ,f (cid:96) ( π , Π ) ≤ . (3.1)We denote the set of all plastically admissible stresses by E sglob . As indicated earlier, theconditions in (3.1) cannot be verified locally by a direct approach, and so an alternativeapproach has to be adopted to determine whether σ ∈ E sglob .Let Π be such that (3.1) holds and div Π exists. Inserting (3.1) into (3.1) , we arriveat the inequality f (cid:96) ( σ D + div Π , Π ) ≤ , (3.2)which is a sufficient condition to be σ ∈ E sglob . In particular, if Π = then (3.2) reducesto f (cid:96) ( σ D , ) ≤ . (3.3)Thus local yield is a sufficient condition for σ to be plastically admissible.To find necessary conditions for σ ∈ E sglob , we develop an alternative approach to thedefinition of plastically admissible stresses using duality theory. The definition based on(3.1) may be interpreted as the static approach while the dual definition will be kinematicin nature. The following derivation is inspired by the techniques of limit analysis in perfectplasticity [29, 3, 19].First, define the auxiliary function ˜ ϕ glob ( σ ) := supremum over all λ ≥ for which thereexists ( π , Π ) such that π − div Π = λ σ D in Ω , (3.4a) Π n = on ∂ Ω F , (3.4b) f (cid:96) ( π , Π ) ≤ . (3.4c)Let E kglob := { σ | ˜ ϕ glob ( σ ) ≥ } . It is readily seen that the inequality ˜ ϕ glob ( σ ) ≥ σ ∈ E sglob and that E sglob ⊂ E kglob . In addition, if ˜ ϕ glob ( σ ) > λ < ˜ ϕ glob ( σ ), i.e., for λ = 1 implying σ ∈ E sglob . On the otherhand, if ˜ ϕ glob ( σ ) = 1 then it is not clear whether σ belongs to E sglob or not. This task is8eyond the scope of this work.From (2.20) the weak form of (3.4a)–(3.4b) reads (cid:90) Ω [ π : q + Π ◦ ∇ q ] dx = λ (cid:90) Ω σ D : q dx ∀ q ∈ W. (3.5)This condition holds if and only if λ = inf q ∈ W (cid:82) Ω σ D : q dx =1 (cid:90) Ω [ π : q + Π ◦ ∇ q ] dx. (3.6)Hence, we arrive at˜ ϕ glob ( σ ) = sup ( π , Π ) ∈ Xf (cid:96) ( π , Π ) ≤ inf q ∈ W (cid:82) Ω σ D : q dx =1 (cid:90) Ω [ π : q + Π ◦ ∇ q ] dx. (3.7)Next, we swop the order of inf and sup in (3.7). In general, we have sup inf ≤ inf sup.Sufficient conditions for sup inf = inf sup are introduced e.g. in [7, Chapter VI]. Forthe present problem, from [7, Proposition VI.2.3 and Remark VI.2.3] it follows that theequality sup inf = inf sup holds if the set { ( π , Π ) ∈ X | f (cid:96) ( π , Π ) ≤ } is bounded.This is the case for the function f (cid:96) introduced in (2.7) and considered in this work.Therefore, we can write˜ ϕ glob ( σ ) = inf q ∈ W (cid:82) Ω σ D : q dx =1 sup ( π , Π ) ∈ Xf (cid:96) ( π , Π ) ≤ (cid:90) Ω [ π : q + Π ◦ ∇ q ] dx. (3.8)Using [7, Proposition IX.2.1], it is possible also to swop the supremum and the integralin (3.8):sup ( π , Π ) ∈ Xf (cid:96) ( π , Π ) ≤ (cid:90) Ω [ π : q + Π ◦ ∇ q ] dx = (cid:90) Ω sup ( π , Π ) ∈ Xf (cid:96) ( π , Π ) ≤ [ π : q + Π ◦ ∇ q ] dx = (cid:90) Ω D (cid:96) ( q , ∇ q ) dx, (3.9)where D (cid:96) is the dissipation introduced in (2.9) or (2.10). Hence,˜ ϕ glob ( σ ) = inf q ∈ W (cid:82) Ω σ D : q dx =1 (cid:90) Ω D (cid:96) ( q , ∇ q ) dx. (3.10)9sing (3.10), the (kinematic) global yield surface in gradient plasticity is defined as follows: E kglob = { σ | ˜ ϕ glob ( σ ) ≥ } = (cid:26) σ | (cid:90) Ω σ D : q dx ≤ (cid:90) Ω D (cid:96) ( q , ∇ q ) dx ∀ q ∈ W (cid:27) . (3.11)If σ ∈ E sglob then the inequality in (3.11) must be satisfied for any q ∈ W . That is, thisinequality is a necessary condition for σ to be plastically admissible. The relation (3.11)has also been obtained in [22, Secton 6.4.1], using arguments based on polar conjugatesof convex positively homogeneous functionals. In this section we assume that the data (that is, ¯ u , ¯ t and b in Section 2) depend linearlyon the time variable t ; that is,¯ u = t ˆ u , ¯ t = t ˆ t , b = t ˆ b , t ∈ (0 , T ) , where ˆ u , ˆ t and ˆ b are given functions defined in Ω. We consider the evolution of thesolution to the problem (2.1)–(2.6), (2.17)–(2.19) over the time interval (0 , T ). Clearly,for sufficiently small t > u = t ˆ u e , ε = t ˆ ε e , σ = t ˆ σ e , ε p = , γ = 0 , π = t ˆ π e , Π = t ˆ Π e . (4.1)The solution components ˆ π e and ˆ Π e need not be uniquely defined as we shall see inSection 6.The elastic threshold t ∗ is the maximal value of t for which the solution to (2.1)–(2.6),(2.17)–(2.19) satisfies (4.1). It can be defined through the plastically admissible stressfields: that is, t ∗ := max (cid:8) t > | ∃ ( π , Π ) ∈ X : t ˆ σ De = π − div Π ,f (cid:96) ( π , Π ) ≤ , Π n = on ∂ Ω F (cid:9) . (4.2)The results of Section 3 can be applied to find lower and upper bounds of t ∗ . With f (cid:96) defined by (2.7), (4.2) can be written t ∗ = max (cid:8) t > | ∃ ( π , Π ) ∈ X : ˆ σ De = π − div Π ,t (cid:112) | π | + (cid:96) − | Π | ≤ Y in Ω , Π n = on ∂ Ω F (cid:9) . (4.3)10he lower bound t ∗ L of t ∗ can then be written in terms of Π and ˆ σ De by eliminating π from (4.3): this gives t ∗ L = Y min Π (cid:13)(cid:13)(cid:13)(cid:113) | ˆ σ De + div Π | + (cid:96) − | Π | (cid:13)(cid:13)(cid:13) L ∞ (Ω) , (4.4)where Π is a third-order tensor-valued function such that Π n = on ∂ Ω F and div Π exists in an appropriate sense. For a sufficiently sharp bound it is necessary to minimizethe denominator in (4.4) with respect to Π . On the other hand, setting Π = , (4.4)reduces to t ∗ L, = Y / (cid:13)(cid:13) | ˆ σ De | (cid:13)(cid:13) L ∞ (Ω) (4.5)(cf. eqn (5.15) in [2] for a similar result in the discrete context). We note that t ∗ L, definesthe elastic threshold in perfect (Hencky) plasticity.The “kinematic” definition of t ∗ reads t ∗ = ˜ ϕ glob ( ˆ σ ) = inf q ∈ W (cid:82) Ω ˆ σ De : q dx =1 (cid:90) Ω D (cid:96) ( q , ∇ q ) dx = inf q ∈ W (cid:82) Ω ˆ σ De : q dx =1 (cid:90) Ω Y (cid:112) | q | + (cid:96) |∇ q | dx. (4.6)It enables us to find the upper bounds of t ∗ from t ∗ U = (cid:90) Ω Y (cid:112) | q | + (cid:96) |∇ q | dx (cid:90) Ω ˆ σ De : q dx , (4.7)for any q ∈ W such that (cid:82) Ω ˆ σ De : q dx > We consider here a subclass of problems for which the stress field ˆ σ e is constant in Ω. Forthis case the formulas (4.2) and (4.6) defining the elastic threshold t ∗ can be simplified.To this end, we consider the following forms of π , Π and q : π = t ˆ σ De ˜ π, Π = t ˆ σ De ⊗ ˜ Π , q = ˆ σ De ˜ q . (4.8)11ere ˜ π, ˜ q ∈ L (Ω), ˜ Π ∈ L (Ω; R ), and [ ˆ σ De ⊗ ˜ Π ] ijk = [ ˆ σ De ] ij [ ˜ Π ] k . Inserting (4.8) into(4.2) and (4.6), we obtain Y | ˆ σ De | ˜ t ∗ ≤ t ∗ ≤ Y | ˆ σ De | ˜ t ∗ , (4.9)where ˜ t ∗ = max (cid:8) t > | ∃ (˜ π, ˜ Π ) ∈ ˜ X : 1 = ˜ π − div ˜ Π ,t (cid:113) ˜ π + (cid:96) − | ˜ Π | ≤ , ˜ Π · n = 0 on ∂ Ω F (cid:9) , (4.10)and ˜ t ∗ = inf ˜ q ∈ ˜ W (cid:82) Ω ˜ q dx =1 (cid:90) Ω (cid:112) ˜ q + (cid:96) |∇ ˜ q | dx, ˜ W = { ˜ q ∈ W , (Ω) | ˜ q = 0 on ∂ Ω H } . (4.11)Here and henceforth W m,p (Ω) denotes the Sobolev space of functions which together withtheir weak derivatives of order ≤ m are in L p (Ω), for integer p ≥
1. Analogously as inSection 3, one can derive the duality between (4.10) and (4.11), implying ˜ t ∗ = ˜ t ∗ =: ˜ t ∗ .Hence, t ∗ = Y | ˆ σ De | ˜ t ∗ , (4.12)where ˜ t ∗ is defined either by (4.10) or (4.11). The lower and upper bounds (4.4) and (4.7)of t ∗ may be rewritten as follows: t ∗ L = Y | ˆ σ De | (cid:13)(cid:13)(cid:13)(cid:113) (1 + div ˜ Π ) + (cid:96) − | ˜ Π | (cid:13)(cid:13)(cid:13) L ∞ (Ω) , ˜ Π ∈ ˜ X Π , (4.13)and t ∗ U = Y | ˆ σ De | (cid:90) Ω (cid:112) ˜ q + (cid:96) |∇ ˜ q | dx (cid:82) Ω ˜ q dx , ˜ q ∈ ˜ W , (4.14)respectively, where˜ X Π = { ˜ Π ∈ L ∞ (Ω; R ) | div ˜ Π ∈ L ∞ (Ω; R ) , ˜ Π · n = 0 on ∂ Ω F } . Penalization methods for estimating the elastic thresh-old
The elastic threshold t ∗ may be estimated according to the results presented in the pre-vious section. Our aim is to find lower and upper bounds t ∗ L and t ∗ U which are sufficientlyclose to t ∗ . Let us recall that the bounds t ∗ L and t ∗ U are defined by appropriate functions Π and q , respectively. In some special cases, it is sufficient to choose Π and q by analyticalformulas. Nevertheless, the usage of numerical methods is more widely applicable.To find numerical solutions, it is necessary to discretize the problems (4.4) and (4.6). Ifwe choose finite-dimensional spaces X h and W h such that X h ⊂ X and W h ⊂ W then anyadmissible Π h ∈ X h and q h ∈ W h define lower and upper bounds of t ∗ , respectively.Next, it is difficult and/or slow to solve the problems (4.4) and (4.6) or their discretecounterparts directly because they contain non-differentiable functionals. Therefore, it isconvenient to transform or modify these problems. To this end, we make use of penal-ization methods. We describe the methods in the context of the simplified problem fromSection 4.1, corresponding to a constant stress field. Returning to equation (4.13) in Section 4.1, consider the following lower bound problem: s ∗ := inf ˜ Π ∈ ˜ X Π I ( ˜ Π ) , I ( ˜ Π ) = (cid:13)(cid:13)(cid:13) (1 + div ˜ Π ) + (cid:96) − | ˜ Π | (cid:13)(cid:13)(cid:13) L ∞ (Ω) , (5.1)where ˜ X Π = { ˜ Π ∈ L ∞ (Ω; R ) | div ˜ Π ∈ L ∞ (Ω; R ) , ˜ Π · n = 0 on ∂ Ω F } . Under the assumptions in Section 4.1, we have: t ∗ = Y | ˆ σ De | √ s ∗ , t ∗ L = Y | ˆ σ De | (cid:113) I ( ˜ Π ) ≤ t ∗ ∀ ˜Π ∈ ˜ X Π . For penalization of the problem (5.1), we replace the L ∞ norm with the L p one, p ≥ (cid:107) f (cid:107) L p (Ω) ≤ | Ω | /p (cid:107) f (cid:107) L ∞ (Ω) , f ∈ L ∞ (Ω) , p ≥ . s ∗ ≈ s /pp | Ω | /p , s p := inf ˜ Π ∈ ˜ X Π I p ( ˜ Π ) , I p ( ˜ Π ) = (cid:90) Ω (cid:104) (1 + div ˜ Π ) + (cid:96) − | ˜ Π | (cid:105) p dx . (5.2)The functional I p is convex and twice differentiable, and finite element approximationsto this problem can be found, for example, with the use of Raviart-Thomas or standardcontinuous piecewise-linear elements, which are conforming in the sense that the corre-sponding discrete space ˜ X Π ,h satisfies ˜ X Π ,h ⊂ ˜ X Π . Then, the discrete solution ˜ Π p,h ∈ ˜ X Π ,h defines the following lower bound of t ∗ : t ∗ ≥ Y | ˆ σ De | (cid:113) I ( ˜ Π p,h ) . (5.3)It is possible also first to discretize (5.1) and then to apply the penalization method tothe discrete problem.Let us note that numerical integration is required to evaluate (5.2). Nevertheless, if wefirst discretize (5.1) using, for example, lowest order Raviart-Thomas (RT0) or continuouspiecewise-linear (P1) elements, then it is not necessary to carry out numerical integration.Indeed, RT0- or P1-based functions ˜ Π h are linear on any simplicial finite element T ∈ R d with vertices V T, , V T, , . . . V T,d +1 . Hence,max x ∈ T (cid:110) (1 + div ˜ Π h ) + (cid:96) − | ˜ Π h | (cid:111) = max i =1 , ,...,d +1 (cid:110) (1 + div ˜ Π h ( V T,i )) + (cid:96) − | ˜ Π h ( V T,i ) | (cid:111) and t ∗ h = Y | ˆ σ De | (cid:113) I ( ˜ Π h ) ≤ t ∗ , I ( ˜ Π h ) = max T ∈T h max i =1 ,...,d +1 (cid:110) (1 + div ˜ Π h ( V T,i )) + (cid:96) − | ˜ Π h ( V T,i ) | (cid:111) . The corresponding penalized functional can be written as I p,h ( ˜ Π h ) = 1 p (cid:88) T ∈T h d +1 (cid:88) i =1 (cid:110) (1 + div ˜ Π h ( V T,i )) + (cid:96) − | ˜ Π h ( V T,i ) | (cid:111) p . For numerical solution of the discrete penalized problem, we combine adaptive contin-uation over p with the Newton method. One can expect that the larger the value of p , the sharper the lower bound (5.3) would be. At the same time, continuation over p is important for the initialization of the Newton method. Within the continuation, we14ncrease the values of p depending on the change of the value I ( ˜ Π p,h ). This approach willbe applied to examples presented in Sections 6.3 and 7.4. Consider the upper bound problem (see(4.14))˜ t ∗ = inf ˜ q ∈ ˜ W (cid:82) Ω ˜ q dx =1 J (˜ q ) , J (˜ q ) := (cid:90) Ω (cid:112) ˜ q + (cid:96) |∇ ˜ q | dx, (5.4)where ˜ W = { ˜ q ∈ W , (Ω) | ˜ q = 0 on ∂ Ω H } . Under the assumptions from Section 4.1, we have: t ∗ = Y | ˆ σ De | ˜ t ∗ , t ∗ ≤ Y | ˆ σ De | J (˜ q ) ∀ ˜ q ∈ ˜ W .
The penalization described below has been successfully used in limit load analysis inperfect plasticity [17, 18, 28, 19]. The local dissipation function, that is, the integrand in(5.4), may be written in the form D (cid:96) (˜ q, ∇ ˜ q ) = (cid:112) ˜ q + (cid:96) |∇ ˜ q | = sup (˜ π, ˜ Π ) ∈ R × R ˜ π + (cid:96) − | ˜ Π | ≤ { ˜ π ˜ q + ˜ Π · ∇ ˜ q } . (5.5)It can be penalized as follows: D α (˜ q, ∇ ˜ q ) = max (˜ π, ˜ Π ) ∈ R × R ˜ π + (cid:96) − | ˜ Π | ≤ (cid:26) ˜ π ˜ q + ˜ Π · ∇ ˜ q − α (˜ π + (cid:96) − | ˜ Π | ) (cid:27) , α > . (5.6)We have: D α (˜ q, ∇ ˜ q ) = α (˜ q + (cid:96) |∇ ˜ q | ) , (cid:112) ˜ q + (cid:96) |∇ ˜ q | ≤ α (cid:112) ˜ q + (cid:96) |∇ ˜ q | − α , (cid:112) ˜ q + (cid:96) |∇ ˜ q | ≥ α , (5.7)with D α ≤ D and D α → D as α → + ∞ . In addition, the function D α is convex, differ-entiable, and its second derivative exists in a generalized sense. The penalized problemreads ˜ t ∗ α = inf ˜ q ∈ ˜ W (cid:82) Ω ˜ q dx =1 J α (˜ q ) , J α (˜ q ) := (cid:90) Ω D α (˜ q, ∇ ˜ q ) dx. (5.8)15s with the penalized lower bound problem, discrete approximations to (5.8) may besought, for example, with the use of conforming piecewise-linear or quadratic finite ele-ments with the corresponding discrete space ˜ W h satisfying ˜ W h ⊂ ˜ W . Then, the discretesolution ˜ q ∗ α,h ∈ ˜ W h defines the following upper bound of t ∗ : t ∗ ≤ Y | ˆ σ De | J (˜ q ∗ α,h ) . (5.9)For a convergence analysis with respect to the discretization parameter h , we refer to[17, 18].For numerical solution of the discrete counterpart to (5.8), we combine an adaptive con-tinuation over α with the semismooth Newton method. One can expect that the largerthe value of α , the sharper the upper bound (5.9) would be. It is possible to study theconvergence ˜ t ∗ α,h → ˜ t ∗ h as α → + ∞ , similarly as in [17, 18]. Continuation over α is alsoimportant for the initialization of the Newton method. Within the continuation, we in-crease α depending on the change of the value J (˜ q ∗ α,h ). This approach will be applied toexamples presented in Sections 6.3 and 7.4. We consider the example introduced in [2, Section 6.2], of a thin plate subjected to biaxialextension. The plate is square with sides of length L and thickness L . Its domain isthus given by Ω = Ω × Ω × Ω , whereΩ = Ω = (0 , L ) , Ω = ( − L / , L / . The macroscopic boundary conditions are as follows (see Figure 6.1): t = t = t = 0 on Ω × ∂ Ω × Ω × (0 , T ) , (6.1a) u = t = t = 0 on ∂ Ω × Ω × Ω × (0 , T ) , (6.1b) u = t = t = 0 on Ω × Ω × { } × (0 , T ) , (6.1c) u = 2 t ˆ u, t = t = 0 on Ω × Ω × { L } × (0 , T ) , (6.1d)16
00 0 1 , ˆ σ e = 2 λ − νν (1 − ν ) ˆ uL ν . (6.2)We see that ˆ ε e , ˆ σ e are spatially constant, so that it is possible to use the results fromSection 4.1. To this end we shall also need the following formulas for ˆ σ De and its norm:ˆ σ De = 2 λ − ν ν (1 − ν ) ˆ uL − ν − − ν
00 0 2 − ν , (6.3) | ˆ σ De | = 2 λ (cid:114)
23 (1 − ν ) √ − ν + ν ν (1 − ν ) ˆ uL . (6.4)Next, we need to define microscopic boundary conditions. To be consistent with the planestress assumptions, we impose the micro-free conditionsΠ ij = 0 on Ω × ∂ Ω × Ω × (0 , T ) . (6.5)17n the remaining parts of the boundary we consider either micro-free or micro-hardboundary conditions; that is,Π ij = 0 on ∂ Ω × Ω × Ω × (0 , T ) , Π ij = 0 on Ω × Ω × ∂ Ω × (0 , T ) , (6.6)or ε p = on [ ∂ Ω × Ω × Ω × (0 , T )] ∪ [Ω × Ω × ∂ Ω × (0 , T )] . (6.7)In the following subsections, we derive lower and upper bounds of the elastic threshold t ∗ based on analytical and numerical approaches. To this end, we use the following inputdata: L = 50 , L = 2 , T = 1 / , λ = 1 . × − , ν = 0 . , ˆ u = 2 / , Y = 1 × − , (cid:96) = 5 . (6.8) t ∗ based on analytical approach If we consider the micro-free boundary condition (6.6), the spaces ˜ X Π and ˜ W correspond-ing to the estimates (4.13) and (4.14) are˜ X Π = { ˜ Π ∈ L ∞ (Ω; R ) | div ˜ Π ∈ L ∞ (Ω; R ) , ˜Π ijk = 0 for x k ∈ ∂ Ω k , k = 1 , , } , ˜ W = W , (Ω) . Choosing ˜ Π = and ˜ q = 1 in (4.13) and (4.14), respectively, we find that the lower andupper bounds coincide and thus t ∗ = t ∗ L = t ∗ U = Y | ˆ σ De | = (cid:114) Y L ν (1 − ν )2ˆ uλ (1 − ν ) √ − ν + ν . = 0 . . (6.9)This coincides also with the solution for the case of conventional plasticity. To completethe elastic solution from (4.1) for any t ∈ [0 , t ∗ ], we set ˆ Π e = and ˆ π e = ˆ σ De .If we consider the micro-hard boundary condition (6.7), the spaces ˜ X Π and ˜ W correspond-ing to the estimates (4.13) and (4.14) are now˜ X Π = { ˜ Π ∈ L ∞ (Ω; R ) | div ˜ Π ∈ L ∞ (Ω; R ) , ˜Π ij = 0 for x ∈ ∂ Ω } , ˜ W = { ˜ q ∈ W , (Ω) | ˜ q = 0 for x i ∈ ∂ Ω i , i = 1 , } . q = 1 (cid:54)∈ ˜ W . Therefore, (6.9) is nowonly a lower bound. To derive a sharper lower bound than t ∗ L, . = 0 . Π in the form (see (4.8))˜Π i ( x ) = a i x i + b i , a i , b i ∈ R , i = 1 , , no sum on i . (6.10)We shall optimize the parameters a i , b i to achieve the largest value of (4.13). To this end,we minimize the functional I ( ˜ Π ) := (cid:13)(cid:13)(cid:13) (1 + div ˜ Π ) + (cid:96) − | ˜ Π | (cid:13)(cid:13)(cid:13) L ∞ (Ω) = max x ∈ Ω (cid:110) (1 + div ˜ Π ) + (cid:96) − | ˜ Π | (cid:111) (6.11)representing the denominator in (4.13). The maximum over x i ∈ Ω i is achieved at oneof the end points of the interval Ω i , i.e., for x i ∈ ∂ Ω i , i = 1 ,
3. To minimize the terms | a i x i + b i | on ∂ Ω i , we set b i = − L a i ( i = 1 , I ( ˜ Π ) = (1 + a + a ) + (cid:18) L (cid:96) (cid:19) ( a + a ) . (6.12)The optimal values of a i are found by minimization of (6.12), leading to the followinglower bound of t ∗ : t ∗ L, = Y | ˆ σ De | (cid:96)L (cid:114) L (cid:96) . = 0 . . (6.13)We see that the elastic threshold t ∗ L, for perfect (Hencky) plasticity is smaller than t ∗ L, ,i.e., then t ∗ . Next, from the presented results, it follows that the microstress Π corre-sponding to the solution to the gradient plasticity problem is not unique. Especially forsufficiently small t >
0, we derived two different solutions corresponding to (4.1): Π = and Π = t ˆ σ De ⊗ ˜ Π , where˜Π ( x ) = L − x L (cid:96) , ˜Π ( x ) = 0 , ˜Π ( x ) = L − x L (cid:96) . Now, we estimate t ∗ from above. We choose ˜ q ( x ) = x ( L − x ) x ( L − x ) ∈ ˜ W . Then(4.14) yields t ∗ U . = 0 . Π and ˜ q through a numerical procedure. t ∗ based on numerical methods We consider the micro-hard boundary conditions (6.7). Under the plane stress assump-tion, the lower and upper bound problems are posed on the domain ˜Ω = (0 , × (0 , t ∗ = Y | ˆ σ De | (cid:114) inf ˜ Π ∈ ˜ X Π I ( ˜ Π ) , I ( ˜ Π ) = (cid:13)(cid:13)(cid:13) (1 + div ˜ Π ) + (cid:96) − | ˜ Π | (cid:13)(cid:13)(cid:13) L ∞ (˜Ω) , (6.14)where ˜ X Π = { ˜ Π ∈ L ∞ ( ˜Ω; R ) | div ˜ Π ∈ L ∞ ( ˜Ω; R ) } , and the upper bound one reducesto t ∗ = Y | ˆ σ De | ˜ t ∗ , ˜ t ∗ = inf ˜ q ∈ ˜ W (cid:82) ˜Ω ˜ q dx =1 (cid:90) ˜Ω (cid:112) ˜ q + (cid:96) |∇ ˜ q | dx, ˜ W = W , ( ˜Ω) . (6.15)We also know that Y / | ˆ σ De | . = 0 . t ∗ based on numerical techniques, we use the penaliza-tion methods presented in Section 5. For the solution of both problems, we use a regularmesh with 401 ×
401 nodes. Problem (6.14) is discretized using RT0 elements, whilepiecewise quadratic P2 elements with 7-point Gauss quadrature are used for solving theproblem (6.15). Both problems are implemented in MATLAB.(a) (b)Figure 2: (a) The norm | ˜ Π | of the solution to (6.14); (b) the corresponding functional I Computed lower and upper bounds depending on the penalization parameters p and α ,20a) (b)Figure 3: (a) The solution ˜ q to (6.15); (b) the corresponding local dissipation D (cid:96) (˜ q )respectively, are depicted in Figure 4. From the largest used values of these parameters,we find the following bounds of t ∗ :0 . ≤ t ∗ ≤ . . We see that both the lower and upper bounds are very sharp. In particular, the lowerbound has been improved in comparison of Section 6.2.The numerical solution to the problem (6.14) arising from the largest value of p is depictedin Figure 2. In particular, the norm of | ˜ Π | is visualized on the left figure, while local valuesdefining the functional I are depicted on the right. This solution does not have a simpleshape, with a corner effect evident. From the right figure, one can see that the local valuesof I are almost constant except at the corner points. Similar results have been obtainedwith the use of P1 elements.The numerical solution to the problem (6.15) corresponding to the largest value of α isdepicted in Figure 3. Also shown is the corresponding local dissipation D (cid:96) (see (5.5)).This figure shows large values of ∇ ˜ q along the edges. Consider a rod of radius R and length L , with 0 < L ≤ ∞ , subject to a uniform torsionwith κ denoting the twist per unit length. The only non-zero displacement, relative to a21a) (b)Figure 4: Lower (a) and upper (b) bounds t ∗ corresponding to the functionals (6.14) and(6.15)cylindrical coordinate system ( r, ϑ, z ) is u ϑ = κzr . The only stress in the elastic region isthen σ ϑz = µκr , where µ is the shear modulus. In accordance with the notation introducedin Section 4, we set κ ( t ) = t ˆ κ , where ˆ κ > σ De = ˆ σ e and | ˆ σ De | = µ ˆ κr . We shall also use the notation ˆ σ ϑz = µ ˆ κr .In the context of size-dependent and microscale mechanical response, this problem hasbeen studied experimentally [10] as well as theoretically and computationally. The study[20] shows the link between strengthening, that is, an increase in yield strength with lengthscale, which characterizes the dissipative problem, the topic of this work. The work [4]studies the energetic problem, in which gradient terms arise through a recoverable energy.In this case the link is between hardening, that is, an increase in slope in the plasticrange, and increase in length scale. In [1] the torsion problem is discussed in the contextof a more elaborate gradient theory that takes account of plastic rotation and dislocationdensity.We assume that Π ϑzr ( r ) and Π rzϑ ( r ) are the only non-zero components of Π . Thenmicroforce balance (2.5) and the yield function (2.7) reduce respectively to σ ϑz = π ϑz − Π ϑzr,r − r − (Π ϑzr + Π rzϑ ) (7.1)and f (cid:96) ( π ϑz , Π ϑzr , Π rzϑ ) = (cid:113) π ϑz + (cid:96) − (Π ϑzr + Π rzϑ ) − Y . (7.2)in (0 , R ). The boundary conditions are a traction-free surface r = R , which is automati-cally satisfied, a microhard condition at r = 0, and a microfree condition at r = R . That22s, u ϑ (0) = 0 , ε pϑz (0) = 0 , Π ϑzr ( R ) = 0 . (7.3)From (4.3), (7.1) and (7.2), the elastic threshold for this problem is given by t ∗ = max (cid:110) t > | ∃ ( π, Π , ¯Π) ∈ X : t ˆ σ = π − Π (cid:48) − r − (Π + ¯Π) , (cid:113) π + (cid:96) − (Π + ¯Π ) ≤ Y in (0 , R ) , Π( R ) = 0 (cid:111) , (7.4)where ˆ σ := ˆ σ ϑz , Π := Π ϑzr and ¯Π := Π rzϑ . If Y is constant, then from (4.4) we have t ∗ = Y min (Π , ¯Π) ∈ X Π max r ∈ [0 ,R ] (cid:112) (ˆ σ + Π (cid:48) + r − (Π + ¯Π)) + (cid:96) − (Π + ¯Π ) , (7.5)where X Π = (cid:8) (Π , ¯Π) ∈ W , ∞ ((0 , R )) × L ∞ ((0 , R )) | Π( R ) = 0 ,r − (Π + ¯Π) ∈ L ∞ ((0 , R )) (cid:9) . (7.6)We may render (7.5) dimensionless by setting˜ r := rR , ˜ (cid:96) = (cid:96)R , ˜Π = Π µ ˆ κR , ˜¯Π = Π µ ˆ κR . Then (7.5) becomes t ∗ = Yµ ˆ κR ˜ t ∗ , ˜ t ∗ = 1 (cid:114) min ( ˜Π , ˜¯Π) ∈ ˜ X Π I ( ˜Π , ˜¯Π) , (7.7)with I ( ˜Π , ˜¯Π) = max ˜ r ∈ [0 , (cid:20)(cid:16) ˜ r + ˜Π (cid:48) (˜ r ) + ˜ r − ( ˜Π(˜ r ) + ˜¯Π(˜ r )) (cid:17) + ˜ (cid:96) − (cid:16) ˜Π (˜ r ) + ˜¯Π (˜ r ) (cid:17)(cid:21) (7.8)and ˜ X Π = (cid:8) ( ˜Π , ˜¯Π) ∈ W , ∞ ((0 , × L ∞ ((0 , | ˜Π(1) = 0 ,r − ( ˜Π + ˜¯Π) ∈ L ∞ ((0 , (cid:9) . (7.9)We see that ˜ t ∗ depends only on ˜ (cid:96) . In (7.7)–(7.9), it is also convenient to introduce23ˆΠ := r − ( ˜Π + ˜¯Π) and to eliminate ˜¯Π in (7.8) by the substitution˜¯Π = r ˜ˆΠ − ˜Π . (7.10)It should be noted that one cannot then use the weighted functional space; see Section7.3.To derive the dual definition of the threshold t ∗ , we introduce a weak form for the mi-croforce balance equation (7.1). Given the homogeneous Dirichlet condition (7.3), wedefine W = { q ∈ W , (0 , R ) | q (0) = 0 } . (7.11)After multiplying (7.1) by qrdr and integrating we obtain t (cid:90) R ˆ σqr dr = (cid:90) R [ πqr + Π( rq (cid:48) ) − ¯Π q ] dr ∀ q ∈ W .
The expression for t ∗ then becomes t ∗ = sup ( π, Π , ¯Π) ∈ X √ π + (cid:96) − (Π + ¯Π ) ≤ Y in (0 ,R ) inf q ∈ W (cid:82) R ˆ σqr dr =1 (cid:90) R [ πqr + Π( rq (cid:48) ) − ¯Π q ] dr = inf q ∈ W (cid:82) R ˆ σqr dr =1 (cid:90) R sup ( π, Π) ∈ X √ π + (cid:96) − (Π + ¯Π ) ≤ Y [ πqr + Π( rq (cid:48) ) − ¯Π q ] dr = inf q ∈ W (cid:82) R ˆ σqr dr =1 (cid:90) R Y (cid:112) q + (cid:96) [( q (cid:48) ) + ( q/r ) ] r dr. (7.12)Substituting ˆ σ = µ ˆ κr , ˜ r := r/R , ˜ (cid:96) = (cid:96)/R , (7.12) can be transformed to t ∗ = YνκR ˜ t ∗ , ˜ t ∗ = inf q ∈ ˜ W (cid:82) ˜ r q d ˜ r =1 (cid:90) (cid:113) q + ˜ (cid:96) [( q (cid:48) ) + ( q/ ˜ r ) ] ˜ r d ˜ r, (7.13)where ˜ W = { q ∈ W , (0 , | q (0) = 0 } . 24 .2 Upper and lower bounds of ˜ t ∗ based on an analytical ap-proach To find an upper bound of ˜ t ∗ we choose q (˜ r ) = 4˜ r , in (7.13). This yields˜ t ∗ ≤ ˜ t ∗ U = 43 (cid:26)(cid:104) (cid:96) (cid:105) / − (2˜ (cid:96) ) / (cid:27) ≤ (cid:104) (cid:96) (cid:105) / . (7.14)If ˜ (cid:96) = 0 then choosing q (˜ r ) = ( n + 1) r n , we obtain t ∗ U → n → + ∞ .To determine a lower bound consider first of all ˜Π = ˜¯Π = 0. This gives, from (7.7),˜ t ∗ ≥ ˜ t ∗ L, = 1 . (7.15)We see that ˜ t ∗ L, = ˜ t ∗ = 1 for ˜ (cid:96) = 0. To obtain a sharper lower bound one can choose˜Π = a ˜ r (˜ r − , ˜¯Π = 0 , ˜ r ∈ [0 , , (7.16)where the parameter a ∈ R is optimized. This bound is depicted and compared againstother estimates in Figure 6. We adopt the following version of the lower bound problem, from (7.8) and (7.10) (tildesymbols are omitted for the sake of simplicity): t ∗ = Yµ ˆ κR ˜ t ∗ , ˜ t ∗ = 1 (cid:114) min (Π , ˆΠ) ∈ ˆ X Π I (Π , ˆΠ) , (7.17)where ˆ X Π = (cid:8) (Π , ˆΠ) ∈ W , ∞ ((0 , × L ∞ ((0 , | Π(1) = 0 (cid:9) and I (Π , ˆΠ) = max r ∈ [0 , (cid:26)(cid:16) r + Π (cid:48) ( r ) + ˆΠ( r )) (cid:17) + ˜ (cid:96) − (cid:20) Π ( r ) + (cid:16) r ˆΠ( r ) − Π( r ) (cid:17) (cid:21)(cid:27) . (7.18)We discretize and penalize this problem as in Section 5.1. Consider a uniform partitionof the interval [0 , r < r < . . . < r N +1 = 1 , r i +1 − r i = h = 1 /N, i = 1 , , . . . , N , X Π :ˆ X Π ,h = (cid:8) (Π h , ˆΠ h ) ∈ C ([0 , × L ∞ ((0 , (cid:12)(cid:12) Π h (1) = 0 , Π h | ( r i ,r i +1 ) ∈ P , ˆΠ h | ( r i ,r i +1 ) ∈ P , i = 1 , , . . . , N (cid:9) . Here P and P are respectively the spaces of polynomials of degree 1, and of constants,on the interval [ r i , r i +1 ]. Clearly, ˆ X Π ,h ⊂ ˆ X Π , which implies t ∗ ≥ t ∗ L,h = Yµ ˆ κR ˜ t ∗ L,h , ˜ t ∗ L,h = 1 (cid:114) min (Π h , ˆΠ h ) ∈ ˆ X Π ,h I (Π h , ˆΠ h ) . (7.19)In addition, we have I (Π h , ˆΠ h ) = max i max r ∈ [ r i ,r i +1 ] (cid:26)(cid:16) r + Π (cid:48) h + ˆΠ h ) (cid:17) + ˜ (cid:96) − (cid:20) Π h + (cid:16) r ˆΠ h − Π h (cid:17) (cid:21)(cid:27) = max i max (cid:110) [ r i + Π (cid:48) h,i + ˆΠ h,i ] + ˜ (cid:96) − [Π h ( r i ) + ( r i ˆΠ h,i − Π h ( r i )) ] , [ r i +1 + Π (cid:48) h,i + ˆΠ h,i ] + ˜ (cid:96) − [Π h ( r i +1 ) + ( r i +1 ˆΠ h,i − Π h ( r i +1 )) ] (cid:111) . (7.20)Here we have used the property that the quadratic function defined in [ r i , r i +1 ] has amaximum either at r i or r i +1 and Π (cid:48) h,i = Π (cid:48) h ( r ) = const , ˆΠ h,i = ˆΠ( r ) = const in ( r i , r i +1 ), i = 1 , , . . . , N .The corresponding penalized functional reads I p (Π h , ˆΠ h ) = 1 p n (cid:88) i =1 (cid:110)(cid:8) [ r i + Π (cid:48) h,i + ˆΠ h,i ] + ˜ (cid:96) − [Π h ( r i ) + ( r i ˆΠ h,i − Π h ( r i )) ] (cid:9) p + (cid:8) [ r i +1 + Π (cid:48) h,i + ˆΠ h,i ] + ˜ (cid:96) − [Π h ( r i +1 ) + ( r i +1 ˆΠ h,i − Π h ( r i +1 )) ] (cid:9) p (cid:111) . (7.21) t ∗ Let us recall the problems (7.13), (7.17) and (7.18), that is,˜ t ∗ = inf q ∈ ˜ W (cid:82) ˜ r q d ˜ r =1 (cid:90) (cid:113) q + ˜ (cid:96) [( q (cid:48) ) + ( q/ ˜ r ) ] ˜ r d ˜ r, ˜ W = { q ∈ H (0 , | q (0) = 0 } , (7.22)26 t ∗ = 1 (cid:115) min (Π , ˆΠ) ∈ ˆ X Π max r ∈ [0 , (cid:26)(cid:16) r + Π (cid:48) ( r ) + ˆΠ( r )) (cid:17) + ˜ (cid:96) − (cid:20) Π ( r ) + (cid:16) r ˆΠ( r ) − Π( r ) (cid:17) (cid:21)(cid:27) , (7.23)where ˆ X Π = (cid:8) (Π , ˆΠ) ∈ W , ∞ ((0 , × L ∞ ((0 , | Π(1) = 0 (cid:9) . To find lower and upperbounds of t ∗ , we use the penalization methods from Section 5.(a) (b)Figure 5: The upper (a) and lower (b) bounds of ˜ t ∗ h depending on p for ˜ (cid:96) = 0 . (cid:96) = 0 . (cid:96) = 1 . q α for various length scales; (b) comparison ofnumerical and analytical bounds of ˜ t ∗ For the upper bound problem (7.22), we use piecewise quadratic (P2) elements with 2-point Gauss quadrature and 1000 elements. By continuation over the penalty parameter α we find that α ≈
10 is sufficiently large for any ˜ (cid:96) ∈ (0 , (cid:96) = 0 . (cid:96) = 0 . (cid:96) = 1 . p ≈
300 is a sufficiently largevalue of the penalization parameter.Numerical solutions of problems (7.22) and (7.23) are depicted in Figure 6 (left) andFigure 7 for three different values of ˜ (cid:96) . Computable majorants of ˜ t ∗ following from Figure5 (a) are 1.2807, 1.9305, and 2.9529 for ˜ (cid:96) = 0 . , . , .
0, respectively. Computableminorants following from Figure 5 (b) are 1.2755, 1.9279, and 2.9493, respectively. We seethat the bounds are sufficiently sharp. Figure 6 (right) compares computable majorants(minorants) with analytical upper and lower bounds presented in Section 7.2.
This work has been concerned with the development of a general approach to determiningbounds on the load parameter corresponding to incipient plastic behaviour, that is, theelastic threshold, for the dissipative model of strain-gradient plasticity. Such an approachhas been motivated by the fact that the yield function, in a generalized associative model,is a function of microstresses and is therefore not a means through which initial yieldmay be determined while following an elastic path. The approach taken in this workis motivated by the theorems of limit analysis, with bounds being determined throughminimization or maximization problems.The theory has been applied to two example problems, both of which have receivedattention in the context of strain-gradient plasticity. For both problems – the first of a28late in plane stress, and subject to a prescribed displacement, and the second concerninga rod in uniform torsion – bounds are determined analytically by choosing candidatefunctions, and numerically by solution of the extremization problems using finite elements.In this latter case the non-differentiable functionals are regularized through the adoptionof penalized approximations.The theory and results presented here shed further light on the behaviour of the dissipativemodel of strain-gradient plasticity. For example, it is evident that the microstresses inthe elastic branch of the loading path are not uniquely defined and their distribution isnon-trivial, even if the Cauchy stress is known analytically.The techniques of limit analysis have proved to be a powerful tool in the study of thismodel. It would be useful to consider further applications of this approach, for exampleto problems of single-crystal plasticity and ensembles of crystals.
References [1] Bardella, L. and Panteghini, A., 2015. Modelling the torsion of thin metal wires bydistortion gradient plasticity.
J. Mech. Phys. Solids
78, 467–492.[2] Carstensen, C., Ebobisse, F., McBride, A.T., Reddy, B.D. and Steinmann, P., 2017.Some properties of the dissipative model of strain-gradient plasticity.
Phil. Mag.
Handbook of Numerical Analysis , Vol IV, Part 2, North-Holland,195–312.[4] Chiricotto, M., Giacomelli, L. and Tomassetti, G., 2016. Torsion in strain-gradientplasticity: energetic scale effects.
SIAM J. Appl. Math.
72, 1169–1191.[5] Chiricotto, M., Giacomelli, L. and Tomassetti, G., 2016. Dissipative scale effects instrain-gradient plasticity: the case of simple shear.
SIAM J. Appl. Math.
76 (2),688–704.[6] Gurtin, M.E. and Anand, L., 2004. A theory of strain-gradient plasticity for isotropic,plastically irrotational materials. Part I: Small deformations.
J. Mech. Phys. Solids
53, 1624–1649. 297] Ekeland, I. and Temam, R., 1974.
Analyse Convexe et Probl`emes Variationnels .Dunod, Gauthier Villars, Paris.[8] Evans, A. G. & Hutchinson, J.W., 2009. A critical assessment of theories of straingradient plasticity.
Acta Mater.
57, 1675–1688.[9] Fleck, N. A. and Hutchinson, J.W., 2001. A reformulation of strain gradient plas-ticity,
J. Mech. Phys. Solids
49, 2245–2271.[10] Fleck, N.A., Muller, G.M., Ashby, M.F. and Hutchinson, J.W., 1994. Strain gradientplasticity: theory and experiments.
Acta Metall. Mater.
42, 475–487.[11] Fleck, N.A., Hutchinson, J.W. and Willis, J.R., 2014. Strain gradient plasticity undernon-proportional loading.
Proc. R. Soc. Lond. A Math. Phys. Eng. Sci.
J. Mech. Phys. Solids
57, 1045–1057.[13] Fleck, N.A. and Willis, J.R., 2015. Strain gradient plasticity: Energetic or dissipa-tive?,
Acta Mech. Sinica
31, 465–472.[14] Gudmundson, P., 2001. A unified treatment of strain gradient plasticity,
J. Mech.Phys. Solids
52, 1379–1406.[15] Gurtin, M. E. and Anand, L., 2004. A theory of strain-gradient plasticity forisotropic, plastically irrotational materials. Part I: Small deformations,
J. Mech.Phys. Solids
53, 1624–1649.[16] Han, W. and Reddy, B. D., 2013.
Plasticity: Mathematical Theory and NumericalAnalysis . Second Edition, Springer, New York and Berlin.[17] Haslinger, J., Repin, S. and Sysala, S., 2016. A reliable incremental method ofcomputing the limit load in deformation plasticity based on compliance: Continuousand discrete setting.
J. Comp. Appl. Math.
Appl.Math.
61, 527–564.[19] Haslinger, J., Repin, S. and Sysala, S., 2019. Inf-sup conditions on convex cones andapplications to limit load analysis.
Math. Mech. Solids
24, 3331–3353.3020] Idiart, M. and Fleck, N., 2010. Size effects in the torsion of thin metal wires.
Modell.Simul. Mater. Sci. Eng.
18, 015009.[21] Martinez-Pa˜neda, E., Niordson, C. F. and Bardella, L., 2016. A finite element frame-work for distortion gradient plasticity with applications to bending of thin foils.
Int.J. Solids Struct.
96, 288–299.[22] McBride, A. T., Reddy, B. D. and Steinmann, P., 2018. Dissipation-consistent mod-elling and classification of extended plasticity formulations.
J. Mech. Phys. Solids
Comp. Meths Appl. Mech. Engng
International journal of solids and structures , 47(1), 100-112.[25] Reddy, B.D., 2011a. The role of dissipation and defect energy in variational formu-lations of problems in strain- gradient plasticity. Part 1: single-crystal plasticity.
Cont. Mech. Thermodyn.
23, 551–572.[26] Reddy, B.D., 2011b. The role of dissipation and defect energy in variational formu-lations of problems in strain- gradient plasticity. Part 2: polycrystalline plasticity.
Cont. Mech. Thermodyn.
23, 527–549.[27] Reddy, B.D., Ebobisse, F. and McBride, A. T., 2008. Well-posedness of a modelof strain gradient plasticity for plastically irrotational materials,
Int. J. Plast.
Comp. Math. Appl.
75, 199–217.[29] Temam R., 1985.