Finite element analysis of cold strip rolling

Int. J. Mach. Tool Des. Res. Vol. 17, pp. 159-168. PergamonPress 1977. Printed in Great Britain

FINITE

ELEMENT

ANALYSIS

OF

COLD

STRIP ROLLING

S. S. R A O * a n d A L O K K U M A R ~

(First received 28 February 1977; in final form 29 June 1977) A b s t r a c t - - T h e cold rolling of strips is analyzed as a plane strain problem using the finite element method. The incremental displacement method is used for the nonlinear analysis of the problem. The stress and strain history of each element can be traced as the material passes through the rolls. The convergence of the method, the flow of material under the rolls, the extent of deformation of the material and the pressure distribution on the rolls are studied in this work. The method is equally applicable for the analysis of rolling of rectangular bars where the spread cannot be neglected.

INTRODUCTION T H E HOMOGENEOUS flow theories of cold rolling due to Von Karman [1], Trinks [2], Tselikov [3], Nadai [4] and Bland and Ford [5] predict many of the parameters needed for the construction of equipment reasonably well. However, these theories are of little value in predicting the stress and strain distributions during and after the deformation due to the presence of non-homogeneity as evidenced by the experimental investigations of MacGregor and Coffin [6] and Hundy and Singer [7]. The slip-line field approach is based on the very feature which produces non-homogeneous deformation and hence is expected to be a powerful analytical tool in metal forming processes. Although accurate slip-line solution for the general problem of rolling is not available [8], Alexander has presented a slip-line field for the hot rolling process [9]. Firbank and Lancaster [10,11] proposed a slip-line solution for the cold rolling of thin strips with a pass reduction of 20~o. However the results predicted by this method differed by about 50~ in the position of neutral point from exit, 1~o in the mean roll pressure and 25~o in the roll torque compared to those predicted by Von Karman's theory [10]. The slip line solutions are thus applicable to situations where reductions are small to permit one to assume that the arc of contact is a straight line. Thus the homogeneous flow and slip-line theories can be seen to be applicable to plane strain problems only. When it comes to the analysis of rolling of rectangular bars, where the spread cannot be neglected, the theories have hardly progressed beyond the empirical stage. The finite element method was originally developed for the analysis of complex structures and is in the process of rapid extension to a wide range of nonstructural problems. This method has been applied to simple cases of metal forming like indentation [12], where studies were made on the development of plastic zone, the load displacement relations, and the stress and strain distributions. More recently the finite element method was applied to the problems of plane strain and axisymmetric extrusion [13], where the effect of frictional coefficient on the spread of plastic zone, the pressure displacement curve and the stress and strain distributions were studied for the initial nonsteady state extrusion using incremental load method of analysis. In this study the cold rolling of strips is analyzed as a plane strain problem using the finite element method. The method is equally applicable to problems where the spread cannot be neglected. The incremental displacement method is used for the nonlinear analysis of the problem. The material to be deformed in cold rolling is subdivided into triangular elements. Initially the material is assumed to be elastic but when it passes through the rolls, it becomes plastic due to deformation. Depending on whether * Assistant Professor, Department of Mechanical Engineering, Indian Institute of Technology, K a n p u r 208016, India. "~Presently Graduate Student, Department of Mechanical Engineering, University of Houston, Houston, Texas, U.S.A. 159

160

S.S. RAO and ALOK KUMAR

the element is in elastic or plastic state, the element stiffness matrix is calculated using the elastic or plastic stress-strain matrix and the linear strain displacement relations. By using the standard procedures of structural analysis, the displacement vector of the continuum and hence the stresses and strains in each element can be computed. This information about the stresses is used for ascertaining whether the element is in elastic or plastic state and is used in subsequent calculations. As incremental displacement method is used, one can trace the stress and the strain history of each element as the material passes through the rolls. The flow of material under the rolls, the extent of deformation of the material and the pressure distribution on the rolls are also studied in this work.

The .finite element method The various steps involved in the finite element analysis can be stated as follows: 1. Subdivide the material under the rolling process into a number of finite elements (idealization). 2. Assume a suitable displacement model within a finite element which satisfies the compatibility conditions with adjacent elements. 3. Derive the elemental stiffness matrices. 4. Assemble elemental stiffness matrices to construct the global stiffness matrix for the whole material. Also formulate the global load vector. 5. Solve the resulting system of algebraic equations (after the incorporation of appropriate boundary conditions) on a high speed digital computer, and obtain the global displacement vector. 6. From the known displacement vector, calculate the stresses and strains in each of the finite elements. In this work, the cold rolling process is treated as a plane strain problem. This is justified for strip or sheet rolling which involves negligible spread compared to the deformation in the other two directions. Due to symmetry between the rolls, only the material lying under one roll is considered for analysis. As the incremental theory of plasticity is adopted in this study, the force-displacement relation can be expressed in incremental form as dF = [k] d6

(1)

where F represents the nodal force vector, [k] the stiffness matrix and 6 the nodal displacement vector. The stiffness matrix [k] is calculated as [k]=

fff[B]T[D][B]d(Vol)

(2)

where [B] is a matrix which relates the nodal displacements to the strains and is called 'strain rate matrix', [B] r is the transpose of [B], and [D] is the stress-strain matrix. The matrix [D] takes care of the state of the element (elastic or plastic). The elastic matrix, [D(el], and the plastic matrix, [D(P~], are derived based on the Hooke's law and the Prandtl-Reuss equations, respectively.

Elastic stiffness matrix, [k I~] By assuming a linear displacement model u ~- a 1 + a 2 x -? a a y ]

v = a4 + a5x + a6y ,(

(3)

within a triangular element (with nodes i j m), the elastic stiffness matrix of the element can be obtained as [k (~)3 = ~ ~ [B3 TED(~)][B3 t dx dy = [B] TED(e)][B] tA dd

(4)

Finite Element Analysis of Cold Strip Rolling

161

where

[B] = 2A Lci

A = ½det 1 1

0

bj

0

bm

0

cg

0

cj

0

Cm

bi

cj

bj

Cm

b,.

xi

y~

xj

yj

Xm

Y,n

(5)

(6)

= area of the triangle i j m, ai

= Xj

Ym -- xm

(7)

Y3'

bi = Yj - Ym,

(S)

Ci ~

(9)

X m -- Xj,

with the other coefficients obtained by a cyclic permutation of the subscripts in the order i, j, m, and [D (e)] =

E (1 + v)(1 - 2v)

1-v v

v 1-v

0

0

0 0

I

- - 1- 2v 2 = matrix relating the elastic stresses and strains.

(10)

Here (xl, Y3, (x j, y j) and (x,,, Ym) indicate the coordinates of the nodes i, j and m respectively, t is the thickness of the element, and the integration in equation (4) is over the area of the triangle. The form of stiffness in equation (4) is now sufficiently explicit for computation. Plastic stiffness matrix, [k (p)]

If an incremental procedure is used, the matrix relating the incremental plastic stresses to strains, for a plane strain problem, can be obtained as (l_z--2-vv [Dtp)] -

a~2-)

symmetric

E

(1 + v)

1 --2v

--~

(11)

where the material is assumed to obey von Mises yield criterion [14], and S=2~2

1+ ~

(12)

with = ~(tTijtT'ij) 1/2

= equivalent stress,

d~" = (~-dE,5. d~5) '/2,

(13)

(14)

and d~ m'

-

d~ p

= slope of the equivalent stress (~) - plastic strain (.[d~p) curve,

M.I I ] R .

173

I)

(15)

162

S.S. RAO and ALOKKUMAR

a'~ and o-'yare deviatoric stresses and a'~y = r~, is the shear stress, G is the shear modulus and E~, e~ and E~, = 7~y are the plastic strains. The total stress and the total strain are obtained by integrating the incremental quantities calculated at each stage. It is to be noted that the elastic compressibility and strain-hardening characteristics of the material are incorporated in the matrix [D¢P)]. Once [D ¢p)] is known, the plastic stiffness matrix of the element ij m can be obtained as

[k tp)] = f f [B]T[D (p)][Bit dx dy = [B]T[D (p)][B].

t.A.

(16)

It can be seen that equation (11) is valid even for elastic perfectly plastic materials for which H' = 0. At each incremental stage of the calculation which traces the expansion of the elastic-plastic interface, the procedure is simply to replace the elastic stressstrain matrix [D ¢e)] by [D ¢p)] for the yielded triangular elements. By comparing equations (10) and (11), it can be seen that the diagonal elements of [D tp)] are definitely less than the corresponding diagonal elements of [Dte)].tThis amounts to an apparent decrease of stiffness or rigidity due to plastic yielding.

Criterion for yielding In this work, the Von-Mises yield criterion is used. For a general state of stress, this criterion postulates that yielding occurs when 2z + "Czx 2 ] = 3k 2 e 2 = l [ ( a x - - ay) 2 + (O'y - - O-z)2 + (G z -- O'x) 2] + 3[Zx2y + "Cy

(17)

where k is a parameter dependent on the amount of prestrain. The value of k is taken as Y/x/3 where Y is the yield stress of the material in uniaxial tension.

Elastic-plastic behaviour of materials The finite element analysis is conducted for two nonlinear models of the material. First, the problem is solved for the elastic-perfectly plastic behaviour of the material for which de u'

-

d6 p

-

0.

(18)

The second model considered is the Ramberg-Osgood model[151 according to which the a - • relation can be expressed as • =

1+ ~

(19)

where e _~ 0.02 corresponds to the usual engineering definition of yielding. Tensile curves for c~ = 0.02 are shown in Fig. 1 for different values of n. The value of H' can be

n=5 2.0 I0 15

oY

3O 1.0

I

I

I

I

1.0

2.0

.3..0

4.0

EE Y

FIG. 1.

Finite Element Analysis of Cold Strip Rolling

163

obtained from equation (19) as =

--

(20)

n@

Equation (19) represents a compact expression between e and a well suited for digital computer applications. Computational

procedure

The computations start with a division of the material into triangular finite elements as shown in Fig. 2. As the strip is rolled through the rolls, the following boundary conditions are applicable: (i) Fx = 0 and Fy = 0 for all the nodes except those which are under the rolls and those lying on the line C D , (ii) Fx = 0 and v = 0 for all the nodes lying on the line C D , (iii) u = 0 and v = vertical displacement of nodes because of the displacement Ax in x-direction = Ay = ~ + R Z

-

[R

2

-

( B ' E ' ) 2 ] 1/2 -

Yi,

(21)

where B ' E ' = C D + [-R(h1 -

h2) -

(h 1 -

h2)2/4-] 1/2 -

( x i -I- A x )

(22)

for all the nodes which are in contact with the rolls. Here (xi, Y0 are the coordinates of node i, which is in contact with the rolls, and C D , as shown in Fig. 2, is the length of the strip considered for the analysis. The derivation of the expression for Ay in equation (21) is given in Appendix A.

/ A

J ~

C

x

+°

B

O

(a) T° /Ay

x

C

C'

~

x

I,'E

D

(b) FIG. 2.

D'

0

164

S.S. RAO and ALOK KUMAR

The process starts from the stress-free condition and the displacement is given to the material in x-direction in an incremental way. By assuming that (i) the arc of contact does not change, i.e. the rolls are rigid, and (ii) the spread is zero, i.e. the problem can be treated as a plane strain problem, the initial non-steady state of cold rolling is analyzed. It can be noted that for each node, either the force or the displacement component is known and hence the governing equations can be expressed as

ull

I

[Kl1]

nxn

II I

i

[K12]

nxN-n

nxl

-t

I-K21]

I

F1 =

11×

F2

I12

[K22]

.

1

(23)

I

N-

nx n I N-

n x N-

n

N -nx

1

N-nxl

where ul and F 2 are the known displacement and force components, respectively. Equations (23) can be solved for the unknown force and displacement components. A flow chart of the computer procedure is given in Fig. 3. The total stiffness matrix l-K] is stored in the banded form and the nodes are numbered in such a way as to obtain the least band width 'B'. A smaller band width is usually obtained simply by numbering the nodes along the shortest dimension of the continuum. The banded form storage of the [ K ] matrix requires only NB locations as compared to about N2/2 locations for storage of the upper or lower triangle of a full symmetric matrix of order 'N'. Also of importance is the savings in execution time that band operations permit. A band-form Gauss elimination solver is used to solve the simultaneous equations that result after the incorporation of the boundary conditions [16]. By counting the basic operations, i.e. one multiplication or division plus one addition or subtraction the solution time can be estimated to be approximately proportional to NB2/2 in the case of band matrix as compared to a value of N3/6 for a full symmetric matrix. Once the nodal displacement vector [ul/uz} is known, one gets the new coordinates of the nodal points. The strains and the stresses in the elements can be obtained from = [B] u

(24)

and ,, = [D] [ B] u.

(25)

Numerical results The material properties taken in this analysis are as follows: Modulus of elasticity of the material rolled = 1.27 x 104 tons/in 2 (17.51 x 101° N/m2); Yield strength of the material in tension = 47.8 tons/in 2 (65.92 x 10 v N / m 2 ) ;

Poisson's ratio = 0.3. The length of the strip analyzed is taken as 0.8 in. (0.02032 m). It is found that the strains and the stresses beyond this length are negligible when the strip is rolled. The convergence of the method with different number of finite elements can be seen from the results given in Table 1. The convergence of these results follows the general trend that is expected with an increase in the number of elements (the results can be seen to be better as the aspect ratio of the triangular elements approaches unity). The idealization involving 255 nodes, 510 degrees-of-freedom and 400 elements gave reasonably good results (the discrepancy between the theoretical and experimental values of the roll force can be reduced by using more number of elements and/or by considering the elasticity of rolls into account. Although the result given by the finite element method (using 400 elements) differs from the experimental value slightly, the finite element

Finite Element Analysis of Cold Strip Rolling

Start Read and print input data Apply a unit displacement ) increment in x direction (Ax) and modify the coordinates of nodes Find the corresponding boundary conditions Calculate element stiffness matrices by using equations (4) and (10) for elements in elastic state and equations (16) and (11) for elements in plastic state

Assemble the total stiffness matrix

Solve the equations [K] Au = AF after applying boundary conditions Calculate stress, strain and displacement increments: Act, Ae and Au Calculate total stress, strain and new coordinates and effective stress in each element a, E, u and ~

!

~ Find element i out ~of elastic elements having maximum of ~ = ~i

Yes

Treat the element , i as yielded and print it has yielded

no

is~ > Y'?

[no

Check whether rolling is completed ~Yes Print the results and stop

FIG. 3. Flow chart of computational procedure.

165

S. S. RAO and ALOK KUMAR

166

TABLE 1. CONVERGENCE STUDY* Number of elements

Number of degrees-offreedom

80

110

240

310

320

410

400

510

Approximate (average) size of an element [(ji x im in Fig. 2(a)] 0.020 (0.000508 0.012 (0.0003048 0.010 (0.000254 0.007 (0.0001778

x x x x x x x x

0.010 in. 0.000254m) 0.007 in. 0.0001778 m) 0.007 in. 0.0001778m) 0.007 in. 0.0001778m)

Experimental value [17]

Roll force 39.42 tons (35.06 x 104N) 23.94 tons (21.30 x 104 N) 18.63 tons (16.58 x 104N) 17.97 tons (15.99 x 104N) 16.56 tons (14.73 x 10'* N)

* The coefficient of friction between rolls and strip is taken as 0.1. Initial thickness of the strip = 0.063 in. (0.0016 m); Final thickness of the strip = 0.0378 in. (0.00096 m); Roll radius = 5 in. (0.1270 m). TABLE 2. ROLL FORCE AND TORQUE REQUIRED TO DEFORM THE MATERIAL Elasticplastic model

Initial strip thickness

Final strip thickness

Element configuration

Elasticperfectly plastic model Elasticperfectly plastic model RambergOsgood model Experimental results [17]

0.063 in. (0.0016 m)

0.0378 in. (0.00096 m)

0.063 in. (0.0016 m)

0.0378 in (0.00096 m)

0.063 in. (0.0016 m)

0.0378 in. (0.00096 m)

0.063 in. (0.0016 m)

0.0378 in. (0.00096 m)

Roll force

Roll torque

1._6 [~] 2 7

18.66 tons (16.61 x 104 N)

2.966 ton-in (670.2 N-m)

1,,_...._.,6 1~ 2 7 1,~._..,6 [ ~ 2 7

18.58 tons (16.53 × 104 N)

2.882 ton-in (651.2 N-m)

17.97 tons (15.99 x 104 N)

2.678 ton-in (605.3 N-m)

16.56 tons (14.73 x 104 N)

2.37 ton-in (535.6 N-m)

method is preferred whenever the detailed deformation, stress and strain pattern inside the strip are needed. The idealization involving 400 elements is used to find the pressure distribution on the rolls. Totally, three analyses have been conducted by taking an initial strip thickness of 0.063 in. (0.0016 m). The final thickness of the strip is taken as 0.0378 in. (0.00096 m) and the radius of rolls as 5 in. (0.1270 m) in all the cases. In the first two cases the material is taken as elasti~perfectly plastic and they differ from one another in the orientation of the elements. In the third case, the orientation of the elements is the same as that of the second case, but the behavioural model of the material is taken as the Ramberg-Osgood model. The results of these analyses are shown in Table 2. Figure 4 shows a comparison of the pressure distribution along the roll face obtained for different orientations of the elements and for different models for the elastic-plastic behaviour of the material. This figure indicates that the peak value of the pressure on the rolls goes up with an increase in the roll force. Figure 5 gives the deformation pattern and flow of material under the rolls and the location of elastic and plastic zones for a strip having 0.063 in. (0.0016 m) initial thickness and for a material behaving according to Ramberg-Osgood model.* From Table 2, it can be seen that the roll force and torque predicted by the present analysis are higher than the experimental values. The reasons for this might be (i) the error involved in the experimental results and (ii) the assumption of 'rigid rolls' made in this work. The elasticity of the rolls * Since the actual deformation pattern of the thin strip cannot be reproduced conveniently, Fig. 5 has been given to indicate the nature of deformation qualitatively.

167

Finite Element Analysis of Cold Strip Rolling I

6

2r~]7

Initial thickness =0.063"(0.0016m)

Final thickness =0.0378'~0.00096m) . . . . Roll radius =5"(0.1270m)

121~]67

80.0

Elastic Perfectly plastic model 72.0

Ramberg-Osgood / model

64.0

E Z

56.0

% x

48.0

tO

40.0

32.0 o

24.0

16.0 O8.0

0.04

0.08

0.12

O J6

0.20

0.24

0.28

0.32

0.56

Length of arc of contact, in. or.O.O254m

F1G. 4.

can be considered to improve the accuracy of the results by using Hitchcock's equation [18]. However, this involves the use of the present analysis program in an iterative loop. The execution time (with rigid rolls assumption) for one finite element analysis is about 14min on an IBM 7044 computer. However, if a subsequent analysis is to be made with changed boundary conditions (geometry and hence the stiffness of the material remains unchanged), the computer time is expected to be less than 14 rain. On the other hand, if the elasticity of the rolls is also considered, the method is expected to converge within 2 or 3 iterations (with usual convergence criterion) in about 28-42 min of computer time. Elastic zone Material comes in

contact

with roils

Plastic zone

/

Material

/ c oof m e srolls out

FIG. 5. CONCLUSION

The finite element method of analyzing cold strip rolling removes many restrictions imposed by other methods and offers a possibility of obtaining the information that the other methods could not provide. The method is applicable for complex geometry and takes into account the actual material properties. It reveals the elastic-plastic boundaries, the development of the plastic zone and the detailed stress and strain distributions. The method is equally applicable for the rolling of thin strips as well as rectangular slabs. However in the case of rolling of rectangular slabs, three dimensional elements like solid tetrahedron [-19] and hexahedron [20] have to be used for the idealization. Here, for the same accuracy, the computer time is expected to be slightly more than that of the plane strain case.

168

S.S. RAO and ALOK KUMAR

A c k n o w l e d g e m e n t s - - T h e authors would like to thank the referees for bringing references [9 11] to their

attention.

[1] [2] [-3] [4] [5] [-6] [7] [8] [9] [-10] [11] [12] [13] [14] [15] [16] [17] [-18] [19] [20]

REFERENCES TH. VON KARMAN, M a t h e m a t i k M e c h a n i k 5. 139 141 (1925). W. TR1NKS, Blast Furnace & Steel Plant 25, 617 619 (1937). A. T. TSELIKOV, Metallurgy 6, 61 76 (1936). A. NADA1, J. Appl. Mech. 6, A54-A62 (1939). D. R. BLAND and H. FORD, Proc. Inst. Mech. Eng. 159, 144 (1948). C. W. MACGREGOR and L, F. COFFIN, J. Appl. Mech. 10, A13 (1943). B. B. HUNDV and A. R. E. SINGER, J. Inst. Metals 83, 401-407 (1954-1955). R. HILL, T h e Mathematical T h e o r y o f Plasticity. Clarendon Press, Oxford (1971). J. M. ALEXANDER,Proc. Inst. Mech. Engrs 169, 1021 (1955). T. C. FIRBANK and P. R. LANCASTER, Int. J. Mech. Sci. 7, 847 852 (1965). T. C. FIRBANK and P. R. LANCASTER Int. J. Mech. Sci. 9, 65-67 (1967). C. H. LEE and S. KOBAYASHI, Int. J. Mech. Sci. 12, 349 370 (1970). K. IWATA, K. OSAKADA and S. FUJINO, J. Enong Ind. 94, 697-703 (1972). Y. YAMADA, N. YOSHIMURA and T. SAKURAI, Int. d. Mech. Sci. 10, 343-354 (1968). W. RAMBERG and W. R. OSGOOD, Description of stress-strain curves by three parameters. NACA TN. No. 902, July (1943). R. D. COOK, Concepts and Applications o f Finite Element Analysis. John Wiley, New York (1974). D. R. BLAND, H. FORD and F. ELLIS, J. Iron Steel Inst. 171, 239-245 (1952). D. R. BLAND and H. FORD, J. Iron Steel Inst. 171, 245 249 (1952). R. H. GALLAGHER, J. PADLOG and P. P. BIJLAARD, J. Am. R o c k e t Soc. 32, 70(~707 (1962). R. W. CLOUGH, Comparison of three dimensional finite elements. Proc. Syrup. the Application o f the Finite Element M e t h o d in Civil Engineering, ASCE, Nov. (1969).

APPENDIX

A

Vertical displacement (Ay) o f a node in contact with the rolls for a displacement o f Ax in x-direction

Let the coordinates of node i be (x~, Yl) and the radius of rolls be R. Let the input and output thicknesses of the strip or the bar rolled be h I and h2 respectively. From Fig. 2(b), CQ = CD + DQ = CD + (R 2 -- OF2)1/2

(A1)

where OE = R - (h 1 - h2)/2.

The new coordinates of node i after a movement of Ax are x~ + Ax and y~ + Ay, and E'P = O P -

OF' = R -

OF'.

Since OF' = [R 2 - (B'E')2] 1/2,

and B'E' = D'Q = CQ - (x i + Ax)

Equation (AI) gives B'E' = CD + [R(hl - h2) - (hi - h2)2/4] 1/2 - (xl + Ax)

(A2)

and Yi + A y = E'P + P Q ~ (h2/2) + R - [R 2 - (B'E')2] 1/2. .'. A y = v = h2/2 -t- R - [R 2 - (B'E')2] 1/2 - Yi

where B'E' is given by equation (A2).

(A3)