# Transient algorithms for heat transfer general developments and approaches for theoretically generating Nth-order time-accurate operators including practically useful second-order forms

код для вставкиСкачатьINTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS ZDZISLAW WI E V CKOWSKI?;? , SUNG-KIE YOUN? AND BYUNG-SIK MOONД Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, 373-1, Kusong-dong, Yusong-gu, Taejon 305-701, Korea SUMMARY A stress-based model of the nite element method is evolved for two-dimensional quasi-static plasticity problems. The self-equilibrating elds of stresses are constructed by means of the Airy stress function, which is approximated by three types of elements: the Bogner?Fox?Schmit rectangle, the Hsieh?Clough? Tocher triangle and its reduced variant. Traction boundary conditions are imposed by the use of the Lagrange multiplier method which gives the possibility of calculation of displacements for boundary points. The concept of multi-point-constraints elements is applied in order to facilitate the application of this technique. The iterative algorithm, analogous to the closest-point-projection method commonly used in the displacement-based nite element model, is proposed for solving non-linear equations for each load increment. Two numerical examples with stress- and displacement-controlled load are considered. The results are compared with those obtained by the displacement model of FEM. Bounds for limit loads are obtained. Copyright ? 1999 John Wiley & Sons, Ltd. KEY WORDS: stress function nite element method; equilibrium model of FEM; plasticity; load limit bounds; variational inequalities; INTRODUCTION The stress-based nite element method utilizing the principle of complementary work or stationary complementary energy is used rather rarely in the engineering practice. The possibility of application of this model of FEM is mentioned only in several nite element textbooks,1?4 written among many others. It seems that the main reasons of the rare use of the model in comparison with the displacement nite element method are two facts: the statically admissible elds of stresses are more dicult to construct than the kinematically admissible elds of displacements (at least in case of plane problems); the stress boundary conditions lead to linear equations when stress functions are used, while?in the displacement formulation?nodal displacements can be prescribed usually. Despite these diculties, the stress-based FE model is attractive since the use of it together with ? Correspondence to: Zdzislaw Wi e Vckowski, Department of Mechanics of Materials, Technical University of L├od├z, Al. Politechniki 6, 93-590 L├od├z, Poland. E-mail: zwi@kmm-lx.p.lodz.pl ? Post-doctoral fellow. Permanently Lecturer, Department of Mechanics of Materials, Technical University of L├ od├z, Al. Politechniki 6, 93-590 L├od├z, Poland ? Associate Professor Д Graduate student Contract=grant sponsor: Korea Science and Engineering Foundation CCC 0029?5981/99/101505?21$17.50 Copyright ? 1999 John Wiley & Sons, Ltd. Received 19 January 1998 Revised 19 June 1998 1506 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON the displacement FE model provides the lower and upper bounds for the strain energy, load limit, etc., and makes possible to estimate the error for the approximate solution in an easy way. As the rst papers concerning the stress approach to FEM, those written by Fraeijs de Veubeke,5 and Fraeijs de Veubeke and Zienkiewicz6 should be cited, where the analogy between plane problems and the bending plate problem has been shown. Using this analogy, Morley7 and Elias8 applied the stress approach to the Kirchho theory of plates. Watwood and Hartz9 solved the plane problems of elasticity. Fraeijs de Veubeke and Hogge10 considered the dual approach to the two-dimensional problem of heat ow using the triangular element with linear shape functions for components of the heat-ux vector. The equilibrium model has been used more frequently in the analysis of torsion of a prismatic bar,1 also in the case of plastic analysis,11; 12 and stationary seepage problems.13; 14 During the last two decades, some newer concepts were introduced. Taylor and Zienkiewicz15 used the penalty method to satisfy the equilibrium equations?this method can be regarded as a mixed approach. Gallagher, Heinrich and Sarigul solved plane problems of elasticity using the Airy stress function and the concept of blending interpolation technique in order to facilitate the satisfaction of the boundary conditions.16; 17 Vallabhan and Azene18 also solved the plane problems, but used Lagrange multipliers to impose the boundary conditions. Other recent works by Gallagher and co-workers19; 20 should also be cited here. It should be noted that the stress boundary conditions can also be satised by the use of the penalty method (see e.g. Reference 1) which does not increase the number of unknown variables?this is a relatively simple but approximate technique. The stress-based formulation of FEM was utilized in the theory of plasticity by Rybicki and Schmit,21 Gallagher and Dhalla,22 and Azene.23 The plane problems are solved by the above-mentioned authors?they used the elastic?plastic constitutive matrix to calculate the stress increments in the iterative procedure applied for each load increment. Belytschko and Hodge24 used the stress approach for the limit analysis of the plane stress problem. The application of complementary variational principles to elastodynamics and elastic buckling was considered by Tabarrok, Sodhi and Simpson in References 25?27. Using dual formulations of FEM, WiVeckowski found lower and upper bounds of eective moduli for an elastic, brous composite material,28; 29 and analysed its eective constitutive relations in the case of elastic?plastic behaviour.30 The application of the stress-based FEM to plane friction problems was outlined by the same author in Reference 31. In the present paper, the Airy stress function is employed to construct the self-equilibrating stress elds. The function is approximated by the Bogner?Fox?Schmit rectangular32 and Hsieh? Clough?Tocher triangular33 elements. The Lagrange multiplier technique is used to full the traction boundary conditions?the concept of multi-point-constraints elements is developed to make this technique easier to implement and use. The proposed iterative procedure for solving the non-linear system of nite element equations for each load increment is analogous to the closestpoint-projection method used commonly in the displacement formulation of FEM. It is shown that the equilibrium model of FEM can be implemented relatively easily in the standard nite element program. SETTING OF THE PROBLEM The problem of equilibrium of an isotropic elastic?plastic body subjected to plane stress or plane strain state is considered. Although these two-dimensional problems are considered, for the sake Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS 1507 of brevity, most of the mathematical formulas will be written for the general three-dimensional case throughout the paper. Let denote an open bounded domain in R3 with regular boundary @ which consists of two parts u and such that u ? =@ ; u ? = ?. Let I denote the time interval, I ?[0; T ]; T ┐0. The solution to the considered quasi-static problem: the elds of displacements, ui (t), and stresses, ij (t), satisfy the following relations: (i) The strain?displacement relations ij = 12 (ui; j + uj; i ) on (1) where ij denotes the tensor of small deformations, ui the displacement vector. (ii) The equilibrium equations ij; j + bi = 0 on (2) where ij is the stress tensor, and bi the vector of body forces. (iii) The Prandtl?Reuss constitutive relations (see e.g. References 34 and 35) ? ?ij = ?eij + ?pij ? ? ?eij = Cijkl ?kl on ? ? p ?ij (ij ? ij )60 ? ? B (3) where the dot indicates the time derivative, ije and ijp are the elastic and plastic parts of the strain tensor, respectively, and Cijkl denotes the tensor of elastic compliance, 1 (4) [(1 + ) ? (1 + ) ] E where is equal to 0 or 1 for the plane stress problem or the plane strain problem, respectively, C = B = {: f() = q ? 0 60} (5) is the set of plastically admissible stresses, where 0 denotes the plastic limit in the uniaxial tension test, and q is the stress invariant dened as follows: q q = 32 sij sij ; sij = ij + pij ; p = ? 13 ii (the Huber?von Mises yield condition is considered). (iv) The boundary conditions ui = Ui on u; ij nj = ti on where Ui and ti are given displacements and tractions on the unit vector outwardly normal to the surface @ . (v) The initial conditions are taken as ui = 0; ij = 0 for t = 0 u and , respectively, and ni is (6) where t is the time variable. In the plane stress problem, the following conditions hold 3i ? 0, 3 ? 0, b3 ? 0, while in case of the plane strain problem 3 ? 0, 3i ? 0, b3 ? 0. It is assumed above that the Latin indices vary from 1 to 3 while the Greek ones from 1 to 2. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1508 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON DUAL VARIATIONAL FORMULATIONS OF THE PROBLEM Let V and W denote the following spaces: V = {v ? [H 1 ( )]3 : vi = U? i (t) on u} 3 W = {v ? [BD( )] : vi = U? i (t) on u} where H 1 ( ) is the Sobolev space (see e.g. Reference 36, and BD( ) is the space of functions of bounded deformation (see References 37 and 38). Let P denote the set P = { ? [L2 (Y )]9 : ij = ji ; (x) ? B(x) a:e: on } where L2 ( ) is the space of square integrable functions (Reference 36), B is the set of plastically admissible stresses dened in equation (5). From the constitutive relations (3), the following inequality can be derived: (?ij ? Cijkl ?kl )(ij ? ij )60 ? ? B After integrating this relation over domain , and using the strain?displacement relations (1) and the denition of set P, we get the inequality Z (Cijkl ?kl ? u?i; j )(ij ? ij ) dx┐0 ? ? P After multiplying the equilibrium equations (2) by (vi ? u?i ) and integrating the result over domain , and then using the Green formula, we obtain the equation expressing the principle of virtual work Z Z Z ij (vi; j ? u?i; j ) dx = bi (vi ? u?i ) dx + ti (vi ? u?i ) ds ?v ? V The problem can be set variationally as follows: Find (; u?) : I ? P О W such that a.e. on I b(; ? ? ) ? (e(u?); ? ) ┐ 0 ? ? P (; e(v ? u?)) = l(v ? u?) (7) ?v ? V (8) (0) = 0 where the following notation is used: Z Cijkl ij kl dx; b(; ) = [e(v)]ij = 1 (vi; j + vj; i ); 2 (9) Z (; ) = ij ij dx Z l(v) = Z bi vi dx + ti vi ds (10) The variational formulation (7)?(9) is the starting point for the formulation of the problem in displacements, in stresses, and in mixed forms. In order to set the problem in stresses, we introduce the following set of statically admissible elds of stresses Y = { ? [L2 ( )]9 : ij = ji ; ij; j + bi = 0 on ; ij nj = ti on Copyright ? 1999 John Wiley & Sons, Ltd. } Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS 1509 The equation of virtual work (8) is satised identically for each ? Y , and the following equation is true Z Z U? i (ij ? ij )nj ds ?u? ? W; ? ? Y u?i; j (ij ? ij ) dy = (e(u?); ? ) ? u Then the considered problem can be formulated in the dual form as follows: Find : I ? K such that a.e. on I b(; ? ? ) ┐ L(U? ; ? ) ? ? K (0) = 0 where (11) (12) Z U? i ij nj ds; L(U? ; ) = K =Y ?P u The quasi-static variational formulation in stresses for the plasticity problem has been given, for example, in References 34 and 38. FINITE ELEMENT SOLUTION Time discretization Let us consider a partition of the time interval I into N subintervals as 0; t1 ; : : : ; tn ; : : : ; tN = T (0Аt1 А и и и Аtn А и и и АtN ). The time derivative is approximated by the backward dierence scheme (t ? n) ? = n n ? n?1 ? tn tn ? tn?1 (13) where n ? (tn ). Substitution of (13) into (11) gives the semi-discrete (with respect to time) formulation: For n = 1; 2; : : : ; N , nd n ? K such that n = n?1 + n n n n (14) n b( ; ? ) ┐ L(U ; ? ) ? ? K 0 =0 (15) (16) n The solution to the inequality (15), , minimizes the functional of complementary energy Z Z 1 1 Cijkl ij kl dx ? Uin ij nj ds (17) () = b(; ) ? L(U n ; ) ? 2 2 u on the set K ? K(tn ). Construction of self-equilibrating stress elds According to the denition of set Y , the elds of stresses?elements of this set?must satisfy the equilibrium equations inside region and on its boundary . In the two-dimensional problems Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1510 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON considered in the paper, the equilibrium equations (2) are satised identically if the stresses are expressed by the Airy stress function, F, as follows: = e e F; + ? (18) where e and are the permutation and Kronecker symbols, respectively, and ? is the potential of volumetric forces such that b = ? ?; : The problem of fullling the stress boundary conditions is described in the Appendix. Space discretization Stresses being the second derivatives of the main variable, the Airy stress function, require the function to be approximated with the help of elements of C 1 class. Several types of elements have been used in previous studies. The Bogner?Fox?Schmit rectangular element (see References 32 and 3) has been used by Gallagher and Dhalla,22 Vallabhan and Azene,18 Gallagher et al.,16 Sarigul and Gallagher,17 and Wieckowski.28; 29; 31 The Bell triangular element (e.g. Reference 3) with 18 degrees of freedom has been used by Azene,23 and Vallabhan and Azene.18 Rybicki and Schmit21 have used the rectangular element with polynomials of the 5th degree as shape functions. Gallagher et al.,16 and Sarigul and Gallagher,17 have used rectangular elements with the so-called ?blending function interpolants? in order to utilize the second derivatives as degrees of freedom at boundary nodes. Watwood and Hartz9 used triangular elements with linearly distributed self-equilibrating stresses. Unlike the above-mentioned works, where nodal parameters related to the Airy stress function are basic unknown variables, they introduced the so-called ?generalized stresses?. Then they applied some concepts of ?building block? in order to obtain equilibrium between the elements. One of these concepts leads to the same results as would be obtained by using the Hsieh?Clough?Tocher macro-element (see References 33 and 3) utilized in the present paper. In the present paper, the Bogner?Fox?Schmit rectangular and Hsieh?Clough?Tocher triangular elements are applied to approximate the Airy stress function. The Bogner?Fox?Schmit element (BFS, see References 32 and 3), shown in Figure 1, has 4 degrees of freedom at each corner, [F @F=@x @F=@y @2 F=@x@y]T . The Airy function is approximated by a polynomial of the order less than 4 with respect to both variables x and y P F= cpq xp yq (19) p; q63 In order to apply the Gauss quadrature, 16 (4 О 4) integration points are introduced. Figure 1. The Bogner?Fox?Schmit element Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS 1511 Figure 2. The Hsieh?Clough?Tocher macro-element Using equations (18) and (19), we obtain normal and tangential stresses distributed linearly and quadratically, respectively, along the element edges. The Hsieh?Clough?Tocher element originally introduced in Reference 33 (see also References 1?3) is a macro-element, subdivided into three triangles (see Figure 2) inside each of them the Airy stress function is a complete cubic polynomial ? on K1 ? F1 = c1 + c2 x + c3 y + и и и + c9 xy2 + c10 y3 2 3 (20) F = F2 = d1 + d2 x + d3 y + и и и + d9 xy + d10 y on K2 ? F3 = e1 + e2 x + e3 y + и и и + e9 xy2 + e10 y3 on K3 The element has 6 nodes: 3 at the corners and 3 in the mid-sides. 12 degrees of freedom are dened: three degrees of freedom, [F @F=@x @F=@y]T , at each corner node, and one, the normal derivative, @F=@n, at each mid-side node. The 30 coecients existing in the polynomial representation of function F (20) are calculated by solving the system of the following equations: the continuity conditions for function F and its rst derivatives, written for corner nodes P1 ; P2 and P3 (18 equations), and the central point, C (6 equations); the continuity conditions for normal derivatives of function F, written for external mid-side nodes P4 ; P5 and P6 , and internal mid-side nodes Q1 ; Q2 and Q3 (6 equations). In the present paper, the reduced HCT element is also utilized (see e.g. References 1 and 3). In this case, the normal derivative at each mid-side node is replaced by the average of the two normal derivatives at the ends of the corresponding side?this leads to reduction of numbers of nodes and degrees of freedom from 6 to 3, and from 12 to 9, respectively. The above system of algebraic equations is established for real positions of nodes. The mapping from the reference space to the real one, used for instance in case of isoparametric elements, is not applied here, because in general it does not preserve angular quantities, which leads to destruction of C 1 -continuity between elements resulting in the loss of equilibrium for tangential stresses. The Hammer numerical integration procedure (see e.g. Reference 39) with 3 points for each sub-triangle is applied for both types of elements: the original (HCT6) and reduced (HCT3) ones. The normal stresses are distributed linearly along the edges of both types of HCT element, while the tangential stresses are linear for the HCT6 element and constant for the HCT3 element. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1512 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON Let us consider the following discrete set of statically admissible elds of stresses, Yh ? Y , Yh = { ? [L2 ( )]4 : ? PK (Ki ) ?Ki ? Th ; = ; ; + b = 0 on ; n = t on } where PK (Ki ) is the space of polynomials dened on element Ki according to the structure of shape functions of the considered elements, Th denotes the triangulation of region . The discrete formulation of the considered problem can be written as follows: For n = 1; 2; : : : ; N , nd hn ? Kh such that hn = hn?1 + hn b(hn ; ? hn ) ┐ L(U n ; ? hn ) ? ? Kh h0 =0 (21) (22) (23) where Kh = Yh ? P. Using matrix notation, we can write the system of nite element equations for the equilibrium model of FEM in the form Z HT U dx ? F = 0 (24) analogous to that obtained for the displacement model, where ? ? @2 N1 @2 N2 и и и ? ? @y2 @y2 ? ? Z 2 2 ? ? @ N1 @ N2 T T ? H n U ds; H = ? F= и и и? ? 2 2 @x ? ? @x u ? ? @2 N 2 @ N2 1 ? иии ? @x@y @x@y ? ? ? x ? Ux n 0 ny U = y ; n = x ; U= 0 ny nx Uy ? ? xy (25) Equation (24) can be derived from (12) or (17) after representation of the Airy stress function by nodal parameters a ? [a1 a2 и и и]T and shape functions N ? [N1 N2 и и и]T , F = Na. Iterative solution For each time increment, tn , the following iterative procedure is employed: (i) Initialization of increments of deformations and stresses (at integration points) ijn(0) = 0; hn(0) = 0 (ii) For the ith iteration, calculation of the successive estimation of hn(i) ? Yh , solving the linear problem b(hn(i) ? hn(i?1) ; ? hn (i) ) = ?(n(i?1) ; ? hn(i) ) + L(U n ; ? hn(i) ) ? ? Yh Copyright ? 1999 John Wiley & Sons, Ltd. (26) Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS 1513 (iii) Calculation of increment of strains, n(i) , for each integration point, from the relations ijn(i) = ijen(i) + ijpn(i) n(i) ijen(i) = Cijkl hkl ( nij if ┐0 pn(i) ij = 0 if 60 where nij = 3 sij 2 q = eijp(i) nij n(i) eijp(i) = ijpn(i?1) + Cijkl (hkl ? [P(hn(i) )]kl ) (27) P is the operator of projection onto the yield surface. (iv) Checking the accuracy of the solution. If the assumed accuracy criterion is satised the calculation process can be stopped for the nth time increment, otherwise the next iteration should be performed starting from point (ii). The matrix of the linear system of nite element equations (26) is constant during the entire calculation process. The operator of stress projection onto the yield surface is dened as follows: (ij ? [P()]ij )Cijkl (kl ? [P()]kl )6(ij ? ij )Cijkl (kl ? kl ) ? ? B (28) and has easy form in case of the plane strain problem 1 0 [P()]ij = kk ij + sij 3 q while in the case of the plane stress problem the constrained minimization problem (28) has to be solved numerically. Following the concept described in Reference 40, the Newton?Raphson technique has been used to nd the minimum point of the augmented function ? 0 ) L(; ) = ( ? )C ( ? )?(q() where is the Lagrange multiplier. In the trivial case when ij = 0, the relation [P()]ij = 0 should be applied, in order to avoid the problem of non-uniqueness of projection onto the yield surface. The procedure for calculation of the increment of plastic strains is illustrated in Figure 3. The procedure is analogous to that commonly applied in the standard displacement nite element method, known as the closest-point-projection method or backward Euler method (see e.g. References 41 and 42). The above procedure and the CPP method give the same results in the case of boundary value problems, where stress elds are homogeneous. Equation (26) has the following representation in matrix notation Z Z T n(i) n(i?1) H C H dx (a ? a )= ? HT ?n(i?1) dx + Fn where terms of the matrix of elastic compliance, C, are dened in equation (4). Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1514 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON Figure 3. Calculation of increment of plastic strains Figure 4. Bending of rectangular plate, dimensions in centimeters NUMERICAL EXAMPLES Two numerical examples are considered: the bending of a rectangular plate, and the tension of a perforated plate. The application of the three described elements: BFS, HCT6 and HCT3 are shown in the case of the rst problem, while in the second one, only the triangular elements: HCT6 and HCT3 are utilized. The results are compared with those obtained by the displacement model of the nite element method with the help of 6-node isoparametric triangular elements, for which all stress components are approximated by linear functions as in the case of the HCT6 element. Bending of rectangular plate The plane strain problem of a rectangular plate is considered. The plate is subjected to the self-equilibrating loading according to Figure 4. Calculations have been made with the following material data (steel): Young?s modulus E = 2 О 1011 Pa, Poisson?s ratio = 0 и3, plastic limit 0 = 2и5 О 108 Pa. To stop the iteration process, the following convergence criterion has been checked: |q(i) ? q(i?1) | 60и01 |q(i) | R where q(i) = HT U(i) dx (see notation in equation (25)) in the case of equilibrium model of the nite element method, and q(i) is the vector of degrees of freedom calculated in the ith Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1515 STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS Table I. Lower and upper bounds of load limit Element size (m) 0и05 0и025 0и0125 Element type pmax (MPa) Dierence (per cent) pmax (MPa) Dierence (per cent) pmax (MPa) Dierence (per cent) BFS HCT6 HCT3 TRI6 102и4 106и8 100и4 117и4 12и777 9и029 14и480 ? 109и2 108и4 106и6 114и0 4и211 4и912 6и491 ? 111и6 111и2 110и0 112и6 0и888 1и243 2и309 ? iteration in the case of displacement-based nite element method, | и | the Euclidean length of a vector. The load increment has been set as p = 1 О 106 Pa, and reduced to p = 0и2 О 106 Pa for the load level close to the limit value. Calculations have been made for three regular meshes with elements having equal horizontal and vertical dimensions: 0и05; 0и025 and 0и0125 m. Because of the symmetry, only the left part of the plate has been analysed. The maximum values of load p, for which the required accuracy of the solution has been achieved, are shown in Table I. The abbreviation TRI6 is used in the table for the 6-node triangular element used in the displacement model of the nite element method. The values obtained by the stress-based model are the lower bounds for the limit load, while those obtained by the standard displacement model of FEM are the upper bounds. The table contains also the relative dierence between the bounds, calculated as follows: kin st ? pmax pmax kin pmax kin st where pmax and pmax denote bounds obtained by the displacement and equilibrium models of FEM, respectively. In the case of the nest mesh the dierences do not exceed 2и4 per cent, which is satisfactory from the engineering viewpoint. The slight dierence between the results obtained for BFS and HCT6 elements is observed. The HCT3 element provides a lower value for pmax due to the fact that the tangential stresses are constant along its boundary, while in the case of HCT6 and BFS elements they are linear and quadratic functions, respectively. The use of the method of Lagrange multiplier gives the possibility to calculate displacements for the boundary where stresses are prescribed (see the Appendix). Assuming that points A and B are xed in the vertical direction, the mean vertical displacement for the segment of boundary L (see Figure 4), located near the symmetry axis, has been calculated, and compared with the kinematically admissible displacement of the corresponding node. The comparison is shown in Figure 5 in the form of load?deection diagrams. The relative dierences between kinematically and statically admissible solutions are also shown in the gure for all elements used: BCT, HCT6 and HCT3. The dierences between displacements calculated by the two methods used are small at the beginning of the load?deection paths, and grow rapidly when the load approaches the limit. The BFS and HCT6 elements give similar results. The HCT3 element produces larger values of displacements, but for the nest mesh, the results are quite satisfactory, as the dierences are smaller than 1 per cent for the initial part of the path. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1516 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON Figure 5. Load?deection diagrams Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS 1517 Figure 6. Plastic zones Table II. Computation time in seconds (Silicon Graphics MIPS R5000 Indy A2) Element type BFS HCT6 HCT3 TRI6 Element size (m) 0и05 0и025 0и0125 5и263 6и406 6и135 1и031 19и438 25и344 19и666 3и609 85и324 102и825 88и375 14и427 The comparison of plastic zones obtained with the help of the three self-equilibrating elements with the kinematically admissible solution is shown in Figure 6 for three levels of load p. The similarity between the results can be observed. It is worth mentioning that the stress-based model appears to require approximately twice the number of iterations and a larger number of integration points than that used in the displacement model. The computation times for the self-equilibrating elements and the TRI6 element are compared in Table II for the rst 30 load increments. The comparison shows that the stress-based model of FEM is more time consuming than the displacement-based one. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1518 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON Figure 7. Geometry of perforated plate, dimensions in millimeters Figure 8. Exemplary meshes Tension of perforated plate The plane stress problem of the perforated plate shown in Figure 7 is considered. Unlike the previous example, the load on the plate is displacement-controlled?it is assumed that the horizontal displacements, U1 = ▒U , on the vertical sides of the plate are given. The considered material data corresponds to aluminium: Young?s modulus E = 7 О 1010 Pa, Poisson?s ratio = 0и2, plastic limit 0 = 2и34 О 108 Pa. Three triangular element meshes, shown in Figure 8, are used in the analysis. The values of maximum element sizes for these meshes are as follows: (a) 3и86 mm, (b) 2и09 mm, (c) 1и12 mm. As the plate is doubly symmetric, only the right upper quarter of its area is analysed. The increment of given displacement is set as U = 2 О 10?6 m. The same accuracy criteria are applied as in the previous example. The results of calculations are shown in the form of P?U paths presented in Figure 9, where P denotes the resultant force calculated for half height of the plate. As in the previous example, the smaller values of forces, P, are obtained for the HCT3 element than for the HCT6 one. The dierences between the load limits, calculated by the two nite element models, are not greater than 1и5 per cent for the nest mesh. The growth of plastic zones is presented in Figure 10, where their shapes are compared for three values of displacement U . Similar plastic zone shapes are obtained for the three types of elements used. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS 1519 Figure 9. P?U paths for perforated plate CONCLUSIONS The equilibrium model of the nite element method is applied to plane, quasi-static problems of plasticity theory. The statically admissible elds of stresses are constructed by means of the Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1520 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON Figure 10. Plastic zones for perforated plate Airy stress function. Three kinds of elements are utilized to approximate the stress function: the Bogner?Fox?Schmit rectangular element, and the Hsieh?Clough?Tocher triangular element in its original and reduced versions. The stress boundary conditions are fullled using the Lagrange multiplier method, allowing the evaluation of the displacements on the boundaries with prescribed tractions. The concept of multipoint-constraints elements is proposed; this leads to an easier implementation of the method as well as the easier use of a computer program. The proposed iterative procedure, for solution of the non-linear system of equations for each load increment, is analogous to the closest-point-projection method routinely utilized in the displacement nite element method. The used procedure and the CPP method give the same results in case of boundary value problems, where stress elds are homogeneous. The constant matrix of the set of algebraic equations is applied in the procedure. The equilibrium model of the nite element method provides a lower bound for the load limit which is important from the design perspective. However, some eort is to be paid?this model is more dicult to implement than its displacement counterpart, because of necessity of using elements of C 1 class and multi-point-constraints imposing the appropriate boundary conditions. The equilibrium model usually need the longer computation time than the displacement one when the similar numbers of elements are applied to the specic problem in both the models of FEM. The concepts, described in the present paper, somewhat soften these two drawbacks, as they make it possible to implement the stress-based model in the standard nite element program without special extensions. The numerical examples show the convergence of the equilibrium model and good agreement with the results obtained by the use of the standard displacement nite element method. Small dierences between the lower and the upper bounds are obtained for load limits. Very good agreement for displacements and shapes of plastic regions are also observed. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1521 STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS Table III. Characteristics of elements Distribution of stresses Element type Normal BFS HCT6 HCT3 Linear Linear Linear Tangential Number of Lagrange multipliers for MPC element Number of d.o.f. in MPC element Quadratic Linear Constant 13 11 9 5 4 3 APPENDIX Imposition of stress boundary conditions Using notation shown in equation (25), we can express the boundary conditions in the form of the matrix equation nHa = t (29) where t = [tx ty ]T is the vector of given tractions. Equation (29) indicates that the multi-point constraints (MPC) problem is to be solved unlike in the displacement formulation of the nite element method, where, for simplest cases, the nodal displacements can be prescribed directly. In the present paper, equation (29) is satised with the help of the method of Lagrange multiplier. This technique leads to increasing the number of unknowns in the nal system of algebraic equations, but enables calculation of displacements for points belonging to boundary , where tractions are given. This follows from the fact that Lagrange multipliers are interpreted as displacements of boundary points when used to impose equation (29). Using (29), we obtain the same distribution of normal stresses and dierent distribution of tangential stresses for the three elements utilized in the paper (see Table III). In the considered nite element model, only the boundary conditions of the forms given in Table III can be satised identically. In other cases, they should be approximated by functions given in the table. In order to facilitate the implementation of the Lagrange multiplier technique as well as handling the boundary conditions by a user of the computer program, the concept of MPC-elements is used. The concept is described for HCT3 element. Using the local coordinate system (n; s) for the edge of the element, shown in Figure 11, the normal and tangential stresses can be written as functions of s, and are dependent on degrees of freedom dened at nodes P1 and P2 only 6 2 (2a1 + la3 ? 2a4 + la6 ) s ? 2 (3a1 + 2la3 ? 3a4 + la6 ) l3 l 1 T (s) = (a2 ? a5 ) l N (s) = (30) where a1 ; : : : ; a6 are terms of the vector of degrees of freedom dened at nodes P1 and P2 aT = [F(P1 ) @F=@n(P1 ) @F=@s(P1 ) F(P2 ) @F=@n(P2 ) @F=@s(P2 )]T If normal and tangential stresses are prescribed by means of constants A, B and C as follows: N (s) = A s + B T (s) = C Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1522 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON Figure 11. MPC-element for HCT3 element it follows from (30) that the boundary conditions are satised for arbitrary s if the following system of equations is fullled: ? ?a ? ? ? 1? 0 6l?2 ?12l?3 0 6l?2 12l?3 ? ? ? ? ? ? a2 ? ? ?A? ? ? ?1 ?2 ?1 ? ? ?6l?2 (31) 0 ?4l 6l 0 ?2l . ? ? . ? = ?B? ? . ? ? C ? ? ?1 ?1 ? ? 0 l 0 0 ?l 0 a6 or briey Ga ? g = 0 (32) The HCT3 element requires to introduce 3 Lagrange multipliers for each side located on in order to satisfy equation (32). Let us dene MPC element corresponding to HCT3 element. Let it consist of three nodes P1 , P2 , and the additional one, P3 , at which the auxiliary three degrees of freedom are dened?3 Lagrange multipliers, associated with equations (31) (see Figure 11). Let us dene the vector of degrees of freedom for this element, ae , and the following element square matrix, ke , and the vector of right-hand-side terms, f e , ? ? F(P1 ) ? ? ? ? ? ? ? ? ? ? ? ? @F ? ? ? ? (P ) 1 ? ? ? ? ? ? @n ? ? 0? ? ? ? ? ? ? ? ? ? ? ? ? ? ? @F 0? ? ? ? ? ? ? ? ? ) (P 1 ? ? ? ? ? ? ? ? @s 0 ? ? ? ? ? ? ? ? F(P ) ? ? ? ? ? ? ? 2 ? ? ?0? ? T 0G e e e ; f = 0 ; k = a = G 0 ? @F ? ? ? ? ? ? 0? ? ? ? (P2 ) ? ? ? ? ? ? ? ? ? @n ? ? ? A? ? ? ? ? ? ? ? ? ? ? ? ? @F ? ? ? ? B ? ? ? ? (P ) ? ? ? 2 ? ? ? ? ? ? ? @s C ? ? ? ? ? ? ? ? 1 ? ? ? ? ? 2 ? ? ? ? ? ? ? 3 Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS 1523 Then treating matrices ke and f e of MPC element, during the assembly process for global matrices K and F; in the same way as exibility matrices and right-hand-side vectors of triangular elements, we impose the boundary conditions (32). Before assembling, matrix ke should be transformed to the global co-ordinate system kge = TT ke T where ? T T= ? 0 0 0 T 0 ? 0 0? ; I ? 1 T = ? 0 0 ? 0 ny ? ; nx 0 nx ?ny ? 1 I = ?0 0 0 1 0 ? 0 0? 1 The element matrices and vectors for MPC elements associated with HCT6 and BFS elements can be dened in a similar way. The numbers of Lagrange multipliers and degrees of freedom of all MPC elements used are summarized in Table III. Lagrange multipliers N0 and T0 (N0 = 2 , T0 = 3 in case of HCT3 element) corresponding with equations written for terms containing s0 in expressions for N (s) and T (s) can be interpreted as Z l Z l uN ds and uT ds 0 0 respectively, so the average normal and tangential displacements for each segment can be calculated from the equations av = uN N0 l and uTav = T0 l The above interpretation of multiplier N0 is utilized in the rst numerical example. The use of Lagrange multiplier technique means that zero terms appear on the diagonal of the element matrix ke . Therefore, the pivoting or the proper ordering of degrees of freedom for the global system of equations should be used. The easiest way is to group all additional nodes, where Lagrange multipliers are degrees of freedom, at the end of the list of nodes, but this measure leads to uneconomical growth of length of prole or half-width of the equations matrix. In the present paper, a simple solution, consisting of two steps, is used. In the rst step, the well-known Cuthill? McKee43 algorithm is applied to renumber only basic nodes related to variables dening the Airy stress function. Then all numbers of the additional nodes are inserted to the list of all nodes in such a way that the number of each additional node is set to n + 1, where n denotes the largest number of a basic node belonging to the MPC element related to the considered additional node. ACKNOWLEDGEMENTS The rst author wishes to thank the Korea Science and Engineering Foundation for the nancial support in the framework of the Post-doctoral Fellowship Program for Foreign Researchers. REFERENCES 1. O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, McGraw-Hill, London, 1989. 2. R. H. Gallagher, Finite Element Fundamentals, Prentice-Hall, Englewood Clis, N.J., 1975. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) 1524 Z. WI E VCKOWSKI, S.-K. YOUN AND B.-S. MOON 3. P. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. 4. M. Kleiber, Incremental Finite Element Modelling in Solid Mechanics, PWN=Ellis Horwood, Warszawa, 1989. 5. B. Fraeijs de Veubeke, ?Displacement and equilibrium models in the nite element method?, in O. C. Zienkiewicz and G. S. Holister (eds.), Stress Analysis, Wiley, New York, 1965. 6. B. Fraeijs de Veubeke and O. C. Zienkiewicz, ?Strain-energy bounds in nite-element analysis by slab analogy?, J. Strain Anal., 2, 265 ?271 (1967). 7. L. S. D. Morley, ?The triangular equilibrium element in the solution of plate bending problems?, Aero Quart., 19, 149?169 (1968). 8. Z. M. Elias, ?Duality in nite element method?, Proc. ASCE, J. Engng. Mech. Div., EM4, 94, 931?946 (1968). 9. V. B. Watwood and B. J. Hartz, ?An equilibrium stress eld model for nite element solution of two-dimensional elastostatic problems?, Int. J. Solids Struct., 4, 857?873 (1968). 10. B. Fraeijs de Veubeke and M. Hogge, ?Dual analysis for heat conduction problems by nite elements?, Int. J. Numer. Meth. Engng., 5, 65?82 (1972). 11. R. Glowinski, J. L. Lions and R. Tremolieres, Analyse Numerique des Inequations Variationelles, Dunod, Paris, 1976. 12. Z. Wi eVckowski, ?Dual nite element analysis for plasticity-friction torsion of composite bar?, Int. J. Numer. Meth. Engng., 38, 1901?1916 (1995). 13. J. J. Connor and C. A. Brebbia, Finite Element Technique for Fluid Flow, Newnes-Butterworths, London, 1976. 14. A. R. S. Ponter, ?The application of dual minimum theorems to the nite element solution of potential problems with special reference to seepage?, Int. J. Numer. Meth. Engng., 4, 85?93 (1972). 15. R. L. Taylor and O. C. Zienkiewicz, ?Complementary energy with penalty functions in nite element analysis?, in R. Glowinski (ed.), Energy Methods in Finite Element Analysis, Wiley, New York, 1979. 16. R. H. Gallagher, J. C. Heinrich and N. Sarigul, ?Complementary energy revisited?, in Atluri, Gallagher and Zienkiewicz (eds.), Hybrid and Mixed Finite Element Methods, Wiley, Chichester, 1983. 17. N. Sarigul and R. H. Gallagher, ?Assumed stress function nite element method: two-dimensional elasticity?, Int. J. Numer. Meth. Engng., 28, 1577?1598 (1989). 18. C. V. G. Vallabhan and M. Azene, ?A nite element model for plane elasticity problems using the complementary energy theorem?, Int. J. Numer. Meth. Engng., 18, 291?309 (1982). 19. R. H. Gallagher, ?Can and should the force method be revisited??, in S. N. Atluri, G. Yagawa and T. A. Cruse (eds.), Proc. Int. Conf. Comp. Engng. Sci:; July 30?August 3; 1995; Hawaii; USA, Springer, Berlin, 1995, pp. 2?7. 20. C. Balakrishna, J. H. Kane and R. H. Gallagher, ?Assumed stress function nite element method: Towards unication?, in S. N. Atluri, G. Yagawa and T. A. Cruse (eds.), Proc. Int. Conf. Comp. Engng. Sci.; July 30?August 3, 1995, Hawaii, USA, Springer, Berlin, 1995, pp. 1578?1583. 21. E. F. Rybicki and L. A. Schmit, ?An incremental complementary energy method of nonlinear stress analysis?, AIAA J., 8, 1805?1812 (1970). 22. R. H. Gallagher and A. K. Dhalla, ?Direct exibility nite element elastoplastic analysis?, Proc. Int. Conf. on Structural Mechanics in Nucl. Reac. Tech., 1972. 23. M. Azene, ?A nite element complementary energy formulation for plane elastoplastic stress analysis?, Ph.D. dissertation, Dept. of Civil Engineering, Texas Tech. Univ., Lubbock, Texas, May 1979. 24. T. Belytschko and P. G. Hodge, ?Plane stress limit analysis by nite elements?, Proc. ASCE, J. Engng. Mech. Div., EM4, 96, 931?944 (1970). 25. B. Tabarrok and D. Sodhi, ?On the generalization of stress function procedure for the dynamic analysis of plates?, Int. J. Numer. Meth. Engng., 5, 532?542 (1972). 26. B. Tabarrok, ?Complementary variational principles in elastodynamics?, Comput. Struct., 19, 239?246 (1984). 27. B. Tabarrok and A. Simpson, ?An equilibrium nite element model for buckling analysis of plates?, Int. J. Numer. Meth. Engng., 11, 1733?1751 (1977). 28. Z. Wi eVckowski, ?Duality in nite element method and its applications to some linear and non-linear problems of mechanics of composite materials?, in Polish, Ph.D. dissertation, Institute of Civil Engineering, Technical University of Lodz, Poland, February 1987. 29. Z. Wi eVckowski, ?Dual nite element methods in mechanics of composite materials?, J. Theoret. Appl. Mech., 33, 233?252 (1995). 30. Z. Wi eVckowski, ?Dual nite element methods in homogenisation for elastic-plastic brous composite material?, in C. K. Choi, C. B. Yun and H. G. Kwak (eds.), Proc. 7th Int. Conf. on Computing in Civil and Building Engineering, Seoul, Korea, 19?21 August 1997, pp. 1051?1056. 31. Z. Wi eVckowski, ?Dual nite element method in friction problems?, J. Theoret. Appl. Mech., 30, 535?543 (1992). 32. F. K. Bogner, R. L. Fox and L. A. Schmit, ?The generation of interelement compatible stiness and mass matrices by the use of interpolation formulas?, Proc. Conf. Matrix Methods in Struct. Mech., October 1956, Wright Patterson A. F. Base, Ohio, 1965. 33. R. W. Clough, J. L. Tocher, ?Element stiness matrices for analysis of plate bending?, Proc. Conf. Matrix Methods in Struct. Mech., October 1956, Wright Patterson A. F. Base, Ohio, 1965. 34. G. Duvaut and J. L. Lions, Les inequations en mecanique et en physique, Dunod, Paris, 1972. 35. R. Hill, ?A comparative study of some variational principles in the theory of plasticity?, J. Appl. Mech., 17, 64?70 (1950). Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999) STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS 1525 36. R. A. Adams, Sobolev Spaces, Academic Press, New York, 1975. 37. R. Temam and G. Strang, ?Functions of Bounded Deformation?, Arch. Rat. Mech. Anal., 75, 7?21 (1980). 38. C. Johnson, ?Mathematical and Numerical Analysis of Some Problems of Plasticity?, in Telega, J. J. (ed.), Functional Analysis Methods in the Theory of Plasticity, Ossolineum, Wrocl aw, 1981, pp. 71?109 (in Polish). 39. G. Dhatt and G. Touzot, The Finite Element Method Displayed, Wiley, Chichester, 1984. 40. A. Samuelsson and M. Froier, ?Finite elements in plasticity?a variational inequality approach?, in J. R. Whiteman (ed.), MAFELAP 1978 (III ), Academic Press, London, 1979. 41. J. C. Simo and R. L. Taylor, ?Consistent tangent operators for rate-independent elastoplasticity?, Comp. Meth. Appl. Mech. Engng., 48, 101?118 (1985). 42. M. A. Criseld, Non-linear Finite Element Analysis of Solids and Structures, Volume 1: Essentials, Wiley, Chichester, 1991. 43. A. George and J. W. Liu, Computer Solution of Large Sparse Positive Denite System, Prentice-Hall, Englewood Clis, N.J., 1981. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)

1/--страниц