- Email: [email protected]
Numerical simulation of lubricated contact in rolling processes
Journal of Materials Processing Technology 125–126 (2002) 405–411
Numerical simulation of lubricated contact in rolling processes Romain Boman, Jean-Philippe Ponthot* Faculte´ des Science Applique´es, De´partment ASMA, LTAS—Milieux Continus et Thermome´canique, Universite´ de Lie`ge, 1, Chemin des Chevreuils, B-4000 Lie`ge-1, Belgium Received 20 December 2001; accepted 8 January 2002
Abstract In this paper, the lubrication problem in numerical simulation of rolling process is presented. In this case, the recent and complex model of Marsault for the solution of the mixed lubrication regime has been implemented and tested. This model requires the use of the finite difference method to work properly. We will discuss the advantages and the difficulties encountered when trying to solve the same problem with the finite element method in a general frame. Finally, a finite element formulation for the solution of the time-dependent Reynolds’ equation coupled with the deformation of the workpiece is proposed. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Rolling; Lubrication; Large deformation; Finite element
1. Introduction The general frame of this paper is in the field of numerical simulation of rolling processes. Numerical studies are generally focused on the development of effective constitutive laws for the formed material. However, when industrial problems are considered, realistic results cannot be obtained even if the material behavior model is accurate: other phenomenon have to be taken into account, like complex frictional contact and lubrication. These boundary conditions are extremely non-linear and their role is always significant when dealing with rolling, deep drawing, extrusion, welding, etc. For this reason, new models, which can take into account the local contact conditions (pressure, relative velocity, lubricant viscosity, temperature, local geometry, etc.), have to be developed. These models are based and rely upon the analytical work of tribologists. The main idea is that four different lubrication regimes can occur. They are determined by the ratio h=s, where h is the lubricant film thickness and s the standard deviation of the height of the asperities. Our goal is to solve a general lubricated contact problem. This is a very challenging task because most of the work performed in the field of lubrication analysis for metal forming has been particularized to a given process. For the moment, an effective and accurate model, which takes *
Corresponding author. Tel.: þ32-4-366-9310; fax: þ32-4-366-9141. E-mail addresses: [email protected] (R. Boman), [email protected] (J.-P. Ponthot).
into account asperities flattening, thermal effects, elastic– plastic materials has been developed by Marsault [4] for the rolling process of thin sheets. This very interesting work which relies on the Runge–Kutta method tries to summarize the current knowledge in lubricated contact analysis. For this reason, we decided to start our researches with this complex model and then subsequently extend it to the finite element method (FEM). This paper is divided into three parts. The first one introduces the four lubrication regimes and the equations to be solved. The method used by Marsault in the case of the rolling process is presented. In the second part, all the difficulties encountered when trying to solve the same problem with the FEM are discussed. Finally, in the third part, we present a method for the resolution of transient hydrodynamic lubricated contacts, which is based on the work of Hu and Liu [3].
2. Lubrication regimes and their modeling The complete resolution of a lubricated contact problem is very demanding task because it implies the resolution of the coupled system composed by the equilibrium equations of the body in contact and the motion of the lubricant in the contact area. The main difficulty comes from the roughness of both surface in contact which can significantly influence the lubricant flow. That is why different lubrication regimes were introduced by Wilson [9]. The current local regime is
0924-0136/02/$ – see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 4 - 0 1 3 6 ( 0 2 ) 0 0 2 9 1 - 1
406
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
determined by the ratio of the lubricant film thickness (h) and the compositepsurface ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi roughness (s). The latter is computed by s ¼ s1 þ s2, where s1 and s2 represent the surface roughness of each surface. This allows to consider the real contact like a contact between a rough surface and a smooth one. 2.1. Hydrodynamic regime This regime can be divided into the thick film and the thin film regimes. This thick film regime is reached when h=s > 10. In that case, both surfaces in contact are well separated and the asperities of the metal can be neglected. Then, the lubricant flow is well described by the classical Reynolds’ equation which can be written for a two-dimensional mechanical problem: @ h3 @pb @ @h uhÞ þ (1) ¼ ð @x 12Z @x @x @t where pb is the lubricant pressure, which equals to the contact pressure in this case, u ¼ ðu1 þ u2 Þ=2 the mean velocity of the two sliding surface and Z the lubricant viscosity. This equation is an expression of the lubricant flow conservation integrated on the thickness of the contact area. When the lubricant film thickness becomes smaller, the asperities start to play a role in the lubricant flow. The thin film regime is an extension of the previous one in which the asperities are not yet in contact but the roughness of the surfaces have to be taken into account. This can be done by the definition of flow factors fx and fs into the Reynolds’ equation, which becomes the average Reynolds’ equation firstly introduced by Patir and Cheng [5]: @ h3 @pb @ @h @ v fx uhÞ þ þ sfs (2) ¼ ð @x @x @t @x 2 12Z @x where v ¼ ðu1 u2 Þ is the relative sliding velocity. In this equation, pb and h are now the mean values of the lubricant pressure and film thickness at the scale of the asperities. Analytical and semi-empirical expressions for the pressure flow factor fx and the shear flow factor fs can be found in the literature (e.g. [5,8,10] among others). 2.2. Boundary regime The boundary regime is the critical situation when nearly the whole pressure is supported by the asperities. In that case, the lubricant film thickness is reduced to a few molecular layers. Complex phenomena, which are not well understood for the moment, play a large role in this regime. For example, the chemistry of the lubricant and the oxide layers on the metal surface, micromechanics and wear will have to be studied in the future. However, a constant coefficient of friction can be used as a first approximation. This regime is observed on the top of the asperities in contact and will be used in the following regime to model the friction between asperities.
2.3. Mixed regime This regime is reached when the contact conditions are between the hydrodynamic lubrication and the boundary regime. In this case, the asperities are in contact and the total pressure (p) is shared between the hydrodynamic pressure of the lubricant (pb ), which is constrained to flow in the asperities valleys, and the asperity contact pressure (pa ). p ¼ Apa þ ð1 AÞpb
(3)
The weighting factor A is the real area of contact. It is the ratio of the actual contact surface, corresponding to the asperities in contact, divided by the nominal area of contact. Since the heights of the asperities decrease during the process, the mean film thickness cannot be directly measured. In fact, the mean curve of the asperities’ heights have to be updated. The thickness measured from this updated curve is the current mean thickness (ht ) compared to the thickness measured from the initial curve (h), which is a mathematical parameter used to compute the value of A and which can be negative if A is close to the unity (see Fig. 1). The total shear stress (t) is decomposed in the same way between the shear stress coming from the lubricant viscosity (tb ) and the shear stress coming from the boundary regime that appears on the top of the asperities in contact (tb ): t ¼ Ata þ ð1 AÞtb
(4)
These equations introduce a new unknown to the system to be solved. Consequently, a new equation, which models the asperity flattening, is added for the computation of A. For the particular case of rolling, Marsault [4] introduces the following equation: dA e_ pll ¼ xx f ðhÞ dx vx Ep
(5)
where e_ pl xx is the plastic strain rate in the rolling direction, l a geometric parameter equals to the half mean distance between two asperities, vx the velocity of the sheet, f ðhÞ the density function of the asperities’ heights and Ep a nondimensional strain rate which has to be computed from an asperity flattening model. The asperity flattening model is a relation between the non-dimensional strain rate Ep , the real area of contact A and the non-dimensional hardness of the asperities Ha ,
Fig. 1. Definition of the actual area of contact (A).
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
which is defined by: pa pb Ha ¼ k0
(6)
where k0 is the yield shear stress of the body in contact. We have considered two asperity flattening models: the first one was obtained by Wilson and Sheu [11], who used an upper bound method in a simple equivalent model of asperity indentation. The final semi-empirical relation can be written in the form Ha ¼ Ha ðEp ; AÞ: Ha ¼
2 f1 Ep þ f2
(7)
f1 ¼ 0:515 þ 0:345A 0:860A2 f2 ¼
1 2:571 A A ln ð1 AÞ
(8) (9)
The second asperity model was proposed by Sutcliffe [7]: 2 Ha ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A ð3:81 4:38AÞEp
(10)
407
elastic. When the first contact between asperities occurs, the hydrodynamic model is no longer valid and the mixed regime equations are used. Then, the pressure in the lubricant is not equal to the total pressure because of the increase of the real area of contact. The boundary between the inlet zone and the work zone corresponds with the beginning of the yielding of the sheet. In the work zone, elastic–plastic equation are used. If the rolling speed is high, the lubricant pressure can increase until it reaches the total pressure. In that case, it can be shown that the system of equations must be rewritten in another form assuming that pb ¼ pa ¼ p for numerical stability reasons. That’s why the first formulation is called ‘low speed mixed regime’ and the second one ‘high speed mixed regime’. The outlet zone starts when the sheet becomes elastic and ends when the total pressure equals to zero. The integration of the different systems of equations are performed with a fourth order Runge–Kutta method with an adaptive step in order to keep the error under a chosen tolerance. 3.1. Iterative solution procedure
3. Runge–Kutta/slab method solution of the lubrication problem The resolution of the rolling process has been proposed by Marsault [4], using the equations described above. In addition, he assumed non-rigid rolls and thermal effects. The equilibrium and the elastic–plastic behavior of the sheet is taken into account by a classical slab method which is well suited for some rolling problems but which is only valid for thin sheets. Since the system of equations is different in each different lubrication regime, the contact area must be divided into a number of zones with unknown boundaries before the computation (see Fig. 2). The resolution starts in the inlet zone, where the lubricant enters between the roll and the sheet. In the first zone, a hydrodynamic regime is observed and the sheet remains
Unfortunately, it is not possible to compute directly the solution because the inlet velocity of the sheet and the lubricant flow are unknown. Nevertheless they can be iteratively computed by a double shooting method. The first loop consists of adjusting the lubricant flow until the lubricant pressure drops back to zero at the end of the outlet zone. In the second loop, the inlet velocity is adjusted until the prescribed normal stress is obtained in the sheet, see Fig. 3 for the complete flowchart. 3.2. Numerical results The parameters of the rolling process studied here can be found in Table 1. The roll is assumed rigid and the material is elastic–plastic with a constant yield stress.
Fig. 2. Total pressure (p), lubricant pressure (pb ) and actual area of contact (A) along the contact line.
408
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
Fig. 5. Pressure fields (Sutcliffe). Fig. 3. Double shooting method (Marsault).
Table 1 Parameters for the rolling simulation Lubricant viscosity, Z0 (Pa s) Pressure coefficient (Barus), gl (Pa1 ) Composite roughness, s (mm) Half mean asperity length, l (mm) Friction coefficient for asperities, ma Inlet sheet thickness, t1 (mm) Outlet sheet thickness, t2 (mm) Young’s modulus, Es (GPa) Poisson ratio, ns Yield stress, s0 (MPa) Roll’s radius (rigid), R (mm)
0.01 108 0.5 30 0.25 1 0.7 70 0.3 200 200
Fig. 4. Pressure fields (Wilson).
Figs. 4 and 5 show the solution obtained with each asperity flattening model (Wilson and Sutcliffe). We see that this model has a strong influence on the lubricant pressure. This is the main difficulty of the model: more sophisticated asperity flattening models have to be studied to match with the experiments. However, in this case, the total pressure and the total shear stress are similar. Fig. 6 shows the total pressure field obtained for different rolling speeds ranging from 1 mm/s to 10 m/s. We see that the maximum of pressure decreases as the rolling speed increases. Finally, Fig. 7 shows the outlet values for the real area of contact A and the lubricant film thickness ht . At low rolling speed, the contact is mainly supported by the asperities. The real area of contact is close to one and the film thickness is very small. At higher speeds, the hydrodynamic effect plays a significant role and the film thickness becomes larger. Consequently, the real area of contact decreases.
Fig. 6. p as a function of the rolling speed.
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
409
for the Runge–Kutta integration, 40 loops for the lubricant flow and 10 loops for the inlet velocity of the sheet. This means 4–5 min of CPU time on a 600 MHz alpha station (for a rigid roll and no thermal effects). If time-dependent problems with deformable tools are considered, this CPU time is multiplied by the number of loops used for the convergence of the tool’s geometry at each time step and by the number of time steps. Finally, it would be better to introduce the resolution of the lubrication model into the Newton–Raphson iterations which are performed at each time step. In this way, it is possible to solve simultaneously the lubricant equations and the mechanical behavior of the body in contact.
Fig. 7. ht and A as a function of the rolling speed.
For these reasons, the complex lubrication model of Marsault cannot be directly used in a finite element code which is not specially devoted to the rolling problem. 4.2. Discretization of the Reynolds’ equation
4. Finite element solution 4.1. Advantages and difficulties In the previous section, we have seen that Marsault’s method coupled with a sophisticated lubrication model can lead to accurate results in the case of the rolling process. However, this model can only be applied to the rolling process and needs to be improved and generalized if we want to compute other material forming processes like e.g. deep drawing. The FEM could be applied to obtain a more general formalism which could be introduced in actual non-linear codes, like METAFOR [6]. That is our goal but a lot of difficulties appear very soon: If transient problems have to be solved as well as stationary problems, we have to deal with an unknown contact area which is a function of time. The discretization of the geometry, and particularly the boundary surface of the studied body leads to a big problem. In fact, due to the discretization of a smooth curve into a broken line, the contact region grows or shrinks discontinuously. The study of the lubrication model of Marsault shows that a good approximation to the geometry is required, especially in the inlet zone, to obtain good results. The geometrical values of the process (like the roll’s radius or the initial sheet thickness in rolling, the punch radius in deep drawing, etc.) cannot appear in the equations of the model. The equations must be solved at the nodes of the mesh. These nodes are fixed to the body in a Lagrangian formulation or can be moved by a remeshing algorithm if an arbitrary Lagrangian–Eulerian (ALE) formulation. In all cases, these nodes move at a different speed than the lubricant one, introducing additional convection terms in the equations. The method developed by Marsault is CPU expensive. In fact, it is common to have 10,000 to 30,000 spatial steps
As depicted in Fig. 8, the proposed method consists of defining finite elements in the contact zone. The number of elements increases or decreases if some nodes enters or leave the contact zone. In our case, the pressure and the velocity fields are supposed to be known (resulting e.g. from a FEM computation) and we try to compute the film thickness along the contact zone. These values lead us to the computation of the shear stress and the tangential force applied on the nodes. In each Newton–Raphson iteration for the resolution of the equilibrium equations, the friction forces are computed according to the current pressure and velocity fields. Then, a numerical tangent stiffness matrix can be derived and used to reach the balance of the external forces and internal forces (see e.g. Ponthot [6] for details). Since the finite element mesh used for the lubrication problem continuously moves during the computation, a supplementary convective term has to be introduced into the Reynolds’ equation: @ h3 @p @ @ V @h @h fx Rq fs þ uM ¼ ðuhÞ þ @x @x @x 2 @t w @x 12Z @x (11)
Fig. 8. Definition of the lubrication finite elements.
410
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
Fig. 9. Initial mesh and geometry of the rolling simulation.
where uM is the velocity of the ALE mesh and w is the ALE reference coordinate associated to this mesh. The deduction of the finite elements equations in a matrix form is a classical procedure and was already performed by Hu and Liu [3] in the stationary case. From the Reynolds’ equation, a weak form can be written using the weighted residual method. An integration by parts on the pressure and mean velocity terms then leads to Z L Z L Z L @dh h3 @p @h @uM fx dx ¼ dx þ h dx dh dh @t 12Z @x @x 0 @x 0 0 Z L Z L @dh @dh V ð uhÞ dx Rq fs dx 2 0 @x 0 @x 2 3L 0 1 6 V h3 @pA7 6 7 þ 6dh@ ur h þ Rq fs fx (12) 7 4 2 12Z @x 5 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} q
0
where L is the current length of the contact zone, dh the test function and q the lubricant flux. This term must be evaluated at each boundary of the contact zone in the case of an outlet zone. Otherwise, in the case of an inlet zone, the film thickness must be prescribed and dh ¼ 0. For the spatial discretization, we use linear shape functions and an SUPG [2] technique is necessary to avoid oscillations due to the convective terms. h¼
N X i¼1
Ni hi ;
dh ¼
N X
Ni dhi
(13)
i¼1
Introducing the spatial discretization of the unknown into the weak form, the following matrix form is obtained: C1 h_ C2 h þ q ¼ Su ur þ S V V Sp p
Table 2 Process parameters and material properties Roll’s radius, R (mm) Half initial strip thickness, t0 (mm) Half final strip thickness, t1 (mm) Strip’s mean roughness, Rq1 (mm) Roll’s mean roughness, Rq1 (mm) Young’s modulus, Es (GPa) Poisson ratio, ns Lubricant viscosity at p0 , Z0 (Pa s) Viscosity/pressure factor, g (Pa1 ) Rolling speed, V (m/min)
50 0.50 0.45 1.0 0.2 200 0.3 0.434 002 108 100
algorithm based on the finite volume method [1], the Godunov method, is used to update the values at the Gauss points. Due to the symmetry, only one-half of the problem is considered. The initial mesh is shown in Fig. 9. It consists of 80 6 Q4P0 finite elements. The process parameters and the material properties for the strip, the rolls and the lubricant are presented in Table 2. The rolls are assumed to be rigid. The following non-linear hardening rule is used: sY ¼ 313:79ð1 þ 0:052ep Þ0:295
(15)
The simulation is performed in two steps. During the first one (103 s), the rotating roll goes down and squeezes the strip until the prescribed reduction (10%) is obtained. During the second one, the roll’s axis is fixed in space and the strip, entrained by the roll’s rotation, moves through the roll’s bite region if the shear stress is large enough until a steady state is obtained. Some interesting results are presented in Figs. 10 and 11. Fig. 10 shows the influence of the roughness of the strip on
(14)
The matrices Su , SV and Sp whose detailed expressions are given in [1] are fluidity matrices and correspond respectively to the film entrainment, the surface sliding and the pressure variation. They are evaluated numerically with a Gaussian quadrature. Time integration procedures, based on a fully implicit algorithm, as well as associated boundary conditions can be found in [1].
5. Numerical example: rolling simulation The following example is a strip rolling simulation, which has been previously presented by Hu and Liu [3]. In order to reduce the size of the problem, the ALE formulation is used. For each time step, a classical Lagrangian step is first performed. Then, the strip is remeshed and a convection
Fig. 10. Lubricant film thickness as a function of the contact angle for different roughnesses of the strip.
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
411
sheets. If more generality is required (especially if the lubrication model has to be introduced into a classical non-linear finite element code), a lot of difficulties appear and solving the time-dependent Reynolds’ equation coupled with the deformation of the workpiece becomes a very difficult task. An extension to the finite element procedure in that sense has been proposed.
References
Fig. 11. Lubricant film thickness as a function of the contact angle for different rolling speeds.
the lubricant film thickness distribution along the contact area. Simulations with Rq1 ¼ 1 and 1.5 mm are compared with the smooth case Rq1 ¼ 0 mm. If the roughness of the sheet increases, the film thickness increases. Fig. 11 shows the influence of the rolling speed. Simulations with V ¼ 100, 200 and 300 m/min have been performed. The higher the roll’s speed, the thicker the lubricant film. In each case, the film thickness gradually decreases from the entrance to the exit of the roll-bite region. These results are very similar to those obtained by Hu and Liu [3].
6. Conclusion In this paper, the lubrication problem in numerical simulation of forming processes has been presented. For the moment, complex lubrication models are available for stationary processes like rolling, but they are limited to thin
[1] R. Boman, J.-P. Ponthot, Numerical simulation of lubricated contact between solids in metal forming processes using the arbitrary Lagrangian–Eulerian formulation, in: Mori (Ed.), Simulation of Materials Processing: Theory, Methods and Applications, Proceedings of NUMIFORM 2001, 2001, pp. 45–53. [2] A.N. Brooks, T.J.R. Hugues, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Eng. 32 (1982) 199–259. [3] W.-K. Hu, W.K. Liu, An ALE hydrodynamic lubrication finite element method with application to strip rolling, Int. J. Num. Meth. Eng. 36 (1993) 855–880. [4] N. Marsault, Modelisation du regime de lubrification mixte en laminage a` froid, Ph.D. Thesis, Ecole Nationale Superieure des Mines de Paris, France, 1998. [5] N. Patir, H.S. Cheng, An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication, J. Lubr. Technol. 100 (1978) 12–17. [6] J.-P. Ponthot, Traitement unifie´ de la me´ canique des milieux continus solides en grandes transformations par la me´ thode des e´ le´ ments finis, Ph.D. Thesis, Universite´ de Lie`ge, Lie`ge, Belgium, 1995. [7] M.P.F. Sutcliffe, Surface asperity deformation in metal forming processes, Int. J. Mech. Sci. 30 (1988) 847–868. [8] J.H. Tripp, Surface roughness effects in hydrodynamic lubrication: the flow factor method, J. Lubr. Technol. 105 (1983) 458–465. [9] W.R.D. Wilson, Friction and lubrication in bulk metal forming process, J. Appl. Met. Work. 1 (1979) 7–19. [10] W.R.D. Wilson, N. Marsault, Partial hydrodynamic lubrication with large fractional contact areas, ASME—J. Tribol. 120 (1998) 1–5. [11] W.R.D. Wilson, S. Sheu, Real area of contact and boundary friction in metal forming, Int. J. Mech. Sci. 30 (1988) 475–489.
Numerical simulation of lubricated contact in rolling processes Romain Boman, Jean-Philippe Ponthot* Faculte´ des Science Applique´es, De´partment ASMA, LTAS—Milieux Continus et Thermome´canique, Universite´ de Lie`ge, 1, Chemin des Chevreuils, B-4000 Lie`ge-1, Belgium Received 20 December 2001; accepted 8 January 2002
Abstract In this paper, the lubrication problem in numerical simulation of rolling process is presented. In this case, the recent and complex model of Marsault for the solution of the mixed lubrication regime has been implemented and tested. This model requires the use of the finite difference method to work properly. We will discuss the advantages and the difficulties encountered when trying to solve the same problem with the finite element method in a general frame. Finally, a finite element formulation for the solution of the time-dependent Reynolds’ equation coupled with the deformation of the workpiece is proposed. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Rolling; Lubrication; Large deformation; Finite element
1. Introduction The general frame of this paper is in the field of numerical simulation of rolling processes. Numerical studies are generally focused on the development of effective constitutive laws for the formed material. However, when industrial problems are considered, realistic results cannot be obtained even if the material behavior model is accurate: other phenomenon have to be taken into account, like complex frictional contact and lubrication. These boundary conditions are extremely non-linear and their role is always significant when dealing with rolling, deep drawing, extrusion, welding, etc. For this reason, new models, which can take into account the local contact conditions (pressure, relative velocity, lubricant viscosity, temperature, local geometry, etc.), have to be developed. These models are based and rely upon the analytical work of tribologists. The main idea is that four different lubrication regimes can occur. They are determined by the ratio h=s, where h is the lubricant film thickness and s the standard deviation of the height of the asperities. Our goal is to solve a general lubricated contact problem. This is a very challenging task because most of the work performed in the field of lubrication analysis for metal forming has been particularized to a given process. For the moment, an effective and accurate model, which takes *
Corresponding author. Tel.: þ32-4-366-9310; fax: þ32-4-366-9141. E-mail addresses: [email protected] (R. Boman), [email protected] (J.-P. Ponthot).
into account asperities flattening, thermal effects, elastic– plastic materials has been developed by Marsault [4] for the rolling process of thin sheets. This very interesting work which relies on the Runge–Kutta method tries to summarize the current knowledge in lubricated contact analysis. For this reason, we decided to start our researches with this complex model and then subsequently extend it to the finite element method (FEM). This paper is divided into three parts. The first one introduces the four lubrication regimes and the equations to be solved. The method used by Marsault in the case of the rolling process is presented. In the second part, all the difficulties encountered when trying to solve the same problem with the FEM are discussed. Finally, in the third part, we present a method for the resolution of transient hydrodynamic lubricated contacts, which is based on the work of Hu and Liu [3].
2. Lubrication regimes and their modeling The complete resolution of a lubricated contact problem is very demanding task because it implies the resolution of the coupled system composed by the equilibrium equations of the body in contact and the motion of the lubricant in the contact area. The main difficulty comes from the roughness of both surface in contact which can significantly influence the lubricant flow. That is why different lubrication regimes were introduced by Wilson [9]. The current local regime is
0924-0136/02/$ – see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 4 - 0 1 3 6 ( 0 2 ) 0 0 2 9 1 - 1
406
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
determined by the ratio of the lubricant film thickness (h) and the compositepsurface ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi roughness (s). The latter is computed by s ¼ s1 þ s2, where s1 and s2 represent the surface roughness of each surface. This allows to consider the real contact like a contact between a rough surface and a smooth one. 2.1. Hydrodynamic regime This regime can be divided into the thick film and the thin film regimes. This thick film regime is reached when h=s > 10. In that case, both surfaces in contact are well separated and the asperities of the metal can be neglected. Then, the lubricant flow is well described by the classical Reynolds’ equation which can be written for a two-dimensional mechanical problem: @ h3 @pb @ @h uhÞ þ (1) ¼ ð @x 12Z @x @x @t where pb is the lubricant pressure, which equals to the contact pressure in this case, u ¼ ðu1 þ u2 Þ=2 the mean velocity of the two sliding surface and Z the lubricant viscosity. This equation is an expression of the lubricant flow conservation integrated on the thickness of the contact area. When the lubricant film thickness becomes smaller, the asperities start to play a role in the lubricant flow. The thin film regime is an extension of the previous one in which the asperities are not yet in contact but the roughness of the surfaces have to be taken into account. This can be done by the definition of flow factors fx and fs into the Reynolds’ equation, which becomes the average Reynolds’ equation firstly introduced by Patir and Cheng [5]: @ h3 @pb @ @h @ v fx uhÞ þ þ sfs (2) ¼ ð @x @x @t @x 2 12Z @x where v ¼ ðu1 u2 Þ is the relative sliding velocity. In this equation, pb and h are now the mean values of the lubricant pressure and film thickness at the scale of the asperities. Analytical and semi-empirical expressions for the pressure flow factor fx and the shear flow factor fs can be found in the literature (e.g. [5,8,10] among others). 2.2. Boundary regime The boundary regime is the critical situation when nearly the whole pressure is supported by the asperities. In that case, the lubricant film thickness is reduced to a few molecular layers. Complex phenomena, which are not well understood for the moment, play a large role in this regime. For example, the chemistry of the lubricant and the oxide layers on the metal surface, micromechanics and wear will have to be studied in the future. However, a constant coefficient of friction can be used as a first approximation. This regime is observed on the top of the asperities in contact and will be used in the following regime to model the friction between asperities.
2.3. Mixed regime This regime is reached when the contact conditions are between the hydrodynamic lubrication and the boundary regime. In this case, the asperities are in contact and the total pressure (p) is shared between the hydrodynamic pressure of the lubricant (pb ), which is constrained to flow in the asperities valleys, and the asperity contact pressure (pa ). p ¼ Apa þ ð1 AÞpb
(3)
The weighting factor A is the real area of contact. It is the ratio of the actual contact surface, corresponding to the asperities in contact, divided by the nominal area of contact. Since the heights of the asperities decrease during the process, the mean film thickness cannot be directly measured. In fact, the mean curve of the asperities’ heights have to be updated. The thickness measured from this updated curve is the current mean thickness (ht ) compared to the thickness measured from the initial curve (h), which is a mathematical parameter used to compute the value of A and which can be negative if A is close to the unity (see Fig. 1). The total shear stress (t) is decomposed in the same way between the shear stress coming from the lubricant viscosity (tb ) and the shear stress coming from the boundary regime that appears on the top of the asperities in contact (tb ): t ¼ Ata þ ð1 AÞtb
(4)
These equations introduce a new unknown to the system to be solved. Consequently, a new equation, which models the asperity flattening, is added for the computation of A. For the particular case of rolling, Marsault [4] introduces the following equation: dA e_ pll ¼ xx f ðhÞ dx vx Ep
(5)
where e_ pl xx is the plastic strain rate in the rolling direction, l a geometric parameter equals to the half mean distance between two asperities, vx the velocity of the sheet, f ðhÞ the density function of the asperities’ heights and Ep a nondimensional strain rate which has to be computed from an asperity flattening model. The asperity flattening model is a relation between the non-dimensional strain rate Ep , the real area of contact A and the non-dimensional hardness of the asperities Ha ,
Fig. 1. Definition of the actual area of contact (A).
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
which is defined by: pa pb Ha ¼ k0
(6)
where k0 is the yield shear stress of the body in contact. We have considered two asperity flattening models: the first one was obtained by Wilson and Sheu [11], who used an upper bound method in a simple equivalent model of asperity indentation. The final semi-empirical relation can be written in the form Ha ¼ Ha ðEp ; AÞ: Ha ¼
2 f1 Ep þ f2
(7)
f1 ¼ 0:515 þ 0:345A 0:860A2 f2 ¼
1 2:571 A A ln ð1 AÞ
(8) (9)
The second asperity model was proposed by Sutcliffe [7]: 2 Ha ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A ð3:81 4:38AÞEp
(10)
407
elastic. When the first contact between asperities occurs, the hydrodynamic model is no longer valid and the mixed regime equations are used. Then, the pressure in the lubricant is not equal to the total pressure because of the increase of the real area of contact. The boundary between the inlet zone and the work zone corresponds with the beginning of the yielding of the sheet. In the work zone, elastic–plastic equation are used. If the rolling speed is high, the lubricant pressure can increase until it reaches the total pressure. In that case, it can be shown that the system of equations must be rewritten in another form assuming that pb ¼ pa ¼ p for numerical stability reasons. That’s why the first formulation is called ‘low speed mixed regime’ and the second one ‘high speed mixed regime’. The outlet zone starts when the sheet becomes elastic and ends when the total pressure equals to zero. The integration of the different systems of equations are performed with a fourth order Runge–Kutta method with an adaptive step in order to keep the error under a chosen tolerance. 3.1. Iterative solution procedure
3. Runge–Kutta/slab method solution of the lubrication problem The resolution of the rolling process has been proposed by Marsault [4], using the equations described above. In addition, he assumed non-rigid rolls and thermal effects. The equilibrium and the elastic–plastic behavior of the sheet is taken into account by a classical slab method which is well suited for some rolling problems but which is only valid for thin sheets. Since the system of equations is different in each different lubrication regime, the contact area must be divided into a number of zones with unknown boundaries before the computation (see Fig. 2). The resolution starts in the inlet zone, where the lubricant enters between the roll and the sheet. In the first zone, a hydrodynamic regime is observed and the sheet remains
Unfortunately, it is not possible to compute directly the solution because the inlet velocity of the sheet and the lubricant flow are unknown. Nevertheless they can be iteratively computed by a double shooting method. The first loop consists of adjusting the lubricant flow until the lubricant pressure drops back to zero at the end of the outlet zone. In the second loop, the inlet velocity is adjusted until the prescribed normal stress is obtained in the sheet, see Fig. 3 for the complete flowchart. 3.2. Numerical results The parameters of the rolling process studied here can be found in Table 1. The roll is assumed rigid and the material is elastic–plastic with a constant yield stress.
Fig. 2. Total pressure (p), lubricant pressure (pb ) and actual area of contact (A) along the contact line.
408
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
Fig. 5. Pressure fields (Sutcliffe). Fig. 3. Double shooting method (Marsault).
Table 1 Parameters for the rolling simulation Lubricant viscosity, Z0 (Pa s) Pressure coefficient (Barus), gl (Pa1 ) Composite roughness, s (mm) Half mean asperity length, l (mm) Friction coefficient for asperities, ma Inlet sheet thickness, t1 (mm) Outlet sheet thickness, t2 (mm) Young’s modulus, Es (GPa) Poisson ratio, ns Yield stress, s0 (MPa) Roll’s radius (rigid), R (mm)
0.01 108 0.5 30 0.25 1 0.7 70 0.3 200 200
Fig. 4. Pressure fields (Wilson).
Figs. 4 and 5 show the solution obtained with each asperity flattening model (Wilson and Sutcliffe). We see that this model has a strong influence on the lubricant pressure. This is the main difficulty of the model: more sophisticated asperity flattening models have to be studied to match with the experiments. However, in this case, the total pressure and the total shear stress are similar. Fig. 6 shows the total pressure field obtained for different rolling speeds ranging from 1 mm/s to 10 m/s. We see that the maximum of pressure decreases as the rolling speed increases. Finally, Fig. 7 shows the outlet values for the real area of contact A and the lubricant film thickness ht . At low rolling speed, the contact is mainly supported by the asperities. The real area of contact is close to one and the film thickness is very small. At higher speeds, the hydrodynamic effect plays a significant role and the film thickness becomes larger. Consequently, the real area of contact decreases.
Fig. 6. p as a function of the rolling speed.
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
409
for the Runge–Kutta integration, 40 loops for the lubricant flow and 10 loops for the inlet velocity of the sheet. This means 4–5 min of CPU time on a 600 MHz alpha station (for a rigid roll and no thermal effects). If time-dependent problems with deformable tools are considered, this CPU time is multiplied by the number of loops used for the convergence of the tool’s geometry at each time step and by the number of time steps. Finally, it would be better to introduce the resolution of the lubrication model into the Newton–Raphson iterations which are performed at each time step. In this way, it is possible to solve simultaneously the lubricant equations and the mechanical behavior of the body in contact.
Fig. 7. ht and A as a function of the rolling speed.
For these reasons, the complex lubrication model of Marsault cannot be directly used in a finite element code which is not specially devoted to the rolling problem. 4.2. Discretization of the Reynolds’ equation
4. Finite element solution 4.1. Advantages and difficulties In the previous section, we have seen that Marsault’s method coupled with a sophisticated lubrication model can lead to accurate results in the case of the rolling process. However, this model can only be applied to the rolling process and needs to be improved and generalized if we want to compute other material forming processes like e.g. deep drawing. The FEM could be applied to obtain a more general formalism which could be introduced in actual non-linear codes, like METAFOR [6]. That is our goal but a lot of difficulties appear very soon: If transient problems have to be solved as well as stationary problems, we have to deal with an unknown contact area which is a function of time. The discretization of the geometry, and particularly the boundary surface of the studied body leads to a big problem. In fact, due to the discretization of a smooth curve into a broken line, the contact region grows or shrinks discontinuously. The study of the lubrication model of Marsault shows that a good approximation to the geometry is required, especially in the inlet zone, to obtain good results. The geometrical values of the process (like the roll’s radius or the initial sheet thickness in rolling, the punch radius in deep drawing, etc.) cannot appear in the equations of the model. The equations must be solved at the nodes of the mesh. These nodes are fixed to the body in a Lagrangian formulation or can be moved by a remeshing algorithm if an arbitrary Lagrangian–Eulerian (ALE) formulation. In all cases, these nodes move at a different speed than the lubricant one, introducing additional convection terms in the equations. The method developed by Marsault is CPU expensive. In fact, it is common to have 10,000 to 30,000 spatial steps
As depicted in Fig. 8, the proposed method consists of defining finite elements in the contact zone. The number of elements increases or decreases if some nodes enters or leave the contact zone. In our case, the pressure and the velocity fields are supposed to be known (resulting e.g. from a FEM computation) and we try to compute the film thickness along the contact zone. These values lead us to the computation of the shear stress and the tangential force applied on the nodes. In each Newton–Raphson iteration for the resolution of the equilibrium equations, the friction forces are computed according to the current pressure and velocity fields. Then, a numerical tangent stiffness matrix can be derived and used to reach the balance of the external forces and internal forces (see e.g. Ponthot [6] for details). Since the finite element mesh used for the lubrication problem continuously moves during the computation, a supplementary convective term has to be introduced into the Reynolds’ equation: @ h3 @p @ @ V @h @h fx Rq fs þ uM ¼ ðuhÞ þ @x @x @x 2 @t w @x 12Z @x (11)
Fig. 8. Definition of the lubrication finite elements.
410
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
Fig. 9. Initial mesh and geometry of the rolling simulation.
where uM is the velocity of the ALE mesh and w is the ALE reference coordinate associated to this mesh. The deduction of the finite elements equations in a matrix form is a classical procedure and was already performed by Hu and Liu [3] in the stationary case. From the Reynolds’ equation, a weak form can be written using the weighted residual method. An integration by parts on the pressure and mean velocity terms then leads to Z L Z L Z L @dh h3 @p @h @uM fx dx ¼ dx þ h dx dh dh @t 12Z @x @x 0 @x 0 0 Z L Z L @dh @dh V ð uhÞ dx Rq fs dx 2 0 @x 0 @x 2 3L 0 1 6 V h3 @pA7 6 7 þ 6dh@ ur h þ Rq fs fx (12) 7 4 2 12Z @x 5 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} q
0
where L is the current length of the contact zone, dh the test function and q the lubricant flux. This term must be evaluated at each boundary of the contact zone in the case of an outlet zone. Otherwise, in the case of an inlet zone, the film thickness must be prescribed and dh ¼ 0. For the spatial discretization, we use linear shape functions and an SUPG [2] technique is necessary to avoid oscillations due to the convective terms. h¼
N X i¼1
Ni hi ;
dh ¼
N X
Ni dhi
(13)
i¼1
Introducing the spatial discretization of the unknown into the weak form, the following matrix form is obtained: C1 h_ C2 h þ q ¼ Su ur þ S V V Sp p
Table 2 Process parameters and material properties Roll’s radius, R (mm) Half initial strip thickness, t0 (mm) Half final strip thickness, t1 (mm) Strip’s mean roughness, Rq1 (mm) Roll’s mean roughness, Rq1 (mm) Young’s modulus, Es (GPa) Poisson ratio, ns Lubricant viscosity at p0 , Z0 (Pa s) Viscosity/pressure factor, g (Pa1 ) Rolling speed, V (m/min)
50 0.50 0.45 1.0 0.2 200 0.3 0.434 002 108 100
algorithm based on the finite volume method [1], the Godunov method, is used to update the values at the Gauss points. Due to the symmetry, only one-half of the problem is considered. The initial mesh is shown in Fig. 9. It consists of 80 6 Q4P0 finite elements. The process parameters and the material properties for the strip, the rolls and the lubricant are presented in Table 2. The rolls are assumed to be rigid. The following non-linear hardening rule is used: sY ¼ 313:79ð1 þ 0:052ep Þ0:295
(15)
The simulation is performed in two steps. During the first one (103 s), the rotating roll goes down and squeezes the strip until the prescribed reduction (10%) is obtained. During the second one, the roll’s axis is fixed in space and the strip, entrained by the roll’s rotation, moves through the roll’s bite region if the shear stress is large enough until a steady state is obtained. Some interesting results are presented in Figs. 10 and 11. Fig. 10 shows the influence of the roughness of the strip on
(14)
The matrices Su , SV and Sp whose detailed expressions are given in [1] are fluidity matrices and correspond respectively to the film entrainment, the surface sliding and the pressure variation. They are evaluated numerically with a Gaussian quadrature. Time integration procedures, based on a fully implicit algorithm, as well as associated boundary conditions can be found in [1].
5. Numerical example: rolling simulation The following example is a strip rolling simulation, which has been previously presented by Hu and Liu [3]. In order to reduce the size of the problem, the ALE formulation is used. For each time step, a classical Lagrangian step is first performed. Then, the strip is remeshed and a convection
Fig. 10. Lubricant film thickness as a function of the contact angle for different roughnesses of the strip.
R. Boman, J.-P. Ponthot / Journal of Materials Processing Technology 125–126 (2002) 405–411
411
sheets. If more generality is required (especially if the lubrication model has to be introduced into a classical non-linear finite element code), a lot of difficulties appear and solving the time-dependent Reynolds’ equation coupled with the deformation of the workpiece becomes a very difficult task. An extension to the finite element procedure in that sense has been proposed.
References
Fig. 11. Lubricant film thickness as a function of the contact angle for different rolling speeds.
the lubricant film thickness distribution along the contact area. Simulations with Rq1 ¼ 1 and 1.5 mm are compared with the smooth case Rq1 ¼ 0 mm. If the roughness of the sheet increases, the film thickness increases. Fig. 11 shows the influence of the rolling speed. Simulations with V ¼ 100, 200 and 300 m/min have been performed. The higher the roll’s speed, the thicker the lubricant film. In each case, the film thickness gradually decreases from the entrance to the exit of the roll-bite region. These results are very similar to those obtained by Hu and Liu [3].
6. Conclusion In this paper, the lubrication problem in numerical simulation of forming processes has been presented. For the moment, complex lubrication models are available for stationary processes like rolling, but they are limited to thin
[1] R. Boman, J.-P. Ponthot, Numerical simulation of lubricated contact between solids in metal forming processes using the arbitrary Lagrangian–Eulerian formulation, in: Mori (Ed.), Simulation of Materials Processing: Theory, Methods and Applications, Proceedings of NUMIFORM 2001, 2001, pp. 45–53. [2] A.N. Brooks, T.J.R. Hugues, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Eng. 32 (1982) 199–259. [3] W.-K. Hu, W.K. Liu, An ALE hydrodynamic lubrication finite element method with application to strip rolling, Int. J. Num. Meth. Eng. 36 (1993) 855–880. [4] N. Marsault, Modelisation du regime de lubrification mixte en laminage a` froid, Ph.D. Thesis, Ecole Nationale Superieure des Mines de Paris, France, 1998. [5] N. Patir, H.S. Cheng, An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication, J. Lubr. Technol. 100 (1978) 12–17. [6] J.-P. Ponthot, Traitement unifie´ de la me´ canique des milieux continus solides en grandes transformations par la me´ thode des e´ le´ ments finis, Ph.D. Thesis, Universite´ de Lie`ge, Lie`ge, Belgium, 1995. [7] M.P.F. Sutcliffe, Surface asperity deformation in metal forming processes, Int. J. Mech. Sci. 30 (1988) 847–868. [8] J.H. Tripp, Surface roughness effects in hydrodynamic lubrication: the flow factor method, J. Lubr. Technol. 105 (1983) 458–465. [9] W.R.D. Wilson, Friction and lubrication in bulk metal forming process, J. Appl. Met. Work. 1 (1979) 7–19. [10] W.R.D. Wilson, N. Marsault, Partial hydrodynamic lubrication with large fractional contact areas, ASME—J. Tribol. 120 (1998) 1–5. [11] W.R.D. Wilson, S. Sheu, Real area of contact and boundary friction in metal forming, Int. J. Mech. Sci. 30 (1988) 475–489.