INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 1909—1922 (1997) AN UPPER BOUND FINITE ELEMENT PROCEDURE FOR SOLVING LARGE PLANE STRAIN DEFORMATION CHUNG-LI HWAN Department of Mechanical Engineering, Feng Chia ºniversity, 100 ¼enhwa Road, Seatwen, ¹aichung, ¹aiwan, R.O.C. SUMMARY A unique and robust upper bound finite element procedure is developed for the analysis of large plastic deformation problems under plane strain condition. It can consistently treat problems with isotropic strain varying materials. It can also effectively solve problems with any initial ‘guessed’ velocity field, even from an random number generator. To explore and demonstrate the capability of this new approach, strip tension and plane strain compression problems are solved. For validation, the computed results are compared with existing analytical or experimental solutions in good agreement. The phenomenon of shear band formation can be simulated and, as expected, is found to develop more distinctly in strain softening materials than in perfectly plastic and strain hardening materials. ( 1997 by John Wiley & Sons, Ltd. KEY WORDS : upper bound; finite element; plane strain; large deformation INTRODUCTION Large plastic deformation problems of metal can generally be solved by two kinds of formulation, known as solid formulation and flow formulation.1 In the analysis of metal-forming problems, elastic strain of a deformed metal can be neglected because it is usually very small as compared to plastic strain. The neglect of elastic strain can simplify the formulation by treating the deforming solid as a non-Newtonian fluid. Formulation based on this simplification is known as flow formulation.2~5 It can provide most important information accurately at a fraction of the computing cost of solid formulation. However, flow formulation approach cannot predict the complete stress—strain history, which can only be obtained by using solid formulation in terms of finite strain elastoplastic incremental analysis.6~10 Classical limit analysis11,12 belonging to the category of flow formulation, it involves two principles which lead to both the lower bound and the upper bound approaches. The former predicts a load which is less than or equal to the exact limit load needed to enable the incipience of plastic flow of metals. The solution from the lower bound approach must satisfy equilibrium equations, a suitable yield criterion, and static boundary conditions. On the other hand, the upper bound approach provides a load which is at least equal to or greater than the exact limit load, and solves for a kinematically admissible velocity field which must satisfy kinematic boundary conditions. Like the classical approach,11,12 modern limit analysis uses a pair of related formulations to bound the exact solutions from above and below. But it is more mature in theory and methodology. It uses an inequality form of constitutive relation and establishes a duality relation that equates the least upper bound to the greatest lower bound.13,14 It applies computational CCC 0029—5981/97/101909—14$17.50 ( 1997 by John Wiley & Sons, Ltd. Received 15 May 1995 Revised 8 July 1996 1910 C.-L. HWAN optimization techniques to approach the corresponding maximum and minimum solutions,15 sometimes simultaneously. The rigorous convex analysis16 and computational optimization techniques help to put modern limit analysis on a solid mathematical foundation. Based on the recent advances in limit analysis mentioned above, static and kinematic solutions for many complex problems17~19 have been obtained with high degree of accuracy and certainty of convergence. A kinematic solution can be interpreted either as a steady large deformation flow solution from an Eulerian view point or an instantaneous velocity field in a Lagrangian co-ordinate system. The latter interpretation enables the integration of velocity field in a small time step to provide a corresponding displacement field, which in turn updates the configuration of the deforming body and the computational grid system. A subsequent limit problem is then solved for this new configuration. This updating process is repeated to form a sequence leading to the solution of a large deformation problem and is hereafter called sequential limit analysis. Since the von Mises yield criterion is employed and expressed in an inequality form, stresses do not appear in this upper bound formulation. As a result, neither complicated stress updating6~8 nor rigid zone treatment1 is needed in sequential limit analysis. The effect of material nonlinearity is incorporated in the analysis by using a yield criterion that varies step-by-step locally with deformation history. Using such a stepwise model, not only strain hardening material but also perfect plastic and strain softening materials can be solved consistently. Another major advantage of sequential limit analysis is the global stability of computation even in the case of increasing deformation under decreasing load, such as necking. Like most flow formulations, sequential limit analysis cannot provide certain information, such as elastic strain, residual stress and spring-back after load removal, which can only be obtained by an elastoplastic incremental analysis. However, sequential limit analysis can provide most important information with reasonable computing cost. Using an inequality to express the von Mises yield criterion, we shall first present the primal (lower bound) formulation for a solid undergoing large plastic deformation under plane strain condition. Through use of the weak equilibrium statement, the general divergence theorem, the constant shear parameter model for interface friction1 and a generalized Hölder inequality,20 we then arrive at the dual (upper bound) formulation. A duality theorem is then stated to equate the greatest lower bound to the least upper bound. Numerical procedure for the upper bound formulation is then presented. Finally, examples related to strip tension and plane strain compression tests are computed for demonstration. LOWER BOUND FORMULATION Lower bound problems arise naturally thus their mathematical formulation is called the primal. Using equilibrium equations, static boundary conditions and a yield criterion, a lower bound problem under the plane strain condition can be posed as follows. P! pnn dS maximize j" 1 subject to + · p"0 in D p · n"0 on ! 0 EpE )k'(e ) V 1 in D INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) (1) ( 1997 by John Wiley & Sons, Ltd. FINITE ELEMENT PROCEDURE FOR LARGE PLANE STRAIN DEFORMATION 1911 where + · denotes the two-dimensional divergence operator, p is the symmetric 2]2 stress matrix function, 0 is the two-dimensional null vector, D is the two-dimensional domain, ! is the traction 0 free boundary and the integral condition j": p dS represents the limit load applied on the !1 nn boundary ! by the applied normal stress p . In this study, j is the tensile load for the tension 1 nn problem or the press load for the compression problem. In equation (1), the von Mises yield criterion is generally modified to take the effect of strain hardening or softening into account. EpE is the norm notation used to denote the von Mises V yield function of stress matrix function p. k is the yield stress in shear determined experimentally. '(e ) is a strain hardening or softening function of e , which is the accumulated equivalent plastic 1 1 strain in a local sense. For a perfectly plastic material, '(e )"1. 1 DUALITY AND UPPER BOUND FORMULATION The exact upper bound formulation under the plane strain condition can be derived through an upper bounding process. Using the virtual work concept, the equilibrium equation in problem (1) is satisfied in a weak sense such that PPD u · (+ · p)dA"0 (2) for all u"Mu(x, y), v(x, y)N3K, where u and v are components of the velocity vector u in both x and y directions, respectively. The set K is the space of all kinematically admissible velocity functions, u, which satisfy the incompressibility and kinematic boundary conditions. After using a general divergence theorem and a sequence of mathematical manipulation,21 equation (2) gives PPDp : eR dA#mk P! D u4 D dS j" (3) & where eR is the symmetric strain rate matrix, u is the tangential component of velocity vector along 4 the frictional boundary ! and 0)m)1 is the lubrication factor commonly used to model the & induced shear stress from friction, q "!mk sign (u ). 4 4 Consequently, the lower bound formulation defined by equation (1) can be restated as maximize j PPDp:eR dA#mkP! Du4 D dS subject to j" & + · u"0 (4) EpE )k'(e ) V 1 kinematic boundary conditions on LD where + · u"Lu/Lx#Lv/Ly"0 denotes the incompressibility condition which is required in the study of metal plasticity and LD is the total boundary of the domain D. ( 1997 by John Wiley & Sons, Ltd. INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) 1912 C.-L. HWAN Using a generalized Hölder inequality,20 Dp : eR D)EpE EeR E , and the von Mises yield criterion V " defined in equation (1), j expressed by equation (3) is bounded such that PPD p : eR dA#mk P! Du4 D dS ) PPD EpEV EeR E" dA#mkP! Du4 D dS )k PPD '(e1)EeR E " dA#mkP! Du4 D dS j" & & & "jM (u) (5) where jM (u) is an upper bound functional and EeR E is the dual norm of the strain rate matrix " function eR , which can be expressed, in terms of its components, as EeR E "J(eR !eR )2#4eR 2 xy " xx yy SA " B A B Lu Lv 2 Lu Lv 2 # ! # Lx Ly Ly Lx (6) through use of the flow rule associated with von Mises yield criterion. Once having obtained the upper bound functional, the upper bound problem can be mathematically defined as the dual one which seeks the smallest value of jM (u) where u3K. The dual problem is a constrained minimization one which can be posed as minimize jM PPD '(e1 )SA Lx!Ly B #A Ly#LxB dx dy #mk P! Du4 D dS Lu subject to jM "k Lv 2 Lu Lv 2 (7) & Lu Lv # "0 Lx Ly kinematic boundary conditions on LD Because the inequality remaining between both the upper bound and the lower bound functionals is sharp (equality inclusive) as shown by equation (5), the upper bound problem yields a weak duality relation to the lower bound problem. The weak duality relation establishes the equivalence of the least upper bound to the greatest lower bound. This is mathematically known as a saddle-point problem which can be stated as maximize j(p)"j*"minimize jM (u) (8) where the exact solution j* can be obtained by either maximizing a lower bound functional j or minimizing an upper bound functional jM . The minimizing of an upper bound functional jM is called INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) ( 1997 by John Wiley & Sons, Ltd. FINITE ELEMENT PROCEDURE FOR LARGE PLANE STRAIN DEFORMATION 1913 the upper bound approach in which the duality condition (equality) can be obtained by a proper choice of K and a minimizer u3K. NUMERICAL IMPLEMENTATION To solve numerically the upper bound problem given by equation (7), the following procedure is employed. First, through the use of the penalty function procedure22 to treat the incompressibility constraint, the constrained minimization problem is transformed into an unconstrained minimization one. Since a standard non-linear programming method is not applicable to solve the unconstrained minimization problem which involves a functional with square root integrand, a combined smoothing successive approximation method4 is then applied to solve the unconstrained minimization problem iteratively by a sequence of quadratic programming problems. A quadratic programming problem in each iteration is then discretized using both the fournode quadrilateral element for the subdomains and the two-node line element for the frictional subboundaries. Finally, by assembling all the element stiffness matrices and element velocity vectors into their corresponding global positions,21 the approximate functional in the nth iteration, jM , becomes a quadratic form in matrix notation as n jM "qT K q n n n n (9) where q is the global velocity vector in the nth iteration and K is the global stiffness matrix which n n is symmetric, banded and positive definite. Having obtained the quadratic functional through use of the finite element discretization, the approximate problem in the nth iteration can be stated as minimize jM n subject to jM "qT K q n n n n (10) kinematic boundary conditions on LD which can be solved for the velocity vector q in terms of the optimality condition. n Since the procedure is iterative, a suitable criterion is needed to terminate the iteration. A convergence criterion used in this study is the error norm of the velocity vectors, which is defined as E "Eq !q E /Eq E . This norm is used to calculate the deviation measured in u n n~1 2 n~1 2 the Euclidean space in terms of the calculated velocity vectors between two successive iterations and required to decrease from iteration to iteration. When it is smaller than 10~5, the solution is regarded as good enough and the iteration stops. The ‘converged’ jM and q are interpreted as the n n limit load, jR, and the corresponding global velocity vector, qR, for the original unconstrained minimization problem. Both geometry and material property for the successive configuration are updated by using the velocity field obtained from the current configuration. Geometry is updated by multiplying the cnodal velocity vector with a small pseudo-time increment, dt, which is a scale factor related to the step size, to obtain a small increment of the displacement vector at each node. The displacement increment is then used to update the current configuration and computational grid system. A sequence of updating leads to large deformation. Material property for the new configuration is updated as follows. The velocity field is first used to calculate the equivalent strain rate, which is defined as eN "J2/3(eR eR )1@2, at each Gaussian ij ij ( 1997 by John Wiley & Sons, Ltd. INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) 1914 C.-L. HWAN integration point. The calculated equivalent strain rate is then multiplied by dt to yield a small increment of equivalent plastic strain. The equivalent plastic strain increment is then used to update the accumulated equivalent plastic strain at each Gaussian integration point, which in turn gives the updated yield stress by an interpolation on a chosen material property curve. To proceed the analysis of large deformation, the pseudo-time increment, dt, must be carefully selected. In general, dt, must be chosen small enough such that dtEMº*NE is much smaller than = ¸, where ¸ is a characteristic length of the problem. This condition validates the small deformation theory. Computationally, the choice of dt may also be subject to other factors, such as the time necessary for the next free node to contact the die or tool surface, the desired maximum strain-increment, and the error in the volume constancy.1 RESULT AND DISCUSSION Strip tension problems Tensile specimens have been known to exhibit localization of plastic flow into neck-shaped regions near the stage of maximum load in tensile tests for many years. Even when great care is taken to achieve uniformity of material properties and symmetry of loading, the formation of necking has been unavoidable. We wish to simulate this well-known necking phenomenon computationally as the first test case for sequential limit analysis. As shown in Figure 1, a subdomain of the strip, D , is chosen as the computational domain abcd with initial dimensions ¸]H and an aspect ratio ¸/H"6. Such a subdomain is discretized by the finite element mesh system with 341 nodes and 300 elements. As the strip elongates, a neck is expected to form at boundary a—d where symmetric boundary conditions and denser grids are specified. Symmetric conditions are also specified for boundary a—b. Shear-free conditions are specified for the right end b—c, where the tensile force, f, is applied and stress-free conditions are specified for boundary c—d. Using a perfectly plastic material in the analysis, an initial test yields the following results. Let 2h denote the neck thickness, d denote the displacement at the right end and k be the yield stress in shear. The thickness ratio, h/H, and the normalized applied load, ""f/4Hk, are plotted against d/¸ in Figure 2 in solid and dash curves, respectively. Since elastic strain is neglected and perfectly plastic material is used, necking begins at the start of the deformation and the load decreases with reduced neck thickness as expected. If the tensile stress in the neck cross-section were uniformly equal to 2k as a state of uniaxial stress would be, we would have f /2h"2k and the curves f/4Hk and h/H in Figure 2 would then coincide. The fact that the load curve is slightly above the thickness curve requires an explanation. The reason is that the biaxial stress states near the surface (point d) and the triaxial stress states near centre (point a) allow the tensile stress in the Figure 1. Schematic diagram of a strip in tension INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) ( 1997 by John Wiley & Sons, Ltd. FINITE ELEMENT PROCEDURE FOR LARGE PLANE STRAIN DEFORMATION 1915 Figure 2. " and h/H as functions of d/¸ for a strip in tension cross-section to be higher than 2k according to the von Mises yield criterion. The tensile stress is maximum at the centre of the neck.23 A sequence of grids and equivalent plastic strain contours for the above problem are shown in Figure 3. The top one is for the undeformed configuration (d/¸"0). The second one denotes a deformed configuration (d/¸"0·025) in which the contour values of equivalent plastic strain are 0·03, 0·06, 0·09, 0·12, 0·15 and 0·18, respectively. The next one depicts the deformed configuration (d/¸"0·05) with contour values 0·15, 0·3, 0·45 and 0·6. The fourth one is another deformed configuration (d/¸"0·075) where strain contour values are 0·25, 0·5, 0·75, 1·0 and 1·25. The last deformed configuration (d/¸"0·1) shows strain contours with values 0·4, 0·8, 1·2 and 1·6. Since the largest contour value in all these configurations is at the centre of the neck and increases with the deformation, it is evident that the phenomenon of strain localization is closely related to the development of necking. To investigate the effect of material property on the development of necking, a mathematical model of strain varying materials, '(e )"p /p "1#Men , is used in the analysis. p is the 1 0 1 : 0 initial equivalent yield stress which is equal to 2k in plane strain condition, p is a subsequent : equivalent yield stress and e is the accumulated equivalent plastic strain. Choosing n equal to 1, 1 such a model describes a class of materials with slope M. The ones with positive slope M correspond to linear hardening materials. On the other hand, those with negative slope M correspond to linear softening materials. Perfectly plastic materials are modelled by choosing M equal to zero. The normalized applied loads, ""f/4Hk, are plotted against d/¸ in Figure 4 for a class of materials with different slope M. It is shown in Figure 4 that the phenomenon of necking develops more quickly for softening materials than for perfectly plastic solids. The development of necking in hardening materials is slower than in perfect plastic solids and the onset of necking may even ( 1997 by John Wiley & Sons, Ltd. INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) 1916 C.-L. HWAN Figure 3. Grids and equivalent plastic strain contours for a strip in tension Figure 4. " as functions of d/¸ for materials with slope M in plane strain tension be retarded due to the effect of strong strain hardening. For the same end displacement, d/¸, softening materials tend to yield at the lower tensile load while hardening materials act to yield at the higher tensile load as compared to perfectly plastic materials. INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) ( 1997 by John Wiley & Sons, Ltd. FINITE ELEMENT PROCEDURE FOR LARGE PLANE STRAIN DEFORMATION 1917 Figure 5. Grids and equivalent plastic strain contours for materials with slopes equal to !1, !0·5, 0, 0·5, 0·75, respectively, in plane strain tension To further investigate the influence of material property on the development of necking, the deformed grids and equivalent plastic strain contours for five kinds of materials at the same end displacement, d/¸"0·0375, are plotted in Figure 5. The top one is for M equal to !1; the second one is for M equal to !0·5; the third one is for M equal to 0 (perfectly plastic materials); the fourth one is for M equal to 0·5; the last one is for M equal to 0·75. Because of strain softening, shear band localization appears distinctly in both softening materials (M"!1 and M"!0·5). On the other hand, the effect of strain hardening tends to retard strain localization and yield a more uniform deformation in the specimen before necking as illustrated by the last two plots of Figure 5. Block compression problems Because necking in a tensile specimen destroys the homogeneous strain/stress field, thus rendering the data useless beyond the point, compression tests are commonly used to determine the flow curve up to high strains. However, there still exist problems associated with compression; tests as well. Short specimens must be used to avoid buckling. The end of the specimen in contact with the press heads must be lubricated to allow free-side expansion. The lubrication is never perfect in reality. Poor lubrication will cause the mid-section of the specimen to bulge. By bulging, the material closer to the centre moves outward faster than the material near the platens. This again destroys uniaxiality. We wish to simulate this well-known bulging phenomenon computationally as the second test case for sequential limit analysis. To facilitate the study of compression problems under the plane strain condition, a block compressed by a pair of flat rigid platens with a compression force per unit thickness, f , is schematically shown in Figure 6. A subdomain of the block, D , is chosen as the computational abcd domain with initial dimensions ¼]H and an aspect ratio ¼/H"1. Such a subdomain is discretized by the finite element mesh system with 289 nodes and 256 elements. Symmetric ( 1997 by John Wiley & Sons, Ltd. INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) 1918 C.-L. HWAN Figure 6. Schematic diagram of a block in compression Figure 7. " as a function of 1!h/H for perfectly plastic materials under plane strain compression boundary conditions are specified at both boundaries a—b and a—d. Stress-free conditions are specified for the boundary b—c. Mixed boundary conditions are specified for boundary c—d by imposing displacement boundary conditions along the travelling direction of press heads and frictional shear stresses along the interface between press heads and the workpiece. Employing the shear parameter friction model, plane strain compression tests are simulated under two extreme lubrication conditions (m"0 for perfect lubrication; m"1 for no lubrication) for perfectly plastic materials. Let the distance between the press heads, 2h (0(h )H), be the control parameter for the compression sequences and k denote the yield stress in shear. The computed results are shown in Figure 7, where the normalized applied load, ""f /4¼k, is plotted as a function of (1!h/H) for the cases m"0 and 1. The results for are verified by the INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) ( 1997 by John Wiley & Sons, Ltd. FINITE ELEMENT PROCEDURE FOR LARGE PLANE STRAIN DEFORMATION 1919 simple exact solution f"4¼Hk/h. As expected, the load—displacement curve for m"1 becomes gradually higher than that for m"0 as deformation continues. The reason is that the effect of friction have gradually caused severe distortion along the diagonals of the block, which requires higher plastic work than that for a uniform deformation without friction. In the simulation of plane strain compression of perfectly plastic materials under no lubrication, the occurrence of severe mesh distortion in the computational domain can cause deterioration of the computational results and even prohibit further analysis at some step. To ensure the quality of our computed results and to continue on the analysis until h/H"0·5, mesh rezoning24 was performed at 1!h/H"0·3. Grids and equivalent plastic strain contours before and after mesh rezoning are shown in Figure 8 for comparison. Grids and equivalent plastic strain contours for m"1 at six different states of deformation (1!h/H"0, 0·1, 0·2, 0·3, 0·4, 0·5) are depicted in Figure 9 for illustration. These plots show that a slight bulge occurs during the deforming process together with a severe distortion along the Figure 8. Grids and equivalent plastic strain contours before and after mesh rezoning Figure 9. Grids and equivalent plastic strain contours for perfectly plastic materials under plane strain compression at 1!h/H"0, 0·1, 0·2, 0·3, 0·4 and 0·5, respectively ( 1997 by John Wiley & Sons, Ltd. INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) 1920 C.-L. HWAN diagonal line a —c. These plots also show that the stress-free edge near the corner of the block progressively folds itself over and comes into contact with the press head due to large rotation. Computationally, the phenomenon of contact cannot be treated in a continuous sense as it must be physically. As a result of using a discrete contact model in computer program, several discontinuities exists in the dotted curve in Figure 7. To investigate the effect of strain hardening on strain localization and the formation of bulging, aluminum is chosen for the simulation of block compression because it exhibits significant strain hardening. The mechanical behaviour of aluminum is expressed by /(e )"p /p "(1#16·4e )0>25, 1 : 0 1 which is a mathematical model that relates equivalent yield stress, p , to the accumulated : equivalent plastic strain, e , for quasi-static compression.25 p is the annealed yield stress of 1 0 106 Mpa. The computed results are shown in Figure 10, where the normalized applied load, ""f/4¼k, is plotted as a function of (1!h/H) for the cases m"0 and 1. Due to the effect of strain hardening, the computed applied loads for aluminum are higher than those for perfectly plastic materials as shown in Figure 7. Grids and equivalent plastic strain contours for aluminum under plane strain compression with no lubrication are depicted in Figure 11 for illustration. These plots show that a more significant bulging and a more uniform strain distribution occurs during the deforming process than in the case of perfectly plastic materials. This is because strain hardening can effectively prevent a concentration of strain from occurring and result in a more uniform strain distribution and a more significant bulging.26 Because of a more uniform strain distribution, the simulation can be performed up to more than 50 per cent deformation without severe mesh distortion. Because of large rotation and friction, the maximum accumulated equivalent plastic strain for each deformed configuration is found near the original corner point c. Figure 10. " as a function of 1!h/H for aluminum under plane strain compression INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) ( 1997 by John Wiley & Sons, Ltd. FINITE ELEMENT PROCEDURE FOR LARGE PLANE STRAIN DEFORMATION 1921 Figure 11. Grids and equivalent plastic strain contours for aluminum under plane strain compression at 1!h/H"0, 0·1, 0·2, 0·3, 0·4 and 0·5, respectively CONCLUSION A new approach, sequential limit analysis, has been presented in this study for solving large plastic deformation problems under plane strain condition. This method is an extension of a modern limit analysis. It is not only mathematically concise in formulation but it is also computationally easy to implement. Both plane strain tension and compression problems have been studied by this method. The formation of shear band localization can be simulated in plane strain tension and found to develop more distinctly in strain softening materials than in both perfectly plastic and strain hardening materials. As a result, necking occurs earlier and develops quicker in strain softening materials. In plane strain compression of perfectly plastic materials, the formation of shear band may eventually lead to severe mesh distortion, which requires the use of mesh rezoning if further computation is required. Because of the effect of strain hardening, no severe mesh distortion occurs in the compression of aluminum and the analysis can be performed up to more than 50 per cent deformation without mesh rezoning. The present analysis is done without the consideration of formability, which is important in practical metal-forming problems and can be evaluated if an adequate ductile fracture model is chosen. REFERENCES 1. S. Kobayashi, S. I. Oh and T. Atlan, Metal Forming and ¹he Finite Element Method, Oxford University Press, New York, 1989. 2. C. H. Lee and S. Kobayashi, ‘New solution to rigid plastic deformation problems using matrix method’, Journal of Engineering for Industry, 95, 865—873 (1973). ( 1997 by John Wiley & Sons, Ltd. INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) 1922 C.-L. HWAN 3. O. C. Zienkiewicz and P. N. Godbole, ‘Flow of plastic and visco-plastic solids with special reference to extrusion and forming processes’, Int. j. numer. methods eng., 8, 3—16 (1974). 4. W. H. Yang, ‘A variational principle for plastic flows’, Proc. 9th º.S. National Congress of Applied Mechanics, Cornell University, ASME, 1982. 5. O. C. Zienkiewicz, ‘Flow formulations for numerical solutions of forming processes’, in J. F. T. Pittman et al. (eds.), Numerical Analysis of Forming Processes, Wiley, New York, 1984. 6. E. H. Lee, ‘Elastic—plastic deformation at finite strain’, J. Appl. Mech., 36, 1—6 (1969). 7. H. D. Hibbit, P. V. Marcal and J. R. Rice, ‘A finite element formulation for problems of large strain and large displacement’, Int. J. Solids Struct., 6, 1069—1086 (1970). 8. R. M. McMeeking and J. R. Rice, ‘Finite element formulation for problems of large elastic—plastic deformation’, Int. J. Solids Struct., 11, 601—616 (1975). 9. A. S. Wifi, ‘An incremental complete solution of the stretch forming and deep drawing of a circular blank using a hemispherical punch’, Int. J. Mech. Sci., 18, 23—31 (1976). 10. E. H. Lee, R. L. Mallett and W. H. Yang, ‘Stress and deformation analysis of the metal extrusion process’, Comput. Methods Appl. Mech. Eng., 10, 339-353 (1977). 11. W. Prager and P. G. Jr. Hodge, ¹heory of Perfectly Plastic Solids, Chapman & Hall, London, 1951. 12. B. M. Fraeijs de Veubeke, ‘Upper and lower bounds in matrix structural analysis’, in B. M. Fraeijs de Veubeke (ed.), Matrix Methods of Structural Analysis, Pergamon Press, Oxford, 1964. 13. R. Temam and F. Demengel, ‘Duality and limit analysis in plasticity’, in W. H. Yang (ed.), ¹opics in Plasticity, AM Press, Ann Arbor, MI, 1991. 14. W. H. Yang, ‘A duality theorem for plastic plates’, Acta Mech., 69, 177—193 (1987). 15. D. G. Luenberger, ¸inear and Nonlinear Programming, Addison-Wesley, Reading, MA, 1984. 16. R. T. Rockaffellar, Convex Analysis, Princeton University Press, Princeton, 1970. 17. H. Huh, ‘Limit analysis in plane stress’, Ph.D. Dissertation, The University of Michigan, Ann Arbor, 1986. 18. K. H. Liu, ‘Limit analysis of plane strain extrusions’, Ph.D. Dissertation, The University of Michigan, Ann Arbor, 1986. 19. T. Tyan, ‘Simulation and optimization for metal forming processes: cutting, rolling and extrusion’, Ph.D. Dissertation, The University of Michigan, Ann Arbor, 1990. 20. W. H. Yang, ‘On generalized Hölder inequality’, Nonlinear Anal., 16, 489—498 (1991). 21. C. L. Hwan, ‘Large plastic deformations by sequential limit analysis: a finite element approach with applications in metal forming’, Ph.D. Dissertation, The University of Michigan, Ann Arbor, 1992. 22. O. C. Zienkiewicz and P. N. Godbole, ‘A penalty function approach to problems of plastic flow of metals with large surface deformations’, J. Strain Anal., 10, 180—185 (1975). 23. H. L. Morrison and O. Richmond, ‘Large deformation of notched perfectly plastic tensile bars’, J. Appl. Mech, 971—977 (1972). 24. Jung-ho Cheng and N. Kikuchi, ‘A mesh rezoning technique for finite element simulation of metal forming processes’, Int. j. numer. methods eng., 23, 219—228 (1986) 25. W. Johnson and P. B. Mellor, Engineering Plasticity, Ellis Horwood, Chichester, U.K., 1983. 26. R. Houlston, G. W. Vickers and D. L. Anderson, ‘A finite element and experimental study of rigid—plastic compression’, Int. j. numer. methods eng., 23, 1407—1437 (1986). . INT. J. NUMER. METHODS ENG., VOL. 40: 1909—1922 (1997) ( 1997 by John Wiley & Sons, Ltd.