viscoplastic element free Galerkin method

Materials Science and Engineering A 479 (2008) 197–212

Bulk metal forming process simulation based on rigid-plastic/viscoplastic element free Galerkin method Ping Lu, Guoqun Zhao ∗ , Yanjin Guan, Xin Wu Mold & Die Engineering Technology Research Center, Shandong University, Jinan 250061, China Received 1 February 2007; received in revised form 14 June 2007; accepted 19 June 2007

Abstract A rigid-plastic/viscoplastic element free Galerkin method is established based on the flow formulation for rigid-plastic/viscoplastic materials to simulate bulk metal forming processes. According to incomplete generalized variational principle, stiffness equations and solution formulas are derived. The transformation method is adopted to impose the essential boundary condition in the local coordinate system. The arctangent frictional model is used to implement the frictional boundary conditions, and the transform matrix from the global coordinate system to local coordinate system is given. The analysis software is developed. The method for dealing with the rigid region is introduced. Pressure projection method is used to solve volumetric locking and pressure oscillation problems. The contact and detachment criterion for the workpiece and die is given to judge the contact state. Numerical examples of bulk metal forming process are analyzed based on rigid-plastic/viscoplastic element free Galerkin method. The material flow, field variable distributions, etc. are analyzed. The numerical analysis results obtained by the method proposed in the paper are in good agreement with those obtained by the finite element method and experiment. The correctness of the presented method is demonstrated. © 2007 Elsevier B.V. All rights reserved. Keywords: Element free Galerkin method; Meshless; Rigid-plastic/viscoplastic; Bulk metal forming; Numerical simulation analysis

1. Introduction Nowadays, numerical simulation plays a more and more important role in metal forming process with the development of numerical analysis and computer technology. Among these numerical simulation methods for metal forming, finite element method (FEM) is used widely. However, because of high mesh-dependence, FEM encounters some difficulties that the precision and efficiency of the simulation degrade when the meshes become severely distorted, especially in handling large deformations. In this case, remeshing that is regarded as an effective way is necessary to complete the whole process simulation. Although remeshing techniques have made some development, remeshing usually leads to not only the loss of the accuracy but also the increase of time-consuming. In addition, the remeshing techniques in the three-dimensional metal forming simulation still need further investigation. In order to avoid the above problems, a new numerical computational method—mesh free (meshless) [1] method was

∗

Corresponding author. Tel.: +86 531 88393238; fax: +86 531 88395811. E-mail address: [email protected] (G. Zhao).

0921-5093/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2007.06.059

proposed. The main characteristic of the method is that the domain of the problem is discretized by a set of nodes (particles). The approximation field function is obtained through information of nodes not meshes. In recent years, the meshless method develops promptly, and many international computational mechanics researchers have devoted themselves to the development of meshless method, and made lots of innovative contributions not only to theoretical researches [2–14] but also engineering applications [15–26]. Many progresses of meshless applications are made in the field of metal forming processes [27–48]. Some researchers have applied meshless method to solve metal forming problems and made certain achievements. Chen et al. [27–29] applied reproducing kernel particle meshless method (RKPM) in metal forming problems such as twodimensional ring compression test and upsetting, in which material was considered as perfectly plastic material. Comparison with experiments demonstrated the performance and effectiveness of the method in metal forming. Bonet and Kulaseqaram [30] addressed corrected smooth particle hydrodynamics (CSPH) method and applied it to metal-forming simulations. The effectiveness of the method in simulating metal forming was demonstrated by the numerical examples: plane strain upsetting,

198

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

plane strain/axisymmetric forging. Xiong et al. [31–35] used reproducing kernel particle meshless method and element free Galerkin (EFG) method for analyzing slightly compressible rigid-plastic materials plane strain rolling. The penalty technique was utilized for enforcing the essential boundary conditions. The main variables to the rolling process were calculated, and the results were in agreement with experimental data. Park [36,37] used the lagrangian reproducing kernel particle method (RKPM) to simulate material processing. A comparative study between elasto-plastic and rigid-plastic RKPM methods was made, and consistent results from elasto-plastic and rigid-plastic simulations for the metal forming application were obtained. Guo et al. [38,39] applied the rigid-plastic point collocation method to the analysis of plane strain forging and backward extrusion processes, and verified the obtained results by comparing with a rigid-plastic finite element solution. Alfaro et al. [40] presented an application of the meshless natural element method in simulating three-dimensional aluminum extrusion, and used ␣shape-based approach to extract the geometry of the domain at each time step, and illustrated the potential of the method using some examples. Kwon et al. [41,42] introduced the least-squares meshfree method for the analysis of elasto-plastic and rigidplastic metal forming. The residuals were represented in a form of first-order differential system using displacement/velocity and stress components as nodal unknowns. The main benefit of the method was that it did not employ structure of extrinsic cells for any purpose. Li et al. [43] utilized a lower order integration scheme and a particle to segment contact algorithm in 3D bulk forming numerical simulation analysis by using reproducing kernel particle method. Liu et al. [44] used the rigid-plastic reproducing kernel particle method for the analysis of three-dimensional bulk metal forming. The results were in good agreement with conventional finite element predictions. Wang et al. [45] presented a parallel meshless method based on the theory of the reproducing kernel particle method for 3D bulk forming process. The parallel contact search algorithm was also presented, and the simulation results demonstrated the efficiency of the parallel reproducing kernel particle method. Li and co-workers [46] employed reproducing kernel particle method for simulating upsetting and rolling process based on the formulation of slightly compressible rigid plastic materials. The meshless method is a node-based numerical method and does not need mesh. From the current researches mentioned above, the meshless method shows obvious merits, such as simple preprocessing and post-processing as well as high computational precision in the numerical simulation of bulk metal forming processes especially large deformation bulk metal forming processes in which severe distortion of meshes often occurs during conventional finite element method analysis. Researchers achieved remarkable progresses on the application of meshless method to metal forming simulations. In the study done by researchers, the main material models are perfectly plastic material, slightly compressible rigid-plastic material, elasto-plastic and rigid-plastic material. And the main meshless method is RKPM. For bulk metal forming processes, such as forging, rolling, extrusion, drawing, etc., the portion of the workpiece undergoing

plastic deformation is much larger than the portion undergoing elastic deformation. Thus, the elastic deformation can be generally neglected. The rigid-plastic/viscoplastic material model that neglects elastic deformation and simplifies the solution of equations is widely used in the simulations of bulk metal forming processes [49]. Element free Galerkin (EFG) has the character of high stability and high accuracy. The paper combines the flow theory and element free Galerkin (EFG) method, and establishes rigid-plastic/viscoplastic element free Galerkin (EFG) method. The related mathematic formulas, as well as key algorithms and technologies are given according to flow formulation for rigid-plastic/viscoplastic materials. Equal channel angular pressing process and cubic billet upsetting process are analyzed by rigid-plastic/viscoplastic EFG method numerically. The related mechanics parameters are calculated and detailed comparisons with the results obtained by finite element method and experiment are made to validate the effectiveness of the method presented in the paper. 2. Rigid-plastic/viscoplastic element free Galerkin method 2.1. Approximation of velocity field Using moving least-squares (MLS) scheme, the velocity approximation uh (x) at any point x in the sub-domain x can be defined by uh (x) =

N 

I (x)dI = T (x)d

(1)

I=1

where dI is “generalized” velocity vector, I (x) is the shape function of node xI and N is the total number of nodes, I is the Ith node. (x) = [1 (x) 2 (x) . . . N (x)]T

(2)

d = [d1 d2 . . . dN ]T

(3)

with

⎡

ΦI (x) ⎢ I (x) = ⎣ 0 0 dI = [ dxI

dyI

0

0

ΦI (x)

0

0

ΦI (x)

T

dzI ]

⎤ ⎥ ⎦

(4)

(5)

Due to the lack of Kronecker delta properties in the element free Galerkin method shape functions, namely ΦI (xJ ) = δIJ and dI = uh (xI ), the essential boundary conditions cannot be imposed directly. In order to exert the boundary conditions directly, the transformation method presented by Chen et al. [47] is adopted to modify shape function. Then the relationship between the nodal value uI and the “generalized” velocity dI is established and expressed as u = T d

(6)

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

where

where

u = [uhx (x1 )

uhy (x1 )

uhz (x1 )

...

uhx (xN )

uhy (xN )

uhz (xN )] ⎡

ΠP =



α 2

Φ1 (x1 )I

⎢ Φ (x )I ⎢ 2 1 =⎢ .. ⎢ ⎣ .

Φ1 (x2 )I Φ2 (x2 )I .. .

... ... .. .

Φ1 (xN )I Φ2 (xN )I ⎥ ⎥ ⎥ .. ⎥ ⎦ .

Πf = −

Φn (x1 )I

Φn (x2 )I

...

Φn (xN )I

SF

(8)

ΠD =

Ω



(7) ⎤

199

(˙εV )2 dΩ

(14)

F¯ i ui dS

(15)

σ¯ ε¯˙ dΩ

for rigid-plastic materials

Ω

ΠD =

Ω

U(˙εij ) =

U(˙εij )dΩ

ε¯˙

(16a)

σ¯ dε¯˙

0

where I is a unit vector. From Eq. (6), the following equation can be derived.

for rigid-viscoplastic materials

d = −T u

Introducing a similar method as in finite element method, the stiffness equations and other equations for meshless solution of metal forming process can be established. Considering the relationship of velocity and strain-rate and using Eq. (10), strain-rate ε˙ can be expressed as

(9)

Substituting Eq. (9) into Eq. (1), we obtain uh (x) =

N 

I (x)uI = T (x)u

(10)

I=1

where (11)

2.2.1. Variational formulation for rigid-plastic/viscoplastic materials When the plastic deformation for the rigid-plastic/viscoplastic material occurs, basic equations [49]: equilibrium equations, compatibility conditions, constitutive equations, yield criterion, incompressibility condition and boundary conditions should be satisfied. By employing penalty function method to constrain incompressibility condition, incomplete generalized variational principle considers that among the admissible velocities which satisfy velocity boundary conditions and compatibility conditions, the actual solution make the total energy rate functional have an minimum value. It is expressed as α 2 ˙ Π= (˙εV ) dΩ − F¯ i ui dS σ¯ ε¯ dΩ + 2 Ω Ω SF for rigid-plastic materials

(12a)

(˙εV ) dΩ − 2

Ω

for rigid-viscoplastic materials

SF

F¯ i ui dS (12b)

where σ¯ denotes the effective stress. ε¯˙ is the effective strain-rate. α is a very large positive constant called penalty constant, a recommended value for α is 105 –106 . ε˙ V represents the volumetric strain-rate, F¯ i is surface tractions. Eqs. (12a) and (12b) can be written as, Π = ΠD + ΠP + Πf

B(x) = [ B1 (x) With

2.2. Meshless formulation for rigid-plastic/viscoplastic models

α Π= U(˙εij )dΩ + 2 Ω

. . . BN (x) ]u = B(x)u

B2 (x)

(17)

where B(x) is strain-rate matrix and expressed as

(x) = (x)−T



ε˙ = [ B1 (x)

(16b)

(13)

⎡

B2 (x)

ΨI,x (x)

⎢ 0 ⎢ ⎢ ⎢ 0 BI (x) = ⎢ ⎢ Ψ (x) ⎢ I,y ⎢ ⎣ 0

. . . BN (x) ] 0

0

ΨI,y (x) 0 ΨI,x (x) ΨI,z (x)

ΨI,z (x)

0

(18)

⎤

⎥ ⎥ ⎥ ΨI,z (x) ⎥ ⎥ 0 ⎥ ⎥ ⎥ ΨI,y (x) ⎦ 0

(19)

ΨI,x (x)

Ψ I,x (x), Ψ I,y (x) and Ψ I,z (x) are the derivatives of transformed shape function of node I with respect to x, y and z, respectively. The matrix of the effective strain-rate ε¯˙ is defined as 1/2 ε¯˙ = (˙εT D˙ε)

(20)

Introducing Eqs. (17)–(20), the following equation is obtained 1/2 1/2 ε¯˙ = (uT BT DBu) = (uT Pu)

(21)

where D is a matrix relating stresses with strain-rates, namely ⎡ ⎤ 2/3 0 0 0 0 0 ⎢ 0 2/3 0 0 0 0 ⎥ ⎥ ⎢ ⎢ ⎥ ⎢ 0 0 2/3 0 0 0 ⎥ ⎥ (22) D=⎢ ⎢ 0 0 0 1/3 0 0 ⎥ ⎢ ⎥ ⎥ ⎢ ⎣ 0 0 0 0 1/3 0 ⎦ 0

0

0

0 ⎡

P11 ⎢P ⎢ 21 P(x) = BT (x)DB(x) ≡ ⎢ ⎢ .. ⎣ . PN1

0

1/3

P12

...

P1N

P22 .. .

... .. .

P2N .. .

PN2

. . . PNN

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

(23)

200

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

⎡

PIJ

PIJxx ⎢ T = BI (x)DBJ (x) ≡ ⎣ PIJyx PIJzx

PIJxy PIJyy

⎤ PIJxz ⎥ PIJyz ⎦

PIJzy

PIJzz

where σ¯ ∂ σ¯ (36a) = − 2 for rigid-plastic materials ∂ε¯˙ ε¯˙ ˙ε¯ ∂ σ¯ σ¯ 1 ∂σ¯ − 2 for rigid-viscoplastic materials (36b) = ∂ε¯˙ ε¯˙ ε¯˙ ∂ε¯˙ ˙ε¯

(24)

Volumetric strain-rate ε˙ V can be expressed as follows ε˙ V = ε˙ x + ε˙ y + ε˙ z

(25)

According to Eq. (35), the second-order variation of D is expressed as follows ⎧ ⎪ ⎪ ⎨ σ¯ ∂ σ¯ 2 T δ ΠD = δuI PIJ dΩ + ⎪ ¯˙ ε¯˙ ε¯˙ Ω ∂ε ⎪ ⎩ Ω

Substituting the corresponding strain-rate components in Eq. (17) into Eq. (25) results in ε˙ v = [

Ψ1,x (x)

Ψ1,y (x) Ψ1,z (x) ... ] u ΨN,x (x) ΨN,y (x) ΨN,z (x)

(26)

ε˙ V = [ CT1

. . . CTN ]u = CT u

CT2

1 × [ PI1 ε¯˙

(27)

where C = [ CT1

T

CT2

. . . CTN ]

T

2.2.3. The variation of P The substitution of Eq. (27) in Eq. (14) yields α 2 ΠP = (CT u) dΩ 2 Ω

(29)

2.2.2. The variation of D The first-order partial derivative of D is defined as follows 

∂ΠD ∂ΠD ∂ΠD T ∂ΠD = (30) ∂uxI ∂uyI ∂uzI ∂uI

(31b)

Ω

PI2

Ω

. . . PIN ]u

(32)

¯ IJ = CI CT C J ⎡ ΨI,x (x)ΨJ,x (x) ⎢ = ⎣ ΨI,y (x)ΨJ,x (x) ΨI,z (x)ΨJ,x (x)

(33)

ΨI,x (x)ΨJ,y (x) ΨI,y (x)ΨJ,y (x)

⎤ ΨI,x (x)ΨJ,z (x) ⎥ ΨI,y (x)ΨJ,z (x) ⎦

ΨI,z (x)ΨJ,y (x)

ΨI,z (x)ΨJ,z (x) (43)



δΠD = δuTI

and the second-order partial derivative and variation of P can be evaluated as ∂ 2 ΠP T ¯ IJ dΩ =α CI CJ dΩ = α (41) C ∂uI ∂uJ Ω Ω ¯ IJ dΩδuJ (42) δ2 ΠP = δuTI α CI CTJ dΩδuJ = δuTI α C where

Substituting Eq. (32) into Eq. (31) results in σ¯ ∂ΠD [ PI1 PI2 . . . PIN ]u dΩ = ¯˙ ∂uI Ω ε Therefore

Ω

∂ ∂ΠD = U(˙εij )dΩ ∂uI ∂uI Ω ∂ε¯˙ = dΩ for rigid-viscoplastic materials σ¯ Ω ∂uI

(38)

By employing a similar treating procedure as that of D , the first-order partial derivative and variation can be evaluated as ∂ΠP =α CI CT u dΩ (39) ∂uI Ω T δΠP = δuI α CI CT u dΩ (40)

By introducing Eqs. (16)–(24) and Eq. (30), the following equations can be obtained. ∂ε¯˙ ∂ΠD σ¯ = dΩ for rigid-plastic materials (31a) ∂uI Ω ∂uI

where 1 ∂ε¯˙ = [ PI1 ε¯˙ ∂uI

(37)

(28)

CI = [ ΨI,x (x) ΨI,y (x) ΨI,z (x) ]



⎫ ⎤ P1J ⎪ ⎪ ⎬ ⎢ ⎥ T ⎢ .. ⎥ . . . PIN ]uu ⎣ . ⎦ dΩ δuJ ⎪ ⎪ ⎭ PNJ ⎡

The above equation can be rewritten as

Ω

σ¯ [ PI1 ε¯˙

PI2

. . . PIN ]u dΩ

(34)

The second-order partial derivative of D can be evaluated as ∂ 2 ΠD = ∂uI ∂uJ

Ω

σ¯ PIJ dΩ + ε¯˙

Ω

∂ σ¯ 1 [ PI1 ∂ε¯˙ ε¯˙ ε¯˙

⎡

⎤ P1J ⎢ . ⎥ ⎥ . . . PIN ]uuT ⎢ ⎣ .. ⎦ dΩ PNJ

(35)

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

2.2.4. Frictional contact conditions In metal forming processes, the frictional conditions along the workpiece–die interface are complex. It is difficult to handle the frictional conditions by a standardized mathematical formula. Commonly, the frictional force can be calculated by some simple formulas. The usually used friction models include Coulomb law f = μp, and the constant friction factor law f = mk, where, μ denotes frictional coefficient, p denotes the pressure on the workpiece-die interface, m is called frictional factor, and k represents shear yield stress. The direction of the frictional stress is opposite to the direction of the relative velocity. There is a point (or a region) along the die–workpiece interface where the velocity of the deforming material relative to the die becomes zero, and the location of this point (or region) depends on the magnitudes of the frictional stress itself. The frictional stress at the neutral point (or region) changes suddenly. The sudden change of the frictional stress near the neutral point can cause difficulties in numerical calculation. The arctangent friction law [49] defined as follows can solve the difficulties effectively. So it is employed in the paper to treat the friction conditions on the workpiece–die interface.

 2 −1 ws f = −mk (44) tan π u0 where ws is the sliding velocity of material relative to the die velocity, u0 is an very small positive constant and recommended as 10−3 to 10−4 . In the three-dimensional problems, establish the tangentnormal local coordinate system s, t, n, which make up of the right hand coordinate system for contact node. And s denotes the first tangential direction, t represents the second tangential direction, and n denotes the normal direction and its positive direction is from workpiece to die, respectively. The first and second tangential direction components of relative slipping velocity on the contact interface for the contact node are defined as us (x) and ut (x). So ws (x) can be defined as  ws (x) = us (x)2 + ut (x)2 (45) And us (x) and ut (x) can be approximated as follows using moving least-squares (MLS) scheme, namely us (x) =

N 

ΨI (x)usI

(46)

ΨI (x)utI

(47)

I=1

ut (x) =

N  I=1

where usI and utI indicate the first and second tangential direction components of relative slipping velocity for node I which is on the frictional contact interface, and can be evaluated as usI = uxI − sT ud

(48)

utI = uyI − tT ud

(49)

where uxI and uyI are the first and second tangential direction components of slipping velocity for node I, ud represents

201

the velocity of die, s is the first tangential direction unit vecT tor on the contact interface expressed as s = ( sx sy sz ) , t is the second tangential direction unit vector expressed as T t = ( tx ty tz ) . The relationship for velocity field between global coordinate system and local coordinate system are given by u = Tu

(50)

where u = [u1 u2 . . . uN ]

T

(51)

with uI = [ uxI and

⎡

⎢ ⎢ T=⎢ ⎢ ⎣ with

T1

⎡

0 .. . 0

tI11

⎢ TI = ⎣ tI21 tI31

uyI 0 T2 .. . 0

uzI ]

T

... ... .. . 0

tI12

(52)

0

⎤

0 ⎥ ⎥ ⎥ .. ⎥ . ⎦ TN

tI13

(53)

⎤

tI22

⎥ tI23 ⎦

tI32

tI33

(54)

In the above matrix, tIjm is the cosine value of angle from direction j of local coordinate system to direction m of global coordinate system. Because T is an orthogonal matrix, it can be obtained u = T−1 u = TT u

(55)

Substituting Eq. (44) into Eq. (15), the following Eq. (56) is obtained

  ws 2 −1 ws Πf = dws dS (56) tan mk π u0 Sc 0 The first-order partial derivative of f in local coordinate system can be evaluated as  ∂Πf us 2 −1 ws = mkΨI tan dS (57) ∂uxI π u w 0 s Sc  ut ∂Πf 2 −1 ws = mkΨI tan dS (58)  ∂uyI π u0 ws Sc and the second-order partial derivative of f can be obtained as follows  u u 2 u0 ∂ 2 Πf = mk ΨI ΨJ s 2s + tan−1   2 2 ∂uxI ∂uxJ π u0 + w s ws Sc  ws ΨI ΨJ us us −1 ws ΨI ΨJ 3 dS × − tan u0 ws u0 ws (59)

202

∂ 2 Πf = ∂uxI ∂uyJ

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

mk Sc

2 π

− tan−1

∂ 2 Πf = ∂uyI ∂uxJ



2 mk π Sc

− tan−1

∂ 2 Πf = ∂uyI ∂uyJ



us ut u0 Ψ Ψ I J w2s + w2s  ws us ut ΨI ΨJ 3 dS u0 ws u20

(60)



u0 u  u ΨI ΨJ t 2s 2 2 ws u0 + w s  ws u u ΨI ΨJ t 3s dS u0 ws

2 mk π Sc

(61)



In the two-dimensional problems, s denotes the tangential direction, n denotes the normal direction and its positive direction is also from workpiece to die. Thus, the following Eqs. can be obtained for two-dimensional problems. ws (x) = us (x) ∂Πf 2 −1 ws dS = mkΨI tan ∂uxI π u0 Sc ∂ 2 Πf 2 u0 = mk Ψ Ψ dS 2 + w2 I J ∂uxI ∂uxJ π u Sc s 0   tI11 tI12 TI = tI21 tI22

namely

∂Πf ∂ 2 Πf

u(n) = f − S+ ∂u∂uT ∂u where



SIJ =

(63)

Ω

σ¯ PIJ dΩ + ε¯˙

1 × [ PI1 ε¯˙

(64)

Ω

∂ ∂ε¯˙

+α

(65)

Ω

(66)

fI = −

Ω



(67)

Introducing Eqs. (34) and (40) into above equation and considering arbitrariness of δu, results in σ¯ ∂Πf [ PI1 . . . PIN ]u dΩ + α CI CT u dΩ + =0 ˙ ¯ ∂uI ε Ω Ω (68) The above equation is called the meshless stiffness equation of rigid-plastic/viscoplastic metal forming problems. The NewRaphson iterative procedure is implemented to solve the above non-linear equation until converged results are obtained. The Eq. (68) can be rewritten as

 

2 ∂Π ∂ Π (n) +

uJ = 0 (69) ∂uI u=u(n−1) ∂uI ∂uJ u=u(n−1)

−α Ω

⎡

. . . PIN ]u(n−1) u(n−1)

T

⎤ P1J ⎢ . ⎥ ⎢ . ⎥ dΩ ⎣ . ⎦ PNJ

¯ IJ dΩ C

σ¯ [ PI1 ε¯˙

(70)

(71)

σ¯ ε¯˙



2.2.5. Stiffness equation Applying incomplete generalized variation principle, make the first-order variation of the total energy rate functional vanish, namely δΠ = δΠD + δΠP + δΠf = 0

⎫ ⎤ P1J ⎪ ⎪ ⎢ . ⎥ ∂ 2 Πf ⎬ (n) ⎥ ⎢ ¯

uJ × ⎣ .. ⎦ dΩ+α CIJ dΩ+ ∂uI ∂uJ ⎪ Ω ⎪ ⎭ PNJ σ¯ =− [ PI1 . . . PIN ]u(n−1) dΩ ˙ ¯ ε Ω ∂Πf −α CI CT u(n−1) dΩ− ∂uI Ω ⎡

ut ut u0 −1 ws Ψ Ψ + tan I J w2s u0 u20 + w2s  ΨI Ψ J ws u u × (62) ΨI ΨJ t 3t dS − tan−1 ws u0 ws



Introducing Eqs. (35) and (41) to Eq. (69), Eq. (69) can be rewritten and re-arranged as ⎧ ⎪ ⎪ ⎨ σ¯ ∂ σ¯ 1 T PIJ dΩ+ [ PI1 . . . PIN ]u(n−1) u(n−1) ˙ ˙ ˙ ˙ ⎪ ¯ ε¯ ε¯ ε¯ Ω ∂ε ⎪ ⎩ Ω

(72)

. . . PIN ]u(n−1) dΩnonumber

CI CT u(n−1) dΩ

(73) (73)

In order to employ frictional contact conditions directly, it is need to assemble the stiffness equation in local coordinate system. Hence, Eq. (71) is rewritten as

∂ 2 Πf ∂Πf (n) (74) TSTT +   T u = Tf − ∂u ∂u ∂u and further simplified as K(n) u

(n)

= F(n)

(75)

where K(n) = TSTT + F(n) = Tf −

∂ 2 Πf ∂u ∂u T

(76)

∂Πf ∂u

The expressions of (57)–(65).

(77) ∂Πf ∂u

and

∂2 Πf ∂u ∂u T

are obtained from Eqs.

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

The integrals of all above equations can be evaluated by means of Gaussian quadrature through the utilization of background cells. 2.3. Treatment of the rigid region In metal forming processes, the rigid region that is characterized by a very small value of effective strain-rate exists, in which the value of the first term in Eq. (12a) or (12b) cannot be calculated when the effective strain-rate approaches zero. This difficulty is removed by the following linear constitutive equation (78) when the effective strain-rate of the Ith node is smaller than ε¯˙ 0 ε˙ ij =

3 ε¯˙ 0  σ 2 σ¯˙ 0 ij

(78)

problems encountered in rigid-plastic/viscoplastic element free Galerkin method. The algorithm is realized by modifying the following penalty function items p of total energy rate functional. α α 2 (˙εV )2 dΩ = (CT u) dΩ (85) Πp = 2 Ω 2 Ω Volumetric strain-rate ε˙ V is mapped onto each integration domain Ωxs which is a lower order space spanned by a set of functions Q = {Q1 (x), Q2 (x),. . .,Qn (x)}. ε˙ sV = {˙εsV 1 , ε˙ sV 2 , . . . , ε˙ sVn }T is obtained by minimizing R(␧˙ sV ) as follows R(˙εsV ) = ||˙εV − Q˙εsV ||2L2 (Ωxs )

Ms ε˙ sV = Zs

Then,

where

δΠD =

δuTI

Ω

σ¯ 0 ε¯˙ 0

ε¯˙ δε¯˙ dΩ

Substituting Eq. (32) into Eq. (79) yields σ¯ 0 ∂ΠD [ PI1 PI2 . . . PIN ]u dΩ = ¯˙ 0 ∂uI Ω ε

(80)

∂R(˙εsV ) =0 ∂ε˙ sV

Then the second-order partial derivative of D and the second-order variation of D are replaced as follows σ¯ 0 ∂ 2 ΠD = PIJ dΩ (83) ¯˙ 0 ∂uI ∂uJ Ω ε σ¯ 0 2 T δ ΠD = δuI PIJ dΩδuJ (84) ¯˙ 0 Ω ε 2.4. Treatment of volumetric locking In order to obtain satisfactory accuracy, a high-order integration rule is often used in meshless method. This often leads to an over-constrained condition and causes volumetric locking and pressure oscillation problems. Although increasing the size of the influence domain can release volumetric locking, it cannot solve the problem of pressure oscillation and deviates from the local fitting character owned by meshless method. In addition, the computational efficiency is reduced as the increase of the size of influence domain. In the paper, a releasing algorithm is established based on pressure projection method that is proposed by Chen [48] to solve volumetric locking and pressure oscillation

(88)



Ms =

The first-order variation of D is expressed as follows σ¯ 0 T [ PI1 PI2 . . . PIN ]u dΩ (82) δΠD = δuI ¯˙ 0 Ω ε

(87)

The following equation can be obtained

Zs = (81)

(86)

where || · ||2L2 (Ωxs ) is least-square norm in the integration domain Ωxs . And

¯ ε, ε¯˙ 0 ), ε¯˙ 0 takes an assigned limiting value. where σ¯ 0 = σ(¯ Thus, in the rigid region the Eq. (31a) or (31b) is replaced by the following equation (79). σ¯ 0 ∂ε¯˙ ∂ΠD ε¯˙ = dΩ (79) ¯˙ 0 ∂uI ∂uI Ω ε

203

QT Q dΩ

(89)

QT ε˙ V dΩ

(90)

Ωxs

Ωxs

and the projected volumetric strain-rate ε˙ ∗s V in the integration domain Ωxs is calculated as follows ε˙ ∗s εsV = QMs−1 Zs V = Q˙

(91)

When the volumetric strain-rate is projected onto a constant field Q = {1, 1,. . .,1}, the following is obtained  T Ωs C u dΩ ∗s ε˙ V = x (92) As where As is the volume of Ωxs . 2.5. Dynamic adjustment of the boundary nodes The contact and detachment state between workpiece boundary and die surface changes continually with the increase of deformation during metal forming processes. Therefore, it is necessary to give the contact and detachment criterion for the workpiece and die. After the velocity field is convergent in one incremental step, it needs to judge whether the free boundary nodes contact die or not in this step, and whether the contact nodes separate from the die or not. The time increment tI for each boundary free node to contact die should be calculated. Where I represents the Ith node. The problem can be solved by treating with the intersection of a ray and the die surface. The starting point coordinate of the ray is the current location coordinate expressed as [xI ,yI ,zI ], and

204

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

the direction of the ray can be expressed by the velocity vector as [uxI − uxd , uyI − uyd , uzI − uzd ], in which [uxI ,uyI ,uzI ] is the velocity components of the node, and [uxd ,uyd ,uzd ] is the velocity components of the die. Assuming the coordinate of the intersecting point is [xc ,yx ,zc ], the value of tI can be calculated as follows 1/2

tI =

[(xc − xI )2 + (yc − yI )2 + (zc − zI )2 ]

1/2

[(uxI − uxd )2 + (uyI − uyd )2 + (uzI − uzd )2 ]

(93)

If tI is smaller than the time increment t, the Ith node is regarded as the contact one and imposed contact restriction. Otherwise, the node is still free. The normal stress σ n is used as the detachment and contact criterion for the nodes that have been contact nodes in the last step. If σ n ≥ 0, the nodes which have contacted with the die surface in the last step separate from the die and become free in the current step. And the restriction should be released. Or else, the nodes are still regarded as the contact ones, and move along the die surface. After the detachment and contact judgment is done, the geometry of the workpiece can be obtained by updating from the previous value by xI (t + t) = xI (t) + uxI t

(94)

yI (t + t) = yI (t) + uyI t

(95)

zI (t + t) = zI (t) + uzI t

(96)

The effective strain of the Ith node can be updated and expressed as follows

Fig. 2. The metal flow patterns by using EFG method at different strokes of punch. (a) Stroke is 0 mm, (b) stroke is 15 mm, (c) stroke is 45 mm and (d) stroke is 75 mm.

a cuboid solid with the length of 80 mm and the cross-section dimension of 10 mm × 10 mm. And the equal channel angular pressing belongs to a plane strain problem. The material is pure aluminum (Al 99.99%), and the material model is regarded as rigid-plasticity. The flow stress of the material is σ¯ = C¯εn

(98)

The schematic diagram of ECAP is illustrated in Fig. 1. Φ is die channel angle and usually ranges from 90◦ to 150◦ , Ψ is die corner angle and ranges from 0◦ to (180◦ - Φ). The billet is

where C = 170 MPa and n = 0.24. In the simulation, the workpiece is discretized by 11 × 81 nodes, friction factor m = 0.12, Φ is 90◦ , Ψ is 37◦ . The speed of the punch is 2 mm/s, the time increment is 0.25 s and the temperature is 20 ◦ C. During the analysis, the Gauss quadrature order is 3 × 3, linear basis function is employed. Fig. 2 gives the metal flow patterns when the strokes of the punch are 0, 15, 45, and 75 mm, respectively. Because of no remeshing, the discrete nodes with the deformation increasing can clearly form flow lines that can reflect the material locus. So the flowing states of material are more visual by using the EFG method.

Fig. 1. The schematic diagram of ECAP process.

Fig. 3. The distributions of effective strain in the workpiece by using EFG method. (a) Stroke is 15 mm, (b) stroke is 45 mm and (c) stroke is 75 mm.

ε¯ I (t + t) = ε¯ I (t) + ε¯˙ I t

(97)

where ε¯ I is the effective strain of the Ith node, and ε¯˙ I is the effective strain-rate. 3. Numerical examples 3.1. Equal channel angular pressing (ECAP)

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

Fig. 4. The distributions of effective stress in the workpiece by using EFG method. (a) Stroke is 15 mm, (b) stroke is 45 mm and (c) stroke is 75 mm.

Fig. 3 shows effective strain contours of the workpiece when the strokes of the punch are 15, 45 and 75 mm. The effective strain is accumulated relying on the shear force when extrusion part goes through the die corner with the increasing of punch stroke to refine grains. It can be observed that the effective strain distributes more uniformly in the horizontal direction except stub bar and two channels intersection regions of the extrusion part. This is because shear deformation only occurs in the two channels intersection zones and it hardly occurs in other regions. In the vertical direction the effective strain of the regions that are close to the upper surface and the center portion is large, the

205

Fig. 5. The distributions of effective strain-rate in the workpiece by using EFG method. (a) Stroke is 15 mm, (b) stroke is 45 mm and (c) stroke is 75 mm.

effective strain decreases with the distance getting smaller to the lower surface. Fig. 4 gives the distribution contours of effective stress for the extrusion part at different strokes 15, 45 and 75 mm, respectively. Obviously, the effective stress in the die corner is more concentrated, and the value and the gradient of the effective stress are large. The effective stress in the inner corner is larger than that in the outer corner. These indicate the workpiece will suffer the large force when it goes through the die corner. Therefore, the appropriate materials should be selected during the die design, and heat treatment process should be carried out to ensure the carrying capacity and strength requirements of the dies.

Fig. 6. Different field variable distributions in the workpiece by using Deform (2D) when the stroke of the punch is 45 mm. (a) Distribution of effective strain, (b) distribution of effective stress and (c) distribution of effective strain-rate.

206

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

Fig. 5 shows effective strain-rate contours of the workpiece when the strokes of the punch are 15, 45 and 75 mm, respectively. The effective strain-rate reduces gradually from the inner corner to outer corner. This explains why the effective strain decreases from the upper surface to the lower surface displayed in Fig. 3. The numerical model is also simulated by the finite element method software Deform (2D) to verify the correctness of the established method. The workpiece is discretized by 800 quadrangle elements with 891 nodes. Other parameters are similar with those during EFG method simulation. During the FEM simulation, mesh distortion occurs when the stroke is 8.5mm, and three times remeshing are needed in order to enable the simulation analysis continue. The effective strain, effective stress and effective strain-rate obtained by Deform (2D) are shown in Fig. 6 when the stroke of the punch is 45 mm. Comparing Fig. 3(b), Fig. 4(b) and Fig. 5(b) with Fig. 6(a)–(c), it finds that the distributions of the mechanics variables are in good agreement. The equal channel angular pressing experiment is carried out. The billet geometry with the dimension 10 mm × 10 mm × 80 mm is similar with the simulation model, and the other process parameters are almost the same as in the numerical simulation. Fig. 7 gives the experimental dies. The forming equipment is 300 kN electronic universal testing machine in the experiment (shown in Fig. 8). Fig. 9(a) and (b) give the geometric shapes of the workpiece obtained by EFG analysis and experiment. The symbol A-B signed in Fig. 9(a) and (b) indicates the shearing direction, and the shearing direction gained from numerical simulation is identical with that gained from actual experiment. The label C indicates the zone in the circle is stub bar portion of the extrusion part, and the geometric shapes obtained from numerical simulation and actual experiment, respectively, are similar. Label D in Fig. 9(a) denotes the metal flowing direction calculated by EFG method, and label D in Fig. 9(b) denotes texture streamline. It is clear that both of the labels have the some trend. It can be seen that the simulation results agree with the experiment results.

Fig. 7. The experimental die of square-section ECAP.

Fig. 8. The equipment for experiments—300 kN testing machine.

3.2. Cubic billet upsetting Upsetting process is a typical metal forming process. A threedimensional cubic billet upsetting process was analyzed by using the EFG method proposed in the paper. The predicted results were compared with those given by commercial FEM softwareDeform (3D). The initial geometrical shape of the billet is shown in Fig. 10. The size of the billet is 20 mm × 20 mm × 20 mm and 11 × 11 × 11 divided nodes are used for the simulation. The upper die moves downward with a velocity of −1 mm/s, while the lower die is stationary. The arctangent friction model was used for modeling the friction behavior on the workpiece–die interface. The constant frictional factor m is selected as 0.3. Due to the symmetry of geometry and deformation, only one eighth of the billet discretized with 6 × 6 × 6 nodes is considered. In the analysis using EFG method, the Gauss quadrature order is 4 × 4 × 4, linear basis function is used, the weight function is cubic spline weight function, and the sphere domains of influence are employed. In order to verify the validity of the predicted solution of the EFG, the problem is simulated by FEM software-Deform (3D) too. The FEM model is discretized by

Fig. 9. The numerical and experimental analysis of the workpiece. (a) The geometric shapes of the workpiece obtained by EFG analysis. (b) The geometric shapes of the workpiece obtained by experiment.

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

207

Table 1 Displacement of maximum bulge at 20% reduction in height

EFG FEM-Deform (3D)

Displacement (mm)

Absolute error (mm)

Relative error (%)

1.50807 1.47342

0.03465

2.35167%

Notes: Absolute error = EFGresult − FEMresult .Relative error = 100%.

Absolute error FEMresult

×

Table 2 Displacement of maximum bulge at 40% reduction in height

EFG FEM-Deform (3D) Fig. 10. The billet shape and the initially discretized meshless model.

999 tetrahedral elements with 233 nodes. The processes from 0 to 40% reduction in height are analyzed, and the incremental step in reduction is 2% of the initial billet height. The billet material is AISI 8620 steel and the temperature is 1100 ◦ C. The relationship between strain and stress is given as m

σ¯ = Cε¯˙

(99)

where C = 75.3 MPa and m = 0.134. Fig. 11 describes the geometrical shapes of the deforming material at different stages of deformation. Fig. 11(a)–(d) give

Displacement (mm)

Absolute error (mm)

Relative error (%)

3.48673 3.47687

0.00986

0.28359

Notes: Absolute error = EFGresult − FEMresult .Relative error = 100%.

Absolute error FEMresult

×

the geometrical shapes of the deforming material obtained by EFG and FEM-Deform (3D) at 20 and 40% reduction in height, respectively. Clearly, the bulge at the stage of 20% reduction in height is not very obvious. Meanwhile, the shape of the formed part obtained from EFG is in good agreement with that obtained from FEM-Deform (3D). When the deformation becomes bigger, the bulge is obvious, and the agreement of the shapes between EFG and FEM is still good. The displacement values of maximum bulge analyzed by EFG and FEM-Deform

Fig. 11. Material flow patterns at different reductions in height. (a) The shape of the deforming body obtained by EFG at 20% reduction in height. (b) The shape of the deforming body obtained by FEM-Deform (3D) at 20% reduction in height. (c) The shape of the deforming body obtained by EFG at 40% reduction in height. (d) The shape of the deforming body obtained by FEM-Deform (3D) at 40% reduction in height.

208

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

Fig. 12. Velocity vector distributions of all nodes at 20 and 40% reduction in height. (a) Velocity vector distributions obtained by EFG at 20% reduction in height. (b) Velocity vector distributions obtained by FEM-Deform (3D) at 20% reduction in height. (c) Velocity vector distributions obtained by EFG at 40% reduction in height. (d) Velocity vector distributions obtained by FEM-Deform (3D) at 40% reduction in height.

(3D), and corresponding absolute errors and relative errors at 20% and 40% reduction in height are shown in Tables 1 and 2, respectively. It can be seen that the displacement values of maximum bulge analyzed by EFG is larger than those analyzed by FEM both at 20 and 40% reduction in height. But the corresponding relative errors of the results calculated by two methods are small, the relative error is 2.35167% at 20% reduction in height and it is only 0.28359% at 40% reduction in height.

Fig. 12 shows the velocity vector plots of all nodes at 20 and 40% reduction in height. It can be seen from these plots that, the velocity fields of nodes evaluated by EFG match very well with those evaluated by FEM. The directions and values of velocity vector and the material flow trends for both EFG and FEM are very close, too. The positions, velocity directions and magnitudes of the characteristic nodes numbered as 1–8 in Fig. 12 are listed in Tables 3–6. The numerical comparisons show that the results of EFG and FEM are in good agreement.

Table 3 The coordinates, velocity directions and values of the characteristic nodes (CN) at 20% reduction in height evaluated by EFG Sequential number of CN

1 2 3 4 5 6 7 8

Coordinate

Velocity value (mm/s)

x

y

z

10.67954 0 0 10.87233 11.50735 0 0 11.23201

0 0 10.68996 10.87127 0 0 11.50807 11.23119

8 8 8 8 0 0 0 0

1.10233 1 1.10577 1.24936 0.83697 0 0.83788 0.99152

Velocity direction α (◦ )

β (◦ )

γ (◦ )

65.1172 90 90 64.8948 0 – 90 44.9681

90 90 64.7355 64.946 90 – 0 45.0319

155.117 180 154.736 143.169 90 – 90 90

Notes: α represents the angle between velocity direction of the node and x direction in global coordinate system; β represents the angle between velocity direction of the node and y direction in global coordinate system; γ represents the angle between velocity direction of the node and z direction in global coordinate system.

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

209

Table 4 The coordinates, velocity directions and values of the characteristic nodes (CN) at 20% reduction in height evaluated by FEM Sequential number of CN

1 2 3 4 5 6 7 8

Coordinate

Velocity value (mm/s)

x

y

z

10.77209 0 0 10.97056 11.47342 0 0 11.24257

0 0 10.75682 10.98967 0 0 11.46725 11.25626

8 8 8 8 0 0 0 0

1.10698 1 1.10275 1.31254 0.82969 0 0.82473 0.97959

Velocity direction α (◦ )

β (◦ )

γ (◦ )

64.6035 90 90 63.1518 0 – 90 45.3089

90 90 65.0696 62.3356 90 – 0 44.6911

154.604 180 155.07 139.63 90 – 90 90

Notes: α represents the angle between velocity direction of the node and x direction in global coordinate system; β represents the angle between velocity direction of the node and y direction in global coordinate system; γ represents the angle between velocity direction of the node and z direction in global coordinate system.

The comparisons of effective strain distributions obtained by using the two methods at 20 and 40% reduction in height are shown in Fig. 13, respectively. For two methods, the value of effective strain in the central upper zone in contact with the die is very small and usually called the dead zone, the material in this zone almost has no deformation. The material

in the middle zone and the corner area has a large effective strain distribution and this zone or area is usually called large deformation zone. The distributions of effective strain calculated by EFG is accordant with that calculated by FEM. The deformation distribution is in good agreement with the practical situation.

Fig. 13. Comparisons of effective strain distributions obtained by EFG and FEM at 20 and 40% reduction in height. (a) Color contours of effective strain obtained by EFG at 20% reduction in height. (b) Color contours of effective strain obtained by FEM-Deform (3D) at 20% reduction in height. (c) Color contours of effective strain obtained by EFG at 40% reduction in height. (d) Color contours of effective strain obtained by FEM-Deform (3D) at 40% reduction in height.

210

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

Table 5 The coordinates, velocity directions and values of the characteristic nodes (CN) at 40% reduction in height evaluated by EFG Sequential number of CN

1 2 3 4 5 6 7 8

Coordinate

Velocity value (mm/s)

x

y

z

11.87739 0 0 12.27201 13.48566 0 0 12.83488

0 0 11.9226 12.26884 0 0 13.48673 12.83913

6 6 6 6 0 0 0 0

Velocity direction α (◦ )

1.21122 1 1.22152 1.52761 1.1349 0 1.13331 1.24863

55.6505 90 90 57.3439 0 – 90 45.843

β (◦ )

γ (◦ )

90 90 54.9503 57.1316 90 – 0 44.157

145.651 180 144.95 130.065 90 – 90 90

Notes: α represents the angle between velocity direction of the node and x direction in global coordinate system; β represents the angle between velocity direction of the node and y direction in global coordinate system; γ represents the angle between velocity direction of the node and z direction in global coordinate system.

Fig. 14 gives the isoline contours of effective strain on the cross-section of x = 3.5 and z = 2.0 of the deforming body obtained by EFG and FEM at the reduction in height of 40%, respectively. Through comparing the isolines, it can be also seen that the change rules of effective strain obtained by the two methods are the same. The effective strains in the central zone of the

deforming material are large, the effective strains in the central upper zone in contact with the die are smallest, and the effective strains in the middle zone close to the surface of the deforming body are small. Fig. 15 shows the load-stroke curves of upsetting process for both EFG and FEM. The curves are in good agreement.

Fig. 14. Isoline contours of effective strain on the different cross-sections given by EFG and FEM. (a) Isoline contours of effective strain on the cross-section of x = 3.5 given by EFG. (b) Isoline contours of effective strain on the cross-section of x = 3.5 given by FEM-Deform (3D). (c) Isoline contours of effective strain on the cross-section of z = 2.0 given by EFG. (d) Isoline contours of effective strain on the cross-section of z = 2.0 given by FEM-Deform (3D).

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

211

Table 6 The coordinates, velocity directions and values of the characteristic nodes (CN) at 40% reduction in height evaluated by FEM Sequential number of CN

1 2 3 4 5 6 7 8

Coordinate

Velocity value (mm/s)

x

y

12.05036 0 0 12.49641 13.47687 0 0 12.91056

0 0 12.01026 12.5749 0 0 13.4566 12.9308

z 6 6 6 6 0 0 0 0

1.28427 1 1.27541 1.67225 1.18324 0 1.17561 1.39875

Velocity direction α (◦ )

β (◦ )

γ (◦ )

51.1374 90 90 56.3474 0 – 90 44.9163

90 90 51.6342 54.6154 90 – 0 45.0837

141.137 180 141.634 126.726 90 – 90 90

Notes: α represents the angle between velocity direction of the node and x direction in global coordinate system; β represents the angle between velocity direction of the node and y direction in global coordinate system; γ represents the angle between velocity direction of the node and z direction in global coordinate system.

is also needed to provide a tool for guiding the process and die design in industrial production. Acknowledgement The research work was supported by National Natural Science Foundation for Distinguished Young scholars of China (Approval No.: 50425517) and National Natural Science Foundation of China (Approval No.: 50575125). References

Fig. 15. Load–stroke curves obtained by EFG and FEM.

4. Conclusions Based on the flow formulation for rigid-plastic/viscoplastic materials, element free Galerkin (EFG) method is applied to the simulation of rigid-plastic/viscoplastic bulk metal forming processes. The EFG mathematic model, the related equations, key algorithms and technologies are proposed in this paper. In order to validate the feasibility and effectiveness of EFG in dealing with rigid-plastic/viscoplastic bulk metal forming problems, equal channel angular pressing and cubic billet upsetting as the examples are simulated numerically. The results such as material flow, forming load and different field variable distributions are calculated, respectively, and the results are in good agreement with those obtained by commercial finite element software and experiment. The meshless method exhibits some advantages for solving bulk metal forming problems since it requires no mesh generation and remeshing in the simulation. For the future lots jobs should be done to promote the application for element free Galerkin (EFG) method in bulk metal forming problems, especially in complex three-dimension metal forming processes. In addition, the computational efficiency should be improved on the basis of ensuring the precision. Reliable and user-oriented EFG based analysis software for bulk metal forming simulation

[1] T. Belytschko, Y. Krongauz, D. Organ, et al., Comput. Methods Appl. Mech. Eng. 139 (1–4) (1996) 3–47. [2] L. Lucy, Astron. J. 8 (12) (1977) 1013–1024. [3] T. Belytschko, Y.Y. Lu, L. Gu, Int. J. Numer. Methods Eng. 37 (2) (1994) 229–256. [4] T. Belytschko, D. Organ, Y. Krongauz, Comput. Mech. 17 (3) (1995) 186–195. [5] W.K. Liu, S. Jun, Y.F. Zhang, Int. J. Numer. Methods Fluids 20 (8–9) (1995) 1081–1106. [6] C.A. Duarte, J.T. Oden, Numer. Methods Partial Different. Eqs. 12 (6) (1996) 673–705. [7] J.T. Oden, C.A.M. Duarte, Comput. Methods Appl. Mech. Eng. 153 (1998) 117–126. [8] S.N. Atluri, T. Zhu, Comput. Mech. 22 (2) (1998) 117–127. [9] T. Zhu, J.D. Zhang, S.N. Atluri, Comput. Mech. 21 (3) (1998) 223–235. [10] N.R. Aluru, Int. J. Numer. Methods Eng. 47 (6) (2000) 1083–1121. [11] G.R. Liu, Y.T. Gu, Struct. Eng. Mech. 11 (2) (2001) 221–236. [12] X. Zhang, X.H. Liu, K.Z. Song, et al., Int. J. Numer. Methods Eng. 51 (9) (2001) 1089–1100. [13] N. Sukumar, B. Moran, A.Y. Semenov, et al., Int. J. Numer. Methods Eng. 50 (1) (2001) 1–27. [14] Y. Cai, H. Zhu, Eng. Anal. Boundary Elem. 28 (6) (2004) 607–613. [15] P.W. Randles, L.D. Libersky, Comput. Methods Appl. Mech. Eng. 139 (1996) 375–408. [16] M.G. Abadi, D.G. Lambas, P.B. Tissera, Int. Astron. Union-Symp. 168 (1996) 577. [17] W.K. Liu, S. Jun, D.T. Sihling, et al., Int. J. Numer. Methods Fluids 24 (12) (1997) 1391–1415. [18] A. Arefmanesh, M. Najafi, H. Abdi, J. Fluids Eng., Trans. ASME 127 (4) (2005) 647–655. [19] T. Belytschko, M. Tabbara, Int. J. Numer. Methods Eng. 39 (6) (1996) 923–938. [20] W.K. Liu, S. Hao, T. Belytschko, et al., Comput. Mater. Sci. 16 (1999) 197–205. [21] T. Belytschko, D. Organ, C. Gerlach, Comput. Methods Appl. Mech. Eng. 187 (3) (2000) 385–399.

212

P. Lu et al. / Materials Science and Engineering A 479 (2008) 197–212

[22] C.A. Duarte, O.N. Hamzeh, T.J. Liszka, et al., Comput. Methods Appl. Mech. Eng. 190 (15–17) (2001) 2227–2262. [23] L. Gao, K. Liu, Y. Liu, CMES—Comput. Model. Eng. Sci. 12 (3) (2006) 181–195. [24] J.S. Chen, C. Pan, C.T. Wu, et al., Comput. Methods Appl. Mech. Eng. 139 (1–4) (1996) 195–227. [25] S. Jun, W.K. Liu, T. Belytschko, Int. J. Numer. Methods Eng. 41 (1) (1998) 137–166. [26] T. Rabczuk, T. Belytschko, Comput. Methods Appl. Mech. Eng. 196 (29–30) (2007) 2777–2799. [27] J.S. Chen, C.M.O.L. Roque, C. Pan, et al., J. Mater. Process. Technol. 80–81 (1998) 642–646. [28] J.S. Chen, C. Pan, C.M.O.L. Roque, et al., Comput. Mech. 22 (3) (1998) 289–307. [29] J.S. Chen, H.P. Wang, Comput. Methods Appl. Mech. Eng. 187 (3) (2000) 441–468. [30] J. Bonet, S. Kulaseqaram, Int. J. Numer. Methods Eng. 47 (6) (2000) 1189–1214. [31] S.W. Xiong, W.K. Liu, J. Cao, et al., Int. J. Mach. Tools Manuf. 43 (2003) 89–102. [32] S.W. Xiong, J.M.C. Rodrigues, P.A.F. Martins, Euro. J. Mech. A/Solids 23 (2004) 77–93. [33] S.W. Xiong, C.S. Li, J.M.C. Rodrigues, et al., Finite Elem. Anal. Des. 41 (6) (2005) 599–614. [34] S. Xiong, W.K. Liu, J. Gao, et al., Comput. Struct. 83 (2005) 574– 587.

[35] S.W. Xiong, P.A.F. Martins, J. Mater. Process. Technol. 177 (1–3) (2006) 49–52. [36] Y.H. Park, J. Mater. Process. Technol. 183 (2–3) (2007) 256–263. [37] Y.H. Park, ASME, PVP, Comput. Technol. Appl. 458 (2003) 231–236. [38] Y.M. Guo, K. Nakanishi, J. Mater. Process. Technol. 140 (1–3) (2003) 19–24. [39] Y.M. Guo, K. Nakanishi, Y. Yokouchi, Adv. Eng. Software 36 (4) (2005) 234–242. [40] I. Alfaro, D. Bel, E. Cueto, et al., Comput. Methods Appl. Mech. Eng. 195 (33–36) (2006) 4269–4286. [41] K.C. Kwon, S.H. Park, S.K. Youn, Int. J. Numer. Methods Eng. 64 (2005) 751–788. [42] K.C. Kwon, S.K. Youn, Int. J. Solids Struct. 43 (25–26) (2006) 7450–7481. [43] G.Y. Li, K. Sidibe, G.R. Liu, Eng. Anal. Boundary Elem. 28 (10) (2004) 1283–1292. [44] Y.H. Liu, J. Chen, S. Yu, et al., Acta Metall. Sin. 19 (5) (2006) 371–378 (English letters). [45] H. Wang, G.Y. Li, X. Han, et al., Adv. Eng. Software 38 (2) (2007) 87–101. [46] Q.L. Cui, C.S. Li, X.H. Liu, et al., Chin. J. Nonferrous Met. 16 (2) (2006) 322–326 (in Chinese). [47] J.S. Chen, C. Pan, C.T. Wu, et al., Comput. Methods Appl. Mech. Eng. 139 (1996) 195–227. [48] J.S. Chen, S. Yoon, H.P. Wang, et al., Comput. Methods Appl. Mech. Eng. 181 (2000) 117–145. [49] S. Kobayashi, S.I. Oh, T. Altan, Metal Forming and the Finite Element Method, Oxford University Press, New York, 1989.