Composite finite elements for rigid-plastic analysis

Journal of Materials Processing Technology, 27 ( 1991 ) 55-71

55

Elsevier

Composite finite elements for rigid-plastic analysis Tatsuhiko Aizawa and Junji Kihara Department o[ Metallurgy, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113, Japan

Industrial Summary In the conventional rigid-plastic finite element analysis, quadrilateral elements or brick-type elements are widely used for two- or three-dimensional simulations in rolling or forging. Through some mathematical studies and numerical tests of errors, the quadrilateral element family, especially the 4-node element with bilinear interpolation for incremental displacement or velocity and constant for pressure, is found to be unsuitable for elasto-plastic or rigid-plastic analyses with large deformations; this type of element often causes large local errors for shear components of strains and stresses, especially in stress concentration problems or in plastic forming problems with local distortion of elements. An alternative family of elements will be proposed and discussed for two- and three-dimensional rigid-plastic simulation in rolling or forging. Since the LBB condition for incompressibility is imposed on the employed velocity interpolations in construction of the adequate finite elements models, the use of conventionaltriangular elements is prohibited in the rigid-plastic plane strain, axisymmetric or three-dimensionalanalyses, where the incompressibilitycondition of plastic strain or strain rate or the volume constancy should be explicitly evaluated by the variational equation or the weak form. Our approach is based on the mixed finite element method where both velocity and pressure are interpolated to satisfy the LBB condition and to become robust for local element distortion. In the present paper, a list of candidate mixed type finite element pairs other than the four-node element will be shown with some comments on their characteristics and features. Special attention will be paid to the composite-type finite element schemes for two- and three-dimensional problems; in particular, the inner-mode composite element family is described with respect to their fundamental behavior in rigid-plasticity. In principle, our developing elements are free from shear/membrane locking and from contamination by zero-energy modes, and stable and robust against the folding behavior of the corner points. Furthermore, the introduced pressure can be excluded from the final finite element equations by the static condensation or the penalty function method. Through the basic numerical tests for upsetting problems, accuracy and convergency will be evaluated for these elements. Finally, some further examples on the complete automatic mesh control will be employed to show the potential of these elements for adaptive mesh control.

1. Introduction

In the conventional rigid-plastic finite element analysis, the quadrilateral elements or the brick-type elements have been, and still are, widely used for 0924-0136/91/$03.50 © 1991--Elsevier Science Publishers B.V.

56 two- or three-dimensional simulations in rolling or forging [ 1,2 ]. Through some mathematical studies and numerical tests of discretization errors, those quadrilateral element families, including the four-node element with bilinear interpolation for incremental displacement or velocity and constant for pressure, are unsuitable for elasto-plastic or rigid-plastic analyses with large deformation: ( 1 ) These elements never satisfy the LBB condition which is a necessary and sufficient condition for finite elements to provide a well-posed and balanced pair of velocity and pressure under the constraints of incompressibility [3-5]; (2) the elements often cause local large errors for shear components of strains and stresses especially in stress concentration problems even for linear analysis [6,7 ]; (3) they suffer from local large distortion and shear/membrane locking behaviors [8]; and (4) with these elements it is difficult to control mesh generation and refinement in a completely automatic manner [9]. Hence, ne{v finite element design is indispensable to further advance the accuracy and reliability of plastic forming simulations. An alternative family of finite elements must be proposed and developed for two- and three-dimensional rigid-plastic simulation in rolling or forging: the bubble-type elements are utilized by Bathe et al. [10] in large-strain/deformation analysis, several candidate elements are compared for accuracy by Zienkiewicz et al. [11] and Yamada [12], and quadratic elements are used in rigid-plastic analysis by Kikuchi et al. [ 13 ]. The authors [ 14,15 ] are particularly concerned with a new type of triangular and tetrahedral elements for this kind of problems, since their logical compatibility makes them preferrable to the other new elements mentioned above. It has to be noted that, due to the above mentioned LBB condition for incompressibility, such conventional triangular elements as the three-node element with linear interpolation for velocity and constant for pressure have been prohibited to be used in rigid-plastic plane strain, axisymmetric and three-dimensional analyses; if the incompressibility condition of plastic strain or strain rate, or the volume constancy, should be explicitly evaluated in the variational equations or in the weak form for those elements, severe locking behaviors could be experienced, such as meshlocking [ 16 ]. Our approach is originally based on the mixed finite element method where both velocity and pressure are appropriately interpolated to satisfy the LBB condition and to become robust for local element distortion. In the present paper, a list of candidate mixed-type finite element pairs other than the fournode element will be shown with some comments on their characteristics and features. Special attention will be paid to the composite-type finite elements for two- and three-dimensional problems: in particular, the inner-mode composite element family is described with respect to their fundamental behaviors in rigid-plasticity. In principle, our developing elements have the following features: (a) they are free from shear/membrane locking behavior and contamination of zero-energy modes; (b) they become robust and stable against

57 the folding behavior of the corner points; (c) the employed pressure terms in these elements can be eliminated from the final finiteelement equations by the staticcondensation or the penalty function method; (d) the element models can be easily generated and controlled by completely automatic handling. Through basic numerical tests for upsetting problems, accuracy and convergency willbe both evaluated for these elements. Finally,some further examples on the completely automatic mesh generation for these elements will be employed to show the compatibility of these elements to element control. 2. Mixed-type variational principle in rigid-plasticity

To understand the theoretical framework of nonlinear analyses, including rigid-plastic analysis, the mixed-type variational principle becomes an essential base where the equilibrium of stresses is taken into account together with the kinematics and the incompressibility condition. Consider a three-dimensional continuum ~ in rigid-plasticity with boundary S ( --S~u Su), as shown in Fig. 1. No body force is applied in Q, and the traction and the velocity are prescribed on So and Su, respectively. In the presence of Su, the rigid part is represented by a small plastically deforming zone with a small tolerance for the equivalent strain rate. Then, the current state of stress, strain rate and velocity field in the plastic body is described by the following mixed type variational equation:

Ti

f ~5~-dV+f P~iidV=~soTi~uidS

(1)

--ftii ~ipdV=O So

x

y

Fig.1.A continuumin rigid-plasticity.

ui

58 where the admissible velocity vi is employed from the subspace V of H 1 satisfying the geometric or essential boundary conditions and the pressure p from the subspace Q of L 2. In the absence of Su, such subsidiary conditions as an equilibrium condition of external and internal forces or an additional velocity condition should be evaluated since Vhas constant ambiguity. In general, both V and Q should be selected in order that the trial function pair (u*, p* ) should be measured for the actual solution (u, p) in the form of

(2)

IIu*--ulIH, + I P * - P l 1~ ~ e

by using the related norms. In other words, there should be some constraints on the selection of V and Q for the weak solution (u*, p*) to converge to the strong solution (u, p) which satisfies the original system of governing equations. This is the LBB condition in the continuum model: infsup

f~ div (u)q dV

>~fl ( > 0 )

(3)

Even in the discretized modeling or the finite element modeling, where (u*, p*) for u*eVh and p*eQh should be constructed by V h c V c H 1 and Qh c Q c L 2, the above condition must be observed to provide the rational finite element solution pair (Uh, Ph) : inf

sup

f a div(u*)q* dV

~*~qh~*~Vh Ilu*llv. lq*l~h

>~fl*

(>0)

(4)

If this condition were violated, both convergency and approximation of the current finite elements should be lost in principle and there should be some possibility of occurrence of locking or zero-energy contamination. Hence, the finite element must be designed so as to satisfy the above LBB condition. It should be noted that this constraint is common to the mixed-type formulations in general; since the incompressibility conditions should be explicitly evaluated in such nonlinear problems as listed in Table 1, the adequate finite element pairs must be selected even in the discretized formulations. When the pressure is approximated by a piecewise constant, say p =Pe, in the selected finite element, the original approximate variational equations for the current element Ve are reduced by the penalty function method:

fv~5~-dVq-fvP~iiidV=f%,~vTiSuidS

-fvdii

8p~ dV=O

First, the second equation of eqn. (5) is replaced by

(5)

59 TABLE 1 Abstract forms in materials nonlinearity Incompressible elasticity

Elasto-plasticity

Rigid-plasticity

Analytical condition

v = 0.5 or v ~ 0.5 Rubber-likematerials

Large deformation when the incremental strains are nearly plastic components

Plane strain condition Axisymmetric condition Three-dimensional state

Variational equation

f aS ~j8%d V + laP 8tii d V = fs~ TiSul dS f,'~ii~P d V= 0

f c~Sij8% d V + f~PSeii dV

f ndSt d V + f,-~P5¢:iidV

Abstract form

= fsoTiSui dS =Jso T,6hi dS f ~ i,Sp d V f a~i,Sp dV=0 1 = - f,~Sp dV K a(u,v)+b*(v,p)=(f,v) VveV b(u,q)- (p,q)=O ¥qeQ

lye(aPe --¢:ii) ~)P d V = 0

(6)

where a is a small positive factor, which is called the penalty factor. Obviously, as a goes to zero, eqn. (6) converges to the original equation of eqn. (5). Since both Pe and 8pe are piecewise constant in an element, Pe can be rewritten by 1 Pe=

~V e

f

Ve

~ii d V

(7)

Substituting eqn. (7) into the first equation of eqn. (5), we have the finite element variational equation to determine the velocity field:

f

ve ~ 8 ~ - d V +

~

1 \(of v o ~ i i d V ) 8 ( f

v~ ~ i i d V ) = f s ~,v~

T~Su~dS

(8)

The pressure field can be obtained by the so calculated strain rates through eqn. (7) in a piecewise manner. It should be noted that the resulting finite element equations from eqn. (8) become nonlinear in vi; accurate time integration together with an iteration scheme and convergence condition should be employed to provide a reliable solution. 3. Element design for rigid-plastic analysis Both for the infinitesimal and the finite rigid-plastic finite element analysis, the incompressibility condition must be adequately evaluated by the employed

60 finite element scheme and the related variational principle. Owing to the mathematical theories for the above mentioned mixed-type variational approach, the LBB condition in general provides a sufficient assessment base in the element design for rigid-plastic analysis. In this situation, two modes of finite element families are recognized as the candidate mixed finite elements which satisfy the LBB condition: [A] continuous-pressure mixed-type finite element schemes; and [B] inner-mode pressure mixed-type finite element schemes. The former set of elements is listed in Fig. 2: the corner nodes share the same nodal velocity and pressure with the neighboring elements. Since the static pressure varies greatly in the local plastic concentration, these types of elements are never suitable for solid analysis with incompressibility. Figure 3 shows various elements where a piecewise constant or linear pressure is as-

Fig. 2. Continuouspressure mixed-typefinite elementschemes.

61

Fig. 3. Inner-mode pressure mixed-type finite element schemes.

sumed independently of each element. The nonconforming element [B1 ] with only middle nodes causes severe shear locking, while [B7] and [B8] are dependent on the mesh alignment. Hence, feasibility tests are necessary to evaluate the adequate elements. In the present study, four elements [B2], [B4], [B5 ] and [B6 ], as shown in Fig. 4, are selected in the feasibility test: [B2 ] has a quadratic interpolation for velocity with constant pressure, [B41 a quadratic function with one cubic function of bubble mode for velocity and a linear pressure, [B5 ] is a composite triangular element, and [B6 ] has a biquadratic function with linear pressure. Besides the LBB condition, the number of degrees of freedom for the velocity field in these elements should be superior to that for the pressure field. In other words, the ratio of the number of degrees of freedom indicates the logical balance of the interpolation pair for both velocity

62 !

I Uh

Quadratic

ph

Constant [B2]

i

I

Quadratic

Composite

Linear

Constant

[84]

:

Quadratic

Constant [B6]

[BS]

(

/

/" (")

/

• /

•

,

/ i

//

~

Fig. 4. Feasible mixed-type finite element schemes. TABLE 2 Ratio of degrees of freedom for velocity and pressure interpolations Type of elements

Ratio of DOF for vi and p

[B1] [B21 [B31 [B4] [B5] [B6] [B7]

3 4 4/3 2 4 8/3 S

[BS]

2

and pressure; i.e. the ratio r of the number of degrees of freedom for the velocity interpolation to the number of constraints or the number of pressure equations. Table 2 shows the ratios r for the above eight elements, from [B1] to [B8 ], when two edges of one rectangular element of paired triangular elements are completely fixed around a single corner. Through empirical observation it is found that: { 1 ) r should be larger than or equal to the number of dimensions (in the two-dimensional situation, r>~ 2), and (2) excess degrees of freedom are indispensable to cope with folding, complex corners or other complex boundary conditions appearing in practice. Now, an adequate benchmark test should be recommended to understand the basic characteristics of the element pairs. The uniaxial compression test is taken to make a fine comparison of rigid-

63 i

V

D

:1

(re=O)

i i

L C (m=0)

2 A

0

D

0

C D

B

A

B

C Fig. 5. Uniaxial compression of work materials to compare several candidate finite element schemes. Quadrilateral

I

Quadratic/Constant

Quadratic/Linear

Composite/Constant

Quadratic/Linear

\l -

\\

Fig. 6. Comparison of numerical behavior in uniaxial compression (bench-mark testing of candidate element schemes).

64 Interpolation

Element Subregion 1 •

! /

\

//

//(~ \

/' /

o 0

/

/

for each s m a l l e r t r i a n g u l a r

'~

,, FD i:

subregion 't

\

Vi = LIVI + Lmv m + Lnvn LI, Lm, Ln: area coordinates

Ill

for each l a r g e r t r i a n g u l a r

subregion

or element #i

\ \ /

~ 1

P=P~ 3

( Before deformation )

( piecewise constant )

( After deformation )

Fig. 7. Elementdesignof candidatecompositetriangularelementschemes. plastic responses. As shown in Fig. 5, a circular disc with radius 2.0 and height 1.0 is compressed by uniform velocity at the fixed end A-O. Nearly the same degrees of freedom are employed in the discretized modeling. As weUknown, the conventional four-node element suffers from severe locking behavior at the corner node O: as shown in Fig. 6, elements at the vicinity of the corner node are forced into a triangular shape and any further deformations are prohibited around this node. Figure 6 shows the deformed shape of elements at 50% reduction. No locking phenomena can be seen in all the results. However, the limit of reduction to which the analysis can be performed becomes the lowest for element [B2]. The bubble elements, [B4] and [B6], are found to be stable in various problems; as predicted through theoretical studies, use of higher interpolations for velocity leads to preferable stabilization at large de-

65

formation. To be noted here is the effect of the corner point: since the velocity changes drastically in the vicinity of this point, some disturbances should occur in each employed element. In case of bubble elements, the obtained free boundary O-D is shaped swinging, while relatively smooth boundaries can be attained by using the composite element [B5 ], due to piecewise linearity in velocity field. 4. Composite triangular/tetrahedral elements

Let us now describe precisely the [B5 ] -type elements or the composite triangular and tetrahedral elements. Figure 7 (a) shows the element subregions and the related interpolation functions: a larger triangular region or an element is subdivided into four smaller subregions. In each smaller subregion, the velocity is interpolated by a linear function, so that the strain rate becomes discontinuous with finite jumps across the subregion boundaries in each element. On the other hand, the pressure is assumed to be piecewise constant in the larger triangle. The features of this element are summarized as follows: (1) Although the initial triangular element is shaped regular, the deformed configuration is allowed to be star-type hexagonal.

\

Element Subregions

Interpolation for each smaller tetrahedral subregion;

v=L~v] + Livi + LkVk + hvl L,, ~, Lk, ~:

¢J

o

2

~7

Volumetric coordinates

j ~ k

for each larger tetrahedral subregion or element #i,

P----Pi ¢D

==

(piecewise constant)

D..

Fig. 8. Element design of candidate composite tetrahedral element schemes.

66

(2) Even if the smaller triangle is distorted as seen in Fig. 7 (b), this element works well until one of the smaller regions completely diminishes (the accuracy of this element is strongly dependent on this type of local distortion ). (3) In other triangular elements, one edge should be in contact with the tool surface in folding, while this kind of contact behavior can be traced by checking the contact conditions of smaller element edges with the tool surface. The above characteristics will be discussed by numerical results. Although some further considerations and tests are still necessary to assure the adaptivity and feasibility of tetrahedral composite elements, the present candidate element is shown in Fig. 8.

\\\\ o

o o

\\\\_

o

2.00

0.00

4.00

0.00

2.00

4.00

2.00

4.00

o

o.

o.

c-J-

o~

I

0

21.00

.00

O

(:3

I

4.00

0.00

....

OI

I

0.00

I

2.00

I

I

I

4.00

Fig. 9 (a). Deformed shape of composite elements in upsetting. Case [A-1 ]: analytical region is triangularized from the upper left to the lower right.

67

5. Numerical examples

The upsetting problem is first taken as a numerical example: a cylindrical billet with dimensions 8.0 X 8.0 is upset with a uniform velocity with the fixed end. Two types of element alignment are employed in simulation: [A] The whole analytical region is triangularized from the upper left to the lower right ( [A-l] ) or from the upper right to the lower left ( [A-2] ); and [B] The triangular alignment is controlled in the near-Union-Jack manner (only the region around the corner is subdivided into elements from upper right to lower left). The deformed shapes with reduction are illustrated for type [A] in Fig. 9(a) and (b); in case of [A-l], an element is collapsed by folding, while a smaller subregion becomes weak at the center in case of [A-2]. In the former case, mesh alignment must be modified at the vicinity of corners; finer elements are recommended at the strain concentration area to cope with the latter affairs• In numerical analysis for case [B ], the volume constancy condition is verified at each reduction: even at the final reduction, the incompressibility condition is satisfied within 0.5% tolerance. As shown in Fig. 9 (c), the folding behavior is accurately traced without locking; even if an element becomes staro o

.qp"

/// //// o

I

0.00

2.00

]

I

0.00

4.00

I

I

2.00

4.00

t~

0 D

C)

I

0.00

2.00

4.00

°0.00

I

2.00

I

I

I

4.00

Fig. 9 (b). Deformed shape of composite elements in upsetting. Case [A-2 ]: analytical region is triangularized from the upper right to the lower left.

68

\\// o

\\\\ I

°0.00

2.00

0.00

2.00

i

4.00

i

I

i

o

4.00

i

I

2.00

0.00

i

I

4.00

i

I

I

I

0,00

2.00

4.00

6.00

0.00

2.00

4.00

6.00

C3

R

o

I

C3

.00

I

2.00

I

I

[

4.00

" 0-

Fig. 9 (c). Deformed shape of composite elements in upsetting. Case [B ]: mixed-type triangulation.

type or some smaller subregions are distorted, the current element is still stable enough to continue the present computation. 6. Discussion and conclusion

Very few researches have been performed to propose and develop more adequate and effective element schemes than the conventional four-node quadrilateral element. This is due to the fact that many practical rolling and forging problems can be dealt with by using those conventional elements with the aid of some tactical algorithms; the use of triangular elements in the transient stages with reduction around the corner node helps the four-node element from locking. Those conventional approaches have been successful in relatively accurate simulations. However, more reliability and accuracy should be required for further advancement of simulation tools which can afford to cope with

69 (a)

(b)

°

-

(c)

°

°

o

.

.

m.

(d)

Fig. 10. Finite element control for composite elements: (a) geometric model; (b) generation of nodes; (c) mesh generation; (d) mesh modeling for composite elements.

prediction of work configuration changes in the actual two- and three-dimensional metal forming and mechanical assessment of physical quantities for materials evaluation. For such purposes, a new family of elements or composite elements is indispensable to create more reliable and accurate simulations. Here to be discussed is the logical compatibility of these composite elements with the finite element control which is necessary for practical simulation. With respect to automatic mesh generation, only a small effort must be added to automatic triangulation to provide a composite mesh alignment for practical simulation. The authors have been concerned with the completely automatic element control in two- and three-dimensional situations on the basis of the computational geometry [17-21]. Through these studies, mesh refinement and derefinement functions can be easily implemented for the composite finite element analysis. Figure 10 illustrates an example of mesh generation of compos-

70

•

Is I

~

i

i

kA---

~L

Fig. 11. Finite element mesh generation of tetrahedral elements in the three-dimensional situation.

ite elements for rotor-materials. Only the initial geometric configuration and the reference element length are needed as input data. In the subsequent stages, both generation of nodes and construction of elements are controlled and managed in a completely automatic manner. Even in three-dimensional situations, the algorithms created for the two-dimensional case can be extended directly; the necessary modifications can be constructed in our developed frame of computational geometry. Figure 11 shows the typical element generation of tetrahedral elements. Rapidly advancing forging technologies strongly require for three-dimensional simulation and the related forgeability evaluation. A reliable element control with adaptive mesh refinement and de-refinement and ALE (Arbitrary

71

Eulerian Lagrangian) is indispensable for the above purpose together with the use of robust and stable finite element schemes like the composite elements [22]. Further discussions of mesh control and ALE will be reported in near future.

References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

18 19 20 21 22

S. Toyoshima et al., Adv. Technol. Plasticity, (1990) 987-992. T. Altan and S.I. Oh, Adv. Technol. Plasticity, (1990) 1779-1787. D.S. Malkus and T.J.R. Hughes, Comput. Methods Appl. Mech. Eng., 15 (1978) 68-81. M. Crouzeix and P.A. Raviart, Rev. Fr. Autom. Inf. Rech. Oper., R-3 {1973) 33-76. J. Boland and R. Nicolaides, SIAM J. Numer. Anal., 20 (1983) 722-731. J. Pitkaranta and R. Sternberg, Report MAT-A222, HUT, 1984. H. Ohtubo et al., RC-89 Report, JSME, 1990. N. Kikuchi, Private Communication, 1990. T. Nakashima, PhD Thesis, University of Tokyo, 1991. T. Sussman and K.J. Bathe, Comput. Struct., 26 (1987) 357-409. O.C. Zienkiewicz, Y.C. Liu and G.C. Huang, Int. J. Numer. Methods Eng., 28 (1989) 21912202. T. Yamada, PhD Thesis, University of Tokyo, 1991. N. Yukawa, N. Kikuchi and A. Tezuka, Adv. Technol. Plasticity, (1990) 1719-1728. T. Aizawa and J. Kihara, Proc. Eur. Conf. on New Advances in Computational Structural Mechanics, 1991, p p . 112-117. T. Aizawa and J. Kihara, IUTAM Conf., Hannover, Germany, August 1991, to be presented. T.J.R. Hughes, The Finite Element Method, North-Holland, Amsterdam, 1988. T. Aizawa, Proc. Int. Conf. Computational Engineering Science, ICES 88, Post-Conference Seminar on Computer-Aided Engineering and Educational Aspects of Computational Mechanics, 1988, pp. 29-48. T. Aizawa, Proc. Annual Meeting of the Computer Graphics and Geometrical Sciences, 1990, pp. V-l-V-8. T. Aizawa, Proc. 5th Int. Conf. on CAD/CAM/Robotics and Future Factories, Springer, Berling, 1991, to be published. T. Aizawa, Computer Aided Innovation of New Materials, Elsevier, to be published. T. Aizawa, Proc. Int. Conf. Computational Engineering Science, ICES 91, August 1991, to be published. T. Aizawa, Proc. Symp. on Computational Methods in Structural Engineering and Related Fields, 1991, to be submitted.