On the formulation of closest-point projection algorithms in elastoplasticity. Part II: globally convergent schemes

REPORT NO. UCB/SEMM-2000/02 STRUCTURAL ENGINEERING, MECHANICS AND MATERIALS BY DEPARTMENT OF CIVIL AND ENVIRONMENTAL ENGINEERING UNIVERSITY OF CALIFORNIA BERKELEY, CALIFORNIA JANUARY 2000 AND ON THE FORMULATION OF CLOSEST-POINT PROJECTION ALGORITHMS IN ELASTOPLASTICITY. PART II: GLOBALLY CONVERGENT SCHEMES. With Appplications to Deviatoric and Pressure-Dependent Plastic Models. F. ARMERO A. PEREZ - FOGUET ' On the Formulation of Closest-Point Projection Algorithms in Elastoplasticity. Part II: Globally Convergent Schemes. With Applications to Deviatoric and Pressure-Dependent Plastic Models by A. P´ erez-Foguet†& F. Armero∗ Structural Engineering, Mechanics and Materials Department of Civil and Environmental Engineering University of California, Berkeley CA 94720, USA Abstract This paper presents the formulation of numerical algorithms for the solution of the closest-point projection equations that appear in typical implementations of return mapping algorithms in elastoplasticity. The main motivation behind this work is to avoid the poor global convergence properties of a straight application of a Newton scheme in the solution of these equations, the socalled Newton–CPPM. The mathematical structure behind the closest-point projection equations identified in Part I of this work delineates clearly different strategies for the successful solution of these equations. In particular, primal and dual closest-point projection algorithms are proposed, in non-augmented and augmented Lagrangian versions for the imposition of the consistency condition. The primal algorithms involve a direct solution of the original closest-point projection equations, whereas the dual schemes involve a two level structure by which the original system of equations is staggered, with the imposition of the consistency condition driving alone the iterative process. Newton schemes in combination with appropriate line search strategies are considered, resulting in the desired asymptotically quadratic local rate of convergence and the sought global convergence character of the iterative schemes. These properties, together with the computational performance of the different schemes, are evaluated through representative numerical examples involving different models of finite strain plasticity. In particular, the avoidance of the large regions of no convergence in the trial state observed in the standard Newton–CPPM is clearly illustrated. KEY WORDS: elastoplasticity; return mapping algorithms; closestpoint projection; globally convergent schemes; augmented Lagrangian. Submitted to International Journal for Numerical Methods in Engineering. †On leave from the Department of Matem`atica Aplicada III, ETSECCPB, UPC, Barcelona, Spain, during the Fall semester 1999. ∗Corresponding author ([email protected]). A. P´erez-Foguet & F. Armero 2 1. Introduction The developments presented in Armero & P´ erez–Foguet [2000] (referred simply as Part I hereafter) have identified a rich mathematical structure behind the closest-point projection equations characteristic of the local integration of elastoplastic models, including their viscoplastic extensions, in both the infinitesimal and finite deformation ranges. Box 1.1 summarizes these equations following the same notation employed in this previous part of this work. More specifically, the updated values (·)n+1 of the stresses and internal variables at the end of a typical time increment [tn,t n+1] in the solution of the global mechanical boundary-value problem are obtained from their known counterparts (·)nfor a given strain increment ∆ε=εn+1 −εn. Associated to this increment, we have the trial values (·)trial n+1 of the different variables obtained under the assumption of incremental elastic response. The case of rate-independent elastoplasticity in a typical plastic corrector step, characterized by the enforcement of the consistency condition fn+1 = 0 from an initial trial solution with ftrial n+1 >0, has been considered in Box 1.1. The discrete forms of the flow rule and hardening laws have been written without considering explicitly the usual assumptions of convexity and associativity required for a rigorous characterization of the aforementioned variational structure. It is in this general context where we present the numerical algorithms proposed in this work, even though the statements on the convergence properties of the different schemes apply rigorously to the associated convex case only. Box 1.1 summarizes also the application of a Newton scheme in the solution of the resulting nonlinear algebraic equations. To this purpose, the equations have been written in residual form and the Jacobian matrix utilized in the iterative process identified. The final scheme, referred to as the Newton–CPPM herein, can be found in Simo & Hughes [1998] and has become a very popular strategy in the numerical solution of elastoplastic problems. As noted in Section 1 of Part I of this work, the asymptotic quadratic rate of convergence of this iterative scheme makes this approach very attractive. However, its limited global convergence properties leads often to a lack of convergence for large ranges of the initial trial state. This situation is related to the strong nonlinearity of the equations which may arise, for example, from the high curvature of the yield function, as it can be found in typical practical applications. The numerical results reported in de Souza Neto et al. [1994], Bi´ cani´ c&Pearce[1996], P´ erez–Foguet et al. [2000a] and herein illustrate these difficulties. These observations identify clearly the goal of this work: the development of efficient globally convergent algorithms for the solution of the closest-point projection equations in elastoplasticity. We focus only on schemes that preserve the local asymptotic rate of convergence characteristic of the original Newton– CPPM. We note also that all the algorithms proposed in this work lead to the exact solution of the original closest-point projection equations. This implies, in particular, that the consistent algorithmic tangent employed in a typical Newton-Raphson solution of the global mechanical boundary-value problem coincides in all cases; see e.g. Simo & Hughes [1998]. Closest-Point Projection Algorithms in Elastoplasticity 3 BOX 1.1. The Newton closest-point projection algorithm for rateindependent elastoplasticity. •For the plastic corrector step (i.e. ftrial n+1 >0), define the residuals rε:= εe n+1 −εe,trial n+1 +∆γmσ(σn+1,qn+1), rα:= −αn+1 +αtrial n+1 +∆γmq(σn+1,qn+1), rf:= f(σn+1,qn+1), for the elastic strain εe n+1 and strain-like hardening variables αn+1, and where σn+1 =∂eψ(εe n+1,αn+1)andqn+1 =−∂αψ(εe n+1,αn+1) the stresses σn+1 and stress-like hardening variables qn+1,forthe stored energy function ψ(εe,α) and yield function f(σ,q). •Let r(x):=   r rα rf  for x:=    εe n+1 −αn+1 ∆γ  . Consider the Newton iterative scheme x(k+1) =x(k)+d(k)for d(k)=−J−1(x(k))r(x(k)), with x(0) =εe,trial n+1 ,−αtrial n+1 ,0Tand the Jacobian J=∇r=  I+∆γ∇m˜ Gm ∇fT˜ G0 , where ∇f=∂σf ∂qf,m=mσ mq,∇m=∂σmσ∂qmσ ∂σmq∂qmq, and ˜ G=∂2 εeεeψ−∂2 εeαψ −∂2 αεeψ∂ 2 ααψ. Remark: The condition ∇f·˜ Gm >0 strictly is assumed for the underlying continuum model. A. P´erez-Foguet & F. Armero 4 An algorithm is globally convergent if all the sequences that it generates converge to a solution of the problem. We refer to the classical global convergence theorem presented in Luenberger [1989] (page 187) for details. In the case of interest herein, this property corresponds to convergence for any initial trial state. As studied in detail in this last reference and other classical references on nonlinear optimization, this desired property can be obtained through the consideration of the appropriate line search scheme under certain technical conditions (usually involving the convexity of the variational problem in optimization and the smoothness of the functions involved). We explore this alternative first, by combining the original Newton scheme with a line search designed specifically for the equations of interest here. In particular, the constrained character of the original closest-point projection equations implied by the non-negative character of the plastic multiplier (∆γ≥0 in Box 1.1) is taken into account explicitly in the iterative process. Appendix I presents the details of the line search schemes developed in this work. We refer to the resulting algorithm as the primal closest-point projection method (the primal– CPPM in short), given its clear relation to the primal variational formulation of the closestpoint projection equations under the usual assumption of convexity and associativity. This first proposed algorithm consists then simply of the Newton scheme with line search. The special considerations required by the constrained character of the governing equations, with the consequent need of special variations in the computer code, motivate the development of alternative techniques. The consideration of the new augmented primal formulations identified in Part I appears then as a clear option, as it is developed next in this paper. We refer to the final scheme as the augmented–primal–CPPM. In addition of involving a regularized unconstrained problem, the augmentation adds also a simple mechanism for improving the computational efficiency of the original primal–CPPM. Numerical experiments identify in this way a range of the regularizing penalty parameter leading to this improvement. It is interesting to observe that the original Newton–CPPM has been often interpreted in the literature as imposing the closest-point projection for each iterate, (·)(k)in the notation of Box 1.1; see Simo & Hughes [1998], page 143 (Figure 3-10 in this reference, in particular). However, a look at the residual form of the equations clearly shows that this interpretation is not valid since, simply, the residuals r(k)corresponding to the flow rule and hardening law do not vanish, in general, during the iteration process, except at the converged value. This observation identifies clearly the role played by the consistency condition and led to the dual variational formulations presented in Part I. In this context, the solution can be obtained simply as the roots of the yield function understood as a scalar function of the plastic multiplier. The evaluation of this “dual consistency function” requires the solution of the discrete forms of the flow rule and hardening law for a fixed plastic multiplier and, hence, unconstrained in nature. These considerations lead naturally to a two-level algorithm that we refer as the dual–CPPM. The consideration of a Newton scheme in combination with a line search in both levels assures the asymptotic quadratic Closest-Point Projection Algorithms in Elastoplasticity 5 local rate of convergence of the algorithm and its globally convergent character, this latter property holding again under the usual technical conditions (i.e. the associated convex case with smooth material relations). The augmented extension of these ideas, referred to as the augmented–dual–CPPM herein, proves also to be a useful tool to add computational efficiency to the final numerical scheme. An outline of the rest of the paper is as follows. Sections 2 and 3 describe the details behind the primal and dual algorithms, respectively, in their non-augmented and augmented forms. Appendix I includes again the details of the line search schemes developed in this work for constrained and unconstrained problems. We focus the developments on the rate-independent problem, with extensions to the simpler viscoplastic problem considered briefly in remarks in the end of each section. The properties and performance of the primal and dual classes of algorithms are evaluated numerically in Sections 4 and 5, respectively. Several representative examples are considered involving purely deviatoric and pressure dependent plastic models (both in an associated framework, with the second one involving non-smooth flow vectors), strain hardening and softening, as well as constant and non-constant elasticities. Finite deformation multiplicative plasticity models are considered in all cases. Appendix II includes a summary of the material models considered. The solution of a finite element problem using the primal–CPPM is shown in Section 6, as an example of practical application. Finally, Section 7 presents several concluding remarks, comparing in particular the different proposed algorithms. 2. Primal Algorithms We develop in this section numerical algorithms based on the primal variational formulations developed in Sections 3 and 4 of Part I of this work for the rate-independent and viscoplastic problems, respectively. More precisely, we address the solution of the corresponding Euler-Lagrange equations or, more generally, their counterparts in the context of non-associated formulations without this variational structure. Extensions involving the corresponding augmented Lagrangian formulations, as developed in Section 5 of Part I, are presented in Section 2.2. We consider the rate-independent elastoplastic problem, with the viscoplastic problem discussed briefly in additional remarks in the end of each section. 2.1. The primal closest-point projection method A simple and efficient choice when trying to improve the convergence properties of the original Newton closest-point projection method for the rate-independent problem summarized in Box 1.1 is the additional consideration of a line search technique, as described in Appendix I. We consider again the residual form of the equations presented in this box, A. P´erez-Foguet & F. Armero 6 namely, r(x):=     εe n+1 −εe,trial n+1 +∆γmσn+1 −αn+1 +αtrial n+1 +∆γmqn+1 fn+1     ,for x:=      εe n+1 −αn+1 ∆γ     ,(2.1) the elastic strains, strain-like internal hardening variables, and the discrete plastic multiplier, respectively, all of them at the end tn+1 of a typical time step [tn,t n+1]. The flow vectors mσn+1 and mqn+1 , and yield function fn+1 in (2.1) are given functions of the conjugate stress variables σn+1 =∂εeψn+1 and qn+1 =−∂αψn+1 for the stored energy function ψn+1 =ψ(εe n+1,αn+1), as considered in Box 1.1. The interest herein is the consideration of a plastic corrector step, detected by the condition ftrial n+1 := f(σtrial n+1 ,qtrial n+1 )>0forthe trial elastic values {εe,trial n+1 ,−αtrial n+1 ,0}, and the corresponding conjugate stress variables σtrial n+1 and qtrial n+1 . The constrained character of the problem arises from the restriction on the last component (say, component nx)ofx, that is, [x]nx:= ∆γ≥0,(2.2) for the plastic multiplier ∆γ. The presence of this constraint requires then the use of the line search scheme described in Section I.2 of Appendix I for unilaterally constrained problems. We denote the final algorithm by the primal closest-point projection method (or, in short, the primal–CPPM). Box 2.1 includes a summary of this scheme for a typical plastic corrector step. The resulting algorithm is then very similar to the original Newton–CPPM, but with the added considerations associated to the line search and constraint detection. In particular, the condition (I.15) in Appendix I for the activation of the constraint (2.2) reads m(k) σn+1 ·(εe,(k) n+1 −εe,trial n+1 )−m(k) qn+1 ·(α(k) n+1 −αtrial n+1 )>0,with ∆γ(k)=0,(2.3) in a generic iteration (k). If the condition (2.3) is satisfied, the modified Jacobian (I.17) presented in Box 2.1 is to be used. We observe that both the trial step and the final solution of the plastic corrector step of interest herein (that is, with ∆γ>0 strictly) do not satisfy this condition, and hence the activation of the constraint by the algorithm near then. To gain a better understanding of the relation (2.3), Figure 2.1 depicts a graphical interpretation of this condition for von Mises and Tresca yield conditions with associated perfect plastic flow and constant elasticities. We note, as an aside, that the constraint (2.2) has not become activated in any of the examples presented in Section 4, involving more complex yield conditions. We conclude that, even though the condition (2.3) needs to be considered from a theoretical point of view, it does not lead to added computational cost in practice. We refer to this section for further details and discussion. As discussed in Appendix I, the global convergence of the algorithm requires the Jacobian to be non-singular and continuously differentiable, with the special considerations Closest-Point Projection Algorithms in Elastoplasticity 7 BOX 2.1. Implementation of the primal closest-point projection method (primal–CPPM). 1. Input data: {εe,trial n+1 ,αtrial n+1 }, with the corresponding {σtrial n+1 ,qtrial n+1 } and ftrial n+1 >0. 2. Initialization: set k=0, x(0) :=      εe,trial n+1 −αtrial n+1 0     ,and r(0) :=      0 0 ftrial n+1     . 3. Compute the Jacobian: J(k):= I+∆γ(k)∇m(k)˜ G(k)m(k) ∇f(k)T˜ G(k)0 where all the quantities in this expression are evaluated at x(k)= {εe(k) n+1,−α(k) n+1,∆γ(k)}Taccording to their definitions in Box 1.1. 4. Compute the direction of advance: aux := 0 IF ∆γ(k)=0THEN aux := m(k) σn+1 ·(εe,(k) n+1 −εe,trial n+1 )− m(k) qn+1 ·(α(k) n+1 −αtrial n+1 ) IF aux ≤0THEN flag := .true. d(k):= −(J(k))−1r(k) ELSE flag := .false. [D(k)]ij :=   (J(k))−1(J(k))−Tij i, j =nxor i=j=nx 0i, j =nxand i=j d(k):= −D(k)(J(k))Tr(k) ENDIF 5. Apply the line search scheme of Box I.2 to the function r(x) by (2.1). INPUT: x(k),r(k),d(k)and flag. OUTPUT: x(k+1) and r(k+1),aswellas(σ,q)(k+1) n+1 from the computation of r(k+1), and other auxiliary quantities used in Item 3. for the Jacobian J(k+1). 6. Check convergence →set (εe,α,σ,q)n+1 =(εe,α,σ,q)(k+1) n+1 and EXIT. 7. Set k←k+1andGO TO 3. A. P´erez-Foguet & F. Armero 8 sn+1 von Mises n+1 s trial n+1 s trial f ( s ) 0 n+1 s (k) n+1 s (k) sn+1 Tresca f ( s ) 0 (k) = 0 n+1 m (k) n+1 m (k) region where the constraint is not activated n+1n+1 m (k)· ( s (k) - s trial ) < 0 n+1n+1n+1 m (k)· ( s (k) - s trial ) < 0 n+1 FIGURE 2.1. Graphical interpretation of the condition (2.3)for the a) von Mises and b) Tresca yield conditions, both with associated perfect plastic flow (m=∇f) and constant elasticities. The region of the deviatoric stress plane (s:= σ−3 i=1[σ]i) where the constraint ∆γ≥0 is not activated, leading to the modified Jacobian in the line search scheme, is shown in gray. on the residually-based descent function presented in Remark I.1 due to the constrained character of the problem. A non-singular Jacobian is implied by the usual assumptions of an associative convex problem (i.e., m=∇fwith ∇2fpositive semi-definite) and strict convexity of the elastic relations (i.e. ˜ Gpositive definite), as considered in Part I of this work. Under these assumptions, the equations (2.1) correspond to the first order necessary and sufficient conditions of the convex primal variational problem considered therein (Proposition 3.1, to be specific). The required smoothness follows from the regularity of the material model, that is, the yield function, and elastic and hardening laws. However, the previous ultimate requirements on the Jacobian can be satisfied in more general situations. As a matter of fact, we consider in Section 4 an example involving strain softening. The performance of the algorithm is also assessed in that section for a case involving a non-differentiable residual due to the consideration of a yield function leading to a non-differentiable flow vector. In both cases, the primal–CPPM shows a dramatic improvement in its global convergence properties over the original Newton–CPPM, while exhibiting locally the desired asymptotic quadratic rate of convergence. Remarks 2.1. 1. Globally convergent primal closest-point projection algorithms can be easily devised for the viscoplastic problem given its unconstrained character; see Section 4 of Part I. Closest-Point Projection Algorithms in Elastoplasticity 15 / c f ( ) (1) (k) = f f ( ) (1) (k) f trial n+1 cf trial n+1 f trial n+1 f c c c = 0c > 0 FIGURE 3.1. Graphical interpretation of the upper level scheme of the dual–CPPM for the non-augmented (c= 0) and augmented (c>0) formulations. for an typical iteration (k) in terms again of the Maucalay brackets ·. The derivative ¯ f in (3.4) is obtained in closed-form using the equations (3.1) and (3.2) as ¯ f=−∇f·˜ G−1+∆γ∇m−1 m.(3.5) At the boundary ∆γ= 0 of the admissible set, this derivative reduces to the value ¯ f∆γ=0 =−∇f·˜ Gm<0 strictly ,(3.6) given the fundamental assumption of the underlying continuum model (see Box 1.1), even in the general non-convex (e.g. with strain softening) and non-associated plastic flow (m=∇f). This important result implies that even though the iteration process (3.4) may lead to the value ∆γ(k)= 0 at some iteration (k), the constraint itself is not activated, that is, the iteration direction (3.4)2leads naturally to ∆γ(k+1) >0 and no modified advancing rules like (I.17) in Appendix I are needed. This situation, together with the aforementioned unconstrained character of the problem in the lower level, avoids altogether the need for the use of complex line search rules as required in the constrained primal formulation presented in Section 2.1. The global convergence of the upper level algorithm defined by the iterative process (3.4) is then assured if the derivative in (3.5) is continuous and does not vanish. The smoothness condition is again implied by the proper smoothness of the material model, that is, the yield surface and elastic and hardening laws. The non-zero condition on the derivative is satisfied under the usual assumptions of convexity and associated plastic A. P´erez-Foguet & F. Armero 16 n+1 , qn+1 f (1) > 0 f n +1 = 0 f trial> 0 n +1 trial, q trial n+1 n+1 (1), q (1) n+1 n+1 Elastic domain f < 0 n+1 (k), q (k) n+1 n+1 f (k) > 0 n+1 n+1 orthogonal projection in the -1 metric FIGURE 3.2. Graphical interpretation of both the dual–CPPM and the augmented–dual–CPPM. The stress variables {σ(k) n+1,q(k) n+1}at each iteration (k) are the closest-point projection onto the current iterate of the yield function ¯ f(∆γ(k))(or ¯ fc(∆λ(k))fortheaugmented algorithm) for the plastic multiplier ∆γ(k)(or ∆λ(k)) obtained by the upper level scheme. A similar figure can be found in the literature but erroneously associated with the original Newton–CPPM of Box 1.1. flow. Under these assumptions, equation (3.5) reveals that the function ¯ f(∆γ) is strictly monotonically decreasing ( ¯ f<0) as obtained in Part I of this work. Figure 3.1 illustrates these results graphically. Under these conditions, the line search scheme considered in this upper level scheme did not need to be activated in the numerical examples reported in Section 5, but its consideration is needed to assure the global convergence of the iterative process in general situations. Therefore, we conclude that the proposed two level scheme based on the dual form of the closest-point projection equations is globally convergent under these assumptions. The numerical examples presented in Section 5 have confirmed this property. Figure 3.2 depicts graphically the idea behind the proposed dual–CPPM for the associated case given by the normality rule m=∇f, and constant elasticities and linear hardening (i.e. ˜ G= constant). Equations (3.2) imply in this case that in a given iteration of the upper level scheme, that is, for each value ∆γ(k), the vector Σ(k) n+1 −Σtrial n+1 in stress space (for Σ={σ,q}) is indeed normal in the metric ˜ G−1to the current value of the yield surface f(k). This property is illustrated in Figure 3.2, showing each iteration value of the stress variables as the closestpoint projection on the elastic domain f(k) n+1 =¯ f(∆γ(k))≤0 for the current iterate ∆γ(k) in the upper level algorithm depicted in Figure 3.1. We note that a similar figure can Closest-Point Projection Algorithms in Elastoplasticity 17 be found in Simo & Hughes [1998], among others, and it is erroneously associated with the original Newton–CPPM. The developments presented in this section clearly identify the approach depicted in this figure with the newly proposed algorithm based on the dual formulation of the governing equations. Remark 3.1. The dual–CPPM described in this section extends to the viscoplastic problem, following the developments of the dual formulation of this problem presented in Proposition 4.3 of Part I. That is, simply replace the scalar equation (3.1) by the relation g(¯ f(∆γ)) −η∆γ=0,(3.7) for a general viscoplastic function g(·) and viscoplastic parameter η≥0. We observe that this is the only needed modification, since the problem (3.2) is not changed, with the lower level algorithm of Box 3.2 applying entirely to the viscoplastic problem as well. Similarly, the same considerations regarding the non-activation of the constraint ∆γ= 0 and, indeed, the same global convergence properties apply under the usual monotonicity assumption of the viscoplastic function g(·). Additional details are omitted. 3.2. The augmented dual closest-point projection method Even though the constrained character of the upper level problem in the dual scheme considered in the previous section was shown not to affect the algorithm itself, the consideration of augmented Lagrangian extensions avoiding altogether this constraint leads to algorithms improving on their numerical performance, as illustrated in the numerical examples presented in Section 5. Following the developments of the dual augmented Lagrangian formulation summarized in Proposition 5.2 of Part I, the plastic consistency function ¯ f(·) in (3.1) is replaced by ¯ fc(∆λ):=f(σc(∆λ),qc(∆λ)) = 0 ,(3.8) with {σc(∆λ),qc(∆λ)}defined by the augmented lower level problem ˜rc(˜x)=0forthe augmented residual ˜rc(˜x)=εe n+1 −εe,trial n+1 +∆λ+cf(σn+1,qn+1)mσn+1 −αn+1 +αtrial n+1 +∆λ+cf(σn+1,qn+1)mqn+1 ,(3.9) for a penalty parameter c≥0 and in terms, again, of the primary unknowns ˜x:= {εe n+1,−αn+1}T. Following the same notation employed in the primal augmented Lagrangian formulations of Section 2.2, we denote the scalar field appearing in these equations by ∆λemphasizing its entirely unconstrained character in front of ∆γ≥0. The final plastic multiplier ∆γis recovered by the simple relation ∆γ=∆λ,(3.10) A. P´erez-Foguet & F. Armero 18 BOX 3.3. Implementation of the upper level of the augmented–dual– CPPM, involving a Newton iterative process in ∆λ. 1. Input data: {εe,trial n+1 ,αtrial n+1 }, with the corresponding {σtrial n+1 ,qtrial n+1 } and ftrial n+1 >0. 2. Initialize: set k=0,(εe,α,σ,q)(0) n+1 =(εe,α,σ,q)trial n+1 ,∆λ(0) =0 and ¯ f(0) c=ftrial n+1 . 3. Compute update direction: δ(∆λ)(k):= ¯ f(k)/∇f(k)·˜ G(k)(I+H(k)˜ G(k))−1m(k), where s(k):= 0,for ∆λ(k)+cf(k)<0, 1,for ∆λ(k)+cf(k)>0, H(k):= ∆λ(k)+cf(k)∇m(k)+cs (k)m(k)⊗∇f(k). 5. Apply the line search scheme of Box I.1 with the scalar residual and descent function rc∆λ(∆λ):= ¯ fc(∆λ),and Mc∆λ(∆λ)=1 2¯ fc(∆λ)2 for x(k)←∆λ(k), and where ¯ fc(∆λ)=f(σc(∆λ),qc(∆λ)) for the values {σc(∆λ),qc(∆λ)}computed by the lower level algorithm in Box 3.4 for a fixed ∆λ. INPUT: ∆λ(k),¯ f(k) cand δ(∆λ)(k). OUTPUT: ∆λ(k+1),¯ f(k+1) c,aswellas(εe,α,σ,q)(k+1) n+1 and other auxiliary quantities used in the update formula in 3., from the evaluation of ¯ f(k+1) cthrough the algorithm in Box 3.4 (lower level). 5. Check convergence →set (εe,α,σ,q)n+1 =(εe,α,σ,q)(k+1) n+1 and EXIT. 6. Set k←k+1andGO TO 3. Closest-Point Projection Algorithms in Elastoplasticity 19 BOX 3.4. Implementation of the lower level of the augmented–dual– CPPM. 1. Input data: the fixed values εe,trial n+1 ,αtrial n+1 and ∆λ(k+1), and the initial values (εe,α)(k) n+1,with(σ,q)(k) n+1. 2. Initialize: set i=0,˜x(0) := εe,(k) n+1 ,−α(k) n+1Tand ˜rc(0) := ˜rc(˜x(0)) by (3.9). 3. Compute the Jacobian: ˜ Jc(i):= I+∆λ(k+1) +cf (i)∇m(i)+cs (i)m(i)⊗∇f(i)˜ G(i), where s(i):= 0,∆λ(k+1) +cf(i)<0, 1,∆λ(k+1) +cf(i)>0. 4. Compute the direction of advance: ˜ dc(i):= −(˜ Jc(i))−1˜rc(i). 5. Apply the line search scheme of Box I.1 based on the function ˜rc(x). INPUT: ˜x(i),˜rc(i)and ˜ dc(i). OUTPUT: ˜x(i+1),and˜rc(i+1),aswellas(σ,q)(i+1) from the computation of ˜rc(i+1), and other auxiliary quantities used in Item 3. for the Jacobian ˜ Jc(i). 6. Check convergence →set (εe,α,σ,q)(k+1) n+1 =(εe,α,σ,q)(i+1) and EXIT. 7. Set i←i+1andGO TO 3. corresponding actually to ∆λ>0 in the plastic corrector step of interest. We observe that the original dual problem (3.1)–(3.2) is recovered for c= 0. The solution of the upper and lower problems defined by (3.8) and (3.9), respectively, is approached again with Newton schemes in combination with the unconstrained line search technique presented in Section I.1. The final scheme is summarized in Boxes 3.3 and 3.4 for the upper and lower level algorithms, respectively. The fully unconstrained character of the problem in the upper level of the algorithm A. P´erez-Foguet & F. Armero 20 leads to the direct consideration of the line search scheme of Box I.1 without the need of the added imposition of the non-negative constraint on ∆γby (3.4)1for the previous constrained dual formulation. The derivative of the augmented plastic consistency function (3.8) can also be obtained in closed-form after using equations (3.9) of the lower level augmented problem as ¯ f c=−∇f·˜ G−1+H−1 m.(3.11) where H:= ∆λ+cf∇m+csm⊗∇f, (3.12) for s:= 0,for ∆λ+cf < 0, 1,for ∆λ+cf > 0.(3.13) Even though no problems have been observed in the actual numerical simulations due to the discontinuity of sat ∆λ+cf = 0 (in fact, values ∆λ+cf < 0 are never reached), we assign the value s= 0 at this point. We can conclude from these results the the same properties of the augmented plastic consistency function ¯ fc(∆λ) as its counterpart ¯ f(∆γ) in the original constrained dual formulation; additional details are omitted. In particular, its strictly monotonically decreasing character (i.e., ¯ fc<0) is concluded under the usual convexity and associativity assumptions. The same globally convergent character of the proposed scheme is concluded in these cases. We also observe that the augmented–dual–CPPM developed in this section has the same graphical interpretation depicted in Figure 3.2 for the non-augmented dual algorithms (i.e c= 0). The upper level algorithm involves, however, a regularized consistency function, a situation that has been depicted in Figure 3.1 for c>0. As argued in the end of Remark 3.2 below and verified in the numerical examples presented in Section 5, this simple modification leads to some improvement in the computational performance of the dual algorithms considered in this section. Remark 3.2. We observe that by considering the initial estimate of the upper level algorithm by ∆λ(0) = 0, the first solution of the lower level problem (3.9) obtained by the scheme in Box 3.4 corresponds to the solution of a viscoplastic problem defined by the viscoplastic parameter c=1/ηand linear viscoplastic model g(f)=<f>. In fact, this problem corresponds to the primal formulation of the viscoplastic problem considered in Remark 2.1 above for this case. The consideration of more general viscoplastic models is easily accomplished through the use of augmented formulations based on the corresponding regularization function g(·) as presented in Remark 5.1.2 (that is, replace <·>by g(·)in the relations of Box 3.4). In this viscoplastic case, the equation (3.10) giving the plastic multiplier ∆γis to be replaced by the viscous relation ∆γ=1 ηg(¯ fc(0)) .(3.14) Closest-Point Projection Algorithms in Elastoplasticity 21 Therefore, the choice c=1/ηleads the final solution of the viscoplastic problem in a single iteration of the upper level algorithm. If the viscoplastic problem is understood as the penalty regularization of the rate-independent problem as η→0, the improvement in the computational performance gained by the consideration of a large regularization parameter cis to be expected when solving the limit rate-independent problem. The usual considerations on the well-conditioning of the resulting problem is to be weighed in the argument, as it is investigated through the numerical examples presented in Section 5below. 4. Numerical Assessment: Primal Algorithms We assess in this section the numerical performance of the new primal algorithms proposed in Section 2. More specifically, it is our interest to evaluate the global convergence properties of these schemes. To this purpose, we study in Section 4.1 the convergence properties of the proposed schemes for a single increment (from, say, tnto tn+1) for different imposed strain increments. Section 4.2 evaluates the overall performance of the schemes in a given imposed path of the deformation gradient, involving different numbers of time increments. Remarks 4.1. 1. The results presented below are expressed in terms of the invariants p∗=−I1(∗) 3,q ∗=3J2(∗)andθ∗=1 3arccos 3√3J3(∗) 2J2(∗)3/2,(4.1) of a generic second order tensor ∗,whereI1(∗), J2(∗)andJ3(∗) denote the first invariant of ∗and the second and third invariants of the deviatoric part of ∗, [dev(∗)]i= [∗]i+p∗, respectively, I1(∗)= 3  i=1 [∗]i,J 2(∗)=1 2 3  i=1 ([dev(∗)]i)2and J3(∗)= 3  i=1 [dev(∗)]i,(4.2) in terms of the three principal points [∗]ifor i=1,2,3. The function θ∗corresponds to the so-called Lode’s angle and varies from 0◦to 60◦. These expressions are used for the Kirchhoff principal stresses ∗=σand the logarithmic principal elastic strains ∗=εe; see Section 2.3 of Part I of this work. 2. All the contour plots presented in this and subsequent sections (see e.g. Figure 4.1) have been obtained with a uniform grid of 80 ×80 sampling values, independently of the ranges of the variables depicted in both axis. A. P´erez-Foguet & F. Armero 22 3. In all the examples presented in this paper, the convergence of the algorithm is detected through the expression e(k) ˆx=||ˆx(k+1) −ˆx(k)||∞ ||ˆx(k+1)||∞≤TOLˆx(4.3) measuring the relative error in iteration (k+ 1) in the maximum norm || · ||∞= maxi|[·]i|, for each component [ ·]idespite they may have different dimensions, in general. Here ˆxrefers the driving variable employed in the particular algorithm under consideration (that is, xgiven by (2.1) in the primal algorithms, ˜xgiven by (3.2) in the lower level of the dual algorithms, and ∆γor ∆λin the upper level of the dual algorithms). The tolerance value of TOLˆx=10 −12 has been employed in all cases. 4.1. Single increment tests We evaluate in this section the convergence properties of the iterative schemes under investigation during a plastic corrector step for an imposed trial elastic state. The initial state is assumed elastic (i.e., vanishing initial values of the plastic internal variables) with vanishing strain. Deviatoric and pressure dependent models are considered in Sections 4.1.1 and 4.1.2, respectively, both in combination of perfect plasticity and Hencky’s hyperelastic law, resulting in a set of constant elasticities in the logarithmic principal elastic strains εe; see Appendix II. The imposed trial state is then characterized in the deviatoric Π-plane by the corresponding Lode angle θεe,trial (= θσtrial ) and the radial measure qεe,trial (= qσtrial /2µfor the shear modulus µ). The pressure parameter pεe,trial (= pσtrial /3κfor the bulk modulus κ) defines completely the imposed trial state in the pressure-dependent models considered in Section 4.1.2. 4.1.1. Deviatoric plastic models We consider first the von Mises–Tresca type yield surfaces described in Section II.2 of Appendix II, defining in terms of the material parameter ma family of deviatoric yield surfaces exhibiting different levels of curvature, as it is of the interest in this work. In particular, the choice m= 0 corresponds to the smooth von Mises circular cylinder, with higher values of mleading to increasing of the curvature of this yield surface at θσ=0,30 (the Tresca non-smooth yield surface is recovered for m→∞). Only the two deviatoric parameters qεe,trial and θεe,trial are needed in this case to define the trial state. Because of the symmetry of the resulting yield functions (see Figure II.1 in Appendix II), Lode angles between 0◦and 30◦need only to be considered. The values of κ= 164.206 kN/mm2and µ=80.1938 kN/mm2are taken for the bulk and shear modulus, respectively, in Hencky’s law (II.2)-(II.3), with the constant yield stress σyo=0.45kN/mm 2. Figure 4.1 compares the performance of the Newton–CPPM and the primal–CPPM in this case. The number of iterations needed for convergence of both schemes for several Closest-Point Projection Algorithms in Elastoplasticity 23 m =5 m =10 m =20 primal{CPPM primal{CPPM primal{CPPM Newton{CPPM Newton{CPPM Newton{CPPM  " e;trial 0 Æ 30 Æ  " e;trial 0 Æ 30 Æ  " e;trial 0 Æ 30 Æ 2 q " e;trial = y o 40 2 q " e;trial = y o 4 0 1 3 5 7 9 11 13 15 17 19 NC FIGURE 4.1. Single increment tests: deviatoric plastic models. Number of iterations needed for convergence by the Newton– CPPM and the primal–CPPM for three different von Mises–Tresca type yield surfaces (m= 5, 10 and 20). “NC” denotes no convergence after more than 100 iterations; extended tests in these regions lead to no convergence after several hundreds of iterations. elastic trial states {qεe,trial ,θ εe,trial }and different values of the material parameter mare depicted in this figure. The regions where no convergence is detected after more than 100 iterations are denoted by “NC”; extended simulations show no convergence after several hundreds of iterations in these regions. The value of qεe,trial varies in the range [0,4· σyo/(2µ)], that is, with corresponding trial stress states up to four times the yield limit σyo. The results depicted in Figure 4.1 reveal that no convergence is detected with the A. P´erez-Foguet & F. Armero 24 m =20 Newton{CPPM primal{CPPM ( j max =1) primal{CPPM  " e;trial 0 Æ 30 Æ  " e;trial 0 Æ 30 Æ  " e;trial 0 Æ 30 Æ 2 q " e;trial = y o 10 0 1 4 8 12 16 20 24 28 32 NC FIGURE 4.2. Single increment tests: deviatoric plastic models (m= 20). Number of iterations needed for convergence by the Newton– CPPM and the primal–CPPM (with a maximum of one curve fitting (jmax = 1), and without a limit on the number of curve fittings) for large excursions outside the elastic domain. standard Newton–CPPM when the solution is close to θεe,trial =0 ◦(and 60◦). As noted above, these values correspond to points where the curvature of the yield surface is higher. Moreover, as mincreases the size of the region of no convergence increases. In fact, the Newton–CPPM becomes basically useless even for moderate values of m. With the primal–CPPM these regions of no convergence are avoided. Convergence is attained for any trial state, with less than 17 iterations everywhere. This improvement is observed mainly in the regions of no convergence of the original Newton–CPPM wherever the Newton– CPPM converges, the primal–CPPM does not reduce the number of iterations sensibly. Note also that for θεe,trial exactly equal to 0◦(and 60◦) the convergence is achieved always in a reduced number of iterations. This very special situation is due to the fact that the gradient of the yield surface has always the same direction at these Lode angles. Both schemes are compared again in Figure 4.2, but with the added consideration of the maximum number of curve fittings in the line search scheme employed in the primal– CPPM, that is, the jmax parameter in Boxes I.1 and I.2 of Appendix I. Larger excursions outside the elastic domain, in the range [0,10 ·σyo/(2µ)] for the radial measure qεe,trial , are also considered. We observe that the primal–CPPM without a limit in the number of curve fittings per iteration avoids completely the regions of no convergence of the original Newton–CPPM. The maximum number of curve fittings needed is only two, though. The small region of no convergence remaining far from the yield surface with a limit of just one Closest-Point Projection Algorithms in Elastoplasticity 31 which involve tensile pressures with this yield surface. In contrast, the primal–CPPM leads to convergence for all the trial states. Still, the overall behavior of the primal–CPPM in the region behind the apex is worse than in the examples of Figure 4.6. There is a high variability in the distribution of the number of iterations, which for some trial states are a relatively large, although less than 50 iterations are needed anywhere. The large number of iterations is again related with to the nondifferentiability of the flow vector, with the associated lack of continuity of the Jacobian at θσ=5 ◦and θσ=55 ◦. The high high sensitivity in the initial trial state can be traced back to the closeness of the iterative process to these points for low values of qσ.The results of Figure 4.7 has been obtained with a maximum number of curve fittings fixed to three. For some trial states, this parameter influences in the convergence results, although the overall behaviour is the same To illustrate better the effect of the non-differentiability of the residual vector and the observed dependence on different algorithmic parameters in the resulting highly sensitive cases, we consider the specific trial state defined by pεe,trial =−3·CMC cot(φ)/(3κ), qεe,trial =2.25 ·√3CMC /(2µ)andθεe,trial =30 ◦in this region of high variability. The primal–CPPM takes 49 iterations with a maximum of three curve fittings (jmax =3),as shown in Figure 4.7. This is the maximum number of iterations observed in this figure. The evolution of the relative error, the yield function, the number of curve fittings, the value of the line search parameter and the Lode angle of the stresses during the iterative process is depicted in Figure 4.8 for a maximum number of curve fittings jmax =5,and 10. In both cases, more curve fittings than this fixed maximum number are needed in a few iterations. These iterations are indicated with a black circle in the plots depicting the required number of curve fittings. Note that both iterative processes (for jmax = 5 and 10) coincide until this point. During these early stages of the iterative process, the iterations are “captured” at the Lode angle θσ=55 ◦, where the Jacobian of the residual is not continuous. The number of curve fittings increases and the iterative increment becomes very small (bounded from below by the backtracking parameter η=0.1 and maximum number of curve fittings; see Box I.1) while the relative value of the yield function remains basically constant and the relative error based on the norm of the increment of the unknowns diminishes. This fact is directly related with the value of the line search parameter. The test with jmax = 10 allows the additional number of curve fittings after this point. After the maximum number of curve fittings jmax is reached, the behavior of the iterative scheme changes dramatically. One curve fitting is needed during some of the subsequent iterations and finally a quadratic rate of convergence is achieved. The iterative process has “escaped” from the discontinuity of the Jacobian. These considerations, together with the discussion presented with Figure 4.2 for smooth cases, allow to conclude that a moderate value of the parameter jmax (≈3−5 seems appropriate) should be used for the maximum number of curve fittings. A. P´erez-Foguet & F. Armero 32 0 2 4 6 8 10 1 1121314151617181 0 15 30 45 60 1 1121314151617181 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1 1121314151617181 1,E-11 1,E-08 1,E-05 1,E-02 1,E+01 1 1121314151617181 0 2 4 6 8 10 1 1121314151617181 0 15 30 45 60 1 1121314151617181 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1 1121314151617181 1,E-11 1,E-08 1,E-05 1,E-02 1,E+01 1 1121314151617181 (k) ex | f / f trial| n+1 ex | f / f trial| n+1 (k) FIGURE 4.8. Single increment tests: pressure-dependent plastic models. Influence of the non-differentiability of the flow vector in the RHMC model for the test given by pεe,trial =−3·CMC cot(φ)/(3κ), qεe,trial =2.25 ·√3CMC/(2µ)andθεe,trial =30 ◦. Evolution of the the relative error e(k) x, the yield function f(k)/ftrial n+1 ,thenumberof curve fittings, the line search parameter α(k)and the stress Lode angle θσ(k)for a maximum of curve fittings jmax = 5 and 10. Closest-Point Projection Algorithms in Elastoplasticity 33 c =0 : 001 =C MC c =0 : 03 =C MC p " e;trial 3  C MC cot(  )  3 0 p " e;trial 3  C MC cot(  )  3 0 q " e;trial 2  p 3 C MC 5 0 -25 -15 -7 -1 1 7 14 21 FIGURE 4.9. Single increment tests: pressure-dependent plastic models. Increment of the number of iterations for the primal– CPPM with different values of the penalty parameter c, with respect to the non-augmented formulation in Figure 4.7 for θεe,trial =30 ◦.. To conclude this section, we consider again the augmented versions of the algorithms. The increment of the number of iterations of the augmented primal–CPPM with respect to the non–augmented version is depicted in Figure 4.9 (elastic trial steps in the meridian plane θεe,trial =30 ◦). The performance of the scheme improves in average with small values of the penalty parameter cwhen compared with the non-augmented formulation. For larger values of c, although a high reduction of the number of iterations is found in some regions, a significant increase is found in another ones (note the results obtained with c=0.03/CMC in the region close to the apex). This effect is accentuated as cincreases, and, as this effect occurs in the region closest to the apex, the method becomes less competitive with large values of c. The effects of the non-differentiability of the residual follow the same pattern in this augmented case as discussed in the previous paragraphs for the dual–CPPM. 4.2. Imposed deformation gradient path We evaluate next the performance of the primal algorithms in more general situations involving Ogden hyperelastic models (i.e., non-constant elasticities in the logarithmic strains) as well as strain hardening and softening relations. To this purpose we consider imposed paths of the deformation gradient. More specifically we consider the deformation gradient history F(t)=b(t)1/3diag(a(t),a(t)−r,a(t)r−1) (4.4) for the pseudotime t, function a(t) (such as a(0) = 1 and a(t)>0) that guides the evolution A. P´erez-Foguet & F. Armero 34 TABLE 4.2. Elastic parameters used with the Ogden hyperelastic model. κ= 164.206 kN/mm2ˆµ1=1.4911 ˆα1=1.3000 µg=80.1938 kN/mm2ˆµ2=0.0028 ˆα2=5.0000 N=3 ˆµ3=−0.0237 ˆα3=−2.0000 of the axial stretch in the first principal direction, and the function b(t) controlling the volumetric response (det F=b(t)). The exponent rin (4.4) controls the desired Lode angle of the elastic trial increment (∆εe,trial =εe,trial n+1 −εn), since r=1 21−√3tan(θ∆εe,trial )(4.5) In particular, the value r=0.5 corresponds to pure tension (θ∆εe,trial =0 ◦) and the value r=−1 to pure compression (θ∆εe,trial =60 ◦). 4.2.1. Deviatoric plastic models We consider the deviatoric plastic models defined by the von Mises-Tresca type yield surfaces described in Section II.2 of Appendix II. The Ogden hyperelastic model defined in (II.4) of Appendix II, with material parameters summarized in Table 4.2 (Miehe [1998]), is considered. Similarly, the hardening/softening potential (II.9) is considered for the saturation exponent δ= 20 and with the saturation stress σy∞=0.6kN/mm 2for the examples with strain hardening and σy∞=0.3kN/mm 2in the examples with strain softening. In all cases, including the examples with perfect plasticity, the initial yield stress is σyo=0.45kN/mm 2. The imposed function of a(t)=1+0.2tfor t∈[0,1] ,(4.6) is considered, with b(t)≡1 leading to a purely isochoric deformation. We consider solutions involving 5, 50 and 500 equal time increments. The resulting Kirchhoff stress path in qσis depicted in Figures 4.10, including details in the small range of the imposed axial stretch. We consider the values of θ∆εe,trial =0 ◦,15 ◦and 30◦, leading to the corresponding exponents rby (4.5). Note the smooth transition from elastic to elastoplastic regime in the curves corresponding to θ∆εe,trial =15 ◦. This is produced by the change of the Lode’s angle of the stresses during the evolution of the simulation. The evolution of qσfor a larger range of the axial stretch and θ∆εe,trial =15 ◦is shown in Figure 4.10.b. The evolution of the Lode angle θσof the Kirchhoff stress is depicted in Figure 4.10.c for the softening case. The solid curves of correspond to the problems solved with 500 equal increments and the Closest-Point Projection Algorithms in Elastoplasticity 35 a) 0,3 0,35 0,4 0,45 0,5 1 1,005 1,01 1,015 1,02 0,3 0,35 0,4 0,45 0,5 1 1,005 1,01 1,015 1,02 0,3 0,35 0,4 0,45 0,5 1 1,005 1,01 1,015 1,02 0o 15o 30o 0o 15o 30o 0o 15o 30o q a(t) q q a(t) a(t) b) c) 0,2 0,3 0,4 0,5 0,6 0,7 1 1,05 1,1 1,15 1,2 0 4 8 12 16 1 1,05 1,1 1,15 1,2 q a(t)a(t) FIGURE 4.10. Imposed deformation gradient paths: deviatoric plastic models (m= 20). a) Details of the evolution of the equivalent Kirchhoff stresses qσin the range a(t)∈[1,1.02], for three prescribed paths (θ∆εe,trial =0 ◦,15 ◦and 30◦) and perfect plasticity, hardening and softening. b) Evolution of qσin the full range a(t)∈[1,1.2] for perfect plasticity, hardening and softening models, and θ∆εe,trial =15 ◦. (solid lines for obtained with 500 time increments, with the stars corresponding to the solution obtained with 5 time increments, both obtained with the primal–CPPM). c) Evolution of the Lode angle θσfor the softening case and θ∆εe,trial =15 ◦. stars to the problems solved with just 5 equal increments. The line search scheme has not been activated in the problem with 500 increments recovering then the Newton–CPPM. It has been activated with the simulation involving 5 time increments only. Still, we observe in Figure 4.10.b a good agreement between both solutions, hence implying that the use of very large elastic trial steps does not imply a loss of accuracy. The convergence results at several pseudotimes tare presented in Figure 4.11 for the hardening and the softening examples, including the evolution of the relative errors e(k) x (see equation (4.3)) in x={εe n+1,−αn+1,∆γ}. The results correspond to the paths solved A. P´erez-Foguet & F. Armero 36 Hardening (5 increments) 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1 3 5 7 9 1113151719 t = 0.2 t = 0.4 t = 0.6 t = 0.8 t = 1.0 Softening (5 increments) 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1357911131517 t = 0.2 t = 0.4 t = 0.6 t = 0.8 t = 1.0 ex ex Hardening (50 increments) 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1357911131517 t = 0.2 t = 0.4 t = 0.6 t = 0.8 t = 1.0 Softening (50 increments) 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1357911131517 t = 0.2 t = 0.4 t = 0.6 t = 0.8 t = 1.0 ex ex Hardening (500 increments) 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1357911131517 t = 0.2 t = 0.4 t = 0.6 t = 0.8 t = 1.0 Softening (500 increments) 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1 3 5 7 9 11 13 15 17 t = 0.2 t = 0.4 t = 0.6 t = 0.8 t = 1.0 ex ex FIGURE 4.11. Imposed deformation gradient paths: deviatoric plastic models (m= 20). Convergence results with the primal–CPPM for several pseudotimes t(hardening and softening, and 5, 50 and 500 increments). Closest-Point Projection Algorithms in Elastoplasticity 37 Softening, 5 increments (t = 0.2) 0 1 2 13579111315 Line-search param. Num. curve fittings Hardening, 5 increments (t = 0.2) 0 1 2 13579111315 Line-search param. Num. curve fittings Softening, 50 increments (t = 0.2) 0 1 2 1 3 5 7 9 11 13 15 Line-search param. Num. curve fittings Hardening, 50 increments (t = 0.2) 0 1 2 13579111315 Line-search param. Num. curve fittings FIGURE 4.12. Imposed deformation gradient paths: deviatoric plastic models (m= 20). Evolution of the line search at t=0.2forthe examples with hardening and softening, and 5 and 50 time increments. with 5, 50 and 500 time increments. In all cases the primal–CPPM converges with the rate of convergence being asymptotically quadratic for each time increment. The evolution of the line search parameter and the number of curve fittings is shown in Figure 4.12, for the runs with 5 and 50 increments, pseudotime t=0.2, and both the hardening and softening problem. As occurred in the examples of Section 4.1, the maximum number of curve fittings is two. Therefore, the computational cost added in each iteration by the considered line search scheme is not significant. We can observe that, although the number of iterations per increment increases as the number of increments reduces, the ratio is clearly favorable to the problems with less time increments. A reduction of the number of increments by a factor 10 implies a reduction of the number of accumulated iterations by a factor of 5, at the least. This reduction is directly related with a reduction of the total computational cost. Therefore, the primal– CPPM is computationally very efficient, including cases with non-constant elasticities. A. P´erez-Foguet & F. Armero 38 0 1 2 3 -2 -1 0 1 2 3 s = -1 2 /(3 ) 0 15 30 45 60 -2 -1 0 1 2 3 s = 0 2 /(3 ) s = 1 2 /(3 ) s = 2 2 /(3 ) s = -1 2 /(3 ) s = 0 2 /(3 ) s = 1 2 /(3 ) s = 2 2 /(3 ) FIGURE 4.13. Imposed deformation gradient paths: pressuredependent plastic models. Evolution of the Kirchhoff stresses in the meridian and the deviatoric planes for the paths defined by r=0and s=−1, 0, 1, 2 ·2µ/(3κ)(withtf=1.58,2.24, 1.58 and 1,respectively), obtained by the primal–CPPM with 10 time increments 4.2.2. Pressure-dependent models We conclude the evaluation of the primal algorithms with the consideration of the rounded hyperbolic Mohr-Coulomb model for imposed path of the deformation gradient. We consider again the path defined by the relation (4.4), with the stretch function a(t)=1+6·10−5tt∈[0,t f],(4.7) for a final time tf(see below) and with the volumetric function b(t)=(a(t))3s,(4.8) allowing to prescribe the parameter p∆εe,trial =−s √r2−r+1 q∆εe,trial √3,(4.9) for the elastic trial strain increment ∆εe,trial in the meridian plane. Hencky’s law is considered in this case which results in a linear relation between the stresses and logarithmic strains. Equation (4.9) is equivalent in this case to s=−r2−r+1 2µ 3κ p∆σtrial n+1 q∆σtrial n+1 /√3.(4.10) where ∆σtrial n+1 =σtrial n+1 −σn. The evolution of the Kirchhoff stresses is depicted in Figure 4.13 for the paths defined by r=0ands=to−1, 0, 1, 2 ·2µ/(3κ)(withtf=1.58,2.24, Closest-Point Projection Algorithms in Elastoplasticity 39 2 increments 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1357911 t = 0.5 t = 1.0 1 increment 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1357911 t = 1.0 ex ex 10 increments 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1357911 t = 0.6 t = 0.7 t = 0.8 t = 0.9 t = 1.0 50 increments 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1357911 t = 0.6 t = 0.7 t = 0.8 t = 0.9 t = 1.0 ex ex FIGURE 4.14. Imposed deformation gradient paths: pressuredependent plastic models. Convergence results for several pseudotimes tof the path s=2·2µ/(3κ)(tf= 1) obtained with the primal– CPPM and 1, 2, 10 and 50 increments. 2 increments (t = 1) 0 1 2 1357911 Line-search param. Num. curve fittings 1 increment (t = 1) 0 1 2 1357911 Line-search param. Num. curve fittings 10 increments (t = 1) 0 1 2 1357911 Line-search param. Num. curve fittings FIGURE 4.15. Imposed deformation gradient paths: pressuredependent plastic models. Evolution of the line search parameter at each last increment of the path s=2·2µ/(3κ)(t=tf= 1) obtained with the primal–CPPM and 1, 2 and 10 increments. A. P´erez-Foguet & F. Armero 40 TABLE 4.3. Imposed deformation gradient paths: pressuredependent plastic models. Accuracy of the final results obtained with the primal–CPPM for the four paths presented in Figure 4.14 with respect to the reference solution computed with 1000 increments. Number incr. Method Approx. θεe n+1 θεe n+1 error qεe n+1 error pεe n+1 error 50 Newton–CPPM 57.7◦0◦0.04% 0.01% 10 primal–CPPM 57.7◦0.005◦0.23% 0.06% 2primal–CPPM 57.0◦0.7◦0.87% 0.40% 1primal–CPPM 56.5◦1.2◦1.14% 0.61% 1.58 and 1,respectively). The stress paths in the meridian and the deviatoric planes are shown. The path defined by r=0ands=2·2µ/(3κ) has been obtained with the primal– CPPM and with 1, 2, 10 and 50 equal time increments. The convergence results for several pseudotimes tare shown in Figure 4.14. Except for the solution obtained with 50 time increments, the line search scheme has been activated in these solutions. The characteristic asymptotic quadratic rate of convergence is observed in all the cases. The values of the line search parameter, αk, and the number of quadratic curve fittings has been depicted in Figure 4.15 for the last time increment (t= 1) of the paths solved with 1, 2 and 10 increments. In all the cases the maximum number of curve fittings is two, and the active constraint condition, equation (2.3), has not been activated. Therefore, the cost per iteration of the primal–CPPM is of the same order than that of Newton–CPPM. In fact, the number of accumulated iterations decreases with the number of increments, hence resulting in a reduced computational cost with the use of the primal–CPPM. The accuracy of the computed final solution is assessed in Table 4.3. The same final deformation gradient has been imposed with 1000 increments, with the final result taken as a reference. Although the results are less accurate as less increments are done, the errors are low (of the same order as the results of Table 4.1) even for only one increment, thus justifying the use of very large trial increments in combination of the primal–CPPM. 5. Numerical Assessment: Dual Algorithms We evaluate in this section the numerical properties of the new dual closest-point projection algorithms presented in Section 3. More specifically, it is our goal to verify their globally convergent character, showing locally an asymptotic quadratic rate of convergence, and to evaluate the computational cost given their two-level structure. For brevity in the Closest-Point Projection Algorithms in Elastoplasticity 47 c =0 : 1 =C MC c =0 : 5 =C MC c =0 c =0 : 01 =C MC p " e;trial 3  C MC cot(  )  3 0 p " e;trial 3  C MC cot(  )  3 0 q " e;trial 2  p 3 C MC 5 0 q " e;trial 2  p 3 C MC 50 110 20 30 40 50 60 NC FIGURE 5.6. Single increment tests: pressure-dependent plastic models. Number of iterations needed for convergence by the dual– CPPM (c= 0) and the augmented–dual–CPPM for different values of the penalty parameter c,andθεe,trial =30 ◦. iterations is lower: with the primal–CPPM it is 49, and only 40 with the dual–CPPM. These maximum values are reached with different trial states. The maximum number of curvefittingsissettojmax = 3 as in the example in Section 4.1.2. The same overall behavior discussed in this section is observed on the dependence of this parameter in the numerical performance of the algorithms in these non-smooth cases. 6. Finite element application In this section, the primal–CPPM is applied to solve the necking of a cylindrical bar, as an example of practical application of the algorithms presented in previous sections. This problem is a well–known benchmark test in large–strain solid mechanics, see Simo A. P´erez-Foguet & F. Armero 48 A B C D E F AC = DF = 26.667 mm CD = BE = 6.413 mm FA = 0.99 CD . FIGURE 6.1. Necking of a cylindrical bar. Problem definition and computational mesh. 0 10 20 30 40 50 60 70 80 01234567 Edge displacement [mm] Edge reaction [kN] von Mises Tresca FIGURE 6.2. Necking of a cylindrical bar. Load–displacement curves obtained with the von Mises model and the Tresca model. & Hughes [1998] for the analysis of this test with the von Mises model, and Miehe [1998] and Peri´ c & de Souza Neto [1999] for the solution with Tresca–like models. The test consists in the uniaxial extension of a cylindrical bar with circular cross– section, with a radius of 6.413 mm and 53.334 mm length. A slight geometric imperfection (1% reduction in radius), see Figure 6.1, induces necking in the central part of the bar. An axisymmetric analysis is carried out with the mesh shown in Figure 6.1. Eight–noded quadrilateral elements with 2 ×2 Gauss–points are used. The Hencky’s hyperelastic law and a Tresca–like plastic model (the von Mises–Tresca Closest-Point Projection Algorithms in Elastoplasticity 49 2 4 6 8 10 12 14 16 18 0 ZOOM ZOOM d = 7 mmd = 5.6 mmd = 2.8 mm d = 4.2 mm FIGURE 6.3. Necking of a cylindrical bar with the Tresca model and the primal–CPPM. Distribution of number of iterations at Gauss-point level for four load increments, top displacements equal to 2.8, 4., 5.6 and 7 mm. model with a shape parameter m= 20, see Section II.2), are used here. The load– displacement curve is depicted in Figure 6.2. The results obtained with the von Mises model are included in the same figure for comparative purposes. Both results agree with those presented by Miehe [1998] and Peri´ c & de Souza Neto [1999]. The load increments used in the simulations are marked on the curves. Load increments of reduced size have been needed in the middle and the end of the test to ensure convergence of the problem with the Tresca model. With the primal–CPPM the convergence at Gauss–point level A. P´erez-Foguet & F. Armero 50 does not limit the size of the load increments, see Figure 6.3 with the distribution of number of iterations at Gauss-point level for four load increments. This is in contrast with the Newton–CPPM, where the size of the load increments is restricted by the poor convergence properties of the algorithm at Gauss–point level, and, therefore, the more demanding closest–point projection problem guides the incremental–iterative solution of the overall finite element problem. 7. Conclusions We have presented in this paper two new families of algorithms for the solution of the closes-point projection equations in elastoplasticity, referred to as primal and dual algorithms. Augmented Lagrangian extensions have been considered in both cases, as well as extensions to the viscoplastic problem. The proposed algorithms have been evaluated in the context of multiplicative models of finite strain isotropic plasticity, based on general models of finite hyperelasticity (Hencky’s law and Ogden models), including examples with perfect plasticity, strain hardening and strain softening. This general setting has been considered without requiring the convexity assumptions necessary for the rigorous theoretical characterization of the variational structure developed in Part I of this work, and which motivated all these developments. Plastic models based on deviatoric and pressure-dependent yield surfaces have been considered, the latter leading in addition to non-differentiable flow vectors and thus allowing the evaluation of the numerical schemes under these conditions. The primal–CPPM consists of the application of a Newton scheme with the appropriate line search scheme applied directly to the primal equations of the problem in the elastic strains, strain-like internal variables and plastic multiplier (that is, the Euler-Lagrange equations in the variational context provided by the assumptions of convexity and associativity). The proposed line search scheme accounts for the constraint of positive plastic multiplier and, up to the existence of spurious accumulation points shown not to appear in common materials models (see Remark I.1.1 in Appendix I for alternative modifications), lead in practice to the globally convergent character of the primal–CPPM. The corresponding augmented extension, the augmented–primal–CPPM, has been proposed avoiding this difficulty. In contrast, the dual algorithms split the closest-point projection equations in two, leading to two-level algorithms: the upper level consisting of the enforcement of the plastic consistency condition through an iteration in the plastic multiplier, with the lower level consisting of the solution of the closest-point projection equations for a fixed plastic multiplier. Both levels are solved with a Newton scheme in combination of a line search scheme, assuring the asymptotic quadratic rate of convergence of both iteration processes. In particular, the dual–CPPM, involving a constrained upper level problem because again of the Closest-Point Projection Algorithms in Elastoplasticity 51 non-negative character of the plastic multiplier, has been shown not to require the activation of special line search schemes, since this constraint cannot be activated. The resulting algorithms have been rigorously proven to be globally convergent under the usual assumptions of convexity. Augmented extensions, referred to as the augmented–dual–CPPM, has still been presented to improve in their numerical performance as summarized below. Based on the numerical results presented in Sections 4, 5 and 6, we conclude that all the newly proposed algorithms show a dramatic improvement over the standard Newton– CPPM. The large regions of no convergence observed with this scheme are avoided entirely, illustrating the improved global convergence properties of the proposed schemes while maintaining locally the desired asymptotic quadratic rate of convergence. Focusing first on the primal–CPPM, we observe that this improvement comes at no extra computational cost for all practical purposes. In fact, the scheme reduces to the Newton–CPPM when the latter does not exhibit difficulties in the convergence. The dual–CPPM has shown to be computationally competitive when the primal–CPPM requires a large number of iterations to converge, namely, the regions where the Newton–CPPM exhibits no convergence. The augmented versions of the algorithms have improved properties with respect to the non–augmented ones. They do not need to enforce explicitly the constraint of nonnegative plastic multipliers, leading to unconstrained problems easily to treat analytically and avoiding altogether complex line search schemes in practice. Moreover, they have been shown to reduce the computational cost of the original algorithms (that is, a reduced number of iterations) when the regularization parameter is chosen in a given range. A dimensionless expression of this range has been identified in the representative numerical simulations presented in this paper, and it covers approximately one order of magnitude. This improvement in computational cost is more significant for the dual–CPPM, resulting in the so-called augmented–dual–CPPM, even though it does not lead to a fully competitive scheme in regions where the primal–CPPM, or even the Newton–CPPM, show no difficulties, as noted above. In conclusion, the local integration of the elastoplastic models, involving commonly used yield surfaces with high curvature and/or complex elastic and hardening/softening laws, lead to highly nonlinear equations where currently available techniques show clear difficulties, even no convergence unless small load increments are considered. Since this limitation, even at a single quadrature point of a typical finite element implementation, hinders directly the solution of the global mechanical boundary-value problem, it is of the main interest to consider robust globally convergent schemes for the solution of the resulting closest-point projection equations. We believe that the new algorithms presented in this work provide efficient and computationally competitive alternatives in this respect. Acknowledgments: Financial support for this research has been provided by the ONR under contract no. N00014-96-1-0818 and the NSF under contract no. CMS-9703000 with UC Berkeley. A. P´erez-Foguet was supported by the Generalitat de Catalunya (grant A. P´erez-Foguet & F. Armero 52 number 1998 BEAI200042) and the Commission for Cultural, Educational and Scientific Exchange between the US and Spain (program of 1999, number 99258). All this support is gratefully acknowledged. References Armero, F. & P´ erez–Foguet, A. [2000], “On the Formulation of Closest-Point Projection Algorithms for Elastoplasticity. Part I: The Variational Structure, ” submitted for publication on Int. J. Num. Meth. Engr. Abbo, A.J. & Sloan, S.W. [1995], “A smooth hyperbolic approximation to the MohrCoulomb yield criterion”, Computers and Structures, 54, 427-441. Bertsekas, D.P. [1982] Constrained Optimization and Lagrange Multiplier Methods, Athena Scientific, Belmont (reprint of 1996). Bi´ cani´ c, N. & Pearce, C.J. [1996], “Computational aspects of a softening plasticity model for plain concrete”, Mechanics of Cohesive and Frictional Materials 1, 75-94. Dennis & Schnabel [1983] Numerical methods for unconstrained optimization and nonlinear equations, Prentice-Hall, New Jersey (reprint of 1996). Luenberger, D.G. [1989] Linear and Nonlinear Programming, Addison-Wesley, Reading. Miehe, C. [1998], “A formulation of finite elastoplasticity based on dual coand contravariant eigenvector triads normalized with respect to a plastic metric”, Comp. Meth. Appl. Mech. Engr., 159, 223-260. P´ erez–Foguet, A., Rodr´ ıguez–Ferran, A. & Huerta, A. [2000a], “Consistent Tangent Matrices for Substepping Schemes”, Comp. Meth. Appl. Mech. Engr. in press. P´ erez–Foguet, A., Rodr´ ıguez–Ferran, A. & Huerta, A. [2000b], “Numerical differentiation for local and global tangent operators in computational plasticity”, Comp. Meth. Appl. Mech. Engr. 189, 277-296. Peri´ c, D. & de Souza Neto, E. A. [1999], “A new computational model for Tresca plasticity at finite strains with an optimal parametrization in the principal space”, Comp. Meth. Appl. Mech. Engr. 171, 463-489. Shultz, G.A., Schnabel, R.B. & Byrd, R.H. [1985], “A family of trust–region–based algorithms for unconstrained minimization with strong global convergence properties”, SIAM J. Numerical Analysis, 22, 47-67. Simo, J.C. & Hughes, T.J.R. [1998] Computational Inelasticity, Springer, New York. de Souza Neto, E.A., Peri´ c, D. & Owen, D.R.J. [1994], “A Model for Elastoplastic Damage at Finite Strains: Algorithmic Issues and Applications”, Engr. Computation Closest-Point Projection Algorithms in Elastoplasticity 53 11, 257-281. A. P´erez-Foguet & F. Armero 54 Appendix I. Line Search Schemes We summarize in this appendix some basic results on the formulation of solution algorithms for nonlinear algebraic systems of equations as they are used in this paper. The interest herein is the development of globally convergent algorithms for the closest-point projection equations of elastoplasticity. To this purpose, Newton algorithms combined with the proper line search schemes are presented for general unconstrained and unilaterally constrained problems in Sections I.1 and I.2, respectively. Complete details on the classical results presented in this appendix section can be found in Bertsekas [1982], Dennis & Schnabel [1983] and Luenberger [1989]. I.1. Unconstrained problems Consider the general unconstrained algebraic system of equations r(x)=0,(I.1) for the unknown vector x∈Rnx, with nx≥1, and the function r:Rnx→Rnx.The general problem (I.1) corresponds in some of the cases considered in this work to the first order necessary conditions of the unconstrained minimization problem min x∈Rnxϕ(x)(I.2) for a function ϕ(x):Rnx→Rcontinuously differentiable with r(x)=∇ϕ(x) (its gradient). In this context, we consider a general algorithm of the form x(k+1) =x(k)+α(k)d(k)k=0,1,... , (I.3) for a given initial value x(0), the update direction d(k)∈Rnxand line search parameter α(k)∈R. We direct our attention to solution algorithms exhibiting an asymptotically quadratic rate of convergence through the consideration of the Newton update direction d(k)=−J(x(k))−1 r(x(k))forJ(x):=∇r(x)∈Rnx×nx,(I.4) so J=∇2ϕin the case given by (I.2). Enough regularity is tacitally assumed for the definition (I.4) to make sense, as it is the invertibility of the Jacobian matrix. The case of a strictly convex objective function ϕ(x), as it appears in some of the variational principles considered in Part I of this work, assures this last condition, given the positive definite character of J(x) in this case. Our main interest is directed to globally convergent extensions of the pure Newton scheme (α(k)≡1), assuring the convergence of the sequence {x(k)}from any initial trial Closest-Point Projection Algorithms in Elastoplasticity 55 value x(0) to a solution of (I.1). The classical Global Convergence Theorem (see Luenberger [1989], page 187) assures this desirable property for bounded sequences {x(k)} if the algorithm (I.3) defines a continuous mapping from x(k)to x(k+1) (closed mapping in the general context if multiple values x(k+1) exist) and a continuous descent function M:Rnx→Rexists (that is, M(x(k+1))<M(x(k)), “≤”ifx(k)is a solution). These two requirements are satisfied by the algorithm (I.3) if the line search parameter α((k)>0is obtained as the minimization of the merit function ˆ M(α):=M(x(k)+αd(k)),(I.5) defined in terms of the descent function itself, and the update direction d(k)defines a descent direction in the sense that ˆ M(0) = ∇M(x(k))·d(k)≤0.(I.6) Practical line search schemes are obtained by considering an approximate solution of the minimization problem of the merit function (I.5). In this case, the global convergence of the final algorithm is preserved if the so-called Goldstein’s conditions (Luenberger [1989], page 214) M(x(k)+α(k)d(k))≤M(x(k))+βα (k)∇M(x(k))·d(k)(I.7) M(x(k)+α(k)d(k))≥M(x(k))+(1−β)α(k)∇M(x(k))·d(k),(I.8) with β∈(0,1/2), are satisfied. Moreover, the value α(k)= 1 is chosen whenever it verifies equations (I.7) and (I.8), assuring at the same time that the asymptotic quadratic rate of convergence of the original Newton update (I.4) is maintained. The first Goldstein condition (I.7) implies the reduction of the merit function, bounding also the line search parameter α(k)>0 from above. The second Goldstein condition (I.8) assures that this parameter is bounded away from zero. Other equivalent conditions can be found in the literature. A descent function for the problem (I.1) can be constructed in terms of the residual r(x)itselfas M(x):=1 2r(x)·r(x),(I.9) for the Euclidean inner product “·”inRnx. The descent property (I.6) follows for this function and the Newton update (I.4) after noting that ∇M(x(k))·d(k)=−r(x(k))·r(x(k))=−2M(x(k))≤0.(I.10) For the minimization problem (I.2), an alternative descent function is given by the simple choice M(x)=ϕ(x) when this function is strictly convex. Given the generality of (I.10), A. P´erez-Foguet & F. Armero 56 BOX I.1. A line search scheme for unconstrained problems. The parameters η=0.1andβ=10 −4are considered in the numerical simulations presented in this paper. 1. Input data: x(k),r(k)and d(k). 2. Initialize: set j=0,α(k) (0) =1, ˆ M(k):= r(k)·r(k)/2 and ˆ M(k):= −2ˆ M(k). 3. Compute the new unknowns, residuals and merit function: x(k+1) (j):= x(k)+α(k) (j)d(k) r(k+1) (j):= r(x(k+1) (j)) ˆ M(k+1) (j):= r(k+1) (j)·r(k+1) (j)/2 4. Check Goldstein’s condition: IF ˆ M(k+1) (j)≤1−2βα (k) (j)ˆ M(k)THEN set x(k+1) =x(k+1) (j)and EXIT. 5. Check for maximum number of quadratic curve fittings: IF j=jmax THEN notify, set x(k+1) =x(k+1) (j)and EXIT. 6. Compute new value of line search parameter: α(k) (j+1) := MAX     ηα (k) (j),−α(k) (j)2ˆ M(k) 2ˆ M(k+1) (j)−ˆ M(k)−α(k) (j)ˆ M(k)     7. Set j←j+1andGO TO 3. we have chosen the function (I.9) in the developments that follow, even if this convexity property holds. Several iterative techniques can be found in the literature for the construction of line search schemes satisfying the conditions (I.7)-(I.8). In this work we consider a curve fitting technique, consisting of a quadratic fit of the merit function ˆ Mfrom the values ˆ M(k),ˆ M(k) and ˆ M(k) (j)at the iteration value α(k) (j), with the new value α(k) (j+1) given by the minimum of the resulting quadratic curve. This process is repeated for a finite number of iterations (j=0,1,...,j max) until the first Goldstein condition (I.7), which reads in this case ˆ M(α(k) (j))=: ˆ M(k) (j)≤1−2βα (k) (j)ˆ M(k),(I.11) Closest-Point Projection Algorithms in Elastoplasticity 63 q y o 1 Tresca (m= ) m= von Mises (m =1) m =1 m =5 m =10 m =20 von Mises - Tresca 1 23 FIGURE II.1. Trace on the deviatoric plane of the von Mises–Tresca type yield surface for different values of m. purpose, we consider two yield functions, one involving a deviatoric plastic model characteristic of metal plasticity and a pressure-dependent yield surface typical in geomechanics applications. Associated plastic evolutions are considered in both cases. 1. Deviatoric models. We consider the von Mises–Tresca type model defined by the general formula (see Miehe [1998]) f(σ,q)=2m−1 2m √33  i=1 [σ]i−[σ]mod(i+1,3)2m1 2m −&2 3(σyo−q),(II.5) for the initial yield limit σyoand the stress-like internal variable q=−∂αψhfor the considered isotropic strain hardening/softening. The material parameter mdetermines the shape of the yield surface in the stress deviatoric plane, recovering the von Mises circle for m= 1 and the Tresca hexagon for m→∞. The elastic limit at θσ=0 ◦and 60◦is the yield stress value qσ=σyofor all values of m. See Figure II.1 for an illustration. High values of the curvature of the yield surface are obtained for high values of the parameter m. For the value m= 20, the difference between the yield stress at θσ=30 ◦between this yield surface and the Tresca hexagon is less than 2%. The combination of the yield surface (II.5) with the regularized elastic models of the previous section leads to closest-point projection equations in the deviatoric stress plane only (i.e. in s:= σ−3 i=1[σ]i), with a purely elastic volumetric response. The final algebraic equations considered in this paper have been implemented in this plane for this case (i.e. the residual associated with the flow rule is written in terms of the deviatoric part of the elastic strains dev[εe]:=εe−3 i=1[εe]i). A. P´erez-Foguet & F. Armero 64 cos ( p tan ( CM C q CMC FIGURE II.2. Trace on the deviatoric and meridian stress planes of the Rounded Hyperbolic Mohr–Coulomb yield surface. 2. Pressure-dependent yield functions. We consider the Rounded Hyperbolic Mohr– Coulomb (RMHC) yield function (Abbo & Sloan [1995]) f(σ)=J2(σ)K2(θσ)+(0.05 CMC cos(φ))2−(pσsin(φ)+CMC cos(φ)) ,(II.6) for the frictional angle φand cohesion CMC , and where the function K(θσ) is defined as K(θσ)=           A1−A2sin(φ)+(B1sin(φ)−B2)cos(3θσ)θσ≤5◦ 3 + sin(φ) 2√3cos(θσ)+1−sin(φ) 2sin(θσ)5 ◦≤θσ≤55◦ A1+A2sin(φ)+(B1sin(φ)+B2)cos(3θσ)θσ≥55◦ (II.7) with A1=cos(25 ◦)+2+√3 3sin(25◦) A2=2+√3 3√3cos(25◦)−1 √3sin(25◦) B1=2√2 3(3 −√3) cos(25◦) B2=2√2 3(√3−1) sin(25◦). (II.8) The RMHC surface (II.6) has been considered in Abbo & Sloan [1995,96] and P´ erez–Foguet et al. [2000a,b] in the context of infinitesimal plasticity. Figure II.2 depicts an illustration of the trace of this yield surface with the meridian and deviatoric planes. The function (II.6) is continuously differentiable for all stress states Closest-Point Projection Algorithms in Elastoplasticity 65 (Abbo & Sloan [1995]), leading to a continuous flow vector everywhere. However, the second derivative of (II.6) in σis not continuous at θσ=5 ◦,θσ=55 ◦,and therefore neither at qσ= 0. This fact, and the high curvature at some the rounded zones, causes the Newton–CPPM not to converge in large regions of the trial state space, in contrast with the new schemes proposed in this work. II.3. Hardening/softening laws The numerical examples presented in Sections 4and 5 consider perfect plasticity as well as isotropic strain hardening and strain softening. These last two cases are considered in combination with the deviatoric yield function (II.5), with the saturation potential ψh(α)=(σy∞−σyo)$α+1 δexp(−δα)% =⇒q=−∂αψh=−(σy∞−σyo)(1−exp(−δα)) ,(II.9) for the initial yield limit σyo, the saturation yield stress σy∞the saturation yield stress and the saturation exponent δ. Perfect plasticity is considered only with the pressuredependent yield surface (II.6).