- Email: [email protected]
A study on residual stresses in rolling
Int. J. Mech. Tools Manuf~:t. Vol. 37, No. 6, pp. 837-853. 1997 © 1997 F.~tevie~"Scieace Ltd All rights reserved. Printed in Glut Bfitlfin 0890--6955~7517.00 + .00
Pergamon
PH: S0890-6955(96)00052-1
A STUDY ON RESIDUAL STRESSES IN ROLLING U. S. DIXITt and P. M. DIXITt:~ (Received 15 November 1995)
Abstract--Determination of residual stresses in rolled material is important for the design of product and process. Different approaches for determining them and difficulties associated with the methods are discussed. A simplified approach to find the longitudinal residual stress (stress in the direction of rolling) is proposed. A parametric study has been carried out to show the influence of process parameters on the residual stress distribution pattern. © 1997 Elsevier Science Ltd.
NOMENCLATURE A b £
E
E"
f
F F~ FP G
h~,h2 IDI2
1 L lp If
N n
nj P r
R"
s_3u So Su t ti
t. tn,ts U,,U2 Uc t' i Vn, Its
v. Wu wp W I . Wvl, Wv2, W$1D WS22, ]4'$12
x~ a a
area of the domain material constant proportionality constant in Equation (19) fourth-order elasticity tensor Young's modulus logarithmic elastic strain coefficient of friction deformation gradient material derivative of deformation gradient elastic part of the deformation gradient plastic part of the deformation gradient shear modulus inlet and exit thicknesses of strip unit vectors along the xm and x2 directions length of a fiber velocity gradient plastic part of 1 common lenl~th of fibers direction of S and material constant component of unit normal vector hydrostatic stress (or pressure) reduction ratio rotation tensor from polar decomposition of F ¢ deviatoric part of stress tensor deviatoric part of objective stress tensor Jaumann rate of S o time component of the traction vector component of the traction rate vector corresponding to 7"0 material derivative of first Piola--Kirchoff stress tensor interracial normal and shear stress components inlet and exit velocities of the strip right stretch tensor from polar decomposition of 1~ component of the velocity vector normal and tangential velocities on interface roll velocity spin tensor plastic part of spin tensor weight functions Cartesian coordinate angle of contact (shown in Fig. I)
a parameter defined by Equation (4)
t ~ of Mechanical Engineering, Indian Institute of Technology, Kanpar 208 016, India ~:Author to whom correspondence should be addressed. 837
838 At
~p
~ij, ~p /.z v
cry
(cr,)o GrU
O'ij, (lr
0
U.S. Dixit and P. M. Dixit time increment equivalent strain infinitesimal strain tensor equivalent strain rate (modulus of ~) modulus of ~' strain rate tensor elastic rate of deformation tensor plastic rate of deformation tensor Levy-Mises coefficient Poisson's ratio flow stress in tension yield stress in tension modulus of Jaumann rate of trij Cauchy stress tensor objective stress tensor 1. INTRODUCTION
Rolling is one of the oldest metal working processes. In the rolling process, a material is subjected to large plastic strains. When the objective is to find roll torque, roll force and normal and shear stress distributions at the roll-strip contact, often rigid-plastic analysis is used, in which elastic strains, being small, are neglected. Elastic effects must be taken into account when the objective is to find plastic boundaries and residual stresses. Residual stresses are the stresses which persist in a material under uniform temperature in the absence of externally applied loads. They are created whenever a member is permanently deformed in a non-uniform manner. The residual stresses are a manifestation of the elastic deformations that the infinitesimal elements of a member must experience in order for them to remain compatible. It is generally believed that a bad shape or poor flatness in the rolled sheet is caused by non-uniform deformation and residual stresses. Additionally, residual stresses have significant effects on further processing of sheet forming. Thus, the distribution of residual stresses across the cross-section of cold rolled sheet plays an important role in rolling technology. Experimental methods for determination of residual stresses have been developed [1, 2]. However, when the process or product is to be optimized, experimental methods do not provide an economical solution. There is, then, a need to develop analytical or numerical methods for elasto-plastic analysis of metal forming processes. Updated Lagrangian formulations have been used and are probably the most common formulations used for elasto-plastic analysis, since the incremental nature of plasticity can be conveniently accommodated in such a formulation [3, 4]. The published results of Updated Lagrangian formulations on rolling differ from experimental data, however; the boundary conditions on the free surface of the strip where the normal and shear stresses must vanish are not fully satisfied. Further, detailed results on residual stress distribution have not been published. Eulerian elasto-viscoplastic formulations have been reported by many authors. Dawson and Thompson [5] have developed "initial stress rate" method. In this approach, the stress and velocity fields associated with the convective terms are assumed known and placed on the right-hand side of the finite element equations for momentum. An iterative procedure is then used during which these terms are corrected until there is no discrepancy between the new stress distribution and the previous one. Variations are performed on pressure and velocity only. Results on residual stress distribution or plastic boundaries have not been reported. Abo-Elkhier et al. [6] also have proposed an Eulerian formulation for solving elasto-plastic problems. For a linearly strain-hardening material, analysis for extrusion process has been conducted and some results on residual stresses are presented. Shimazaki and Thompson [7] used the mixed formulation to solve for velocity, pressure and stress simultaneously. In their method, they de,coupled the momentum and continuity equations from the constitutive equation. The former set was used to obtain velocity and pressure with stress assumed constant. The latter set was used to obtain stress with velocity and pressure held constant. Iterations between these two sets were done until convergence. Thompson and Berman [8] employed the procedure developed by Shimazaki and Thomp-
A study on residual stresses in rolling
839
son [7] for the analysis of plane strain rolling. For a perfectly plastic material, analysis was done for cold rolling. However, the rate of convergence was slow and results presented are without getting full convergence. Effective stress contours are plotted, but individual stress components have not been presented. Thompson and Yu [9] have presented a flow formulation to solve the rate equilibrium equation, in their steady state form, for velocity and stress. The authors have used the formulation for the analysis of strip drawing for a perfectly plastic material. It is noted that the rate of convergence deteriorates as the inelastic response becomes large. Effective stress contours are plotted without getting full convergence. No result is presented for residual stress distribution. Maniatty et al. [10] have proposed an Eulerian elasto-viscoplastic formulation in which the kinematics consists of a multiplicative decomposition of the deformation gradient into elastic and plastic components. Hart's [11] constitutive model, which contains an internal scalar variable, has been used for solving rolling and extrusion problems. The system of equations is highly ill-conditioned and the algorithm lacks stability as the material becomes less rate sensitive. An improved algorithm for solving the problem has been presented by Maniatty [12] which makes use of discontinuous pressure field and allows the unknown nodal pressure to be eliminated from the formulation. Only longitudinal stress variation for a typical case is presented in these references [10, 12]. Results for other components of residual stresses and hydrostatic pressure are not presented. Malinowski and Lenard [13] have proposed an elasto-plastic finite element scheme, using the Prandtl-Reuss constitutive relations, in an Eulerian-updated Lagrangian reference frame. The problem is solved in two stages: first, the velocity field is determined by minimizing the total power of deformation and then the stresses are obtained by path line integration of the Jaumann stress rate. Path line integration, though an efficient method, gives unrealistic results owing to non-satisfaction of equilibrium condition. A less efficient approach which involves minimization of a functional under certain constraints is used. Two types of shape functions (parabolic and Hermitian) are used and the results differ significantly from each other. In the proposed material model, Poisson's ratio is assumed to be a continuous material function defined both for elastic and elastoplastic behavior. Malinowski [14] used the methodology developed by Malinowski and Lenard [13] for determining residual stress distribution in plane strain cold rolling, in which roll deformation was also taken into account. An additional functional is minimized iteratively for the roll deformation. Only limited results axe presented and computational difficulties are not discussed. In this paper, an attempt is made to analyze elasto-plastic rolling problem by different approaches. The difficulties encountered are discussed. It is found that there are many unsolved difficulties in getting a complete solution of the elasto-plastic problem of rolling. A simplified method to determine the more dominating longitudinal residual stresses is proposed, which may prove to be an economical analytical method. 2. MIXED FORMULATION
A mixed formulation to solve for velocity, pressure and deviatoric stress simultaneously, is presented here. As the bulk modulus of most metals is large, introduction of compressibility does not make much difference in results. Therefore, elastic deformation is also considered incompressible. The formulation is similar (but not exactly the same) to that presented in [7, 8]. The governing equations for the slow, steady-state and incompressible flow of metals during forming process are: Continuity:~a = 0
3p
(1)
3So
Equilibrium:- ~xi + ~jxj = 0
(2)
a 1 C°nstitutive:to = ~ So + fG Y;o
(3)
840
u.s. Dixit and P. M. Dixit
where {~ for elastic region ot =
for plastic region
(4)
l [~v, Ovjl ' ~iJ"~ 2 ~Xj + ~Xi/
(5)
o bS o Sij -~ Vk OX~k -- SikWkj-- ajkWki,
(6)
w'J=:
(7)
xj/
Here p is the pressure or hydrostatic part of the Cauchy stress tensor o"v, Sij is the deviatoric part of trU, ¢ij is the strain rate or rate of deformation tensor, ~ is the LevyMises coefficient, G is the shear modulus of elasticity, S~j is the Jaumann rate of Sij, wgj is the spin tensor, and v~ is the velocity vector. For materials yielding according to the von Mises criterion,/z is given by =
(8)
where equivalent strain rate ~ (or the modulus of ~ij) is defined by
=~
¢;j,~ij
(9)
and the flow stress Cry is approximated as try = (O'y)O(1 +~)n
(10)
Here, (try)o is the yield stress of material and b and n are material dependent coefficients determined from experiments. The equivalent strain ~ is obtained by the time integration of equivalent strain rate as given below. t
= f ~ dt
(11)
0
Strictly speaking, plastic strain rate should be used in Equations (8), (9) and (11). However, here total strain rate is used, which does not differ much from the plastic strain rate. The domain (or the control volume) is shown in Fig. 1. The conditions on the boundaries AF and DE are: vl = U1, v2 = 0 on AF; vl = U2, v2 = 0 on DE
(12)
where U1 and U2 are, respectively, the inlet and exit velocities related by U~ =/-/2(1-r)
(13)
A study on residual stresses in rolling
841
Xt,li
L
3hi/2 it= 0
h=O
I
3hi/1
. ,
.t/1 1..otnI
~
Vt = 0
tl - 0
I VI =lJ2 r
hill vt-ui I V2=0
D
a o
L
i
F
h[/l A. =Xlll I
(PIoal of simmiry) v t = 0 , tl=0 Fig. 1. The domain and the boundary conditions.
and r is the reduction ratio given by r-
h I -h 2
hi
(14)
Here h~ and h2 are, respectively, the inlet and exit thicknesses of the strip. Actually, in a rolling problem inlet and exit velocities are not specified; instead roll velocity is prescribed. In the present formulation, the roll velocity VR corresponding to the specified value of U2 is obtained as the velocity at the neutral point (a point where shear stress changes its sign). The velocity field for the prescribed VR Can then be obtained by multiplying the whole velocity field by an appropriate scaling factor so as to make the velocity of neutral point equal to the roll velocity. The neutral point is obtained from the rigid plastic analysis of Dixit and Dixit [15]. The procedure to take into account roll deformation is also described there. The boundaries AB and CD are free surfaces, therefore on these boundaries both of the stress components (ti = tren j) are zero. However, it is possible to use the condition v2 = 0 in the x2 direction on the surface except near the entry to roll gap. Therefore, in the present work, the following boundary conditions are used, tl = 0, v2 = 0 on AB and CD (except two nodes near entry to roll gap)
(15)
t~ = O, t2 = 0 on the two nodes near entry to roll gap
(16)
The boundary FE lies on the plane of symmetry. The symmetry implies the following conditions: l 1 ~--"0 , 122 ~-" 0 o n F E
(17)
On the roll-strip interface (boundary BC), the normal component of velocity must be zero. Thus, v. = 0 on BC
(18)
The second boundary condition is provided by the friction law,
I/,I = cl/.I
(19)
Here Wanheim and Bay's friction model [16] is used and the expression for c is given in [15].
842
U.S. Dixit and P. M. Dixit
Let vl, 1'2, p, $11, S22 and S~2 be the functions which satisfy all the essential boundary conditions exactly. Then v~, v2, p, S~1, $22 and $12 constitute a weak solution if the following integral equation is satisfied:
f
(20)
[ll + I2 +13 + I4 + I5 +16]dA = 0
A
where (21)
11 = (611 + 622)Wp
( ap aSl,+ aS14w axz / ~'
12 = - ~xl + ~
ap
(22)
as22 aS,qw v
13 = - ~x2 + ~
+ 3x,]
a
z
(23)
1 o )
14 = 6~1- ~ Sl,- ~ S,1 Ws,,
(24)
a 1 o) Is = 622- ~ $22- ~-~ $22 Ws22
(25)
~
1 o )
16= 612- ~ S~2- ~ S12 Ws,2
(26)
where Wp, w,,, Wv2,Ws~,, Wsz:, Ws~=are the weight functions which satisfy the homogeneous versions of the boundary conditions, and A represents the area of the domain. It is to be noted that in [5, 7, 8], the constitutive relation is substituted in the equilibrium equation to express the deviatoric stress as a function of velocity gradient. Such a formulation generates elements of matrix where /x appears in the numerator. The/z becomes large as the plastic response becomes small. It approaches infinity for regions of pure elasticity. Here the deviatoric stress term in the equilibrium equation is retained as it is. Integrating the second and third parts of Equation (20) and substituting the finite element approximations for velocity, pressure and stress into it leads to matrix equation of the form
(27) //
[0]
|
[K,d [K,s(v3lJ[,{S}J
/
/
[,{0}J
Here {P}, {V} and {S} contain the nodal values of pressure, velocity components and deviatoric stress components, respectively. The vector {F} contains the nodal forces corresponding to the specified traction on the boundary. These boundary terms arise from integration by parts of the terms 12 and/3. This formulation contains tL in the denominator only. The system of Equation (27) is highly ill-conditioned and solving them even by Householder method [17] gives an unreasonable velocity field. The decoupled system of equations is given by
[Kvp] [Kv~].Jt{V}J
l{F(vi,Sij)}J [Ks,(v3]{S] =-[K~]{V].
(28)
A study on residual stresses in rolling
843
Here all the submatrices except [Kvv] are the same as those in Equation (27). The subma-
trix [Kv,,] arises when the constitutive equation is substituted in the equilibrium equation by treating the whole domain as plastic, i.e. by taking ot = 1. The right-hand side vector {F(S;j)} now contains an additional term compared with {F} of Equation (27). This arises from the Jaumann stress rate term of the constitutive equation. A typical steel is selected for study. The material properties are as given below: E = 210 GPa, v = 0.3, (O'y)O = 324 MPa, b = 0.052, n = 0.295
(29)
E and v being Young's modulus and Poisson's ratio, respectively. The system was solved for the following set of process parameters:
R/hl =65, hi = 1 mm, r = 2 4 % , f = 0.08.
(30)
The value of/x in the inlet and exit regions becomes very high. As a result, the equations become highly ill-conditioned. Absurd values of pressures are obtained and solution diverges. Therefore, the value of/.~ has to be restricted to some maximum value in order to obtain convergence. This is found by numerical experiments, starting from some high value o f / x and going on reducing it, until a reasonable pressure field is obtained, i.e. pressure values are not drastically different from those obtained from rigid-plastic analysis. Solving the system, it is observed that stresses in the plastic zone are very close to those obtained from the rigid-plastic analysis of Dixit and Dixit [15], maximum variation being of the order of 5%. Roll torque and roll force are 5.15 kN-m/m and 3.08 MN/m, respectively, which are less than those obtained by the rigid plastic analysis by 0.2% and 2.4%, respectively. This indicates that elastic effects are small. However, stresses in the inlet and exit zones are not predicted accurately. In the inlet zone, equivalent deviatoric stress should reach the yield stress from zero. In the exit zone, the equivalent yield stress should go on decreasing up to some residual value. Such behavior is not obtained here. This may be due, firstly, to the assumption that the whole domain is plastic, made while substituting the constitutive equation in the momentum equation and, secondly, due to restricting the value of/~ in the inlet and exit zones. 3. METHODOF MULTIPLICATIVEDECOMPOSITIONOF THE DEFORMATIONGRADmNT In the mixed formulation of Section 2, the decomposition of total rate of deformation tensor (~)* into the elastic and plastic parts was considered to be additive. However, in the present formulation, it is the deformation gradient tensor F which is decomposed, and the decomposition into the elastic and plastic parts is considered to be multiplicative [18]. The plastic deformation of course is assumed to be incompressible. Thus, F = FeFp, det F v = 1
(31)
where ICe and FP are the elastic and plastic deformation gradients, respectively. As in the previous section, the elastic deformation is also taken to be incompressible. This implies det F" = det F = 1
(32)
Further, the plastic part (W p) of spin tensor is assumed to be zero, as the introduction of an arbitrary skew symmetric component of the plastic velocity gradient tensor does not affect the Cauchy stress for material response as long as the remaining constitutive equations are form invariant [19]. All these assumptions lead to the following decomposition of the rate of deformation tensor:
*In this section, tensor quantities are represented by boldface letters.
844
u.s. Dixit and P. M. Dixit = ~ + R¢&PRCr
(33)
where U is the elastic rate of deformation tensor, ~P is the plastic rate of deformation tensor and R e is the rotation tensor resulting from the polar decomposition of 1~: F e = lieu ¢
(34)
Here U ~ is called the right stretch tensor. The constitutive equation for the plastic deformation is g = 2~P
(35)
where /x is the Levy-Mises coefficient defined by Equation (8) and S is the deviatoric part of an objective stress tensor dr which is related to the Cauchy stress tensor ar by the relation dr = (det F)Rero'Re-R~ro'R e
(36)
Substitution of Equation (33) and the deviatoric part of Equation (36) into Equation (35) leads to the following relationship between S and ~P: S = 2~,(&-~ e)
(37)
which is the same as appearing in Equation (3). The elastic constitutive equation, however, is not expressed in the rate form using the Jaumann rate as done in Equation (3). Instead, the logarithmic strain E ~ = In U e
(38)
and the objective deviatoric stress tensor S are used. Thus, = 2GE e
(39)
Assuming the whole domain to be plastic, substituting the constitutive Equation (37) in equilibrium Equation (2) and applying the finite element procedure to continuity Equation (1) and the equilibrium equations, we get the following systems:
[K~p] [K~d/ t {V}J =
(40)
where the submatrices [Kpv], [Kv,,] and [Kv~] are the same as in Equation (28). However, the right-hand side vector {F(~")} now has a different expression from {F(vi, Su)} of Equation (28). The vector contains the nodal forces corresponding to the boundary terms arising from integration by parts of the equilibrium equation, and an additional term associated with 21,~e of Equation (37). The system Equation (40) is solved iteratively by the Householder method to find the pressure and velocity fields. In each iteration, the elastic rate of deformation tensor ~* has to be determined to compute the right-hand side vector. It is determined by integrating the equation ~'=LF
(41)
along the stream-lines where L is the velocity gradient tensor. The time integration procedure proposed by Weber and Anand [20], which involves decomposition of F into F" and In', is used here. Let Foe and F8 denote the values of 1~ and F p at the initial time.
A study on residual stresses in rolling
845
Starting with FI] = 1 and a suitable choice for ~0 in the (n+l)th time increment, the method involves the following steps. (The quantities at (n+l)th time increment are denoted by subscript (n+l). The trial quantities are denoted by an asterisk.) 1. Compute a trial elastic deformation gradient: Fne*+ 1 "~"F~ + 1FPn- 1
(42)
2. Compute the corresponding trial (logarithmic) elastic strain: Ee + l = In U,~*+ l
(43)
3. Compute the trial (deviatoric objective) stress: S ; +1 =
2GE~*+,
(44)
4. Compute the modulus and direction of the trial stress:
~'~,+ 1 =
-*
-*
*
S. + 1.S. + l, N. + 1 =
Sn+ 1
:z--
fin+ 1
(45)
It can be shown that the actual direction N,+~ is the same as NT, + 5. Check the region in which the point lies: If ~,~+ 1<(cry)o, the point lies in the elastic region and Fen ÷ l ~ - - F n + l , E• en + l ~ ~. +1
(46)
If #* + l > (%)0, the point lies in the plastic region. 6. If the point lies in the plastic region, compute the modulus ({P + 1) of 6P + l: ({Pn + t) = (~*, + t - C r y ) / ( 3 G A t ) ,
(47)
where At is the time increment. Note that modulus of a rate of deformation tensor is defined by Equation (9). Equation (47) follows from a finite difference approximation of Equation (41) and the yield condition. 7. Compute the plastic rate of deformation tensor:
EP + 1 = ~
~=p+
IN. + 1
(48)
Note that in view of Equation (35), the directions of S. ÷ ~ and ~.P+ ~ are the same. (Their moduli, however, will be related by a relation d'. + 1 = 3P~¢-P.+ 1 which is similar to Equation (8).) 8. Compute the elastic rate of the deformation tensor: er l ~ + l = ~ . + 1 - R ~ + l~P+ 1R.+
(49)
9. For the (n+2)th increment, ~ + 1 is needed. It is given by a finite difference approximation of Equation (41): FP~+ 1 = exp(At~P + 0F~,
(50)
The values of material properties and process parameters are taken from Equation (29) and Equation (30). As in Section 2, here also the value o f / z has to be restricted. Figure
846
U.S. Dixit and P. M. Dixit
2 shows the equivalent stress contours where the equivalent stress values go on increasing from zero along the rolling direction. This is in accordance with expected behavior. Stresses in the plastic zone are almost the same as those obtained from the method of Section 2 and the rigid-plastic analysis. Owing to this, the roll force and roll torque obtained from this method are the same as those obtained from the method of Section 2 and the rigid-plastic analysis. After the plastic zone, stresses should unload to residual values. However, here the equivalent stress values in the exit region seem to be unusually high. Further, it is observed that longitudinal stresses are not self-equilibrating and stresses in the thickness direction are far away from zero. Note that in the method of Section 2, stresses were correct only in the plastic zone. However, this method gives the expected stress field in the inlet region also. Thus, this method has an advantage over the earlier method. However, residual stresses are not predicted satisfactorily even by this method. Roll force and roll torque obtained by this method are the same as those obtained by the method of Section 2. 4. RATE FORMULATION
The elasto-plastic formulations described in the previous two sections give the correct stresses in plastic zone. However, as stated earlier, the numerical procedures do not converge in inlet and exit regions where the elastic response is large, unless the Levy-Mises coefficient/z is modified appropriately. On the other hand, the rate formulation proposed by Thompson and Yu [9] converges rapidly for purely elastic behavior, but its rate of convergence deteriorates as the plastic zone increases in size. Therefore, in this paper, the rate formulation is applied only to the exit zone. Assuming the second plastic boundary to be a straight line (line CG in Fig. 1), the exit zone is represented as the region CDEG in Fig. 1. In this formulation, the equilibrium equation is represented not in terms of the Cauchy stress tensor, as done in the earlier formulations, but in terms of the time rate (or material derivative) of the first Piola-Kirchoff stress tensor: ~(T ij
~V i
30 = vk axk
3xk
~V k
+
(51)
' ij
Thus, the equilibrium equation is
axj
_ o
(52)
The rate form of the constitutive equation has already been presented in Section 2. The
r =24 % , R/hi = 65 f = 0.08 /4
1.2 ,4 1.6
I.80
1.77
.3~ .2~
/I.0 t
/Ll9
/~I.09
.1~
.05~
/~1.04
.01~ .89 Fig. 2. Equivalent stress contours.
.94
"
A study on residual stresses in rolling
847
same is used here but with two modifications. Firstly, the compressibility has been introduced. Secondly, now the response is purely elastic, Then the constitutive equation takes the form
~o-vk ~
--0"~kW~j--0"j~Wk~= Cok~kl
(53)
where Cijkt = 2 G ~ik~.it + (1 --v2V~--) ~i/Skt]
(54)
is the fourth-order elasticity tensor. Here, ~'0 is the Jaumann stress rate tensor. Note that, the rate of deformation tensor ~k~and spin tensor w o have been defined earlier by Equations (5) and (7). Applying finite element procedure to Equation (52) and Equation (53), we get
[Kv(0"o)]{V}= {rv}
(55)
[K,,(v,)l{y,} = {ro~vi)}
(56)
Here, {V} and {E} contain the nodal values of velocity and stress components, respectively. The vector {Fv} contains the nodal forces corresponding to the boundary terms arising from integration by parts of the equilibrium equation. The vector {Fo(vl)} contains the terms associated with Cok~k ~. In this formulation, the natural boundary conditions are also expressed in the rate form. For the problem under consideration, the boundary conditions are as follows (see Fig. 1): The velocity components v~ and v2 prescribed on the entry surface (CG)
(57)
tl = 0 and t2 = 0 on the top free surface (CD)
(58)
~ = 0 and v2 = 0 on the line of symmetry (GE)
(59)
t~ = 0 and t2 = 0 on the exit boundary (DE)
(60)
Initial conditions need to be prescribed on the stress components 0-~, 0"22 and 0"12 on the entry surface of the control volume (boundary CG). The values of v~ and 0"o needed on surface CG can be computed either from the decoupled equations of the mixed formulation or from Maniatty's method. It was observed that they give almost identical results. The program was run for the material and process parameters given by Equations (29) and (30) except for the coefficient of friction, which was taken as 0.14. Figure 3 shows the thicknesswise variation of the longitudinal residual stress. It varies from a negative value at the center of the strip to a positive value at the top of the strip. Fibers near the center move smoothly, but friction at the roll-strip interface restrains those fibers from moving as fast. Therefore, tensile stress is observed at the top of the strip, while compressive stress is found at the center. It is observed that residual stresses in the thickness direction are not very close to zero, contrary to expectations in rolling of thin strips. This may be due to inaccuracy in the prescribed velocity component in the thickness direction (as these velocities are small), and improper modeling of the free surface profile. The program was run for a number of different process parameters. It was noted that in some cases longitudinal stresses were not self-equilibrating and, therefore, may not represent the c o r rect values. It seems that accuracy of results depends upon the accuracy of prescribed
848
U.S. Dixit and P. M. Dixit
0.4 0.3 L
r = 24%,
R/hi =65 f = 0.14"
0.2L 0.1
t A
0.0
0
I:~ -o.I
i:~ -0.2 -0.3 -0.4
-0.5 "~'"0.0 0.1 0.2 0.3 0.4 O.S 0.6 0.7' 0.8 0.9 X~l(hal2)
1.0
Fig. 3. Distribution of longitudinal residual stress across the thickness (rate formulation).
velocities and stresses at the end of plastic zone, and the shape of the second plastic boundary. 5. SIMPLIFIED APPROACH TO FINDING LONGITUDINAL RESIDUAL STRESS
In this approach, first the plastic deformation of fibers on different streamlines is estimated using Maniatty's method, described in Section 3. Then the longitudinal elastic strain retained by the fibers after unloading is determined from the compatibility requirement. Finally, the longitudinal residual stress is computed using a simplified elastic constitutive equation, Let l be the length of an infinitesimal fiber at unloading, which was along x~ direction in the undeformed configuration and whose undeformed length was unity. Let lp be the part of I which is due to plastic deformation. Then lp is given by lp = ~/FPT~ ¢,d,
(61)
where ~1 is the unit vector along the xl direction, and the plastic deformation gradient F p at the point of unloading can be determined using the method of Section 3. Note that streamlines in the rolling process are such that the fiber will be along the xrdirection even at the point of unloading. Because of the irthomogeneity of the plastic deformation, lp will be different for different streamlines. This implies that, during unloading, the elastic strains will not drop to zero values, as it would result in incompatible deformation. Therefore, every fiber must retain a part of the elastic strain which will ensure that all fibers have the same length after unloading. Let If be that common length. Then assuming the elastic strains to be small, the longitudinal elastic strain (i.e. elastic strain in the xrdirection) retained by fibers can be expressed as
I~11 =
l¢-lp Ip
(62)
A study on residual stresses in rolling
849
Now we make some simplifying assumptions. Since the top and bottom surfaces are stress free and the strip thickness is very small, we assume that tr22 (residual stress in the XE-direction) is zero. Since this is a plane-strain problem in the x~ - x2 plane, e33 is zero. Thu~, the elastic stress-strain relations become 1
kS (
(63)
(64)
Eliminating o'33 from Equation (63) and Equation (64) and substituting for e~j from Equation (62), we get the following expression for the longitudinal residual stress: E (/f--lp) trll -- (1 --V2) lp
(65)
Note that If is still to be determined. This can be done from the requirement that the net longitudinal force should be zero. Using symmetry, this condition can be stated as h2/2
f o'11 dxz=O
(66)
0
which, in view of Equation (65), leads to h2 h2/2
If=
2
I ~l
(67) dx2
0
So in this approach, first Maniatty's method is used to calculate lp at unloading for various streamlines using Equation (61). Then, Equation (67) is used to compute lf. Finally, trll on various streamlines is determined with the help of Equation (65). The above procedure was applied to calculate the longitudinal residual stress using the material properties given by Equation (29). Figure 4 shows the thicknesswise variation of the longitudinal residual stress for the same set of process parameters as used in Fig. 3. Comparison of these two figures shows that, although there is an overall qualitative agreement between the two distributions, the values and trends (especially near the center and at the top) differ. Figure 4 shows that residual stresses are higher near the top fiber. Further, the stress gradients are also large. This behavior is consistent with the condition of high friction. Figure 3 (which has been obtained by rate formulation of Section 4) does not exhibit these characteristics. Instead, it exhibits high gradients and high residual stress values at the center. Thus, the trend of residual stress distribution obtained by the rate formulation does not seem to be right. This could be due to numerical inaccuracy as stated in Section 4. Therefore, the simplified method has been used to investigate the effects of various parameters on the residual stress distribution. Figure 5 depicts the residual stress distribution for a reduced coefficient of friction. With reduced friction, deformation becomes more homogeneous, thereby reducing the values of residual stresses in the vicinity of top fibers. Figure 6 shows the residual stress distribution without hardening. Hardening does not have much influence on the pattern of residual stress distribution. Figure 7 shows residual stress distribution for a smaller
850
U.S. Dixit and P. M. Dixit
0.70 R/hi = 65
r-24%,
0.60
f =0.14
0.50
0.::5( 0.20
f
o
0.I0
" 0.00
-o.lo
lO'~ -0.20 -0.30 -0.40
- 0.50 -0.60 - 0.7(
0.1
0.2
0.3
0.4
0.5 0.6
0.7 0.8
0.9
1.0
x2/(h2/2) Fig. 4. Distribution of longitudinal residual stress across the thickness (simplified approach).
0.5t I
0.4-
r = 2 4 % , R/hi =65 f = 0.08
0.:5=0.2l
O.I-
b ~ O.Ob~
-0.I -
-0.2-
,,o..-o . , b ~ j -0.4
0.1
0.2
0.3
0.4
0.5
0.6
xz/(hz/2)
0.7
---""
0.8
0.9
Ii0
Fig. 5. Distribution of longitudinal residual stress across the thickness for low friction.
A study on residual stresses in rolling
0.50
851
r =24%,R/hi =65 f = 0.08
0.40 0.3.0 0.20 0.10
t0 x
0.00 -0.10 - 0.21:3
-0.3(
-0.40
•
0.1
0.2
0.3
0.4
0.5
0.6
(x2/hz/2)
"-"
0.7
0.8
0.9
1.0
Fig. 6. Distribution of longitudinal residual stress across the thickness (without hardening).
r = 8 % , R / h I =65
f = 0.08
0.4
a31.0.2
t
0.1
~
0.0
[:~ -0.1 -0.2
• -0.3 -0.4
-'-0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
LO
X=l(hel2)"" Fig. 7. Distribution of longitudinal residual stress across the thickness for low reduction.
852
U.S. Dixit and P. M. Dixit
reduction. Comparison of Fig. 5 and Fig. 7 shows that the pattern of residual stress distribution does not change much with reduction. Figure 8 shows the effect of R/h, on the residual stress distribution. Process parameters are the same as those in Fig. 5 except for R/h,, which is now higher. Increasing R/h, makes deformation more homogeneous. Therefore, the residual stress values are lower now compared with those in Fig. 5. 6. CONCLUSIONS
Different approaches were tried for solving the residual stress problem in plane strain cold rolling. It is observed that the mixed formulation in which pressure, velocity and deviatoric stress are taken as variables, does not yield good results owing to ill-conditioning. When the constitutive equation is decoupled from the equilibrium equation, and the whole domain is treated as plastic, the equations can be solved provided the Levy-Mises coefficient/x is modified appropriately in the inlet and exit regions. In Maniatty's method also, the whole domain is treated as plastic and/z is modified in the inlet and exit regions. Both these methods (the method of decoupled equations and Maniatty's method) give good results in the plastic region. Thus, these methods can be used to estimate the elastic effects in the plastic region. However, comparison of the stresses in the plastic region and the roll force and roll torque obtained by these methods with the corresponding rigidplastic results shows that the elastic effects are small. As far as the stress field in the exit region is concerned, both the methods give unsatisfactory results. This is because the elasto-plastic constitutive equation is used even after unloading, instead of the elastic constitutive equation. Maniatty's method, however, has one advantage over the method of decoupled equations in that it gives good results in the inlet region. Rate formulation applied to the exit region seems to give good results for the distribution of longitudinal residual stress (o-~,) if the values of stresses and velocities are known accurately at the end of the plastic zone. As far as distribution of o22 is concerned, there is an additional difficulty of modeling the free surface of the exit region. A simplified approach based on the fact that residual stresses are caused by non-uniform plastic defor0.4 /1 r - 2 4 % , 0.3
R/hi
=
130
f = 0.08
L
0.2
0.1i i
o
0.0
x
-o.i
-0.2 -0.:3
0.1
0.2
0.3
0.4
0.5
0.6
0.7' 0.8
0.9
1.0
xz / ( h = / 2 ) - - ~ Fig. 8. Distribution of longitudinal residual stress across the thickness for high R/hL.
A study on residual stresses in rolling
853
marion is proposed in this work, and a parametric study is conducted. Results seem to be reasonable qualitatively. Experimental or rigorous analysis is still needed to verify quantitative accuracy. The simplified method may prove to be an economical method to find the longitudinal residual stress for optimizing process and product. REFERENCES [1] K. Feja, V. Hauk, W.K. Krug and L. Pintshovius, Residual stress evaluation of a cold-rolled steel strip using X-rays and a layer removal technique, Mater. Sci. Engng 92, 13-21 (1987). [2] T. Wu, C.S. Hartley, X.M. Wang and C.T. Tsai, Residual stress distribution in cold rolled brass sheet, J. Mat. Proc. Tech. 45, 111-116 (1994). [3] C. Liu, P. Hartley, C.E.N. Sturgess and G.W. Rowe, Elastic-plastic finite element modeling of cold rolling of strip, Int. J. Mech. Sci. 27, 531-541 (1985). [4] C. Liu, P. Hartley, C.E.N. Sturgess and G.W. Rowe, Simulation of the cold rolling of strip using an elasticplastic finite element technique, Int. J. Mech. Sci. 27, 829-839 (1985). [5] P.R. Dawson and E.G. Thompson, Finite element analysis of steady-state elasto-visco-plastic flow by the initial stress rate method, Int. J. Numer. Methods Engng 12, 47-57 (1978). [6] M. Abo-Elkhier, G.A'E. Oravas and M.A. Dokainish, A consistent Eulerian formulation for large deformation analysis with reference to metal-extrusion process, Int. J. Non-Linear Mech. 23, 37-52 (1988). [7] Y. Shimazaki and E.G. Thompson, Elasto visco-plastic flow with special attention to boundary conditions, Int. J. Numer. Methods Engng 17, 97-112 (1981). [8] E. G. Thompson and H. M. Berman, Steady-state analysis of elasto-viscoplastic flow during rolling, Numerical Analysis of Forming Process (edited by J. F. T. Pittman et al.), 269-283. Wiley, New York (1984). [9] E.G. Thompson and S. Yu, A flow formulation for rate equilibrium equations, Int. J. Numer. Methods Engng 30, 1619-1632 (1990). [10] A.M. Maniatty, P.R. Dawson and G.G. Weber, An Eulerian elasto-viscoplastic formulation for steady-state forming processes, Int. J. Mech. Sci. 33, 361-377 (1991). [11] E.W. Hart, Constitutive relation for inelastic deformation of metals, J. Engng Mat. Tech. 98, 193-202 (1976). [12] A.M. Maniatty, Predicting residual stresses in steady-state forming processes, Computing Systems Engng 5, 171-177 (1994). [13] Z. Malinowski and J.G. Lenard, Experimental substantiation of an elastoplastic finite element scheme for flat rolling, Comput. Methods Appl. Mech. Engng 104, 1-17 (1993). [14] Z. Malinowski, Elastoplastic finite element solution to the stress problem in the plane-strain cold rolling process, Metallurgy Foundry Engng 19, 323-336 (1993). [15] U. S. Dixit and P. M. Dixit, A finite element analysis of flat rolling and application of fuzzy set theory. Int. J. Mach. Tools Manufact., 36 947-969 (1996). [16] T. Wanheim and N. Bay, A model for friction in metal forming processes, Ann. CIRP 27, 189-194 (1978). [17] K. J. Bathe, Finite Element Procedures in Engineering Analysis, Chap. 8, Sec. 8.3.2. Prentice Hall of India, New Delhi (1982). [18] E.H. Lee, Elastic-plastic deformation at finite strain, ASME J. Appl. Mech. 36, 1-6 (1969). [19] M.C. Boyce, G.G. Weber and D.M. Parks, On the kinematics of finite strain plasticity, J. Mech. Phys. Solids 37, 647-665 (1989). [20] G. Weber and L. Anand, Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids, Comput. Methods Appl. Mech. Engng 79, 173-176 (1990).
Pergamon
PH: S0890-6955(96)00052-1
A STUDY ON RESIDUAL STRESSES IN ROLLING U. S. DIXITt and P. M. DIXITt:~ (Received 15 November 1995)
Abstract--Determination of residual stresses in rolled material is important for the design of product and process. Different approaches for determining them and difficulties associated with the methods are discussed. A simplified approach to find the longitudinal residual stress (stress in the direction of rolling) is proposed. A parametric study has been carried out to show the influence of process parameters on the residual stress distribution pattern. © 1997 Elsevier Science Ltd.
NOMENCLATURE A b £
E
E"
f
F F~ FP G
h~,h2 IDI2
1 L lp If
N n
nj P r
R"
s_3u So Su t ti
t. tn,ts U,,U2 Uc t' i Vn, Its
v. Wu wp W I . Wvl, Wv2, W$1D WS22, ]4'$12
x~ a a
area of the domain material constant proportionality constant in Equation (19) fourth-order elasticity tensor Young's modulus logarithmic elastic strain coefficient of friction deformation gradient material derivative of deformation gradient elastic part of the deformation gradient plastic part of the deformation gradient shear modulus inlet and exit thicknesses of strip unit vectors along the xm and x2 directions length of a fiber velocity gradient plastic part of 1 common lenl~th of fibers direction of S and material constant component of unit normal vector hydrostatic stress (or pressure) reduction ratio rotation tensor from polar decomposition of F ¢ deviatoric part of stress tensor deviatoric part of objective stress tensor Jaumann rate of S o time component of the traction vector component of the traction rate vector corresponding to 7"0 material derivative of first Piola--Kirchoff stress tensor interracial normal and shear stress components inlet and exit velocities of the strip right stretch tensor from polar decomposition of 1~ component of the velocity vector normal and tangential velocities on interface roll velocity spin tensor plastic part of spin tensor weight functions Cartesian coordinate angle of contact (shown in Fig. I)
a parameter defined by Equation (4)
t ~ of Mechanical Engineering, Indian Institute of Technology, Kanpar 208 016, India ~:Author to whom correspondence should be addressed. 837
838 At
~p
~ij, ~p /.z v
cry
(cr,)o GrU
O'ij, (lr
0
U.S. Dixit and P. M. Dixit time increment equivalent strain infinitesimal strain tensor equivalent strain rate (modulus of ~) modulus of ~' strain rate tensor elastic rate of deformation tensor plastic rate of deformation tensor Levy-Mises coefficient Poisson's ratio flow stress in tension yield stress in tension modulus of Jaumann rate of trij Cauchy stress tensor objective stress tensor 1. INTRODUCTION
Rolling is one of the oldest metal working processes. In the rolling process, a material is subjected to large plastic strains. When the objective is to find roll torque, roll force and normal and shear stress distributions at the roll-strip contact, often rigid-plastic analysis is used, in which elastic strains, being small, are neglected. Elastic effects must be taken into account when the objective is to find plastic boundaries and residual stresses. Residual stresses are the stresses which persist in a material under uniform temperature in the absence of externally applied loads. They are created whenever a member is permanently deformed in a non-uniform manner. The residual stresses are a manifestation of the elastic deformations that the infinitesimal elements of a member must experience in order for them to remain compatible. It is generally believed that a bad shape or poor flatness in the rolled sheet is caused by non-uniform deformation and residual stresses. Additionally, residual stresses have significant effects on further processing of sheet forming. Thus, the distribution of residual stresses across the cross-section of cold rolled sheet plays an important role in rolling technology. Experimental methods for determination of residual stresses have been developed [1, 2]. However, when the process or product is to be optimized, experimental methods do not provide an economical solution. There is, then, a need to develop analytical or numerical methods for elasto-plastic analysis of metal forming processes. Updated Lagrangian formulations have been used and are probably the most common formulations used for elasto-plastic analysis, since the incremental nature of plasticity can be conveniently accommodated in such a formulation [3, 4]. The published results of Updated Lagrangian formulations on rolling differ from experimental data, however; the boundary conditions on the free surface of the strip where the normal and shear stresses must vanish are not fully satisfied. Further, detailed results on residual stress distribution have not been published. Eulerian elasto-viscoplastic formulations have been reported by many authors. Dawson and Thompson [5] have developed "initial stress rate" method. In this approach, the stress and velocity fields associated with the convective terms are assumed known and placed on the right-hand side of the finite element equations for momentum. An iterative procedure is then used during which these terms are corrected until there is no discrepancy between the new stress distribution and the previous one. Variations are performed on pressure and velocity only. Results on residual stress distribution or plastic boundaries have not been reported. Abo-Elkhier et al. [6] also have proposed an Eulerian formulation for solving elasto-plastic problems. For a linearly strain-hardening material, analysis for extrusion process has been conducted and some results on residual stresses are presented. Shimazaki and Thompson [7] used the mixed formulation to solve for velocity, pressure and stress simultaneously. In their method, they de,coupled the momentum and continuity equations from the constitutive equation. The former set was used to obtain velocity and pressure with stress assumed constant. The latter set was used to obtain stress with velocity and pressure held constant. Iterations between these two sets were done until convergence. Thompson and Berman [8] employed the procedure developed by Shimazaki and Thomp-
A study on residual stresses in rolling
839
son [7] for the analysis of plane strain rolling. For a perfectly plastic material, analysis was done for cold rolling. However, the rate of convergence was slow and results presented are without getting full convergence. Effective stress contours are plotted, but individual stress components have not been presented. Thompson and Yu [9] have presented a flow formulation to solve the rate equilibrium equation, in their steady state form, for velocity and stress. The authors have used the formulation for the analysis of strip drawing for a perfectly plastic material. It is noted that the rate of convergence deteriorates as the inelastic response becomes large. Effective stress contours are plotted without getting full convergence. No result is presented for residual stress distribution. Maniatty et al. [10] have proposed an Eulerian elasto-viscoplastic formulation in which the kinematics consists of a multiplicative decomposition of the deformation gradient into elastic and plastic components. Hart's [11] constitutive model, which contains an internal scalar variable, has been used for solving rolling and extrusion problems. The system of equations is highly ill-conditioned and the algorithm lacks stability as the material becomes less rate sensitive. An improved algorithm for solving the problem has been presented by Maniatty [12] which makes use of discontinuous pressure field and allows the unknown nodal pressure to be eliminated from the formulation. Only longitudinal stress variation for a typical case is presented in these references [10, 12]. Results for other components of residual stresses and hydrostatic pressure are not presented. Malinowski and Lenard [13] have proposed an elasto-plastic finite element scheme, using the Prandtl-Reuss constitutive relations, in an Eulerian-updated Lagrangian reference frame. The problem is solved in two stages: first, the velocity field is determined by minimizing the total power of deformation and then the stresses are obtained by path line integration of the Jaumann stress rate. Path line integration, though an efficient method, gives unrealistic results owing to non-satisfaction of equilibrium condition. A less efficient approach which involves minimization of a functional under certain constraints is used. Two types of shape functions (parabolic and Hermitian) are used and the results differ significantly from each other. In the proposed material model, Poisson's ratio is assumed to be a continuous material function defined both for elastic and elastoplastic behavior. Malinowski [14] used the methodology developed by Malinowski and Lenard [13] for determining residual stress distribution in plane strain cold rolling, in which roll deformation was also taken into account. An additional functional is minimized iteratively for the roll deformation. Only limited results axe presented and computational difficulties are not discussed. In this paper, an attempt is made to analyze elasto-plastic rolling problem by different approaches. The difficulties encountered are discussed. It is found that there are many unsolved difficulties in getting a complete solution of the elasto-plastic problem of rolling. A simplified method to determine the more dominating longitudinal residual stresses is proposed, which may prove to be an economical analytical method. 2. MIXED FORMULATION
A mixed formulation to solve for velocity, pressure and deviatoric stress simultaneously, is presented here. As the bulk modulus of most metals is large, introduction of compressibility does not make much difference in results. Therefore, elastic deformation is also considered incompressible. The formulation is similar (but not exactly the same) to that presented in [7, 8]. The governing equations for the slow, steady-state and incompressible flow of metals during forming process are: Continuity:~a = 0
3p
(1)
3So
Equilibrium:- ~xi + ~jxj = 0
(2)
a 1 C°nstitutive:to = ~ So + fG Y;o
(3)
840
u.s. Dixit and P. M. Dixit
where {~ for elastic region ot =
for plastic region
(4)
l [~v, Ovjl ' ~iJ"~ 2 ~Xj + ~Xi/
(5)
o bS o Sij -~ Vk OX~k -- SikWkj-- ajkWki,
(6)
w'J=:
(7)
xj/
Here p is the pressure or hydrostatic part of the Cauchy stress tensor o"v, Sij is the deviatoric part of trU, ¢ij is the strain rate or rate of deformation tensor, ~ is the LevyMises coefficient, G is the shear modulus of elasticity, S~j is the Jaumann rate of Sij, wgj is the spin tensor, and v~ is the velocity vector. For materials yielding according to the von Mises criterion,/z is given by =
(8)
where equivalent strain rate ~ (or the modulus of ~ij) is defined by
=~
¢;j,~ij
(9)
and the flow stress Cry is approximated as try = (O'y)O(1 +~)n
(10)
Here, (try)o is the yield stress of material and b and n are material dependent coefficients determined from experiments. The equivalent strain ~ is obtained by the time integration of equivalent strain rate as given below. t
= f ~ dt
(11)
0
Strictly speaking, plastic strain rate should be used in Equations (8), (9) and (11). However, here total strain rate is used, which does not differ much from the plastic strain rate. The domain (or the control volume) is shown in Fig. 1. The conditions on the boundaries AF and DE are: vl = U1, v2 = 0 on AF; vl = U2, v2 = 0 on DE
(12)
where U1 and U2 are, respectively, the inlet and exit velocities related by U~ =/-/2(1-r)
(13)
A study on residual stresses in rolling
841
Xt,li
L
3hi/2 it= 0
h=O
I
3hi/1
. ,
.t/1 1..otnI
~
Vt = 0
tl - 0
I VI =lJ2 r
hill vt-ui I V2=0
D
a o
L
i
F
h[/l A. =Xlll I
(PIoal of simmiry) v t = 0 , tl=0 Fig. 1. The domain and the boundary conditions.
and r is the reduction ratio given by r-
h I -h 2
hi
(14)
Here h~ and h2 are, respectively, the inlet and exit thicknesses of the strip. Actually, in a rolling problem inlet and exit velocities are not specified; instead roll velocity is prescribed. In the present formulation, the roll velocity VR corresponding to the specified value of U2 is obtained as the velocity at the neutral point (a point where shear stress changes its sign). The velocity field for the prescribed VR Can then be obtained by multiplying the whole velocity field by an appropriate scaling factor so as to make the velocity of neutral point equal to the roll velocity. The neutral point is obtained from the rigid plastic analysis of Dixit and Dixit [15]. The procedure to take into account roll deformation is also described there. The boundaries AB and CD are free surfaces, therefore on these boundaries both of the stress components (ti = tren j) are zero. However, it is possible to use the condition v2 = 0 in the x2 direction on the surface except near the entry to roll gap. Therefore, in the present work, the following boundary conditions are used, tl = 0, v2 = 0 on AB and CD (except two nodes near entry to roll gap)
(15)
t~ = O, t2 = 0 on the two nodes near entry to roll gap
(16)
The boundary FE lies on the plane of symmetry. The symmetry implies the following conditions: l 1 ~--"0 , 122 ~-" 0 o n F E
(17)
On the roll-strip interface (boundary BC), the normal component of velocity must be zero. Thus, v. = 0 on BC
(18)
The second boundary condition is provided by the friction law,
I/,I = cl/.I
(19)
Here Wanheim and Bay's friction model [16] is used and the expression for c is given in [15].
842
U.S. Dixit and P. M. Dixit
Let vl, 1'2, p, $11, S22 and S~2 be the functions which satisfy all the essential boundary conditions exactly. Then v~, v2, p, S~1, $22 and $12 constitute a weak solution if the following integral equation is satisfied:
f
(20)
[ll + I2 +13 + I4 + I5 +16]dA = 0
A
where (21)
11 = (611 + 622)Wp
( ap aSl,+ aS14w axz / ~'
12 = - ~xl + ~
ap
(22)
as22 aS,qw v
13 = - ~x2 + ~
+ 3x,]
a
z
(23)
1 o )
14 = 6~1- ~ Sl,- ~ S,1 Ws,,
(24)
a 1 o) Is = 622- ~ $22- ~-~ $22 Ws22
(25)
~
1 o )
16= 612- ~ S~2- ~ S12 Ws,2
(26)
where Wp, w,,, Wv2,Ws~,, Wsz:, Ws~=are the weight functions which satisfy the homogeneous versions of the boundary conditions, and A represents the area of the domain. It is to be noted that in [5, 7, 8], the constitutive relation is substituted in the equilibrium equation to express the deviatoric stress as a function of velocity gradient. Such a formulation generates elements of matrix where /x appears in the numerator. The/z becomes large as the plastic response becomes small. It approaches infinity for regions of pure elasticity. Here the deviatoric stress term in the equilibrium equation is retained as it is. Integrating the second and third parts of Equation (20) and substituting the finite element approximations for velocity, pressure and stress into it leads to matrix equation of the form
(27) //
[0]
|
[K,d [K,s(v3lJ[,{S}J
/
/
[,{0}J
Here {P}, {V} and {S} contain the nodal values of pressure, velocity components and deviatoric stress components, respectively. The vector {F} contains the nodal forces corresponding to the specified traction on the boundary. These boundary terms arise from integration by parts of the terms 12 and/3. This formulation contains tL in the denominator only. The system of Equation (27) is highly ill-conditioned and solving them even by Householder method [17] gives an unreasonable velocity field. The decoupled system of equations is given by
[Kvp] [Kv~].Jt{V}J
l{F(vi,Sij)}J [Ks,(v3]{S] =-[K~]{V].
(28)
A study on residual stresses in rolling
843
Here all the submatrices except [Kvv] are the same as those in Equation (27). The subma-
trix [Kv,,] arises when the constitutive equation is substituted in the equilibrium equation by treating the whole domain as plastic, i.e. by taking ot = 1. The right-hand side vector {F(S;j)} now contains an additional term compared with {F} of Equation (27). This arises from the Jaumann stress rate term of the constitutive equation. A typical steel is selected for study. The material properties are as given below: E = 210 GPa, v = 0.3, (O'y)O = 324 MPa, b = 0.052, n = 0.295
(29)
E and v being Young's modulus and Poisson's ratio, respectively. The system was solved for the following set of process parameters:
R/hl =65, hi = 1 mm, r = 2 4 % , f = 0.08.
(30)
The value of/x in the inlet and exit regions becomes very high. As a result, the equations become highly ill-conditioned. Absurd values of pressures are obtained and solution diverges. Therefore, the value of/.~ has to be restricted to some maximum value in order to obtain convergence. This is found by numerical experiments, starting from some high value o f / x and going on reducing it, until a reasonable pressure field is obtained, i.e. pressure values are not drastically different from those obtained from rigid-plastic analysis. Solving the system, it is observed that stresses in the plastic zone are very close to those obtained from the rigid-plastic analysis of Dixit and Dixit [15], maximum variation being of the order of 5%. Roll torque and roll force are 5.15 kN-m/m and 3.08 MN/m, respectively, which are less than those obtained by the rigid plastic analysis by 0.2% and 2.4%, respectively. This indicates that elastic effects are small. However, stresses in the inlet and exit zones are not predicted accurately. In the inlet zone, equivalent deviatoric stress should reach the yield stress from zero. In the exit zone, the equivalent yield stress should go on decreasing up to some residual value. Such behavior is not obtained here. This may be due, firstly, to the assumption that the whole domain is plastic, made while substituting the constitutive equation in the momentum equation and, secondly, due to restricting the value of/~ in the inlet and exit zones. 3. METHODOF MULTIPLICATIVEDECOMPOSITIONOF THE DEFORMATIONGRADmNT In the mixed formulation of Section 2, the decomposition of total rate of deformation tensor (~)* into the elastic and plastic parts was considered to be additive. However, in the present formulation, it is the deformation gradient tensor F which is decomposed, and the decomposition into the elastic and plastic parts is considered to be multiplicative [18]. The plastic deformation of course is assumed to be incompressible. Thus, F = FeFp, det F v = 1
(31)
where ICe and FP are the elastic and plastic deformation gradients, respectively. As in the previous section, the elastic deformation is also taken to be incompressible. This implies det F" = det F = 1
(32)
Further, the plastic part (W p) of spin tensor is assumed to be zero, as the introduction of an arbitrary skew symmetric component of the plastic velocity gradient tensor does not affect the Cauchy stress for material response as long as the remaining constitutive equations are form invariant [19]. All these assumptions lead to the following decomposition of the rate of deformation tensor:
*In this section, tensor quantities are represented by boldface letters.
844
u.s. Dixit and P. M. Dixit = ~ + R¢&PRCr
(33)
where U is the elastic rate of deformation tensor, ~P is the plastic rate of deformation tensor and R e is the rotation tensor resulting from the polar decomposition of 1~: F e = lieu ¢
(34)
Here U ~ is called the right stretch tensor. The constitutive equation for the plastic deformation is g = 2~P
(35)
where /x is the Levy-Mises coefficient defined by Equation (8) and S is the deviatoric part of an objective stress tensor dr which is related to the Cauchy stress tensor ar by the relation dr = (det F)Rero'Re-R~ro'R e
(36)
Substitution of Equation (33) and the deviatoric part of Equation (36) into Equation (35) leads to the following relationship between S and ~P: S = 2~,(&-~ e)
(37)
which is the same as appearing in Equation (3). The elastic constitutive equation, however, is not expressed in the rate form using the Jaumann rate as done in Equation (3). Instead, the logarithmic strain E ~ = In U e
(38)
and the objective deviatoric stress tensor S are used. Thus, = 2GE e
(39)
Assuming the whole domain to be plastic, substituting the constitutive Equation (37) in equilibrium Equation (2) and applying the finite element procedure to continuity Equation (1) and the equilibrium equations, we get the following systems:
[K~p] [K~d/ t {V}J =
(40)
where the submatrices [Kpv], [Kv,,] and [Kv~] are the same as in Equation (28). However, the right-hand side vector {F(~")} now has a different expression from {F(vi, Su)} of Equation (28). The vector contains the nodal forces corresponding to the boundary terms arising from integration by parts of the equilibrium equation, and an additional term associated with 21,~e of Equation (37). The system Equation (40) is solved iteratively by the Householder method to find the pressure and velocity fields. In each iteration, the elastic rate of deformation tensor ~* has to be determined to compute the right-hand side vector. It is determined by integrating the equation ~'=LF
(41)
along the stream-lines where L is the velocity gradient tensor. The time integration procedure proposed by Weber and Anand [20], which involves decomposition of F into F" and In', is used here. Let Foe and F8 denote the values of 1~ and F p at the initial time.
A study on residual stresses in rolling
845
Starting with FI] = 1 and a suitable choice for ~0 in the (n+l)th time increment, the method involves the following steps. (The quantities at (n+l)th time increment are denoted by subscript (n+l). The trial quantities are denoted by an asterisk.) 1. Compute a trial elastic deformation gradient: Fne*+ 1 "~"F~ + 1FPn- 1
(42)
2. Compute the corresponding trial (logarithmic) elastic strain: Ee + l = In U,~*+ l
(43)
3. Compute the trial (deviatoric objective) stress: S ; +1 =
2GE~*+,
(44)
4. Compute the modulus and direction of the trial stress:
~'~,+ 1 =
-*
-*
*
S. + 1.S. + l, N. + 1 =
Sn+ 1
:z--
fin+ 1
(45)
It can be shown that the actual direction N,+~ is the same as NT, + 5. Check the region in which the point lies: If ~,~+ 1<(cry)o, the point lies in the elastic region and Fen ÷ l ~ - - F n + l , E• en + l ~ ~. +1
(46)
If #* + l > (%)0, the point lies in the plastic region. 6. If the point lies in the plastic region, compute the modulus ({P + 1) of 6P + l: ({Pn + t) = (~*, + t - C r y ) / ( 3 G A t ) ,
(47)
where At is the time increment. Note that modulus of a rate of deformation tensor is defined by Equation (9). Equation (47) follows from a finite difference approximation of Equation (41) and the yield condition. 7. Compute the plastic rate of deformation tensor:
EP + 1 = ~
~=p+
IN. + 1
(48)
Note that in view of Equation (35), the directions of S. ÷ ~ and ~.P+ ~ are the same. (Their moduli, however, will be related by a relation d'. + 1 = 3P~¢-P.+ 1 which is similar to Equation (8).) 8. Compute the elastic rate of the deformation tensor: er l ~ + l = ~ . + 1 - R ~ + l~P+ 1R.+
(49)
9. For the (n+2)th increment, ~ + 1 is needed. It is given by a finite difference approximation of Equation (41): FP~+ 1 = exp(At~P + 0F~,
(50)
The values of material properties and process parameters are taken from Equation (29) and Equation (30). As in Section 2, here also the value o f / z has to be restricted. Figure
846
U.S. Dixit and P. M. Dixit
2 shows the equivalent stress contours where the equivalent stress values go on increasing from zero along the rolling direction. This is in accordance with expected behavior. Stresses in the plastic zone are almost the same as those obtained from the method of Section 2 and the rigid-plastic analysis. Owing to this, the roll force and roll torque obtained from this method are the same as those obtained from the method of Section 2 and the rigid-plastic analysis. After the plastic zone, stresses should unload to residual values. However, here the equivalent stress values in the exit region seem to be unusually high. Further, it is observed that longitudinal stresses are not self-equilibrating and stresses in the thickness direction are far away from zero. Note that in the method of Section 2, stresses were correct only in the plastic zone. However, this method gives the expected stress field in the inlet region also. Thus, this method has an advantage over the earlier method. However, residual stresses are not predicted satisfactorily even by this method. Roll force and roll torque obtained by this method are the same as those obtained by the method of Section 2. 4. RATE FORMULATION
The elasto-plastic formulations described in the previous two sections give the correct stresses in plastic zone. However, as stated earlier, the numerical procedures do not converge in inlet and exit regions where the elastic response is large, unless the Levy-Mises coefficient/z is modified appropriately. On the other hand, the rate formulation proposed by Thompson and Yu [9] converges rapidly for purely elastic behavior, but its rate of convergence deteriorates as the plastic zone increases in size. Therefore, in this paper, the rate formulation is applied only to the exit zone. Assuming the second plastic boundary to be a straight line (line CG in Fig. 1), the exit zone is represented as the region CDEG in Fig. 1. In this formulation, the equilibrium equation is represented not in terms of the Cauchy stress tensor, as done in the earlier formulations, but in terms of the time rate (or material derivative) of the first Piola-Kirchoff stress tensor: ~(T ij
~V i
30 = vk axk
3xk
~V k
+
(51)
' ij
Thus, the equilibrium equation is
axj
_ o
(52)
The rate form of the constitutive equation has already been presented in Section 2. The
r =24 % , R/hi = 65 f = 0.08 /4
1.2 ,4 1.6
I.80
1.77
.3~ .2~
/I.0 t
/Ll9
/~I.09
.1~
.05~
/~1.04
.01~ .89 Fig. 2. Equivalent stress contours.
.94
"
A study on residual stresses in rolling
847
same is used here but with two modifications. Firstly, the compressibility has been introduced. Secondly, now the response is purely elastic, Then the constitutive equation takes the form
~o-vk ~
--0"~kW~j--0"j~Wk~= Cok~kl
(53)
where Cijkt = 2 G ~ik~.it + (1 --v2V~--) ~i/Skt]
(54)
is the fourth-order elasticity tensor. Here, ~'0 is the Jaumann stress rate tensor. Note that, the rate of deformation tensor ~k~and spin tensor w o have been defined earlier by Equations (5) and (7). Applying finite element procedure to Equation (52) and Equation (53), we get
[Kv(0"o)]{V}= {rv}
(55)
[K,,(v,)l{y,} = {ro~vi)}
(56)
Here, {V} and {E} contain the nodal values of velocity and stress components, respectively. The vector {Fv} contains the nodal forces corresponding to the boundary terms arising from integration by parts of the equilibrium equation. The vector {Fo(vl)} contains the terms associated with Cok~k ~. In this formulation, the natural boundary conditions are also expressed in the rate form. For the problem under consideration, the boundary conditions are as follows (see Fig. 1): The velocity components v~ and v2 prescribed on the entry surface (CG)
(57)
tl = 0 and t2 = 0 on the top free surface (CD)
(58)
~ = 0 and v2 = 0 on the line of symmetry (GE)
(59)
t~ = 0 and t2 = 0 on the exit boundary (DE)
(60)
Initial conditions need to be prescribed on the stress components 0-~, 0"22 and 0"12 on the entry surface of the control volume (boundary CG). The values of v~ and 0"o needed on surface CG can be computed either from the decoupled equations of the mixed formulation or from Maniatty's method. It was observed that they give almost identical results. The program was run for the material and process parameters given by Equations (29) and (30) except for the coefficient of friction, which was taken as 0.14. Figure 3 shows the thicknesswise variation of the longitudinal residual stress. It varies from a negative value at the center of the strip to a positive value at the top of the strip. Fibers near the center move smoothly, but friction at the roll-strip interface restrains those fibers from moving as fast. Therefore, tensile stress is observed at the top of the strip, while compressive stress is found at the center. It is observed that residual stresses in the thickness direction are not very close to zero, contrary to expectations in rolling of thin strips. This may be due to inaccuracy in the prescribed velocity component in the thickness direction (as these velocities are small), and improper modeling of the free surface profile. The program was run for a number of different process parameters. It was noted that in some cases longitudinal stresses were not self-equilibrating and, therefore, may not represent the c o r rect values. It seems that accuracy of results depends upon the accuracy of prescribed
848
U.S. Dixit and P. M. Dixit
0.4 0.3 L
r = 24%,
R/hi =65 f = 0.14"
0.2L 0.1
t A
0.0
0
I:~ -o.I
i:~ -0.2 -0.3 -0.4
-0.5 "~'"0.0 0.1 0.2 0.3 0.4 O.S 0.6 0.7' 0.8 0.9 X~l(hal2)
1.0
Fig. 3. Distribution of longitudinal residual stress across the thickness (rate formulation).
velocities and stresses at the end of plastic zone, and the shape of the second plastic boundary. 5. SIMPLIFIED APPROACH TO FINDING LONGITUDINAL RESIDUAL STRESS
In this approach, first the plastic deformation of fibers on different streamlines is estimated using Maniatty's method, described in Section 3. Then the longitudinal elastic strain retained by the fibers after unloading is determined from the compatibility requirement. Finally, the longitudinal residual stress is computed using a simplified elastic constitutive equation, Let l be the length of an infinitesimal fiber at unloading, which was along x~ direction in the undeformed configuration and whose undeformed length was unity. Let lp be the part of I which is due to plastic deformation. Then lp is given by lp = ~/FPT~ ¢,d,
(61)
where ~1 is the unit vector along the xl direction, and the plastic deformation gradient F p at the point of unloading can be determined using the method of Section 3. Note that streamlines in the rolling process are such that the fiber will be along the xrdirection even at the point of unloading. Because of the irthomogeneity of the plastic deformation, lp will be different for different streamlines. This implies that, during unloading, the elastic strains will not drop to zero values, as it would result in incompatible deformation. Therefore, every fiber must retain a part of the elastic strain which will ensure that all fibers have the same length after unloading. Let If be that common length. Then assuming the elastic strains to be small, the longitudinal elastic strain (i.e. elastic strain in the xrdirection) retained by fibers can be expressed as
I~11 =
l¢-lp Ip
(62)
A study on residual stresses in rolling
849
Now we make some simplifying assumptions. Since the top and bottom surfaces are stress free and the strip thickness is very small, we assume that tr22 (residual stress in the XE-direction) is zero. Since this is a plane-strain problem in the x~ - x2 plane, e33 is zero. Thu~, the elastic stress-strain relations become 1
kS (
(63)
(64)
Eliminating o'33 from Equation (63) and Equation (64) and substituting for e~j from Equation (62), we get the following expression for the longitudinal residual stress: E (/f--lp) trll -- (1 --V2) lp
(65)
Note that If is still to be determined. This can be done from the requirement that the net longitudinal force should be zero. Using symmetry, this condition can be stated as h2/2
f o'11 dxz=O
(66)
0
which, in view of Equation (65), leads to h2 h2/2
If=
2
I ~l
(67) dx2
0
So in this approach, first Maniatty's method is used to calculate lp at unloading for various streamlines using Equation (61). Then, Equation (67) is used to compute lf. Finally, trll on various streamlines is determined with the help of Equation (65). The above procedure was applied to calculate the longitudinal residual stress using the material properties given by Equation (29). Figure 4 shows the thicknesswise variation of the longitudinal residual stress for the same set of process parameters as used in Fig. 3. Comparison of these two figures shows that, although there is an overall qualitative agreement between the two distributions, the values and trends (especially near the center and at the top) differ. Figure 4 shows that residual stresses are higher near the top fiber. Further, the stress gradients are also large. This behavior is consistent with the condition of high friction. Figure 3 (which has been obtained by rate formulation of Section 4) does not exhibit these characteristics. Instead, it exhibits high gradients and high residual stress values at the center. Thus, the trend of residual stress distribution obtained by the rate formulation does not seem to be right. This could be due to numerical inaccuracy as stated in Section 4. Therefore, the simplified method has been used to investigate the effects of various parameters on the residual stress distribution. Figure 5 depicts the residual stress distribution for a reduced coefficient of friction. With reduced friction, deformation becomes more homogeneous, thereby reducing the values of residual stresses in the vicinity of top fibers. Figure 6 shows the residual stress distribution without hardening. Hardening does not have much influence on the pattern of residual stress distribution. Figure 7 shows residual stress distribution for a smaller
850
U.S. Dixit and P. M. Dixit
0.70 R/hi = 65
r-24%,
0.60
f =0.14
0.50
0.::5( 0.20
f
o
0.I0
" 0.00
-o.lo
lO'~ -0.20 -0.30 -0.40
- 0.50 -0.60 - 0.7(
0.1
0.2
0.3
0.4
0.5 0.6
0.7 0.8
0.9
1.0
x2/(h2/2) Fig. 4. Distribution of longitudinal residual stress across the thickness (simplified approach).
0.5t I
0.4-
r = 2 4 % , R/hi =65 f = 0.08
0.:5=0.2l
O.I-
b ~ O.Ob~
-0.I -
-0.2-
,,o..-o . , b ~ j -0.4
0.1
0.2
0.3
0.4
0.5
0.6
xz/(hz/2)
0.7
---""
0.8
0.9
Ii0
Fig. 5. Distribution of longitudinal residual stress across the thickness for low friction.
A study on residual stresses in rolling
0.50
851
r =24%,R/hi =65 f = 0.08
0.40 0.3.0 0.20 0.10
t0 x
0.00 -0.10 - 0.21:3
-0.3(
-0.40
•
0.1
0.2
0.3
0.4
0.5
0.6
(x2/hz/2)
"-"
0.7
0.8
0.9
1.0
Fig. 6. Distribution of longitudinal residual stress across the thickness (without hardening).
r = 8 % , R / h I =65
f = 0.08
0.4
a31.0.2
t
0.1
~
0.0
[:~ -0.1 -0.2
• -0.3 -0.4
-'-0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
LO
X=l(hel2)"" Fig. 7. Distribution of longitudinal residual stress across the thickness for low reduction.
852
U.S. Dixit and P. M. Dixit
reduction. Comparison of Fig. 5 and Fig. 7 shows that the pattern of residual stress distribution does not change much with reduction. Figure 8 shows the effect of R/h, on the residual stress distribution. Process parameters are the same as those in Fig. 5 except for R/h,, which is now higher. Increasing R/h, makes deformation more homogeneous. Therefore, the residual stress values are lower now compared with those in Fig. 5. 6. CONCLUSIONS
Different approaches were tried for solving the residual stress problem in plane strain cold rolling. It is observed that the mixed formulation in which pressure, velocity and deviatoric stress are taken as variables, does not yield good results owing to ill-conditioning. When the constitutive equation is decoupled from the equilibrium equation, and the whole domain is treated as plastic, the equations can be solved provided the Levy-Mises coefficient/x is modified appropriately in the inlet and exit regions. In Maniatty's method also, the whole domain is treated as plastic and/z is modified in the inlet and exit regions. Both these methods (the method of decoupled equations and Maniatty's method) give good results in the plastic region. Thus, these methods can be used to estimate the elastic effects in the plastic region. However, comparison of the stresses in the plastic region and the roll force and roll torque obtained by these methods with the corresponding rigidplastic results shows that the elastic effects are small. As far as the stress field in the exit region is concerned, both the methods give unsatisfactory results. This is because the elasto-plastic constitutive equation is used even after unloading, instead of the elastic constitutive equation. Maniatty's method, however, has one advantage over the method of decoupled equations in that it gives good results in the inlet region. Rate formulation applied to the exit region seems to give good results for the distribution of longitudinal residual stress (o-~,) if the values of stresses and velocities are known accurately at the end of the plastic zone. As far as distribution of o22 is concerned, there is an additional difficulty of modeling the free surface of the exit region. A simplified approach based on the fact that residual stresses are caused by non-uniform plastic defor0.4 /1 r - 2 4 % , 0.3
R/hi
=
130
f = 0.08
L
0.2
0.1i i
o
0.0
x
-o.i
-0.2 -0.:3
0.1
0.2
0.3
0.4
0.5
0.6
0.7' 0.8
0.9
1.0
xz / ( h = / 2 ) - - ~ Fig. 8. Distribution of longitudinal residual stress across the thickness for high R/hL.
A study on residual stresses in rolling
853
marion is proposed in this work, and a parametric study is conducted. Results seem to be reasonable qualitatively. Experimental or rigorous analysis is still needed to verify quantitative accuracy. The simplified method may prove to be an economical method to find the longitudinal residual stress for optimizing process and product. REFERENCES [1] K. Feja, V. Hauk, W.K. Krug and L. Pintshovius, Residual stress evaluation of a cold-rolled steel strip using X-rays and a layer removal technique, Mater. Sci. Engng 92, 13-21 (1987). [2] T. Wu, C.S. Hartley, X.M. Wang and C.T. Tsai, Residual stress distribution in cold rolled brass sheet, J. Mat. Proc. Tech. 45, 111-116 (1994). [3] C. Liu, P. Hartley, C.E.N. Sturgess and G.W. Rowe, Elastic-plastic finite element modeling of cold rolling of strip, Int. J. Mech. Sci. 27, 531-541 (1985). [4] C. Liu, P. Hartley, C.E.N. Sturgess and G.W. Rowe, Simulation of the cold rolling of strip using an elasticplastic finite element technique, Int. J. Mech. Sci. 27, 829-839 (1985). [5] P.R. Dawson and E.G. Thompson, Finite element analysis of steady-state elasto-visco-plastic flow by the initial stress rate method, Int. J. Numer. Methods Engng 12, 47-57 (1978). [6] M. Abo-Elkhier, G.A'E. Oravas and M.A. Dokainish, A consistent Eulerian formulation for large deformation analysis with reference to metal-extrusion process, Int. J. Non-Linear Mech. 23, 37-52 (1988). [7] Y. Shimazaki and E.G. Thompson, Elasto visco-plastic flow with special attention to boundary conditions, Int. J. Numer. Methods Engng 17, 97-112 (1981). [8] E. G. Thompson and H. M. Berman, Steady-state analysis of elasto-viscoplastic flow during rolling, Numerical Analysis of Forming Process (edited by J. F. T. Pittman et al.), 269-283. Wiley, New York (1984). [9] E.G. Thompson and S. Yu, A flow formulation for rate equilibrium equations, Int. J. Numer. Methods Engng 30, 1619-1632 (1990). [10] A.M. Maniatty, P.R. Dawson and G.G. Weber, An Eulerian elasto-viscoplastic formulation for steady-state forming processes, Int. J. Mech. Sci. 33, 361-377 (1991). [11] E.W. Hart, Constitutive relation for inelastic deformation of metals, J. Engng Mat. Tech. 98, 193-202 (1976). [12] A.M. Maniatty, Predicting residual stresses in steady-state forming processes, Computing Systems Engng 5, 171-177 (1994). [13] Z. Malinowski and J.G. Lenard, Experimental substantiation of an elastoplastic finite element scheme for flat rolling, Comput. Methods Appl. Mech. Engng 104, 1-17 (1993). [14] Z. Malinowski, Elastoplastic finite element solution to the stress problem in the plane-strain cold rolling process, Metallurgy Foundry Engng 19, 323-336 (1993). [15] U. S. Dixit and P. M. Dixit, A finite element analysis of flat rolling and application of fuzzy set theory. Int. J. Mach. Tools Manufact., 36 947-969 (1996). [16] T. Wanheim and N. Bay, A model for friction in metal forming processes, Ann. CIRP 27, 189-194 (1978). [17] K. J. Bathe, Finite Element Procedures in Engineering Analysis, Chap. 8, Sec. 8.3.2. Prentice Hall of India, New Delhi (1982). [18] E.H. Lee, Elastic-plastic deformation at finite strain, ASME J. Appl. Mech. 36, 1-6 (1969). [19] M.C. Boyce, G.G. Weber and D.M. Parks, On the kinematics of finite strain plasticity, J. Mech. Phys. Solids 37, 647-665 (1989). [20] G. Weber and L. Anand, Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids, Comput. Methods Appl. Mech. Engng 79, 173-176 (1990).