Modeling boundary friction of coated sheets in sheet metal forming

In forming processes, friction is a local phenomenon influenced by the contact conditions at the tool-sheet metal interface. A multi-scale friction model applicable for coated sheets is developed for the boundary lubrication regime which accounts for the physical behavior of coating and measured surface topographies of sheet and tools. The contact patches and therefore the real area of contact is determined at the tool-sheet metal interface for different contact loading conditions. A single asperity micro-scale ploughing model is adapted at each contact patch to determine the friction force from which the overall coefficient of friction is determined. The friction model is validated using different sets of lab-scale friction tests and cup drawing experiments on zinc coated (GI)


Introduction
The coefficient of friction is one of the key parameters in FE simulations of forming processes which can influence the prediction accuracy of these analyses. In sheet metal forming, different friction regimes may occur depending on the lubricant amount and surface topographies of the tool and sheet metal. Commonly, the boundary lubrication regime prevails in sheet metal forming applications in particular when low amount of lubricant is used. In this regime, the sheet surface is lubricated but the lubricant amount is often not sufficient to fill up the valleys in the sheet surface. Therefore, normal and tangential contact loads are carried only by the asperities on the sheet. Hence, an accurate determination of the real area of contact is necessary for a reliable friction behavior prediction. This study focuses on the modeling and prediction of friction in the boundary lubrication regime. In this regime, the coefficient of friction depends on the micro-scale asperity interactions between the contacting surfaces. During forming processes, the coefficient of friction varies due to the spatial and time dependent evolution of the contact conditions at the tool-sheet interface. The topography of the contacting surfaces, contact loads and material behavior influence the contact mechanics. Traditionally, a constant coefficient of (Coulomb) friction is used in FE simulations. However, it has been observed through various experimental studies [1][2][3] that the coefficient of friction is not constant but depends on the contact conditions at tool-sheet interface in forming processes.
Challen and Oxley [4,5], using the plane strain slip line analysis developed an analytical expression for the coefficient of friction accounting for ploughing and adhesion effects for the case of a single wedge shaped rigid asperity ploughing through a softer material. Expressions for the coefficient of friction are developed for the cutting, wear and ploughing regimes based on the attack angle of the asperity relative to the sliding direction and the material behavior of the substrate. Ma et al. [6] proposed a multi-asperity macro-scale friction model for the aluminum extrusion process by adapting a wedge shaped single asperity model [4] to determine the friction force at each tool-sheet metal contact patch from which the overall coefficient of friction is determined. Similarly, Hol et al. [1] developed a multi-scale friction model for uncoated sheets in deep drawing applications by adapting the coefficient of friction expressions developed by Challen and Oxley [4]. The single asperity model developed by Challen and Oxley [4] is valid for a 2D wedge shaped asperity ploughing through the substrate. Hence adapting such model for ploughing of arbitrary shaped 3D asperities requires model calibration using friction experiments [1]. Furthermore, such 2D models cannot consider the effect of the asperity geometrical orientation relative to the sliding direction on friction forces [7]. Mishra et al. [7] recently proposed an analytical model for an ellipsoidal and elliptical paraboloid shaped asperities ploughing through a softer substrate. In this model, the asperity is assumed to be rigid and the substrate is considered to have ideal plastic behavior. The model also considers the asperity orientation relative to the sliding direction. The total shear force is determined as the contribution from ploughing through the substrate material and shearing of the lubricant boundary layer formed on the lubricated sheet surface. This new single asperity model is more realistic and considered more accurate for a macro-scale friction model.
In addition to the accuracy of the single asperity ploughing model, the reliability of the friction model also depends on the prediction accuracy of the real area of contact at the contacting surfaces. In the forming process, the sheet metal is subjected to normal load, sliding and stretching. Asperity deformation under such loads is determined by the flattening models. Greenwood and Williamson [8] presented a model (GW) to determine the asperity deformation under the normal load using an elastic Hertz contact theory. Later, the GW model was extended to consider asperities of ellipsoidal shapes, elastic-plastic behavior and volume conservation [9][10][11]. However, these models are valid for lightly loaded contacts with small real areas of contact and do not account for strain hardening. Westeneng [3] developed a contact model for a fully plastic contact between an uncoated rough surface and a rigid flat tool. In this model, material hardening effects are considered but experiments are required to calibrate the model. The current research is aimed at modeling the friction for coated steel (GI) sheets. In the literature, there are few contact models available for the coated sheet [12][13][14][15][16]. A semi-analytical normal load contact model proposed by Shisode et al. [16] accounts for the coating and substrate materials, coating thickness and measured surface topography of the rough surface. The model also considers material hardening.
When sliding is present, the real area of contact further increases. According to Tabor [17], the real area must increase to accommodate the additional tangential load as shown by experimental studies [1,[18][19][20][21][22][23][24]. The increase in real area of contact due to sliding can be determined by an analytical expression proposed by Tabor [17]. However, sliding experiments are required to determine the model parameter called shear factor [1].
The surface roughness and therefore the real area of contact evolves as the sheet metal is strained. In the case of free surface straining, surface roughening occurs due to the grain rotations [25,26]. On the contrary, a normal load on the sheet surface combined with a bulk strain substantially increases the asperity deformation and real area of contact as confirmed by several experimental studies [3,[27][28][29][30][31]. In this case, macroscopic plastic deformation of the underlying bulk material reduces the effective hardness leading to easier deformation of asperities which increases the real area of contact. Wilson and Sheu [30] proposed a 2D model to investigate the wedge shaped asperities flattening in the presence of normal load and bulk strain using an upper bound method.
Similarly, Sutcliffe [29] used a slip-line analysis in a plane strain situation to predict asperity deformation. However, these 2D models are too simplistic to account for realistic surface topography such as arbitrary shaped asperities and their distribution. Westeneng [3] adapted these 2D models to derive a model for arbitrary shaped asperities and height distribution assuming an ideal plastic material behavior. Fig. 1 shows the friction mechanism in sheet metal forming. The goal of this paper is to develop a multi-scale friction model for the boundary lubrication regime. In the current model, the physical behavior of the coating and its influence on the friction coefficient is considered as well.
The paper is organized as follows. In section 2, normal load and combined normal load and bulk strain flattening models are summarized. The sliding load flattening model is described in section 3.4. In section 3, the single asperity ploughing model is presented followed by multi-asperity macro-scale friction model in section 3.2. Two types of validation experiments are performed: lab-scale friction experiments (Section 4) and component level cup drawing experiments (Section 5). The implementation of the friction model in a FE code for forming simulations is described in Section 5.

Asperity flattening models
The first step in the friction model is to estimate the real area of contact and therefore identifying the local contact locations at the tool-sheet interface subjected to contact loads. For this purpose, flattening models are required to determine the workpiece asperity deformation. In this section, the flattening models used to determine the asperity deformation under a normal load and combined normal load and bulk strain are presented for the coated sheet.
The asperity flattening due to sliding is covered separately in Section 3.4. The sliding model requires to determine the tangential load at a given normal load to estimate the real area of contact. Therefore, the coefficient of friction at the interface is needed to determine the tangential load. However, the coefficient of friction depends on the real area of contact. Therefore, the real area of contact and the coefficient of friction are solved iteratively. Due to this, the sliding model is summarized after the description of the multi-asperity friction model.

Asperity flattening due to normal load
Asperities on the sheet metal surface are plastically deformed or flattened as the hard and smooth surface of the tool comes into contact with the comparatively soft and rough surface of the workpiece. This leads to a real area of contact which is much smaller than the apparent (nominal) contact area. To determine the real area of contact, different meso-scale models have been proposed [12][13][14][15][16]. The model proposed by Mulki and Mizuno [14] is limited to ideal plastic wedge shaped asperities. Similarly the models proposed by Chen et al. [12] and Chang et al. [13] are limited to spherical asperities.
In the current study, the semi-analytical normal load flattening model developed by Shisode et al. [16] for coated sheets is used. In this study, the same material set is used as in Refs. [16]. The model is able to capture the mechanical behavior of asperities of various sizes and shapes. In this model, results from the series of FE analyses performed on asperities of various sizes, shapes and required coating thickness are used to estimate the contact pressure carried by asperities of the measured rough surface. The real area of contact between the flat rigid tool and rough surface is determined by solving equations of force balance and volume conservation. The results of this model and its comparison with the experiments are presented in Section 4.2.

Asperity flattening due to combined normal load and bulk strain
Combined normal load and bulk strain can increase the real area of contact significantly compared to normal load alone. Westeneng [3] adapted the 2D models [29,30] applicable for ideal plastic asperities to describe flattening of a general 3D rough surface subjected to combined normal load and bulk strain. In the current study, Westeneng's model [3] is used to determine the real area of contact on GI sheet under a combined normal load and bulk strain.
An analytical relation (Equation (1)) is developed by Westeneng [3] to determine the real area of contact α ε as a function of uniaxial strain ε in the sheet. The tool is assumed rigid and smooth which is plastically deforming the soft asperities modeled as ideally plastic material.
where, φ(z) is surface height distribution function of the undeformed surface. d ε is the position of tool from the mean plane of the rough surface and u ε is the uniform rise of non-contacting asperities. φ(z) is determined by fitting a function through the height distribution data of the measured surface. In this study, the height distribution data is fitted by a B-spline function [32] using 10 cubic splines. l m is the mean half asperity spacing and W is the velocity parameter which is defined as, v d and v u are downward velocity of contacting asperities and upward velocity of non-contacting asperities respectively. Westeneng [3] used the results of slip-line analyses performed by Sutcliffe [29] on ideal plastic wedge shaped asperities under a combined normal load and bulk strain to derive the relation for velocity parameter as follows: where, H eff is the effective hardness of the workpiece material defined as H eff = Pnom αε and τ w is the shear strength of the workpiece. τ w = σy0 ̅̅ 3 √ following the von Mises criterion under pure shear. H = Bσ y0 is the initial hardness, σ y0 is the yield strength of the sheet metal (zinc coating for GI sheet) and B is the hardness factor equal to 2.8 for metals as determined by Tabor [33]. The mean half asperity spacing l m in Equation (1) is determined as [3], where, ρ α is the asperity density (number of asperities in contact with the tool per unit area).
Equation (1) is solved using the Euler method to determine α ε at strain ε and nominal pressure P nom . The initial values of α 0 , d 0 and u 0 at nominal pressure P nom and ε = 0 are determined using the normal load flattening model. The asperity density is still an unknown parameter. The asperity density as determined from the measured rough surface is highly dependent on the measurement resolution. Therefore, ρ α for a random surface is described by a relation proposed by Hol et al. [1] as follows: where, the constants c 1 = 1.49 10 3 c 2 = − 0.458, c 3 = 3.14 10 2 and. c 4 = − 0.0460

Estimation of coefficient of friction
Having determined the real area of contact between sheet metal and tool using the flattening models, this section presents a new friction model for sheet metal forming applications of coated sheets in the boundary lubrication regime. In the boundary lubrication regime, the lubricant forms an interfacial boundary layer on the sheet surface due to the chemical and physical adsorption. Therefore, the friction force in the boundary lubrication regime is the result of shearing of the boundary layer and ploughing of the sheet surface by tool asperities. To account for these two effects, the single asperity ploughing model proposed by Mishra [7] is adapted.

Single asperity model
Mishra et al. [7] proposed a model to determine the ploughing forces on a rigid ellipsoidal and elliptical paraboloid shaped asperities ploughing through an ideal plastic substrate. The single asperity model was validated by dedicated single asperity ploughing experiments [34]. The ploughing forces are due to the plastic deformation of the substrate and shearing of the lubricant boundary layer at the interface. Fig. 2 shows the schematic of the single asperity ploughing model. During sliding, only the frontal portion of the asperity (N-S in Fig. 2 right) is in contact with the substrate. The contact surface of the asperity varies with the ploughing depth d and orientation angle β. In their study, the ploughing depth d is determined from the applied normal load F asp and hardness H of the ideal plastic substrate. The contact surface is divided into triangular elements. The forces are calculated on the elemental areas dA. The contact pressure acts normal to dA with the unit normal vector n, whereas the shear stress at the interface acts in the tangential (or the plastic flow) direction t to dA. The total friction force dF due to ploughing and shearing of boundary layer acting over the contact area dA can be written as [7], τ BL is the lubricant boundary layer shear strength. The elemental forces are integrated over the contact area to determine the total shear force from which the coefficient of friction is determined. τ BL is determined experimentally for a specific combination of lubricant and the sheet material. Mishra et al. [35] performed boundary layer shearing experiments on DX56-GI sheets with different types of lubricants. The details of the experiments can be found in Refs. [35]. The boundary layer shear strength is determined at different nominal pressures and described as, τ BL = c 0 P m 0 [35] where c 0 , m are model parameters and P 0 = (P nom /α) is the contact pressure based on the fractional real area (α) of contact. In this study, DX56-GI sheet material lubricated with Fuchs Anticorit PLS100T lubricant is used for which c 0 = 7.34 and m = 0.78 [35]. The model accounts for the change in magnitude and direction of the forces due to ploughing and shearing of the boundary layer based on the orienation, size and ploughing depth of the asperity. The results of the model are shown in Fig. 3 where the influence of the orienation angle and ellipticity ratio (ratio of major and minor axes) on the coefficient of friction is demonstrated. For this sensitivity study, a rigid asperity is ploughing through DX56-GI sheet (lubricated with Fuchs Anticorit PLS100T). The hardness H of the sheet is equal to 210 MPa (H = Bσ y0 ), based on the yield strength of zinc coating σ y0 (75 MPa) [36]. For the example shown in Fig. 3a, an ellipsoidal asperity with ellipticity ratio of 2, semi-major axis of 200μm in the base area is used to demonstrate the effect of asperity orientation. Fig. 3b shows the effect of ellipticity ratio on friction coefficient for an asperity with orientation angle of zero. The applied normal load is 2 N. The orientation angle β is defined as the angle between the major axis of the asperity base and sliding direction (Fig. 2).
During ploughing, the substrate material piles-up in front of the asperity while it is displaced on either side of the asperity. The material pile-up in front of the asperity increases the resistance to ploughing which increases the total ploughing force. The model does not consider this effect. Mishra et al. [34] performed numerical simulations of single ellipsoidal asperities ploughing through an ideal plastic substrate to determine the coefficient of friction. Details of the simulations can be found in Ref. [37]. The results from the simulations are compared with the results of the analytical ploughing model. Consistent material properties, asperity geometry and normal load are used in the simulations and the analytical model. Fig. 4 shows the results of the simulations performed using the material point method (MPM) [34] and the analytical model. The asperity is rigid with different ellipticity ratios ploughing through an ideally plastic steel substrate at different orientation angles, asperity normal load (F asp ) is 1 N, substrate hardness is 450 MPa and boundary layer shear strength (τ BL ) is 0.9H. Mishra et al. [34] concluded that the analytical model and simulation results correlate well at lower ellipticity ratios and orientation angles but deviates as the ellipticity ratio increases. Fig. 4a shows that for an asperity with ellipticity ratio of 2, the results of the analytical model are in good agreement with simulations but deviates at ellipticity ratio of 6.2. The deviations in the results are mainly coming from the effects of material pile-up which are explicitly considered in the MPM simulations. It is noticed that the pile-up effect is more pronounced for ellipsoids of higher ellipticity ratios and orientation angles [34].
Comparion of coefficient of friction between the simulations and the  As the orientation angle increases, the frontal projected area of the asperity increases. This increases the amount of displaced material and the average pile-up height. Hence a pile-up factor is introduced in the analytical model to include the rise in ploughing depth due to the material pile-up. The proposed pile-up factor ψ is determined as follows.   Where, d' is the corrected ploughing depth, d is the initial ploughing depth, R my is the semi-projected length of the asperity perpendicular to the sliding direction and R min is the semi-minor axis of the asperity (see Fig. 2 right). Equation (9) shows that with increase in orientation angle, (R my /R min ) increases leading to higher ploughing depth (d') which results in a higher coefficient of friction as predicted by the numerical study. Similarly, at lower ellipticity ratios and orientation angles, d' ≈ d for which the model results and simulations correlate well. Therefore a geometrical pile-up factor, ψ = 0.3 is determined so as to calibrate the ploughing depth in the analytical model as shown in Fig. 4b. ψ = 0.3 is used in the current study with DX54-GI sheet material.

Multi-asperity friction model
The single asperity ploughing model described in the previous section is adapted to determine the coefficient of friction at macro-scale. The schematic of the macro-scale friction model is shown in Fig. 5. A measured surface height data of the tool and workpiece are used. The asperities deformation and flattening of the workpiece surface is determined using the flattening models. The tool is assumed to be rigid and flat during this step. The surface points of the deformed surface which are in contact with the tool and adjacent to each other form contact patches on the workpiece surface as shown in Fig. 6. Once the contact patches are determined, an elliptical paraboloid shaped asperity is fitted through the penetrated tool height data at each contact patch (Section 3.3 describes the method). The single asperity ploughing model is then applied for each fitted asperity. The initial ploughing depth d at the contact patch is set equal to the height h of the elliptical paraboloid. The ploughing depth d' is determined by Equation (9) using the geometrical parameters of the fitted elliptical paraboloid (asperity height, base area, ellipticity ratio and orientation angle). The friction force contribution from all the contact patches is determined based on the single asperity model from which the total friction force is calculated. Finally, the overall coefficient of friction μ is determined as follows: where, n is the total number of contact patches.

Fitting the elliptical parabola to tool asperities
A procedure described by de Rooij et al. [38] is followed to determine the geometry and orientation of the fitted elliptical paraboloid as shown in Fig. 6. The base of the elliptical paraboloid is approximated by an ellipse with area equal to the area of the contact patch. The lengths of semi-major (R maj ), semi-minor (R min ) axes and orientation angle (β) are determined by having the same normalized second central moments as the contact patch.
To determine the height h of the elliptical paraboloid, height data of the tool is required. Surface height data of the tool is overlaid on the deformed rough surface. The tool height data is lowered such that the total bearing area between the tool and the workpiece is equal to the  Fig. 8. Schematic of rotating friction tester (RFT).
total area of the contact patches. The height h at each contact patch is determined by equating the penetrated volume V pen and the volume of the elliptical paraboloid as follows:

Updating the real area of contact due to sliding
The real area of contact grows as the tool and sheet metal slides relative to each other. The first analytical expression to determine the real area of contact due to sliding was proposed by Tabor [17]. Tabor's analysis is based on understanding that the asperities which have been plastically deformed due to normal load should grow due to the additional tangential load from sliding. Junctions are formed at the contacting points which grow due to the adhesion and ploughing of the tool asperities through the substrate. This is also referred as the junction growth (Fig. 7) [17]. Assuming a perfectly plastic flow at the yielded contact points, the mean pressure at contact locations reduces to accommodate the additional shear stress to maintain a constant von Mises stress. The analytical expression proposed by Tabor [17] is as follows: where, A 0 is the real area of contact due to normal load only and A s is the real area of contact due to sliding, μ is the coefficient of friction and k α is the shear factor constant [32]. During sliding, only the frontal part of the penetrating tool asperity comes in contact with the substrate (see Fig. 7 right). Therefore, the real area of contact must grow by a factor approximately equal to 2 in order to satisfy the force equilibrium. However, due to the additional tangential load, A s > 2A 0 . The theoretical value of k α for simple 2D ideal plastic lightly loaded contacts is 3.0 [17]. However, there is no analytical solution available for 3D contact situations, therefore friction experiments are required to determine k α . The value of k α is expected to be higher for 3D contacts compared to the 2D cases. In this study, the shear factor k α is determined using sliding experiments on GI sheet samples as presented in Section 4.2.

Friction experiments
Normal load and friction experiments are performed on two different GI sheets at different nominal pressures. The real area of contact measured from normal load experiments is compared with the predictions of the flattening model. The friction experiments are used to determine the shear factor k α required in the sliding model. Furthermore, the coefficient of friction measured from the friction experiments are compared with the predictions of the friction model. Two different samples of DX56-GI (CR4-GI) grade with distinct surface roughness (1.0 and 1.3mum) are used. Both samples are textured by electro-discharged textured (EDT) rolls. The sheet thickness is 0.7 mm and the mean coating thickness is 7mum for both samples. Table 1 shows the parameters of the experiments. The samples are lubricated with a small amount of lubricant (0.6 g/m 2 ). The lubricant amount is not enough to completely fill the surface void volume which ensures that the boundary lubrication regime applies in the friction experiments. Fig. 8 shows the schematic of the rotating friction tester (RFT) developed at TATA Steel Europe for normal load and sliding experiments. The set-up consists of a stationary punch and a rotating sample holder. The punch includes 3 highly polished flat notches (S a =27 nm) which are located equidistant from the center of the punch. For the normal load experiments, the punch is moved slowly downward by a hydraulic actuator and pressed against the fixed sample to reach the desired nominal pressure. The surface height measurement of the sample under the notch is performed using confocal microscopy before and after the loading to determine the real area of contact.
For the friction experiments, first the required nominal pressure is applied and afterwards the sample holder is rotated at a predefined sliding speed and rotation angle. The resultant torque due to sliding is measured from which the coefficient of friction is determined.

Material properties
The flow curves of DX56 steel substrate [16] and zinc coating [36] are shown in Fig. 9. The yield strength of DX56 steel substrate is 153 MPa and for zinc coating is 75 MPa.

Real area of contact measurement
The real area of contact of the deformed sample is determined by analyzing the measured surface height data using confocal microscopy. Surfaces of deformed and undeformed samples are measured on the same area of 2 × 2 mm at a pixel resolution of 1.3mum [16]. The surface height distribution function is plotted in Fig. 10a which is used to determine the real area of contact as described in Ref. [16]. Fig. 10 shows the surface topographies and height distribution functions of deformed and undeformed surfaces of sample 1 (S a =1.0mum) at 20 MPa after normal load only and after the friction experiments using the setup depicted in Fig. 8. The measured real areas of contact for samples 1 and 2 at different nominal pressures are shown in Fig. 11.
The normal load flattening model [16] presented in Section 2.1 is used to determine the real area of contact for samples 1 and 2. The results (Fig. 11) shows that the model can predict the real area of contact due to normal load very well. The real area of contact after sliding is used to determine the shear factor k α .

Determine shear factor k α
The real areas of contact measured after sliding are used to determine the shear factor k α in Equation (12). The real areas of contact due to normal load A 0 as used in Equation (12) comes from the normal load flattening model. To determine A s , the shear factor k α and the coefficient of friction μ are required. k α is determined by minimizing the difference between the measured real areas of contact due to sliding and the real areas of contact A s determined by Equation (12)(see Fig. 11). The iterative scheme is shown in Fig. 12. The shear factor k α is optimized for sample 1 and used for sample 2 to check the results. The optimized value Fig. 11. Experiments vs. model: sample 1 and 2 real area of contact due to normal load and sliding. of the shear factor k α from this study is 20 for which the cumulative root mean square (RMS) error in real area of contact due to sliding is <1.5% for both samples.

Measurement of the friction coefficient and comparison with model predictions
The coefficient of friction is determined using the measured torque during sliding and the applied nominal pressure. The coefficient of friction vs. the sliding angle is plotted at different nominal pressures in Fig. 14a. The coefficient of friction shows a decreasing trend with increase in nominal pressure which is in agreement with previous observations [1,3]. The coefficient of friction during the start is high due to the transition from static friction to dynamic friction after which it stabilizes. For comparison with the friction model, the coefficient of friction averaged over the stable part of the curve (40-80deg) is used. Fig. 13 describes the implementation of the friction model. The coefficient of friction determined from the experiments and the model for samples 1 and 2 at different nominal pressures are shown in Fig. 14b. The comparison between the experiments and the model shows that the model predicts the magnitude and the trend of the coefficient of friction very well. Note that only the shear factor k α is fitted from the experimental results to calibrate the sliding load flattening model. The coefficient of friction is the result of the friction model using a calibrated value of shear factor k α .

Application to deep drawing processes
The boundary friction model is implemented for FE simulations of deep drawing processes. Cup drawing experiments are performed on GI sheet to validate the boundary friction model for metal forming processes. The implementation of the boundary friction model for forming simulations is discussed followed by the comparison results between the experiments and the simulations.

Cup drawing experiments
The process parameters of cup drawing experiments are shown in Table 2. The cup is drawn using a deep drawing press available at TATA Steel Europe. Fig. 15 shows the schematic of the set-up and the product.  The experiments are performed at 3 different blank holder forces: 10, 30 and 50 kN. The blank is lubricated on either side with a lubricant amount of 0.3 g/m 2 . The lubricant amount is so small that it falls well within the boundary lubrication regime. The surface roughness and topography are different at different locations of the tooling. Therefore, the surface height measurements are performed at key locations on the tools that will be used in the friction model. Fig. 15 shows the surface roughness of the blank and tools at different locations. To check repeatability, experiments are repeated 3 times at each blank holder force. Fig. 17 shows the punch force vs. punch displacement at different blank holder forces.

Model implementation to forming simulations
During the forming processes, the sheet surface is subjected to different types of loads such as normal load, sliding and stretching which varies at each time increment. This means that the sheet surface evolves continuously during the process. In the FE simulation, contact pressure and equivalent plastic strain at each node of the blank mesh for each time increment can be determined. With these inputs and using the flattening models it is possible to determine the local contact conditions (asperity deformation and real area of contact). Furthermore, the coefficient of friction at each node can be determined using the multiasperity friction model discussed in Section 3.2. However, running the friction model at each node of the blank mesh for each time increment is computationally expensive. Therefore, an offline look-up table is generated which consists of the coefficient of friction and real area of contact for the range of nominal pressures and strains. The look-up table is supplied to the simulation as an input. During the simulation, the nominal contact pressure and the equivalent plastic strain at each node of the blank mesh are extracted from the FE contact to determine the real area of contact for a combined normal load and bulk strain. The coefficient of friction at each contacting node is determined using the lookup table and the nodal contact loads. In this way, the coefficient of friction at each contacting node is applied in the forming simulation. It should be noted that, the look-up table is generated for each tool surface to account for different tool topographies at different parts of the die, punch and blank holder.
The blank is meshed with 3 node triangular Kirchhoff shell elements. Due to symmetry, one quarter of the cup is modeled. Tools are assumed to be rigid. The Vegter yield model [39] is used to describe the yield surface of the sheet metal. The hardening behavior is described by the Bergström-Van Liempt relation [40]. Table 3 shows the parameters of the yield surface and Table 4 shows the parameters used in the hardening relation. Simulations are performed at blank holder forces of 10, 30 and 50 kN. Fig. 16 shows results of the simulation for punch side and die side of the blank at a punch displacement of 32 mm and blank holder force of 30 kN. On the punch side of the blank, the contact pressure is maximum at the punch corner region. The real area of contact is high at  the punch corner region due to a high contact pressure. In the blank holder region, the real area of contact is high due to a high equivalent strain. The coefficient of friction in the punch corner region is high compared to the blank holder region. This is due to the higher tool roughness at the punch corner region compared to the blank holder region (see Fig. 15). On the die side of the blank, the contact pressure and the equivalent strain are maximum at the die corner region resulting in a large real area of contact. Moreover, within the die corner region, the coefficient of friction is high at low contact pressure and low at high contact pressure.

Comparison of experiments and FE simulations
Punch force vs. punch displacement, real area of contact and sheet      shown in Fig. 17b. The results show that punch forces deviate from the experiments at higher blank holder forces.
As an additional validation case, the real area of contact is measured on the die side of the cup at different locations. The results are compared with the simulations as shown in Table 5. The real area of contact correlates well in the blank holder region (R1, D1) at blank holder forces of 30 kN and 50 kN. Elastic recovery of deformed asperities on the measured surfaces is one of the reasons of deviation but it is expected to be low as the plastic deformation is dominant. However, the results are not favorable at a blank holder force of 10 kN; the large difference between measured and predicted real area of contact at D1 (30%) is mainly due to wrinkles formed at the cup edges. Furthermore, the real area of contact measurements at R1 are not reliable due to wrinkling at 10 kN blank holder. On the die corner region (location R2 and D2), the simulations over predicts the real area of contact. But the difference reduces at higher blank holder forces. Overall, the predicted real area of contact at blank holder forces of 30 kN and 50 kN correlate well with the measurements.
Draw-in of the blank is measured and compared with the simulations as well. The material behavior and coefficient of friction influence drawin. For instance, higher the coefficient friction higher the resistance to sliding which results in less draw-in. Fig. 19 compares the measured and FE predicted draw-in values at blank holder forces of 10 kN and 50 kN. The results show a good correlation between the simulations and the measurements. The measured draw-in also shows the earing phenomenon in the rolling and transverse directions due to the material anisotropy which is predicted well by the simulations too.

Conclusions
A new multi-scale boundary friction model for coated sheet is developed. The model combines the recently developed flattening models and single asperity ploughing model to determine the macroscopic friction forces. It considers the influence of contact loads on the evolution of the sheet surface topography and friction. Measured workpiece and tool surfaces are used as the input to determine the coefficient of friction. The contact patches between the tool and the workpiece at normal load, sliding and combined normal load and bulk strain are determined using the flattening models. A single asperity ploughing model is applied at each contact patch from which the coefficient of friction is determined. The model uses the novel elliptical paraboloid single asperity ploughing model to determine the shear forces due to the ploughing of substrate and shearing of the interfacial lubricant boundary layers. The friction model is validated using labscale friction experiments. The model is used in FE analyses of deep drawing applications. Offline friction look-up tables are generated for the measured workpiece and tool surfaces as a function of nominal contact pressure and strain. The coefficient of friction at each node is determined based on the tool surface in contact, the nodal contact pressure and the nodal equivalent strain. The increase in simulation time using the new friction model is less than 2% compared to using constant coefficient of (Coulomb) friction. The cup drawing experiments are performed on GI sheet at different blank holder forces. The results from the experiments and simulations show a good agreement for punch force vs. punch displacement, blank draw-in and real area of contact. This demonstrates the applicability of the developed model for industrial scale forming simulations.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.