Adaptive ALE finite elements with particular reference to external work rate on frictional interface

Computer Methods in Applied Mechanics and Engineeiing 93 11991) 189-216 North-Holland

Adaptive A L E finite elements with particular reference to external work rate on frictional interface* Wing K a m Liu, Jiun-Shyan Chen ~, Ted Belytschko and Yi Fei Z h a n g 1Vorthwestern University, Department of Mechanical Engineering, Evanston, Illinois, 60208, USA Received January 1990 Revised manuscript received July 1990 Simulation methods for material forming processes based on arbitrary Lagrangian-Eulerian (ALE) finite elements and adaptive mesh deployment techniques are developed. Special emphasis is on the ALE formulation with respect to the external virtual work rate on the friction interface. An extension of the solution algorithm from the explicit scheme to Newton iteration is also given. Successful strategies for the use of ALE finite element models in nonlinear problems are explored. The implementation of ALE finite element modeling to rolling material forming processes with friction and the fusion of ALE methods with adaptive finite elements is explored. Applications to metal rolling problems are then given.

I. Introduction W h e n a material is formed by processes such as forging, cutting, sheet pressing and molding, it experiences large unrecoverable deformation. This deformation leads to the d e v e l o p m e n t of zones of high strain concentration and, consequently, the onset of internal (or surface) cracks. This strain localization is the cause of m a n y defects, such as wrinkling and buckling, puckering, necking and tearing. In general, the formation of these defects cannot be identified easily during the design of the forming process. In solids, the strain localization is often called shear bands. Shear bands, which can be described as regions of inhomogeneity in which intense shear strains are developed, are formed when the material flows at different strain rates (or velocities) during the forming process. A n accurate forming simulation should provide information on the forming force, strain and stress distributions; the undesirable processes which lead to localization deformations might then be avoided. The deformation which occurs in forming processes exhibits characteristics of both solid and fluid flow. T h e r e f o r e , either a fluid mechanics approach or a solid mechanics a p p r o a c h are

* The support of this research by NSF Grant EET-8806347 and the National Center for Supercomputing Applications at Urbana-Champaign is gratefully acknowledged. Presently at GENCORP, Research Division, 2990 Gilchrist Road, Akron, OH, 44305, USA. 0045-7825/91/$03.50 © 1991 Elsevier Science Publishers B.V. All rights reserved

190

W.K. Liu et al., Adaptive A L E finite elements

applicable and indeed may be used interchangeably to model the same physical process. In the development of finite element methods for forming simulations based on a :,olid mechanics approach, the Lagrangian description has been employed extensively for moderate deformation and moderate rotation analyses. This method is popular due to the fact that constitutive equations for most solids, ranging from common alloys to highly exotic metals, are historydependent. In Lagrangian meshes, convective effects are absent since the material points coincide with the finite element mesh quadrature points throughout the deformation, so that history dependent materials present no difficulties. However, for large flow simulations, the Lagrangian description is not satisfactory because large material distortions leads to severe element distortions and entanglement. At the present time, interactive mesh rezoning techniques have been introduced into Lagrangian codes to partially alleviate this problem. However, rezones require stress and strain interpolations to initiate the new mesh. Furthermore, a complete computer simulation might require many rezones and stress and strain interpolations. The programming effort is enormous and awkward. We do not follow along this line, since we believe that rezones and interpolation of stresses and strains can be incorporated automatically within the finite element analysis. The mesh distortion problem can be overcome by using a fluid mechanics approach, the Eulerian description. An Eulerian description is convenient in modeling a fixed region in space for situations which may involve large flows, large distortions and mixing of materials. However, convective effects arise because of the relative motion between the material and the fixed mesh, and these introduce diffusion in the material history and other numerical difficulties. Furthermore, the material interfaces and boundaries which may move through the mesh require sophisticated mathematical mappings between the stationary and moving boundaries. It is thus clear that material forming processes are not well-suited to either Lagrangian or Eulerian meshes because of their intrinsic limitations. To circumvent these difficulties, the arbitrary Lagrangian-Eu!erian (ALE) formulation for path dependent materials has been developed. In this description, the marriage of the Lagrangian and Eulerian formulations is effected through the intr,oduction of a mesh motion which is represented by an additional set of independent variables. The use of ALE meshes is well established in the areas of hydrodynamics and fluidstructure interaction. This concept was first proposed in [1] under the name 'coupled Eulerian-Lagrangian', a finite difference scheme for two-dimensional hydrodynamics problems with moving fluid boundaries. Similar techniques have been extended to compressible flow problems in [2]. Later, in [3], the ALE method was employed in conjunction with an implicit formulation for the solution of two-dimensional flows. A version of this computational technique for three-dimensional flows contained within arbitrarily shaped moving boundaries is reported in [4, 5]. The ALE method was introduced into the finite element method in [6, 7] in response to the need for nonlinear simulation techniques for nuclear safety analyses. The advantages of the ALE method in fluid-structure interaction are apparent since the fluid domain can be treated by the ALE formulation and the structural domain can be handled by the usual Lagrangian description. In the above articles, the effort was primarily directed toward inviscid compressible fluids, while in [8] a finite element procedure for viscous incompressible flows and free

W.K. Liu et al., Adaptive A L E finite elements

191

surface flows is presented in conjunction with a general kinematical theory for the ALE description. The capability of ALE method to handle an expanding gas bubble immersed in a fluid has been demonstrated in [9, 10]. The mesh motions for this problem are prescribed according to a simple algorithm in which the boundary nodes are considered Lagrangian and the mesh velocities for the intermediate nodes are linearly interpolated between the boundary velocities. A general ALE formulation was later developed by Liu and Gvildys [11] and Huerta and Liu [12], where the general theoretical framework is established, and the transient algorithm and computer implementations are detailed and applied to three-dimensional liquid storage tanks. Recently, the ALE concept was applied to contact problems between flexible structures in [13], which subdivide the displacement into Lagrangian and Eulerian parts. The slip compatibility conditions are met by making the Lagrangian displacements common to elements on either side of the interface and using separate tangential Eulerian displacements. This concept is extended in [14], for quasi-static solid mechanics and to elastic fracture mechanics. The A L E literature has mainly been directed toward linear-path independent materials like Hookean solids and fluids. The stress states for these materials are solely determined by the instantaneous displacement or velocity fields. When the ALE method is applied to materials with memory, the state variables for an element are affected by the migration of material points which may carry different stress and strain histories. This difficuity arises because the A L E mesh, just like an Eulerian mesh, does not model the same set of material points, through the evolution of the problem. A similar situation is encountered in [15], where the Eulerian description is employed for plastic forming analysis. In this paper, the difficulty is side-stepped by interpolating the stress histories at each increment step. A related numerical technique for the interpolation of the stresses is to use an approximation to the conjugate basis of the stresses [16]. More recently, Benson [17] has implemented an operating splitting approach to the ALE formulation; some example problems are presented to illustrate the strengths and weaknesses of the formulation. A first attempt to develop a general formulation for materials with memory for nonlinear ALE finite element analyses is reported in [18, 19], this study was in the framework of explicit time integration schemes where no linearization is needed. The discrete conservation laws, the constitutive equations, and the equations of state for path-dependent materials are formulated for an ALE finite element method; both the geometrical and material nonlinearities were included. In this paper, the linearized forms of the equations are presented so that the methods can be extended to Newton solution methods. These types of solvers are advantageous for quasistatic problems such as metal forming, where it is difficult to achieve convergence with vector-implicit methods, such as conjugate gradients. The ALE formulation with respect to the external work rate on the friction interface is also emphasized. In the next section, fundamental equations for equilibrium, constitutive law and friction are reviewed. In Section 3, the ALE surface update equation is discussed. The referential rate expressions and Petrov-Galerkin finite element formulations are derived in Section 4. Moving mesh techniques are discussed in Section 5. In Section 6, the computer implementation aspects are given. Numerical examples are presented in Section 7 to illustrate the effectiveness of the method. Finally, discussions and conclusions are given in Section 8.

W.K. Liu et al., Adaptive A L E finite elements

192

2. A L E f o r m u l a t i o n and f u n d a m e n t a l e q u a t i o n s 2.1. N o t a t i o n a n d p r e l i m i n a r i e s

The material coordinate X, which is attached to a material particle, can be treated as a 'marker' of a material particle when it is in motion. One can describe the material motion by a function which gives the position of the particle at different times, i.e., by using x = 4~(X, t ) ,

(1)

where x is the spatial coordinate and ~ is the material motion function or the mapping function between X and x. The motion of the finite element grid need not coincide with ~ ; it can be controlled by the observer (the finite element user) in accordance with the aser's judgment. To describe the governing equations for a finite element mesh which moves independent of the material, a coordinate system, associated with the finite element nodes is needed; it is called the referential coordinate X. The finite element motion in space is then described by x : ~(X, t ) ,

(2)

where ~ is the mesh motion or the mapping between x and X. When the material coordinate, X, and the referential coordinate, X, are taken to be initially coincident with the spatial coordinate x = ~(X, 0), then X and ,!" can also be recognized as the original material and mesh configurations. The material and mesh displacements u and ~i are respectively defined as u = x-

X,

ti = x -

X-

(3), (4)

Throughout this paper, standard indicial notation is adopted; lower case subscripts denote the components of a tensor and repeated subscripts imply a summation over the number of space dimensions (NSD). The 'spatial derivative', i.e., derivative with respect to x, is designated by a comma followed by a subscript; t denotes time. In material balance laws, the material time derivative of a function f is required, which is usually defined as a function of x and t. In an ALE formulation, since X is our computational reference coordinate, it is convenient to express these quantities in terms of X and t. The time derivatives with the material (X) and referential (X) coordinates fixed are denoted by a superposed dot and prime, respectively. The relationship between these two can be established using the chain rule (sce [8]): f = f ' + wk Of/Ogk ,

(5)

where wk = i(k ,

f'=

Of/Otll~ 1 •

(6a,b)

Replacing f by xi, (5) becomes c i = w~ OX,/dXk ,

(7)

W.K. Liu et al., Adaptive A L E finite elements

193

where c, = v , - 6,,

v, = ~ , ,

6, = x;.

(8), (9), (10)

The terms c~, v, and 6, are the convectiw: velocity, material velocity and mesh velocity, respectively. Finally, substituting (7) into (5) and applying the chain rule, the following relationship between the material time derivative and the referential time derivative is obtained:

(11)

f'=f-cif.,.

The term c i f , represents the convective effect due to the relative motion between the material and the mesh. REMARK 1. The quantity wi defined in (6) is recognized as the particle velocity viewed from the referential coordinate. Therefore (7) indicates that the relative velocity between the material and the mesh can be interpreted as the material velocity viewed from the referential coordinate with the relative deformation Ox,/OXk taken into account. 2.2. Basic equations Equilibrium. A quasi-static approach in which the inertia is neglected is appropriate for forming processes such as the low-speed rolling. With the same reason, the body force can be neglected, and the equilibrium equation becomes

r~i.i = 0,

(12)

where ru is the. Cauchy stress tensor. Constitutive equation. In this type of problem the velocities are more important than the total displacements and the rate form for constitutive law is therefore employed. The constitutive equation is therefore written as a rate type of equation in the form [20] (13)

7'u = --CkZ, i.k + C,jklV(k,, ) + S,j~,vlk., I .

In this equa63n, --Ck%.k is the transport of the stress, C,jk~vtk.~ ) and SiyktVlk.t I are the pure deformation and the rotational parts, respectively, for the rate of change of the Cauchy stress. Cij~t describes the characteristics of the material which ,=an be either solid or fluid, elastic or elastic-plastic. If the Jaumann stress rate is used, c C, ikl = C,jkl ,

(14)

c where C,jkt is the co-rotational constitutive tensor. If the Truesdell stress rate is employed,

C,ikt = C~ik, + C~k , ,

where C~, = --T,iak, + ½(~'a6ik + ~,6~k + ~',k6/, + ~'jkaa), (15), (16)

t

and C,jk~ is defined by

W.K. Liu et al., Adaptive A L E finite elements

194

V ~ I = CiiklV(k.i , ),

(17)

•

t where VTq are the components of a frame invariant incremental Cauchy stress tensor, with the current state taken as the referential state. O t h e r stress rate constitutive equations can also be used by defining C~k~. The term Sukwikj ~ is the rotational part for the rate of change of the Cauchy stress, and the fourth order tensor Suk~ is given by

s,,,, =

+ *j,¢k --

(18)

8j,-

It is noted that this fourth-tensor is the same for all rate type constitutive equations. 2.3. Interface condition When two bodies are in contact, friction forces develop due to the relative motion of the two bodies. In order to employ the A L E approach to problems with interfaces such as rolling simulations, the referential rate of the frictional traction has to be developed. An interface element technique or a friction model m e t h o d is usually used to handle the friction between two materials. This technique involves inserting an extra layer of elements on the interface where friction acts. The stiffness of each of these elements is multiplied by a factor m which is determined by the yield stress. The second method involves introducing a friction model, where the friction force is determined by the relative motion between the material and the interface and is applied directly on the elements. Two studies of the latter method are reviewed, the first is the classic friction model, and the second is from O d e n and Martins [21] who offer a regularized friction model based on the empirical data. Classic interface condition During the rolling process, the material may be divided into non-contact and contact regions. When the surface is in contact, the velocity continuity in the normal direction has to be satisfied. N

N

v = v~ ,

(19)

N

where v N and vs are the material and boundary normal velocities on the interface, respectively. In the tangential direction of the interface, the contact can be further classified into 'slip' and 'stick' categories. Since the Coulomb friction is assumed, the slip-stick condition can be determined by slip

if I/T] ~ ~;r]!NI,

stick

if [tt[ < 6~[tn],

(20), (21)

in which 6 r is the static friction coefficient, and l N and t T are the normal and tangential components of the surface traction on the interface, respectively, t N and t v can be written as t N = t . n = ('gjnj)n i

and

t t = t . n t = (.g~nj)nVi,

(22), (23)

where t is the surface traction on the interface, n is the surface unit normal, and n v is the unit

W . K . Liu et al., Adaptive A L E finite elements

195

vector in the direction of relative tangential velocity, i.e., /,/T =

6T/II TII,

6T

=

T

~.)T __ l o B ,

(24)

where l0 T , l0TB and ~ T are the material, boundary and relative material velocities, respectively. The mathematical formulation of the stick and slip conditions is as follows: (1) In the stick region, velocities are prescribed to be the same as the boundary velocities. VN =

N, 10B

vT =

T" loB

(25)

(2) In the slip region, velocity and friction force are prescribed in the normal and tangential directions, respectively. I/~N

N = l0 B ,

tT

= -Glt

N

In T •

(26), (27)

Note that t N < 0 on the contact region; therefore, (27) can be rewritten as (28)

tT = cdNn x ,

where c r is the kinetic friction coefficient. In summary, velocities are prescribed in the stick region and in the normal directior~ of tht. slip region. The union of these two is classified as an essential boundary and is denoted a s if'Ixg. However, in the tangential direction of the slip region, a friction force is applied. This is classified as a natural boundary and is denoted as F~xt. O d e n a n d .Martins m o d e l

The major computational drawback of the classic model is the discontinuity of Coulomb's friction law at zero sliding velocity. Oden and Martins [21] proposed a constitutive law which takes into account the deformability of the interface asperities. They considered an asperity covered surface contacting a perfectly fiat solid. The highest asperity is initially a distance g (the initial gap) from the rigid surface. They defined the approach of the surfaces as

A

=

0,

tiN
tiNg,

bt

-N

(29), (30)

~>g,

where tin is the relative normal displacement. Based on experimental data, Oden and Martins then define the constitutive properties of the interface by a power law relation between the normal stre ~ t N and the surface approach A by t N = -t~A p

,

(31)

in which a and p are material properties independent of A. According to experimental observation, a preliminary microdisplacement occurs before gross-sliding. Therefore, the distinction between static and kinetic coefficients of friction need cot be considered and the stick condition L,; removed by regularizing the frictional stress as

196

W.K. Liu et al., Adaptive A L E finite elements

tT

_[3AJ(2-llv

T - vglll~)(v T - v~)l~ (v ~ - v~)lllv T- o~II

v if llv ~ - ~BII <~ ~,

(32)

if llv T - vBll > ~,

(33)

v where v T and v B are the material and boundary tangential velocities, t T is the frictional traction, 13 and q are material parameters, and e is the regularization parameter. Consequently, the kinematic continuity conditions which have to be imposed on the interface are replaced by the need to modify the global force vector by adding in the appropriate friction and normal force terms.

3. Adaptive ALE surface update equation Since automatic mesh generation can be viewed as a boundary value problem, the initialization of the nodes on the boundary is crucial to the location of the internal nodes. The Lagrangian type boundary, where the boundaries move with the particle velocities, are simple and accurate for tr~,e free surface updating. The accurate mapping of the moving boundaries, which are the material surfaces, is an important task in A L E description, see e.g. [12]. Consider a material particle and a finite element mesh are mapped to the same position on a surface in the space at time instant t. This mapping can be described by x = x(t, X) = x(t, X),

(34a)

where X and X are represented by

x = x ( a , , ,/2),

x = x ( a , , ,LO,

(34b,c)

where (d I, d2) and (&l, &2) are the surface parameters in the material and referential coordinate systems, respectively. Equation (34a) can be rewritten as x = x(t, X) = x(t,

X(ff,, ~2)) ~ aT(t, ~,, 82),

(35a)

x = x(t, X) = x(t,

X(a,, a,)) - ,(t, a,, & ) .

(35b)

The motion of the mesh with respect to the material can be defined by (36a,b) So, we have .r(t, 8,, ~2) = a~(t, &,(t, a,, c~2), &2(t, a , , if2))"

(37)

The material velocity is defined by aa7 V='~-

=

_ = 0--7 lat,a.,I

a& i at

= t3 + ~/Oei .

(38)

197

W . K . Liu et al., Adaptive A L E finite elements

Note that (bJ~/0~) (i = 1, 2) are the surface tangent vectors. By taking an inner product of (38) with the surface normal n, we have o.

n = ft. n + n.

utt i

Since n-~=O 0%

(39)

-g77~-. a i "

fori=l,2,

(40)

(39) becomes (using (38)) c. n = 0.

(41)

The physical interpretation of (41) is that no normal convective velocity occurs across the surface if the surface particles remain on the surface.

EXAMPLE 1. If an Eulerian description is used in x~ direction, then 6~ = 0, and t3z is an unknown to be solved by (41): (ol-Oi)nl+(o,--Oz)nz=O -

and

6, = -nl- v, + u~ -

n,

- '

(42), (43)

and the mesh position of the surface can then be obtained by time integration of (43).

4. T h e rate f o r m o f the e q u i l i b r i u m e q u a t i o n for p r o b l e m s w i t h frictional b o u n d a r y

4.1. The weak form To discretize the system ¢~f partial differential equations, (12), a weak form of the equilibrium equation (12) is used:

fR ~Uirij.j dR x = O,

(44)

where Rx is the spatial domain and ~u~ is the test function. By applying the divergence theorem, a virtual work statement is obtained:

fn6ui.jr# dR, = fr; 6uit~ dFx or ~W i"' = 6W

TM

,

(45), (46)

x

w h e r e t i is the surface traction and Ftx is the traction boundary. Note that the traction

boundazy in a rolling problem, F'x, contains two parts: the free surface F~ (traction free) and the interface/'ixt (friction). We consider a rate form of the problem which can be obtained by taking the time derivative of the virtual work statement (eq. (45)). The referential time derivative of ~W i"' [19] is

W.K. Liu et al., Adaptive ALE finite elements

198

( a w t ) i n l = fR algi'j['l'fiJ -I- TijOl:,k -- "l'ikOj,k] d R

a

x

=fR,.~)Ui,jOuklOk,ldRxq'-fR~3Ui,jTiiklOk,ldRx (47)

+ fR 8Ui'i(TikCj'k-- 7iiCk'k -- Ckru'k) dR, , in which

CiTk,

Dqkt

for Truesdell r a t e ,

{ Cijkl - Ciikl

(48a)

for Jaumann rate

and (48b)

Tijkl = "l'il6ik .

Note that the last three terms in (47) are the transport effect which arises from the A L E description. The tangential stiffness matrix can be obtained when the shape functions are introduced. The ~W TM in the present work is due to friction. Both classic and Oden-Martins friction models are to be investigated. Classic model. The expression of ~W TM for friction is obtained by substituting (28) into the right-hand side of (45): ~WC"==£

-.Irt

8uitidF~.=Crfr,?~ui%,,n,,n,,,nT dr,,

(49)

where F.{.t is the slip region of interface, n is the unit surface normal and n v is the unit tangent which is defined in (24). To obtain the referential time derivative of ~W TM, Nanson's Law (see [22]) has to be introduced:

n,, dI i= JN~ ~aXk dF~ ,

(50)

where n and N are the surface normals in the x and X coordinate systems, respectively, and

J

=

de,t (Oxi/Ox/).

(51)

By substituting (50) into (49), ~W TM can be expressed as follows: SW TM = c r .,,~Suir,,,,,N k ~OXk n,,,nTj dr,,,

(52)

and (SW') TM can then be obtained: (~W')ext - Cr

f,

"h ~ U iN

{

OXk

t

{ OXk ~

T.

OXk

, T

7'ran -~xn l~'nn i J q- ~'mnk -~X n ) Yl'nl:li J "{- "l'm',~'~Xn n m l:li J

ox~ n,.(nT),j +

OXk n,, T :, ] dr~

(53)

W.K. Liu et al., Adaptive ALE finite elements

199

Two identities are now employed:

(aX'k)_ Ox,,

aXk adp Oxp Ox,, and

J' = Jd~. k .

(54), (55)

Equations (54) and (55) can be substituted into (53) to yield

[ ,

OX~ n,,,n:J- r,.,, OXp OXk Ox. Odp nmn~J

( ~ w ' ) e x t = Cr "~' ~uiNk rm" ~

OXk n,.ni , T J + r,.. -~x,, OXk nm(nT)%r~?rm,, ~-xf OX~' I,,,,,.,T^v . . # ] d r . , + rm. -~x,

(56)

which can be transformed back to the spatial domain as i T + r,,,,,(vk.kn,.n ^ T~ i,- ' t (BW') T M = c~ ~/ .,. gu~ [ z,,,,,n,,n.,n~ n,,,n i + nm(n ti ) , )n,,

-

Odp ] r,,,, ~ n,,npn~ OE~.

(57)

The time derivative of the normal and tangential unit vectors can be derived as follows. Let s be the curvilinear coordinate along the path of the mesh motion on the inteHace, then

dn x

_ d n t ds

1

ds txl

1

(58)

and

n,

=

dn ixl

- -

~-

dn d~ l,l c~s

=

(1 b - -

r

- -

1 n T ) ~T p

(59)

- -

'

where p is the radius of curvature, r is the radius of torsion, d t is the mesh speed along the path v T = d~tt I~1 =

116TII

(60)

and b is called the bi-normal to the curve which is defined by b = n+ × n.

(61)

R E M A R K 2. d t can also be recognized as .~he magnitude of the mesh ~tangential speed on the interface. Equations (58) and (59) can be rewritten as f.r (niT)' = P n , ,

d'r d'r n~ = --r b~ -- P niT.

(62), (63)

200

W. K. Liu et cd., Adaptive ALE finite elements

By substituting

(62) and (63) into (57) we get

(~~‘)err = c, I, =C

r

I

,.;,

W&;,

+ w&,,

sUi{[(Dm”k~

+

+

Tm”k,)Uk.,

TmnBrnnil

d’x

+

(7mkc,,.k

+ ~,,J,,:,,il dT,

-

‘,,&k.k

-

CkT,,,,,.k.)b,,n,,n~

(64)

7

*T

..T

Bmni = 5

b,n,n;

+ ;

REMARK

‘ninBmni

- r,,&&m~,~~

3. By introducing can be simplified as

(n,rz; - n;fin;)nt,.

the traction

(65)

equilibrium

equation

on the interface,

the term

-T

=-c,lL~

‘mtl Bn,ni

and the derivatives

mnnmnn7 I, I

P

(661

are shown in Appendix

Oden and Martins model.

A.

The surface traction,

t, on the contact interface can be decom-

posed into two parts: t = tN + tT

or

ti = t” + t: ,,

(670)

where t” = --aAPn. I 3

t; = -pA”cp,

)

(6%

(69)

where (2 - IIUTII/c)UT/.s,

if /Iti’ll s

E ,

UT/ Ilv”ll ,

if IItiT(I >

E ,

(70)

A =

where E is the regularization interface: -T V

=:u T -_u

parameter

and VT IS . the relative

tangential

velocity

on the

T B’

(71)

The SweX’can then be written as awe”’ = _ 5rl’

Gui(aAPni + PA”&) dT, .

The surface transformation rule, (50), can be simplified by taking the dot product with both sides of the equation: dT,=JN,$r,dT,. (73) can be substituted

14

on

(73)

” Equation

(72)

into (72) to yield

W.K. Liu et al., Adaptive ALE finite elements aW

TM

201

= al,V~" + a W ~ / ' ,

(74)

where

OXki u= ' Fx a W eN~ t = - fr ~"auit~A'JN k ~x

and

ext

1"

0Xk

~WT = - J r " ~ u i f A q ~ i J N k - ~ x . , n "

~

alx"

"(75),(76) By following the same procedures as in the previous derivation and applying the orthogonality conditions for n, n T and b, (~W') T M can then be obtained: (~W') T M = (~W~) ext + (SW~) ~' ,

(77)

where (~W~)e"t = - f r " ~uia[(PAp-lA' + AP6*'k)n' - APvk'ink] dE" =

-frl' Buia[(PAP-lA'+ A P v k ' k ) n i

-- APvk'inl< -- A P C k ' k n i

+ AOck.?~k] d r x ,

(78)

(~W~r)e~t= - fr ]: ~u~fl[qA q-l A , dp~+ Aqdp,' q- Aq~piOk,k -- A q ~ i O k , l n k n l ]

dF~

= -jr~, rau,~[qAq-~A'¢h, + Aq4,~ + A q4~iVk.k -- Aq4~ivk.tnknt -- Aqqbick.~. + A~4~iCk.tnknt] d/'..,.

(79)

and 1 -

n

' =

4~

(6v)

' +

-re

[(n~/r)6

,

2 -

ni 6T

if 6 "r <~ e , '

(80)

if 6 "r> e ,

in which -T

v

= II~TII

(81)

and r is the radius of torsion.

4,2. Rate form o f the equilibrium equation The rate form of the equilibrium equation can be obtained by performing integration by parts on (BW') i"~ and (BW') ext to yield f g ~Ui[T'ij "1- "rijUk. k -- ~.ikl~j.k],j dRx x

= Jr'; au,(Ir',, + r,,~k., - r,k~,.~ln, + ~} dC + ~r au,tr',, + ~,,~k.~ - r,~O,.kln, drx, (82)

202

W.K. Liu et al.. Adaptive ALE finite elemenL,

where = - c r { [ ( D m , k , + Tm,k,)Vk., + ("r,,kC,,.k -- %,Ck. k -- ck~-.k)]n,,n,nV~

+ ~',,,Bm,i}

for the classic model, ]i =: a [ ( P A p - I A ' + APvk.k)ni -- APvk.ink -- APck.kni + At'ck.ink] + fl[(qAq-lA'dgi + A't~b~) + Aq¢iVk.k -- Aq4apk.tnkn, -- Aqcb~ck.k

for the Oden-lV[artins model,

+ Aq
(83)

where F~ denotes the free surface boundary. Equation (64) leads to a rate form of equilibrium: (~'ii

+

"r~J'Sk.k -- r~kdJ.k).J = 0

in R x ,

(84)

the friction boundary condition

+

-

on r]'

=

(85)

and the free surface boundary (traction free) (~-'~j+ r,jOk.k - ~'ikOi.k)nj = 0

on F : ,

(86)

where n v, g, ~bi, B,,,, are defined in Section 4.1. The rate form of the equilibrium equation in the ALE description derived in the previous section can be restated as foflows. Given c~ (or a, fl, p, q), v~:(x, t), t b - t ( x , t), T(x, t), determine v(x, t) and ~-'(x, t) such that they satisfy a. equilibrium: eq. (84),

(87)

b. traction boundary condition: eq. (85)

and

eq. (86),

(88), (89)

c. prescribed displacement boundary condition: v(x, t) = vg(x, t)

on F~g ,

(90)

d. mesh motion: (1) Let T be the mesh motion function in i and j directions

6~=T'(x,t), a = i , j

on F~

and the mesh motion in k direction is governed by the surface update equation

(91)

W.K. Liu et al., Adaptive A L E finite elements

1

1

6 k = - - ( n , , c , , ) + v k = - - [n,~(v,~ - T ' ) ] + v k , nk tlk

203

a=i,j

on E~,

(92)

(2) mesh motion inside the domain 6 = X ' = [ ~ b ' ( x , 0] -I 4.3.

ALE

in R x.

(93)

stress update procedure

In (13), the term Ck"ril,k represents the convective effect due to the relative motion between the material and the mesh. In nonlinear finite element problems, lower order interpolation functions are preferred. The element stresses are u~ually C-1. The ambiguity of the stress derivatives at the element interface renders the calculation of the spatial derivatives of stress a difficult task. To remedy this situation, (13) is reformulated by introducing a stress-velocity product (94)

Yqk = L j c , .

Chain rule on (94) and substitute it in (13) (i.e., the term Ckrq.k) yields t

(95)

Tij q- Yqk.k -- Ck,k Tij -~- CijklO(K.t) + SijktUlkd l •

where the stress-velocity product and the initial stress at time t = 0 %(X, 0) = %"

(96)

is given. In order to apply the Petrov-Galerkin finite element method to (94) and (95), these equations must be transformed into equivalent weak forms. Here, the Petrov-Galerkin formulation is used for the stress-velocity product term (94), while the Galerkin method is applied to the Cauchy stress update (95). In a finite element computation, the stresses are evaluated at M quadrature points. Therefore, within a finite element, M subdomains are defined so that no two subdomains overlap and each subdomain, designated by R~, contains the quadrature points x 8. Associated with R B, a stress interpolation function N~ is assigned and its value is prescribed to be unity only at the quadrature points x = xB; that is, for the element, M

7"q=

Z

NRrqB

and

N~(x=xa)=l

no sum on B .

(97),(98)

B=I

The stress weighting function in R~, A = 1 . . . . . that for each element

M is chosen to be the Dirac delta function so

M

~zq= ~'~ N ~ ' # A

and

N A=~(x-x,),

A=l

....

,m.

(99),(100)

A=I

Substitution of these functions into the constitutive equation represents a mathematical

W.K. Liu et al., A ~/~.?~..':~" .ALE finite eLments

204

requirement that the residual of the weak fi:~.,r~ varfi¢".j,:,~~.~ ~zac?~collocation quadrature point. Because the collocation point is coincident with the q~ .~d:~ture point, the algebraic equations resulting from the Galerkin formulation of (95) are dependent only on the stress history associated with its host points. Since the Dirac delta function has the important property that

fR F(X)~(X-- XR) dR,~ = F(XB) B

no sum on B ,

001)

where F(x) is a continuous function and Rx is the element spatial region, all the matrices and vectors resulting from (95) can be easily worked out without numerical integration. To complete the Galerkin formulation of (95), the stress-velocity product is interpolated by a C O function so that NEN

Yiik---- Z

B=I

(102)

NByijke,

where N R and NEN are C ° functions and the number of nodes per element, respectively. The same interpolant is used for the convective and material velocity. The integral equation associated with (95) is

-- ~e fR ~ ~'l";{CijklV'k")dP'l'kjV'i'k] "~ TkiV[j'k]} d g x = o

(103)

where tl~e spatial domain, R~, is discretized into element subdomains, R~, and Ze symbolizes the sum over all elements. For the stress-velocity product, (95), the streamline upwind/Petrov-Galerkm formulation requires a discontinuous weighting function of the form (104)

~Yijk : ~tOijk + oPijk ,

where btOUkis continuous in R x and the shape functions are the same as those given in (102) (i.e., Ns), and ~Pu~ is the discontinuous streamline upwind perturbation; ~Puk is assumed smooth in the element interior. The variational equation resulting from (94) is

~.(r~t°ukYukdRx--~IR~t°Uk~#c~.dR~+~fR~PukYukdR~ -- ~e fR: ~Puk~'uCkdRx=O"

(105)

where ~tOUkand ~Puk are chosen to be NEN

~tOijk = aZ= l N A ~£Oijka ,

~Puk=

F~ c, X0(~t°uk) ' -0- " - - ' ~

(106a,b)

W.K. Liu et al., Adaptive ALE finite elements

205

F~ may be defined using a spatial criterion [23]: F~ = (otec~h e +

/211cll 2

(107a)

or a temporal criterion [23] F~= ½At,

(107b)

with c~ = e~ • c ,

cn = e~ • c ;

a s = ± ½ for c e ~ 0 ;

a~ = ± ~ for c~ ~ O, (108),(109),110)

and e e, h e, en and ho being the unit vector and element length in the ~ and 9 direction, respectively. By substituting the above test and trail functions into (103) and (105), the resulting element matrix equations are My=Ns

and

M~'=z-GVy+Ds,

(111),(112)

where M, N, M ~, z, D, z are constructed from the continuous and discontinuous shape functions. Since 8pijk weighs only on the element interior, it does not affect the traction boundary. In our quasi-static equilibrium equation, the inertia term is ignored and the external forces are due to surface tractions, so that the standard Galerkin formulation is applicable to the equilibrium equation.

5. Adaptive moving mesh techniques 5.1. Background One complicating aspect in the application of the finite element method to a problem like rolling processes is the extreme mesh distortion. One class of mesh rezoning methods is the 'rubber band' stretching coordinate technique, in which the interior nodes are moved proportionally to the free surface or interfaces. This concept is first used by Murray and Landis [24] and is quite useful for problems having z primary boundary motion in one direction. Back extrusion is a problem suitable for this method. Another class of the rezoning schemes involves solving a differential equation together with well-arranged boundary nodes as the boundary conditions. This is discussed in the next section and will be applied to the numerical examples. All the mesh generation techniques which have been investigated are aimed at one objective: to provide a mesh with as little distortion as possible. However, one should also pay attention to the loss of self-adjointness of the ALE formulations from transport. A good mesh rezoning may develop a well-shaped mesh but bring in a significant transport effect. We have here used x "~W= x °td + ~(x tri - x°l~) ,

(113)

206

W.K. Liu et al., Adaptive ALE finite elements

where a=e where

'

b-

e=

II~lf'

(114),(115),(116)

v - (x"' - x°'~)/at,

jfold and x '~ are the mesh positions before and after remesh.

R E M A R K 4. From (4), two extreme cases result: (Xtri -- X°td)/At = V ,

I1~11 = 0 ,

~ = 1,

(X~i - X°'d)/At = O,

II~ll =: v ,

x "eW = x °'d .

X new =

(117)

X 'ri

and (118)

5.2. Automatic mesh generation Automatic mesh generation has been studied for m~my years, see e.g. [25]. The Laplace equation scheme based on the theory of mapping the node space (I, J) into real space (x, y) by solving the Laplace differential equation is the most commonly used one. The objective of this boundary value problem is to find the following functions: x(l, J)

and

y(!, J) ,

such that they satisfy the following equations:

OZy

O2X 02X L2(x) = ~

+ -g~ = o,

O2y

15(y) = T ? +

=0

in R]~',

(119)

with the boundary conditions x(I, J) = £(I, J ) ,

y(l, J) =.~(1, J)

on F~',

(120)

where R m and F m are the domain and boundary of the remeshed region and L 2 is the second-order differential operator. Another useful mesh generation scheme consists of solving a fourth order differential equation [26]:

0 4X L n ( x ) - 0210~ = 0,

0 4y L4(y) = ~

= 0.

(121)

Equations (119) and (121) can be solved by the finite difference method with a Gauss-Seidel iteration scheme.

6. Computer implementation aspects If ( ),, and ( ),,÷i denote the matrices at time t,, = n At and t,,+t = (n + I) At, respectively, where At is the time increment, the flowchart of the computational procedure for adaptive

W.K. Liu et al., Adaptive A L E finite elements

207

A L E finite elements for rolling problems is as follows: 1. 2. 3. 4.

5. 6.

7.

8.

9.

10. 11.

Initialization. d,, v., 6o = O, T = 0. d, is the initial material displacement. Time step loop T <-- T + At, if T > T .... stop. Determine contact nodes according to the new position of the roller. Decide slip or stick for each interface node: if [t+l > it is a slip node and t + = c~[tN[ sign (v + - v +) (eqs. (27) and (29)); N "1" = T if It+l < crtr l it is a stick node and v N:- v~3 and v B (eq. (26)). Reconstruct the system of equations based on the stick-slip interface condition. Predictor phase ( d e n S e by a superscript '~'): for stick nodes Ad,, + i = (AdB),, + ~ , for slip nodes Ad,, + i = (Ad ~),, + i, i = v,, At , for other nodes Ad,,+ t = v,, A t . If an A L E mesh generation is required, on the contact surface, the stick nodes are Lagrangian; whereas the slip node~ are equally spaced between the two adjacent stick nodes. Automatic mesh generation in the A L E region. Determine the position of the interior nodes by solving Laplace equations (119) and (120) or the fourth order equations (121). D e n o t e x ''~d = X,,, and x ' ~ = solutions of the new mesh, x "~W is then defined according to (113). Update Cauchy stress, back stress, yield radius and plastic strain by considering the convective effect (i.e. c - v , , - (x .... - x"ld)/At). y = I Us,, (eq. (111)) ( s ' ) , + | = ( M ' ) - | [ z - G~y + Ds,,] (eq. (112)), s,,+| = s,, + ( S r ),,+| ~t. Iteration loop, i is the iteration counter and for i = 0 all,+ | = d,,; and Adl, + | = Ad,,+ ~ or Ad,,+ ~ according to step 6. Construct the local coordinate system and form rotational matrix R for contact nodes, that is,

cr[tN[,

K ' n ~--RVK'"R

for contact region,

where K t"n is the tangential stiffness matrix arising from (47), and ~-- means replaced by. t2. F o r m f i"' ~ q d f TM resulting from (45) on the spatial coordinate system. For contact region:

13. Solve the equilibrium equation K'"" Bd =fc~t _ f i , t . 14. Update:

i+1 = Adi,+l + B d , Ad,,+l

di+: ;+l l , .+1 =d,,-¢ Ad,,+

;+1 = A./;+~ v;,.~ -,,+~ ~At.

W.K. Liu et al., Adaptive ALE finite elements

21)8

15, C o n v e r g e n c e check, if c o n v e r g e d go to step 2. 16. E n d of iteration loop. 17. E n d o f time step loop.

7. Applications to rolling problems 7.1. Rolling simulation-comparison between A L E and Lagrangian methods A rolling simulation is m a d e to e x a m i n e the p e r f o r m a n c e of the A L E f o r m u l a t i o n . T h e d i m e n s i o n s of the rigid rolls and the workpiece are described in Fig. 1. D u e to symmetry~ only half of the workpiece is mode!e,~ by 600 IPS2 e l e m e n t s [19]. T h e rolling process is s e p a r a t e d by two stages. D u r i n g the f i - , ;tage, the rigid roll is i n d e n t e d into the w o r k p i e c e with indentation speed % = 1 and . . . . ~iing speed w = 0.5 f r o m t = 0.0 t h r o u g h t = 0.8. A t the e n d

r L = 10.0 in H = 2.0 in R = 3.0 in ~' = 1.0 * f(t) in/sec

L ea = 0.5 in/sec p = 2.5,105. g(t) psi E = 1.0.108 psi ET= 1.0, 106 psi

10.~_

l.Ot/'

#

08

-~ Cr = 0.2 v= 0.3 "l;y= 3.0* 106psi isotropic hardenin~

[

l,O

t

'2 l.OT( '~ 0.6

'~ 0.4 02 7ig. I. Prob!2m

~

i

t

6

8

i0 ~"

.~em,mtfor the example o~ Section 7.1: rolling simulat r~

W.K. Liu et al.. Adaptive A L E finite elements

209

of this stage, the rolls are indented into the workpiece by 0.8 in, which is 40% of the original depth. In the second stage, the indentation speed vt, is decreased monotonically to zero. The final depth of the indentation is 0.99 in, and the remaining depth of the workpiece is 0.495% of the original depth. The indentation depth d t' versus time t is shown in Fig. 1. The workpiece is assumed to be plane strain in the z direction and the material properties are shown in Fig. 1. Following the Oden and Martins friction model, the normal pressure and the frictional traction are evaluated by

tx= - a A P n

and

tr = -~Aqd#

.

The data used in [21] are

p=q

and

/3=c,a,

and this is equivalent to the classic model with friction coefficient c r. In our computation, we take c~ = 0.2. The deformed configurati,ans of the Lagrangian analysis are shown in Figs. 2 and 3. The time incremem At is taken to be 0.02 and the computation is terminated at t = 10.0. Since the material yields under the contact region, the meshes are totally distorted on the top layer of

LAGRANGIAN

t=0.2

LAGRANGIAN

t=2.0

iiiiiiiiiii!iiiiiiiiiiii!iiiiiiiiii iiii: t=0.4

t -- 4.0,

iiiiiiiiiiiiiiii:: t=0.8

iiiili:ii iiiiiiiii ili!!iiiiiii

~=6.0

"

t = 1.2

t = 8.0

l = 1.6

t = 10.0

,

iiiiii::iiii:i::iiii:iiiii!!!!!!Z! ! Fig. 2. D e f o r m e d c o n f i g u r a t i o n o f L a g r a n g i a n c a l c u l a tion for N S T E P = 1 0 - 8 0 .

Fig. 3. D e f o r m e d c o n f i g u r a t i o n o f L a g r a n g i a n c a l c u l a tion for N S T E P = 1 0 0 - 5 0 0 .

21tl

W.K. Liu et al., Adaptive A L E finite elements

the workpiece. The final axial length of the workpiece which is measured along the center line is increased by 67.7%. In A L E analysis, the meshes inside the A L E region (the darker lines) are arranged by solving the Laplace equations (119). The time increment is taken to be the same as in the Lagrangian analysis and the computation is also terminated at t = 10.0. The deformed configurations of the A L E analysis (Figs. 4 and 5) show a smoother free surface than the Lagrangian analysis. The final axial length of the workpieces is increased by 70.3%. The effective stresses along the center line of the workpiece are reported when the left-hand side of the workpiece is 3 inches away from the center of the roll in the x-direction and the comparisons are shown in Fig. 6. Both calculations induce stress oscillations. This is probably due to the stick-slip condition in the contact region. In A L E analysis, the effective stress has a significant decrease around x = 0 which is the position where material leaves the interface. This result is more reasonable compared to the Lagrangian analysis where the significant decrease of the residual stress is around x = 1.5. The yield radius distribution of the ~'wo contributions agree well with each other and the comparisons are shown in Fig. 7. In Fig. 8, the effective stress distributions along several cross sections are represented. ~

7.2. Plate rolling For a more realistic geometry of rolling of a metal workpiece with a large roller as shown in Fig. 9 (plane-strain is used), the Lagrangian method failed to converge to a solution after many trials. Therefore, only A L E solutions are presented. ALE

t = 0.2

,GLE

"x x

t = 2.0

/ /

0.4

t=4.0

t=0.8

t=6.0

I =

K

t = 1.2

I=8.0 "x

I = 1.6

/

t = 10.0 •

-8 -7

/

-6 -5 -4

-3 -2 -I

fl

1

2

3

4

5

6

7

8

Fig. 4. D e f o r m e d configuration of A L E calculation for N S T E P = 10-80.

r

-~-~ -;-;-~-'3-~-; ~ i ~ 3' ~ ~ i ~ ~ d ~;~'~6 6fi Fig. 5. D e f o r m e d configuration of A L E c a l c ul a t i on for N S T E P = 100-500.

W.K. Liu et al.. Adaptive A L E finite elements"

211

effective stress along the center line 3, OOe+6

3.(~e+6

I

•~

I

i f I i s

2.00c +6

.ca

,

I

AU.:'

I

1.00¢+6

/

! I I

. /\

(),(IllC +(1

-3

-'~

-I

0

I

t

v

i

,

,

i

,

i

i

,

,

I

2

3

4

5

6

7

~

9

10

11

12

13

14

x-eoord. Fig. 6. Resiclual effective stresses comparison between Lagrangian and ALE calculations for the example o f Section 7.1.

The half thickness of the workpieces is H = 1.0, the roller radius is R = 400. The smallest half gap height between the two rollers (symmetric boundary condition has becp :~'pplied in Fig. 9) is set to 0.5 H. Therefore, the workpiece is to be reduced by 50% in ~,hid:ncss. The length of the workpiece is chosen to be L = 10. The material properties of the pla~e are also depicted in Fig. 9. The rolling process is done by applying an initial force P (see Fig. 9) and setting a constant angular velocity W = 0.00173. Since the rolling I:rocess is symmetric between the two rollers, only the upper half of the workpiece is modeled by a 10 x 60 IPS2 elements. For presentation of the large deformation process due to rolling, the y coordinates (i.e., in the yield radius along the center line Seth "¢+6 6¢+6 5c ~6 4 e+6"

2¢+6 ......

ALE

i

LAG

[

~e*6 '

t (]e+O

-3

,

,

-2

-1

0

1

2

3

4

5

6

7

8

0

1(}

11

12

13

14

x-coord.

Fig. 7. Yield radius comparisons between Lagrangian and ALE calculation for the example of Section 7 . 1 .

W.K. Lit+ et al., Adaptive A L E finite elements

212

effective stress along cross section

4.00c+6

3.00e+6

"~ 2.llOe+6 -

~

A B line C line D line

1.00e+6 "

line

().00c+0

+1.8

-1.6

-t.4

-l.2

-1.0

-0.8

-0.6

-0.4

4) 2

{>.0

y-coord.

°+°t ii.

J J,

-1.0

--

-2.0

line A

line B

line C

line D

Fig. 8. Effective stress along the cross section for the example of Section '7.1.

plate thickness direction) and the y and x coordinates are multiplied by a factor of 1.9 and 0.5, respectively, in Fig. 9 as welt as the subsequent plots (Figs. 10-12). The workpiece is dragged underneath the roller by a force P, up to a time t = 0.0114. This initial push of the workpiece is necessary since there is not enough frictional force induced by the rolling process prior to t = 0.0114. Because of the severe mesh distortion, adaptive ALE finite elements as described in the previous sections are required to roll the workpiece out. The variable time step At, adaptive

a

~"'"--(-~'~O

(_) C) (~

t)

()

~1o.o

~

•

H--I.O

"

~-1.~ I:0.101~

6.000

.ooo .ooo

1.a4o

T

(x

1

10--:)

F i g . 9. Problem statement for the e~ample of Section 7.2: plate rolling.

W.K. Liu et al., Adaptive A L E finite elements

213 ~:~.~

J

........................

f....

i-

[

~

m

~

TM ..

2 ...........................

~

J

i

r~49~

I

Fig. 10. Deformed configuration of ALE calculation

at selected time.

Fig. 11. Streamline of the deformed workpiece at selected time.

A L E mesh generation, and the average number of iterations for convergence (a tolerance on the normalized square root displacement norm of 0.0001 is employed), are tabulated in Table 1. The adaptive A L E mesh movement at selected times of the rolling process is presented in Fig. 10. The small circle attached to the roller represents the movement of the roller. As can be seen, the workpiece is r,olled by friction only. Steady-state deformation is roughly started at Table 1 Computational parameters Time

ALE generation

T

every m steps

0-0.043 0.043 - 5.62 5.62-11.98 11,98-29.32 29. 32- 49.5 4

10 10 5 5 1

At 1.43 × 1.43 x 1.43 × 7.17 × 0.143

Average number of iterations 10 ~ 10 -2 10 2 i0 --~

2.0 6.4 7.6 3,3 4.0

W.K. Liu et al., Adaptive A L E finite elements

214

L__~_J

w

.........................

11:~21

i

~::

T1~--~t,98

. . . .

L _.

%,iI

'I]~:29.32

ST~ 20: 5~4~1~7

Fig. 12. Effective stresses of the d e f o r m e d workpiece at selected time.

Fig. 13. Values of selected effective stresses.

t = 24.98 and continues up to approximately T = 43.77. At the beginning, the roller moves at the same speed of the workpiece. Once the workpiece is underneath the roller, the workpiece moves faster thon the roller. Nonsteady-state deformation then follows and continues to the end of the process, T = 49.54. The shapes of the lead and trail ends look opposite (see the A L E meshes), this implies that the velocity profile across the thickness changes during the rolling process. This can be seen in the streamline plots depicted in Fig. 11. From the plots, it is also shown that the velocity at the contact area is greater than that of the middle of the workpiece (symmetry line) until time T = 32.21, and then it gradually reverses. The speed across the thickness may vary as much as twenty-five percent ( I m a x i m u m l - Iminimum])/ [maximum[). The effective stress contours are shown in Fig. 12~ high stress is observed at the top-front of the workpiece. The values of selected effective stresses are shown in Fig. 13. It is also seen that the effective stress is highest at the rolled surface due to friction, and that stress decreases non-uniformly toward the middle of the workpiece. From the same figure, it is noted that the steady and uniform stress distribution during the steady-state deformation under the roller, T = 24.98 to roughly T = 43.77.

W.K. Liu et al., Adaptive A L E finite element~"

215

8. Conclusions Three aspects of ALE formulations have been described: the development of ALE mesh procedures in which the ALE mesh is deployed dynamically in subdomains to concentrate elements and nodes where needed; the development of efficient and accurate ALE stress and frictional traction update procedures; and the development of rapidly convergent, computationally-efficient ALE equations of motion in adaptive ALE meshes in material forming simulation. The application of adaptive ALE finite elements to metal rolling simulation has been presented. It is shown that for small rollers, both Lagrangian and ALE meshes are feasible and the agreement between the two is quite good. For realistic roller sizes, we were not able to complete a simulation with a Lagrangian mesh. However, the ALE mesh performed quite well, although some tricks are needed to initiate the passage of the workpiece.

Appendix A The term r,,,,,B,,,,,i can be simplified as follows: Rewrite Principle of Virtual Work: . , gui.jrii d R , = c r

(A.1)

"~' gui(r,,,,nmn,,n ~ ) d F , .

By using integration of parts on (A.1). we have the strong form equations: ~,., = 0

an R , ,

riin / = Crr,,,,,n,,tt,,tlT~

(A.2)

o n rl~ t .

Therefore, 7,,,,,b,,,,,, can be written as Tr

, T

Z,,,,,1~,,,,,i = c r r k l n k n / n , , l n , , , n i

/ . T\,~

(A.3)

+ n,,,,d, ) 1 -

Note that T m= Flmll

(A.4)

0

and substitute (63) into (A.3) to yield T [ 15~ T, nnOmn ' =

CrTkinknlllml, T 0

--Cr ~

Tl'lnknlrli

b,,, T

•

__

~T P

T

)

T

n m n, = - q

r)T T

T

T

rk'nkn'n"n'"n'

T

(A.5)

References [1] W F. Noh. CEL: A time-dependent two-space-dimensional coup!ed Eulerian-Lagrangian code. in: B. Adler. S. Fernbach and M. Rotenberg. eds.. Methods in Computational Physics. Vol. 3 (Academic Press. New York. 1964).

216

W.K. Liu et al., Adaptive A L E finite elements

[2] J.G. Trulio, Theory and structure of the AFTON codes, Air Force Weapons Laboratory, AFWL-TR-66-19, 1966. [3] C.W. Hirt, A.A. Amsden and J.L. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 227-253. [4] W.E. Pracht, Calculating three-dimensional fluid flows at all flow speeds with an Eulerian-Lagrangian computing mesh, J. Comput. Phys. 17 (1975) 132-159. [5] L.R. Stein, R.A. Gentry and C.W. Hirt, Computational Simulation of tran~dent blast loading on threedimensional structures, Comput. Methods Appl. Mech. Engrg. 11 (1977) 57-74. [6] J. Donea, P. Fasoli-Stella and S. Guiliani, Lagrangian and Eulerian finite element techniques for tlansient fluid-structure interaction problems, Paper BI/2, Trans. 4th SMiRT Conf., San Francisco, 15-19 August 1977. [7] T. Belytschko and J.M. Kennedy, Computer methods for subassembly simulation, Nuclear Engrg. Des. 49 (19781 17-38. [8] T.J.R. Hughes, W.K, Liu and T.K, Zimmermann, Lagrangian-Eulerian finite element formulation for incompressible viscous flows, U.S,-Japan Seminar on Interdisciplinary Finite Element Analysis, Cornetl University, Ithaca, NY, 7-11 August 1978; also in Comput. Methods Appl. Mech. Engrg. 29 (1981) 329-349. [9] T. Belytschko, D.P. Flanagan and J.M. Kennedy. Finite element methods with user-controlled meshes for fluid-structure interaction, Comput. Methods Appl. Mech. Engrg. 33 (1982) 669-688. [101 I". Belytschko, An overview of semidiscretiza~'ion and time intcgration procedures, in: T. Belytschko and T,J.R. Hughes, eds., Computational Methods for Transient Analysis (North-Holland, Amsterdam, 1983) t --63.

[11] W.K. Liu and J. Gvildys, Fluid-structure interact?on of tanks with an eccentric core barrel, Comput. Methods Appl. Mech. Engrg. 58 (1986) 51--77. [12] A. Huerta and W.K. Liu, Viscous flow with large free surface motion, Comput. Methods Apph Mech. Engrg. 69 (1988) 277-324. [13] R.B. Haber and J.F. Abel, Contact-slip analysis using mixed displacements, ASCE J. Engrg. Mech. 109 (1983) 411-429. [14] R.B. Haber and H.M. Koh, Explicit ext,ressions ,'or energy release rates using virtual crack extensions, Internat. J. Numer. Methods Engrg. 21 (~985) 3111-315. [15] K.A. Derbalian, E.H. Lee, R.L. Mallett 'rod R.M. tClcMeeking, Finite e,ement metal forming analysis with spacially fixed mesh, in: H. Armen and R F. Jones, Jl., eds., Application,, of Numerical Methods to Forming Process, AMD-Vol. 28, ASME Winter A,nnual Meeting, San Francisco, CA, 10-15 December 1978) 39-47. [16] J.T. Oden, Finite Elements of Nonlineal' Continua (McGraw-Hill, New York, 1972). [17] D.J. Benson, An efficient, accurate, simple ALE meth~d for nonlinear finite element programs, Department of Applied Mechanics, University of California, San D!cgo, La Jolla, C::lifornia. preprint, 1987. [18] W.K. Liu, T. Belytschko and H. Chang, An arbitra,'y Lagrangian--Eulerian finite element method for path-dependent materials, Comput. Methods Appl. Meck~. Engrg. 58 (1986) 227-245. [19] W.K. Liu, H. Chang, J.S. Chen and T. Belytschko, ALE Petrov-Galerkin finite elements for nonlinear continua, Comput. Methods Apph Mech. Engrg. 68 (1988~. 259-310. [20] ~.K. Liu, E.S. Law, D. Lam and T. Belytschko, Resul,-~ant-stress degenerated-shell element, Comput. Methods Appl. Mech. Engrg. 55 11986) 259-300. [21] J.T. Oden and J.A.C. Martins, Models and computational methods for dynamic friction phenomena, Comput. Methods Appl. Mech. Engrg. 52 (19841 527-634. [22] L.E. Malvern, Introduction to the Mechanics of a Continuous Media (Prentice-Hall, Englewood Cliffs. NJ, 1969). [23] T.J.R. Hughes and T.E. Tezduyar, Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible Euler equations, Comput. Methods Api;,I. Mech. Engrg. 45 (1984) 217-284. [24] W.D. Murray and F. Landis, Numerical and machine solutions of transient heat-conduction problems involving melting or freezing, Part I, ASME J. Heat Transfer (1959) 106--112. [25] A.M. Winslow, Equipotential zoning of two-dimensional meshes, University of California, Lawrence Radiation Laboratory Report, UCRL-7312, 1963. [26] H.E Wang and H.S. Lee, in: C. Tucker, ed., Fundamentals of Computer Modeling for Polymer Processing (Hanser, Publications. Munich, 1989) Chapter 8.