Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 48, 2, pp. 381-395, Warsaw 2010 RADIAL RETURN METHOD APPLIED IN THICK-WALLED CYLINDER ANALYSIS Grzegorz Widłak Cracow University of Technology, Institute of Machine Design, Cracow, Poland e-mail: widlak@mech.pk.edu.pl The radial returnmethod for elasto-plastic constitutive equations in the special case of analysis of a thick-walled cylinder is presented. Sensitivity study concerning the integration step sizewas performed analytically for the specified case. This simplemodel was also utilized in the adjustment study of nonlinear solution strategies as a preliminary phase of FEM investigations of high pressure reactors. The robustness and step-size insensitivity of this integrationmethodwas proved and efficient solution strategies were chosen. Key words: radial return method, rate-independent plasticity, thick- walled cylinders 1. Introduction Application of the finite element method to nonlinear structural problems is connectedwith iterative calculations of nodal displacementswhich satisfy a set of simultaneous algebraic equations of the standard linear form. An iterative operation is performed for every load level n K i ndu i n = r i n r i n =pn−f(u i n) (1.1) where Kin is the current stiffnessmatrix, while r i n indicates unbalance betwe- en the external load vector pn and a vector of internal forces dependent on the current nodal displacements uin. Subscripts refer to the n-th load level while superscripts to the i-th equilibrium iteration occurring on a specified load level. Moreover, the small step incremental approach is essential to ob- tain satisfactory results when the constitutive law relating stress and strain increments is path dependent. 382 G. Widłak Selection of the solution strategy and definition of convergence criteria for the iterative process depends on specified loading conditions andmaterial characteristics. Before FEM study of complex structures, it is good practice to perform some adjustments based on the same material models but sim- pler geometry and loading conditions. Another essential preparatory step is sensitivity study of the constitutive integration algorithm in which constitu- tive equations can be enforced in a process with different increments of some imposed displacement history. The main field of the author’s investigation is an elastic-plastic analysis of a cyclically loaded thick-walled reactor with strong stress concentrators, in which shakedown effects can occur (Zieliński and Widłak, 2007; Widłak and Zieliński, 2008). In such analysis, the size of the computational model made the respective adjustment study impossible. Therefore, a simple problem of the internally loaded cylinder with plain strain conditions has been chosen. The problem of elastic-plastic analysis of thick-walled cylinders is well known and widely presented in fundamental plasticity texts (Życzkowski, 1981).Most of themgive direct solution forms.However, these analytical solu- tions are based on theHencky-Iliuszyn simple process assumptions, which are not realistic in context of significant plastic deformations in contradistinction to flow theories requiring iterative approaches. Furthermore, some assump- tions present in these solutions neglect certain significant aspects for specific analyses, e.g. strain hardening during cyclic loadings. The problem of thick walled cylinder analysis using incremental plasticity can be considered as the benchmark problem for procedures used in non-linear solution strategies as well as for algorithms used in the integration of constitutive relations. 2. Elastic-plastic constitutive equations. Radial return method The finite element procedures used in structural plasticity problems require specific integration methods for constitutive and flow equations. These me- thods are aimed at solving the main problem in computational plasticity, in which the deformation increment is assumed to be given, while the aim is to accurately and effectively compute the corresponding stress state. If the obta- ined state results in residual (1.1)2 larger than the prescribed tolerance, then a newequilibrium iteration is generated, and the next increment calculated.The secondary problem, which is related to constitutive relations, is calculation of the elasto-plastic stress-strainmatrix alongwith certain strategies as ameans to finding the solution. Radial return method applied... 383 The radial-returnmethod (SimoandTaylor, 1986;Hughes, 1984; Simoand Hughes, 1998) is themostwidely used and effective algorithm for this purpose. When striving for benchmark objectives, it suffices to consider this algorithm with simplified plasticity models. An essential idea concerning the problem, which assumes an elastic-perfectly plastic casewith theHuber-vonMises yield surface and the associated flow rule, is outlined below. Generally, the following set of governing equations has to be satisfied (Ba- the, 1995; Belytschko et al., 2000) σ̇ij =C E ijrs(ε̇rs− ε̇ pl rs) (2.1) in which σ̇ij represents stress rates, CEijrs are components of the elastic con- stitutive tensor and ε̇rs, ε̇plrs are components of total and plastic strain rates, respectively. The Huber-vonMises yield function can be stated in the following form f = 3 2 (sij ·sij)−Y 2 (2.2) where sij are deviatoric stress components and Y denotes the yield stress. This function gives the yield condition f = 0, which must hold throughout the plastic response. Assuming the active process, the plastic strain rates are given by the Prandtl-Reuss flow rule ε̇ pl ij = λ̇ ∂f ∂σij (2.3) where λ̇ is related to the effective plastic strain rate and the flow direction is indicated as normal to the von Mises cylinder. Equations (2.1) and (2.3) are further reduced by enforcing the consistency condition f = 0, which forms a nonlinear equation on λ̇ and leads to the following rate constitutive equations σ̇ij =C ep ijrsε̇rs (2.4) where Cepijrs are components of the tensor which is often referred as a matrix of the continuum elasto-plastic tangent moduli. In fact, in FEM computations, we have no longer continuum problems described by the previous constitutive equations, therefore, the incremental models as the return mapping algorithm are employed. In the following, an outline of thismethod is presentedwhich further is illustrated in theparticular case of analysis of a thickwalled cylinder.Moreover, the elasto-plastic tangent 384 G. Widłak operator, consistentwith the radial-returnalgorithm, isderived in the specified case. In the development of this algorithm, it is beneficial to present geometrical interpretation of equation (2.1). It is observed that in a plastic process the stress σ̇ is an orthogonal projection of the trial vector σ̇tr onto the tangent plane of the yield surface. This is caused by normality of the involved plastic strain rate vector ε̇pl (Fig.1). Fig. 1. Geometric interpretation of the constitutive relation for elastic-perfectly-plastic cace with the Huber-vonMises yield surface The first task of the algorithm is to calculate the elastically induced trial stress σtrn for a given strain increment dε. The trial stress is examined to confirmwhether it is inside or outside the yield surface. If it falls within or on the yield surface, the process is recognized as an elastic one, and the update stress is simply equal to the trial stress σ tr n =σn−1+C E dε (2.5) where CE is an elastic constitutive matrix (see (4.1)). If the trial stress lies outside the yield surface, it is necessary to bring it back accordingly to the consistency condition. In the radial return process, σn is defined as an intersection of the line connecting the trial stress to the center of the yield surface with this surface (Fig.2). The algorithm is convergent and exact in the limit when dε → 0. Mo- reover, the algorithm is recognized as being much more accurate than more complex procedures, and works satisfactorily even for large increments resul- ting from a small number of the applied load steps. The sensitivity study of this phenomenon was analytically performed for the thick walled cylinder problem in Section 7, and the algorithm robustness was proved. Radial return method applied... 385 Fig. 2. Illustration of the radial-return algorithm 3. Advanced nonlinear solution strategies To proceed further, it is worthwhile to review the role played by different solution strategies in the overall nonlinear analysis process. A chosen strategy has significant influence on the generated strain increment size and the rate of solution convergence. In equation (1.1)1, the matrix K i n takes on the role of a stiffness matrix in linear analysis, but now it relates small changes of load to small changes of displacements. This global stiffness matrix is computed or modified for a given deformation state, then inverted and applied to the nodal out-of-balance force vector to obtain the displacement increments. The general objective of each iteration in equation (1.1)1 can be represented as out-of-balance force reduction r i+1 n =pn−f(u i+1 n )=0 (3.1) The form of the stiffness matrix, which is used at subsequent iterations, de- pends on the assumed solutionmethod. Various schemes have been developed to solve nonlinear problems. In this section, a group of Newton methods are discussed,which are themost commonalgorithms used inFEMcomputations. In the most often used Newton-Raphson procedure, equation (3.1) can be approximated to the first order as r i+1 n ∼= rin+ (∂r ∂u )i n duin =0 (3.2) The stiffnessmatrix, corresponding to the tangent direction, can be presented in this case as (KT) i n = (∂f ∂u )i n =− (∂r ∂u )i n (3.3) 386 G. Widłak This matrix is assembled in the actual FEM procedures using element matri- ces, which are determined on the basis of constitutive elasto-plastic tangent operators (KeT) i n = ∫ e B ⊤(Cep)inBdVe (3.4) where B matrix defines relation between the strain vector and nodal displa- cements, and (Cep)in is the constitutive matrix on the n-th step and i-th iteration defined in a tensoric form by equation (2.4) and calculated using the returnmapping algorithm. Finally, equation (3.2) can be rewritten as follows (KT) i ndu i n = r i n (3.5) contributing to actual displacements u i+1 n =u i n+du i n (3.6) The N-R process is distinguished by very rapid convergence in which the asymptotic rate is quadratic. The main drawbacks of this process are con- nected with computation of the tangent stiffness matrix and factorization at each iteration. An obviousmodification of this procedure involves retention of the original tangent stiffness matrix, calculated and factorized at the first iteration to current load levels. In this case, however, the rate of convergence is much lower – linear asymptotic instead of quadratic. Other methods eliminating the tangent stiffness formation at every itera- tion utilize secant directions, which can be computed taking advantage of the known previous displacement increment and out of balance force reduction on this increment. A secant slope can be used in such a way that (KS) i+1 n du i n = r i n−r i+1 n (3.7) This is called the quasi-Newton equation. It does not uniquely define (KS)i+1n , but places a restriction on its definition. By using this approach, one may obtain convergence much more rapidly than in the modified N-R process. Up to now, several different procedures of determining (KS)i+1n have be- en proposed. All of them use some form of updating previously determined matrices or their inverses. For example, one of them, preserving the matrix symmetry and positive definiteness, is the BFGS (Broyden, Fletcher, Gold- farb, Shanno) update (Dennis and More, 1977). This method is effectively used with recomputation and factorization of the tangent matrix after some Radial return method applied... 387 number of secant matrix updates, which can, however, sometimes lead to in- stability when a higher number of updates is involved. It will be investigated in Section 6. The common feature of iterative techniques discussed above is the resul- tant displacement in the vector form duin, the purpose of which is to reduce ri+1n to zero. This is not always easy to achieve in several iterations and de- pends on loading conditions, material characteristics and the chosen solution strategy. To get a solutionwhich represents faster convergence, amodification of equation (3.6) can be introduced as (Crisfield, 1991) u i+1 n =u i n+s i ndu i n (3.8) where sin is a scalar search parameter. This, so called, line search technique is an important feature ofmost numerical techniques of optimization and can be used in a wide range of iterative solution procedures. For the simple iterative procedures outlined above, the scalar sin would be set to unity. In equation (3.8), the scalar sin becomes a variable, which can be determined for a scalar nonlinear system by evaluating ri+1n for various values of u i+1 n . For multi- degree-of-freedom systems such an approach is rather complex. It requires the introductionof a scalar parameter S(sin) of residual forces,which canbe stored while selecting sin, when r i+1 n =pn−f(u i+1 n ) (3.9) has a zero component in the direction of duin. This yields S(sin)= du i n ·r i+1 n =0 (3.10) which is a nonlinear scalar equation for sin. During one iteration several evalu- ation of this function shouldbeperformed.Matthies andStrang (1979) suggest accepting sin if |S(s i n)|< 0.5|S(0)|. With the line search technique, acceleration of convergence is remarkable, especially when applied to themodifiedNewton or quasi Newtonmethods. In many commercial systems, it is even enclosed into the standardN-Rprocedure giving significant reduction to the number of iterations (Hallquist, 2006). 4. Implicit integration algorithm for thick walled cylinders As written previously, the cylindrical high-pressure vessels are the main in- terest of the author. Therefore, a thick-walled cylinder has been chosen to 388 G. Widłak illustrate the problempresented earlier in general form.The plane strain state has been assumed, i.e. ε̇z =0. In this case, in the plane strain conditions the system of equations (2.1) changes to      σ̇r σ̇θ σ̇z      =CE      ε̇r − ε̇ pl r ε̇θ − ε̇ pl θ −ε̇plz      (4.1) where C E = 2G 1−2ν    (1−ν) ν ν ν (1−ν) ν ν ν (1−ν)    G= E 2(1+ν) Yield function (2.2) can be defined as f = 1 2 [(σr −σθ) 2+(σθ−σz) 2+(σz −σr) 2]−Y 2 (4.2) where σr, σθ, σz are the radial, hoop, axial stress, respectively and Y is the yield stress. In the particular case, the flow rule takes the form      ε̇plr ε̇ pl θ ε̇plz      = λ̇      2σr−σθ−σz 2σθ−σr−σz 2σz −σr−σθ      (4.3) where λ̇ is related to the effective plastic strain rate ε̇pleq as follows λ̇= ε̇pleq 2Y [ 1 MPa · s ] (4.4) The effective plastic strain rate is defined as ε̇pleq = √ 2 3 [(ε̇plr )2+(ε̇ pl θ )2+(ε̇plz )2] (4.5) For the numerical implementation and compatibility with the Newton type algorithms, it is required to represent equation (4.1) in the incremental form. There are several techniques of calculation of the stress states σn located on theyield surface.Herein, adirectmethodcorresponding to the returnmapping algorithm is presented. Equations (4.1) and (4.3) can be developed to Radial return method applied... 389      σr σθ σz      (n) =      σr σθ σz      (n−1) +CE      dεr−dε pl r dεθ −dε pl θ −dεplz      = (4.6) =      σr σθ σz      (tr) − dεpleq 2Y C E      2σr −σθ−σz 2σθ −σr−σz 2σz −σr−σθ      (n) where      σr σθ σz      (tr) =      σr σθ σz      (n−1) +CE      dεr dεθ 0      dλ= dεpleq 2Y were introduced. By taking advantage of the known mean plastic strain increments (dεplm =0), we can state equations (4.6) as      σr σθ σz      (n) =      σr σθ σz      (tr) −2Gdλ      2σr −σθ−σz 2σθ −σr−σz 2σz −σr−σθ      (n) (4.7) In an implicit algorithm, in contradistinction to the explicit algorithm, the right hand side of equations (4.7) is computed in the final state, which is indicated by superscript n. Using (4.7), we can calculate σn as σ(n)r = 2Gdλ(σ (tr) r +σ (tr) θ +σ (tr) z )+σ (tr) r 6Gdλ+1 σ (n) θ = 2Gdλ(σ (tr) r +σ (tr) θ +σ (tr) z )+σ (tr) θ 6Gdλ+1 (4.8) σ(n)z = 2Gdλ(σ (tr) r +σ (tr) θ +σ (tr) z )+σ (tr) z 6Gdλ+1 These valuesmust satisfy the yield condition. By introducing them into (4.2), we obtain a nonlinear equation in dλ which is dependent on dεpleq (σ (tr) r −σ (tr) θ )2+(σ (tr) θ −σ (tr) z )2+(σ (tr) z −σ (tr) r )2 2(6Gdλ+1)2 −Y 2 =0 (4.9) Once dλ is determined, all terms of the vector σn can be calculated by sub- stituting specified values into (4.8). 390 G. Widłak 5. Elastic-plastic tangent stiffness matrix for plane strain When the presented algorithm is included into theNewton type solution stra- tegies, the tangent stiffnessmatrix is needed. This can be obtained employing equation (3.4) which uses the tangent stress-strain operators. It is well pro- ved (Simo and Taylor, 1985; Ristinmaa, 1994; Hughes and Belytschko, 2008) that significant loss in the rate of convergence occurs when the elasto-plastic continuum tangent matrix is used in place of the tangent consistently derived from the integration algorithm. The consistent operator is defined as C ep n = ∂σ(σn−1,εn−1,εpleqn−1,ε−εn−1) ∂ε ∣ ∣ ∣ ∣ ε=εn (5.1) and can be obtained in practice by differentiation of stresses (4.8)with respect to time and application of condition (4.9). This enables us to find the rate of dεpleq as a function of the deformation rate. The elasto-plastic matrix finally takes the following form C ep n = 2G ( 6G dε pl eq 2Y +1 ) (1−2ν)    A+(1−ν) A+ν A+ν A+(1−ν) A+ν A+ν    + (5.2) − 1 2Y ( 1− dεpleq Y ∂Y ∂(dεpleq) ) 1 ∂f ∂(dε pl eq)    d1d1 d1d2 d2d1 d2d2 d3d1 d3d2    where subsequently A= Gdεpleq Y (1+ν) d1 = 2G(2σr −σθ−σz) 6G dε pl eq 2Y +1 d2 = 2G(2σθ −σr−σz) 6G dε pl eq 2Y +1 d3 = 2G(2σz −σr−σθ) 6G dε pl eq 2Y +1 The above matrix is of the same kind as the one presented by Simo and Taylor (1985),which incontradistinction to thecontinuummatrixused in (2.4) preserves quadratic rate of convergence in conjunction with the full Newton method. Radial return method applied... 391 6. Adjustment study of the solution strategy Theproblemof thick-walled cylinder analysis in the elasto-plastic range is con- sidered as the benchmark test for procedures used in the iterative nonlinear solution as well as for the radial return algorithm. The dimensions and inter- nal pressure of the investigated cylinder correspond to geometry and loading conditions of actual high-pressure reactors, i.e. internal radius a = 150mm, external radius b = 300mm and pressure p = 150MPa. An elastic-perfectly plasticmaterial with the yield stress Y =200MPawas assumed for the tests. The nonlinear iterative strategies presented in the previous section were employed by the author in the solution of the problem listed above, which was comparatively solved with modern FE codes – Ansys and LS-Dyna. The full Newton method and quasi Newton with the BFGS updates method pro- ved respectively, the fastest convergence rate. However, the second one was more effective because it needed calculation and factorization of the tangent stiffness matrix only at the first iteration of every load increment. The most effective solution procedure was the BFGS method in conjunction with line search technique. The displacement incrementwas scaled using the line search as long as the scalar parameter of residual forces S(sin) satisfied the condition |S(sin)|< 0.5|S(0)|. This condition was sufficient in obtaining convergence in 4 iterationswithoutmatrix recalculations (Fig.3).The rate of convergencewas almost quadratic, characteristic for the fullNewtonmethod, and the computa- tional timewas proved to be three times shorter. Therefore, the quasi Newton method with the BFGS updates was chosen in the detailed investigations of the thick walled reactors with stress concentrators. Fig. 3. Iterations in the quasi Newtonmethod with the line search 392 G. Widłak There are some ways of determining the nonlinear solution convergence. Themost popular are criteria based on specified norms and certain threshold values. The procedure continues to calculate equilibrium iterations until the convergence tolerance is satisfied. Different values of tolerances are suggested for various problems in the commercial nonlinearFE systems.The influence of the convergence criteria on the calculated hoop stress is shown inFig.4.Other stress components present a very similar course, therefore, the author resigned from presenting these plots. It suffices to use εd = 1% relative tolerance for the displacement and residual force norms. The energy relative convergence tolerancemustbeset to εc =0.1% inorder toobtainaccurate results.This last value is approximately one order smaller than this suggested by nonlinear FE systems (Hallquist, 2006). It is best to use a double-check on the convergence, e.g. simultaneousmonitoring of the displacement and energy norms. This last approach was utilized in the investigation of the reactor. Fig. 4. Convergent hoop stresses for various tolerance norms (LS-Dyna) 7. Radial return algorithm study The presented return-mapping integration algorithm, in conjunction with the quasiNewton solution strategy, is stable for awide rangeof nonlinearproblems that involve small deformations. The accuracy considerations are limited to investigations of the load increment size. In the examination of this algori- thm, certain deformations were applied, from which the trial stress state was analytically calculated. With the use of Eq. (4.9), dλ was determined and σn components were calculated from (4.8). The solution for a large number Radial return method applied... 393 of 10000 steps was treated as the reference deformation history. This was fur- ther divided into different numbers of increments and used in the study of the algorithm convergence. Even in the one-step calculations, a very good accuracy was achieved (Fig.5). The stress components differed approximately by 1MPa from the reference level. The inaccuracy associated with the choice of convergence cri- teria had a greater influence on the obtained results in the overall solution process (see Fig.4). Hence, this investigation proved robustness and step-size insensitivity of the radial return algorithm. Fig. 5. Sensitivity study of the radial return algorithm used in an analytical way (Eqs. (4.8)-(4.9)) 8. Conclusions The radial return algorithm was described and illustrated with the special case of analysis of the thick-walled cylinder in the plane strain conditions. This simplemodel was also utilized in the adjustment study of nonlinear solu- tion strategies as a preliminary phase of FEM investigations of high pressure reactors. The return mapping algorithm, applied in the Huber-von Mises as- sociated plasticity, is a robust and step-size insensitive integration method. 394 G. Widłak A closer scrutiny is suggested for a choice of the solution strategy and set- ting of parameters of the nonlinear solution convergence, because theymainly influence the efficiency and accuracy of the nonlinear solution process. Acknowledgement Iwish to thankProfessorAndrzejP.Zieliński,PhD.Eng. forhis essential guidance and for revising mymanuscript. This work was supported by PolishMinistry of Science and Higher Education – grant No. 1353/T02/2007/32. References 1. Bathe K.J., 1995,Finite Element Procedures, Prentice-Hall 2. Belytschko T., Liu W.K., Morgan B., 2000, Nonlinear Finite Elements for Continua and Structures, J.Wiley & Sons, NewYork 3. Crisfield M.A., 1991, Non-linear Finite Element Analysis of Solids and Structures, Vol. 1: Essentials, J.Wiley & Sons, NewYork 4. Dennis J.E., More J., 1977,Quasi-Newtonmethods, motivation and theory, SIAM Review, 19, 1, 46-89 5. Hallquist J., 2006,LS-DynaTheoryManual, LivermoreSoftwareTechnology Corporation 6. Hughes T., 1984, Numerical implementation of constitutive models: rate- independent deviatoric plasticity,Theoretical Foundation for Large-Scale Com- putations for Nonlinear Material Behaviour, Martinus Nijhoff Publishers, The Netherlands 7. Hughes T., Belytschko T., 2008, A Course of Nonlinear Finite Element Analysis, Paris 8. Jetteur P., 1986, Implicit integration algorithm for elastoplasticity in plane stress analysis,Engineering Computations, 3, 3, 251-253 9. Matthies H., Strang G., 1979, The solution of nonlinear finite element equations, International Journal for Numerical Methods in Engineering, 14, 11, 1613-1626 10. Ristinmaa M., 1994, Consistent stiffness matrix in FE calculations of elasto- plastic bodies,Computers and Structures, 53, 1, 99-103 11. Simo J.C., Taylor R.L., 1985, Consistent tangent operators for rate- independent elastoplasticity,ComputerMethods in AppliedMechanics and En- gineering, 48, 101-118 Radial return method applied... 395 12. Simo J.C., Taylor R.L., 1986, A returnmapping algorithms for plane stress elastoplasticity, International Journal for Numerical Methods in Engineering, 22, 3, 649-670 13. Simo J.C., Hughes T., 1998, Computational Inelasticity, Springer-Verlag, NewYork 14. Widłak G., Zieliński A.P., 2008, Local shakedown analysis of reactors sub- ject to pressure and thermal loads,Proc. 8thWorld Congress onComputational Mechanics 15. Zieliński A.P., Widłak G., 2007, Local shakedown analysis in regions of holes in high-pressure vessels,Proc. IX International Conference on Computa- tional Plasticity 16. Życzkowski M., 1981,Combined Loadings in the Theory of Plasticity, PWN- Polish Scientific Publishers,Warsaw Metoda rzutowania powrotnego w analizie stanu naprężeń cylindra grubościennego Streszczenie W pracy zaprezentowanometodę rzutowania powrotnego dla szczególnego przy- padku analizy sprężysto-plastycznej cylindra grubościennego. Wykorzystując anali- tyczne zależności, zbadano dokładność tej metody dla różnych wielkości przyrostów odkształceń. Zagadnienie analizy sprężysto-plastycznej cylindra grubościennego zo- stało również wykorzystane do oceny metod nieliniowego rozwiązywania zagadnień MES, które następnie będą użyte w badaniu reaktorówwysokociśnieniowych.Wpra- cy wykazano dużą dokładność analizowanego algorytmu. Manuscript received April 20, 2009; accepted for print September 24, 2009