Halfspace Simulation of Rough Surface Contact in Metal Forming

Friction has a significant influence on the product properties and tool lifetime in metal forming processes. Therefore it frequently limits the possibilities in the product design. A detailed understanding of friction is necessary to overcome this obstacle. Numerical simulations as well as experimental work are required to gain more knowledge on the processes in the contact interface. The multiscale character of surface roughness makes a fine discretisation necessary, which leads to a huge numerical effort. Consequently alternative approaches to a Finite-Element model are of great interest. Halfspace models have shown their potential in elastic calculations and in plastic simulations with a simplified modelling of plastic deformation. In order to advance the modelling of rough contact in metal forming an enhanced plasticity model is applied within the halfspace framework.


INTRODUCTION
In the contact of rough surfaces only the peaks of surface roughness are in contact for moderate loads.Thus the real contact area A real is smaller than the apparent contact area A o , which is the total area of the contact surface.The primary causes for friction forces are adhesion within A real and ploughing of surface roughness or hard particles through the softer contact partner [1].Hard particles do not occur in the contact interface in metal forming.For ground tool surfaces the ploughing of hard surface peaks over a weaker contact partner can be precluded, since the tool is usually much harder than the surface structured workpiece.Therefore adhesive forces are assumed to be the major contribution to dry friction in metal forming.An upper limit for the shear stress in the contact interface is given by the shear limit of the surface and the ratio between A real and A o .Due to the complex composition of the surface its shear limit is difficult to determine, but an upper limit is given by the shear limit of the underlying material.The important effect of lubrication shall be disregarded for the moment and the focus shall lie on the accurate determination of the real contact area.Determining the real contact area is an important matter in contact mechanics since a long time.A linear relation between the normal force and A real is expected for many contact situations, since Amontons laws state a linear relation between normal force and friction force and an independence of friction from the apparent contact area.Archard found such a linear dependency for a hierarchy of elastic contacts [2].But the assumption of purely elastic contact does not hold for the contact of a rough surfaces in initial contact.Bowden and Tabor argue that plastic surface deformation is induced by very high local stresses caused by small real contact areas [1] [3].They suggest a linear relationship between the mean pressure p mean and A real .
The maximum surface pressure the softer contact partner can withstand is denoted as the surface hardness H.Both the elastic contact and the plastic contact, as described by Bowden and Tabor, were implemented in halfspace models by many authors.Kalker did pioneering work concerning the elastic halfspace model [4].Coupling between normal and tangential contact was studied by Willner [5].The plastic surface deformation based on the surface hardness model was studied for instance by Tian and Bhushan [6] and Willner [7] [8].A more realistic plasticity model taking into account the three-dimensional stress field underneath the surface was presented by Jacq et al. [9] and was subsequently enhanced by Wang and Keer [10], Nélias et al. [11] [12] [13] and other authors.This contact model will be applied to the contact of rough surfaces in metal forming in this article.

Halfspace contact problem
In the halfspace contact model the surface is discretised into N x × N y rectangular patches of the size 2 a × 2 b.The gap between the contact surfaces is defined by the roughness height h ij , the maximum roughness height h max , the surface deformation u ij and the normal approach u o of the surfaces.
In order to account for the roughness of both contact partners the surface height h ij of the rough surface in the numerical model is the superposition of the surface roughness of both contact partners.When solving the contact problem the nonadhesion and nonpenetration conditions must be fulfilled.These complementarity constraints are defined for surface points inside and outside the contact domain I c as follows.
An interdependence between the elastic surface displacement u e at a point (x k , y k ) and the surface pressure p l at a load point (x l , y l ) is given by the Boussinesq solution.
The combined elastic modulus E * enables the modelling of the contact of two bodies with different material parameters.In purely elastic contact the elastic surface deformation equals the total surface deformation u.
For the case of constant pressure within the surface patches the integral depends only on the scalar pressure value p l and the distance between the centre of the loaded patch (x l , y l ) and the displaced point (x k , y k ).This allows the use of a compliance matrix K of the size , which contains all necessary influence coefficients.
The closed for of this integral was first published by Love [14].
Here the distance between the centre of the loaded patch (x l , y l ) and the displaced point (x k , y k ) is denoted as x and y respectively.Now the total displacement at a surface point can be calculated by the summation over the product of the surface pressures p l at all the M surface segments in contact and the corresponding entries of K.
The balance equation for the surface load is defined as follows.

Solution of the contact problem
The contact problem can be formulated as a field problem by introducing the pressure and displacement vectors p and u.The length of the vectors equals the number M of surface segments in the potential contact area.
The corresponding assembled compliance matrix A is then defined as follows.
The complete contact problem consisting of the gap function ( 2), the compliance relation ( 9) and the balance equation ( 10) can be reformulated in compact form as follows [15].
The subscripts c and n identify surface segments of contact and noncontact respectively.The vectors q c and q n contain exclusively ones.The vectors h c and h n contain the surface height h ij reduced by the maximum surface height h max .The noncontact region can be excluded from the system of equations since it contains no information neither about the pressure nor the normal approach.
The problem can be split into two subproblems by the following substitutions.
Each of these subproblems has to be solved.The solutions for the normal approach and the pressure in the contact zone can then be found by back substitution of b and c.
The presented algorithm allows the calculation of the pressure distribution for a specified mean pressure simultaneously by solving only two systems of equations without an outer loop searching for the corresponding normal approach.In order to fulfil the complementarity constraints of the contact problem (3)(4) the members of the contact domain I c have to be identified correctly via an Active Set Strategy, as illustrated in figure (2).It excludes surface segments with negative pressure from the contact domain and includes segments with a negative gap.Before turning towards the actual solution strategy, the storage requirement of  the system of equations shall be discussed.The size of A is M × M which tends towards 1 shows the required memory for storing the vectors and matrices of the contact problem in full contact.Obviously the memory necessary for storing the vectors p and u and the compliance matrix K are small for relevant contact problem sizes, while the assembled compliance matrix A is in the range of gigabytes or even terabytes.Consequently it should be avoided to store the assembled compliance matrix for large contact problems.Due to its symmetry and its positive definiteness the system of equations can be solved with a variety of efficient solvers.On the part of the direct solvers there is the Cholesky decomposition, while the Gauss-Seidel method itself, the method of successive over-relaxation (SOR), which is a variant of the Gauss-Seidel method, and the conjugate gradient (CG) method represent iterative solvers.The Gauss-Seidel method allows enforcing the nonadhesion condition directly in the solver, which makes an exterior time-consuming Active-Set loop unnecessary.Though the SOR method is preferred by many authors, because it offers the same option and is faster with suitable relaxation factors [6] [7].Due to its high efficiency the CG solver is used to solve the systems of equations ( 16) in the presented halfspace algorithm.In the following b stands for the right sides of the systems of equations q c and h c , while x stands for the unknowns b and c.For simplicity A cc is denoted as A. Starting with an arbitrary initial solution x (0) , the search direction d (0) for Beginning with these initial values the step size α, the solution, the residuum and the search direction can be determined for every iteration κ.
The assembled compliance matrix A is only required in the products A x (0) and A d (κ) .These products can be replaced by the convolutions K * x (0) and K * d (κ) of the matrices x (0) and d (κ) , which are of the size N x × N y and contain the entries of x (0) and d (κ) at the respective positions, with the elementary compliance matrix K.This reformulation of the CG algorithm was proposed by Allwood [15], it avoids the storage of the large matrix A and increases the performance of the CG solver by utilising the convolution theorem.

Plastic halfspace contact problem
Pressure loads on the halfspace surface lead to elastic stresses within the halfspace.This stress field is calculated from the pressure field on the surface via influence functions.Plastic strains are calculated where the equivalent stress exceeds the yield strength of the material.These residual strains in turn lead to residual stresses in the bulk and to a residual surface displacement, which are once more calculated via influence functions.Now the total stress is the superposition of the elastic stresses and the residual stresses.Consequently the yield condition has to be evaluated for the superposed stresses in the iteration loop for the plastic strains.In the convergence test the yield function and the changes in plastic strain are evaluated.A relaxation of plastic strains is performed in order to increase the stability of the algorithm.When convergence in the plasticity loop is obtained the elastic contact is recalculated with the surface updated with the residual surface displacement.Once more a relaxation is performed in order to improve the convergence behaviour of the algorithm.Both the plasticity and the contact loop are illustrated in figure 3.

Elastic stresses
The elastic strains in the bulk due to a uniform pressure on a surface segment can be expressed by derivatives of the Boussinesq potential functions [16].The stress components due to the surface load are obtained by inserting derivatives of these strains into Hooke's law.The Boussinesq potential functions are defined as follows.
The displacements can be expressed as derivatives of the potential functions.
Hooke's law states the stresses as follows.
Finally the stresses in the bulk can be expressed as combinations of derivatives of the Boussinesq potential functions: In this manner a tensor T ij which relates the six stress components in the bulk to the surface pressure in normal direction is obtained.

Residual stresses
Residual strains at an arbitrary position (x, y, z) in a halfspace due to a strained cuboid centred at a point (x o , y o , z o ) can by determined by the superposition of three contributions.The first contribution are the stresses due to plastic strains p = ( p 11 , p 22 , p 33 , p 12 , p 23 , p 31 ) T in a cuboid in infinite space as derived by Chiu [17].This solution is superposed with the solution for a mirror cuboid at the opposite position relative to the halfspace surface with plastic strains pm = ( p 11 , p 22 , p 33 , p 12 , − p 23 , − p 31 ) T .The combination of these stresses leaves the surface free from shear stresses, but there remains a virtual normal pressure σ surf on the surface, which is induced by the stress solutions for infinite space.The final solution is obtained by subtracting the stresses in the bulk due to these virtual normal pressures [18].The calculation of the stresses due to the virtual normal pressures is identical to the calculation of the stresses due to surface pressures as outlined before.The tensor B ijkl relates the six strain components in the strained cuboids to the six stress components in the evaluated cuboid.The virtual normal stresses can be calculated as follows.
A detailed description of this procedure and some benchmarks for the influence functions can be found in Chiu's papers [17] [18] and in [19].Figure 5a shows the pressure distribution under the Hertzian contact of a rigid ball of radius 40 mm on a smooth deformable plane.The analytical solution is in good agreement with the result of the elastic halfspace simulation.Also the elastic-plastic Finite Element simulation and the halfspace simulation are in good agreement.They show a strongly reduced maximum surface pressure at the centre of contact.On the other hand the contact radius is increased so that the same normal load is supported by the contact.The residual surface displacement at the symmetry axis is plotted in Figure 5b.It can be seen that there is a large indent at the center of the contact and material is build up at the border of the contact due to plastic deformation.Again the Finite Element and halfspace results agree very well.Young's modulus is set to 205 000 M P a and Poisson's number is 0.34 in both the analytical and numerical calculation.The yield strength is defined by a linear hardening law depending on the equivalent plastic strain γ as follows.
The equivalent stress at the final stage of contact and the equivalent plastic strain are illustrated in Figure 6.The equivalent stress is close to the yield strength in a roughly spherical zone underneath the whole contact region.This highly stressed zone reaches the surface at the border of the contact.At the centre of the contact there is a small layer with small equivalent stress.Knowing that there is high normal pressure at this location, it can be concluded that there is high hydrostatic stress.The zone of comparatively high equivalent plastic strain is confined to a much smaller zone.Its radius is similar, but it is restricted to low depth.

Contact of rough surfaces
Figure 7 show a rough surface in initial and in deformed state after being loaded with a mean surface pressure of 210 M P a.A distinct reduction of the highest surface peaks can be observed.The size of the real contact area compared to the apparent contact area is shown in figure 8.It can be seen that the elastic and the elastic-plastic simulation generate the same result until a normal pressure of about 70 M P a, where the onset of plastic flow occurs.For higher contact loads there is an increase of the real contact area of about 10 % in the elasticplastic simulation compared to the purely elastic calculation.

Conclusion
A simulation model based on the pioneering work of Jacq [9] has been applied to the contact of rough surfaces.The simulation model was made more efficient by the application of a fast solver with small memory requirements.However the plasticity model is not yet stable for arbitrary contact loads.This challenge has to be faced in order to make the model applicable for the investigation of friction in sheet-bulk metal forming.Furthermore the model has to account for tangential loads in the contact.Last but not least the effect of lubrication has to be modelled in order to provide a fully applicable simulation model for the contact in sheet-bulk metal forming.

Figure 1 .
Figure 1.Surface in contact (a) and discretisation (b) in the halfspace model.

Figure 5 .
Figure 5. Surface pressure (a) and residual surface deformation (b) under Hertzian contact.

Figure 6 .
Figure 6.Equivalent stress (a) and equivalent plastic strain (b) under Hertzian contact.

Figure 7 .
Figure 7. Initial surface (a) and deformed surface after contact (b).

Figure 8 .
Figure 8. Load-area curve in elastic-plastic and purely elastic contact.

Table 1 .
Memory requirements for the contact problem It equals the residuum r (0) for the initial solution.