Semi-analytical contact model to determine the flattening behavior of coated sheets under normal load

Friction influences the formability in sheet metal forming processes. It depends on the local contact condition between tool and sheet metal. Therefore, precise estimation of real area of contact is the first step for an accurate prediction of friction. In this study, a multi-scale contact model is developed to predict deformation of asperities on a rough uncoated and coated surfaces under normal load. The model accounts for the coating thickness and material behavior of coating and substrate. Finite element simulations are performed to determine the contact pressure of single asperities of different sizes. These are used to determine the real area of contact. The model is validated relative to the experiments performed on uncoated and zinc coated steel sheets.


Introduction
The real area of contact between two surfaces is determined by the persistence of the highest asperities to deformation mainly at the softer surface. Therefore mechanical properties as well as geometry of the asperities or texture of the contacting surfaces govern the real area of contact. Furthermore, it is generally assumed that real area of contact is solely controlled by plastic deformation of the asperities while a perfectly-plastic material model is utilized.
Reliable modeling of friction in sheet metal forming simulations is vital for the accuracy of the formability analyses [1]. Friction in sheet metal forming is a local phenomenon which depends on continuously evolving contact conditions during the forming process. The real contact area is influenced by local contact pressure, surface topographies of the sheet metal as well as the forming tool, their material behavior, lubricant condition and temperature in the sheet metal-tool contact [1][2][3] However, the first step to have a reliable prediction of friction coefficient is to accurately estimate real area of contact at sheet metal-tool interface.
The real area of contact is determined either numerically using finite element (FE) analyses or by analytical contact models. These two approaches can be understood as direct and indirect approach respectively. In the indirect approach, the rough surface is described by a statistical function which leads to an efficient modeling and analytical solution of the surface. Greenwood and Williamson [4] for example, used the elastic Hertz contact theory (GW model) for the contact between rigid flat tool surface and elastic rough surface. The rough surface asperities are described by a summit height-based stochastic parameter. The model assumed spherical asperities near their tips with all asperities having the same radii. Following the GW model, Bush et al. [5] extended the elastic contact model by incorporating the variations in the radius, ellipticity ratio and the height. Furthermore, Chang et al. [6] proposed an elasto-plastic contact model. The model implements elastic-plastic deformation in the asperity, based on the critical indentation depth derived by equating the hardness of the material and maximum contact pressure derived from the Hertz contact theory on the spherical asperity. The model also takes into account the volume conservation in the asperity. The major shortcoming with this approach is that the contact load at the onset of plastic deformation is not continuous. According to their formulation, the contact pressure is allowed to increase at the onset of yielding. Another such example of elasto-plastic contact models is proposed by Zhao et al. [7]. Though these statistical models are simple and computationally effective, these models are still limited to lightly loaded contacts as it does not account for material hardening, asperity interactions and accurate asperity geometry. Furthermore, the use of summit height-based statistical description is ambiguous due to its dependency on the scanning resolution of the surface. Westeneng [8] developed a fully plastic statistically based contact model where the rough surface is modeled as measured height data represented by bars. The tool surface is assumed to be rigid and perfectly flat. Energy and volume conservation equations are solved to determine the indentation and rise of asperities. The model takes into account the hardening effects of the material but requires experimental results to calibrate the model. Furthermore, the model assumes a simple and constant relation between the material hardness and material yield strength H ¼ Bσ y , where for metals B�2.8 [9]. The constant parameter which relates the hardness to yield strength is derived using hardness tests. However, this relation holds for approximately ideal-plastic materials [9]. Moreover, the constant relation is not always true and can vary based on the type of hardness test, indenter shape and size. Furthermore, there is no conclusive finding on the similarity between an indentation situation in the hardness test (rigid indenter on a polished surface) and a flat rigid tool upsetting soft asperities on a sheet surface.
In the direct approach, using e.g. the FE technique, rough surface can be modeled explicitly and solved for the assumed nominal pressures [10,11]. Full scale FE modeling of rough surface can consider all the material and geometrical complexities such as material hardening, arbitrary shape of the asperities, their interactions and bulk straining of the substrate to predict real area of contact. However, the major drawback of the full scale FE method is the excessive complexity in modeling the rough surfaces at the asperity level and the computational time.
Nowadays, coated steels are widely used in the automotive industry to improve the formability, corrosion resistance and paint appearance of the sheet metals. Although a coating would not directly alter the mechanical behavior of steels, its frictional behavior may change the strain path during the forming process and ultimately affect the formability of the sheet metal. Therefore, in contact modeling, it is necessary to take the micro-mechanical behavior of the coating into account. Though there are various analytical and semi-analytical contact models available for uncoated sheet, very few have been extended to coated sheet (12)(13)(14)(15). The major challenge in modeling the contact mechanics of the coated system is to capture the contribution from coating and substrate for the given load. This is done by defining the effective hardness of the coated system for which there are many empirical relations developed based on the hardness tests on coated sheets [12,[16][17][18][19]. For example, empirical relations proposed by Halling et al. [19] are function of coating and substrate hardness, coating thickness and radius of the asperity. However, hardness results are derived from a typical indentation tests with spherical indenter on polished coated sheet. Moreover, the hardness result from the test depends on various parameters such as the type and size of indenter, load and indentation depth. Nevertheless, determination of effective hardness makes it possible for coated systems to adopt the model developed for uncoated systems. Chang et al. [13] extended an elastic-plastic contact model [6] to a coated surface. Chang et al. [13] proposed a hardness relation for plastically deforming asperity with surface coating. However, the model is limited to spherical asperity which is far from the asperity geometry found on the rough surface. Also, the model does not account for strain hardening. Furthermore, based on plane strain upper bound method, Mulki and Mizuno [14] proposed an ideal plastic contact model for GI sheet. However, the model makes simplifying assumptions about the material behavior, asperity shape and their interaction. Also, the model is not suitable for a general 3D rough surface with asperities of random shapes and distribution. Chen et al. [12] has proposed load independent effective hardness relations for a coated sphere. They used FE analyses to develop the relations for different coating thicknesses under the realistic contact situation of a smooth rigid tool flattening the spherical asperities. The proposed relations are more reliable as it also considers the material hardening but is limited to spherical asperity. Recently, Shisode et al. [15] proposed a normal load flattening model for coated sheet which is an extension of the Westeneng's model. The model allows to input the actual measured height data for coating and substrate material. The model also considers the material hardening effects. However, the limitations of Westeneng's model still exists in this model.
The goal of this study is to develop a normal load contact model for a generalized rough surface which describes more accurate mechanical behavior of the asperity and effect of coating-substrate material combination. For this purpose, the model uses a series of FE analyses of smooth ellipsoidal asperities of various sizes and aspect ratios to determine the contact pressure of a single asperity. Later, this database of contact pressure is used to determine the contact pressure of individual asperities on the rough surface. In this way, the deformation of the rough surface and thus real area of contact can be calculated.
The paper is organized as follows: In section 2, a semi-analytical model is proposed for predicting real area of contact, followed by its implementation procedure. In the next section, experimental results of the normal load flattening experiments are presented for samples of different steel grades, zinc coating thickness and surface roughness. Finally, the model results are compared with the experimental results followed by a discussion and conclusions.

Flattening model
The resistance of the asperity on the rough surface against deformation depends on its material, shape and size. Hence it is necessary to include the effects of asperity geometry on its flattening behavior. One way to account for this shape and size dependency is to determine the contact pressure of each asperity under a flattening condition by a virtual test using the FE method. However, the practical challenge with this approach is the necessity to perform large number of simulations with various sizes and shapes of asperities. Chen et al. [12] considered a spherical asperity and determined its effective hardness or contact pressure under a normal load flattening. In the current model, asperities on the rough surface are approximated by an ellipsoid which provides the flexibility to better approximate the shape of asperities. FE simulations are performed on the ellipsoidal geometry by varying its size and aspect ratio (height to radius ratio) in a certain range. Contact pressure of each asperity under the flattening situation is determined and the contact pressure database is built as a function of the geometrical parameters of the asperity. The contact pressure is determined as the ratio of load carried by asperity and its contact area. Due to material hardening and evolution of the asperity shape with flattening, the contact pressure evolves with upsetting height. Hence the contact pressure is determined at different levels of flattening.
Key assumptions 1. Tool is assumed to be rigid and smooth. These are routine assumptions [6,8] in normal load flattening models. In sheet metal forming processes, tool hardness is normally much higher than the workpiece hardness, moreover, tool surface roughness is generally one order lower than the sheet metal roughness. To accord with these assumptions, in this study, a sapphire notch is used as a tool in the experiments. 2. In the analytical flattening model, displaced material from the deformed asperities is redistributed uniformly in the noncontacting area [8,20]. 3. Frictionless slip contact condition is assumed between a single asperity and the rigid tool. This is on the basis that during normal loading of the surface, there is a limited sliding between tool and a single asperity and therefore friction is not playing a significant role. 4. In FE analyses, the elasto-plastic response of the material is described by a J 2 isotropic plasticity theory. 5. Coating thickness is assumed uniform and perfectly bonded to substrate. This is also confirmed by measuring coating thickness at different location on the sheet surface using Fischer FMP40 Dualscope coating thickness measurement probe. The coating thickness varies by � �1.5 μm.

Single asperity FE model
A contact pressure database for different materials and coating thicknesses is built for ellipsoids of different sizes and aspect ratios. The base radius and height are varied from 10 to 100 μm and 0.8-5 μm respectively. The range is decided by carefully observing the dimensions of the fitted ellipsoids on the samples used in this study at different tool planes. The base radii of 10, 20, 50, 70 and 100 μm and heights of 0.8, 1.5, 3 and 5 μm are used to form 20 axisymmetric simulations for each substrate material and coating thickness combination corresponding to the samples used in the experimental validation (see Table 2).

Material properties
True stress vs. true strain flow curves are available from uniaxial tensile tests of all steel substrate materials as shown in Fig. 1. Flow curves from the tensile tests are extrapolated for higher values of true strain using biaxial tensile test data [21]. Material behavior of the zinc coating (see equation (1)) derived using nano-indentation tests on various zinc grains [22] is used in this study. There, the elasto-plastic properties from the load-indentation response are derived for the zinc coating using an analytical algorithm by Dao et al. [23].
where, σ y is the flow stress, ε p is the plastic strain, n ¼ 0.14, E is the Young's modulus and σ 0 is the initial yield strength. The Young's modulus, Poisson's ratio and initial yield strength of zinc coating [22,24] and steel substrates corresponding to materials used in the experiments are shown in Table 1. Fig. 2 shows the models of asperity for uncoated and coated specimens. For coated samples, the coating is modeled as a mean coating thickness. A commercial FE solver (MSC Marc 2016) was used for nonlinear, elastic-plastic, axisymmetric FE simulations of ellipsoidal geometries. In the simulations, axisymmetric 4 node quadrilateral elements with full integration and constant dilatation are used. The FE mesh is divided into different mesh density zones. A finer mesh density is used in the asperity compared to the other parts of the model. Sufficient model extent in the lateral direction and in the thickness direction (> 3r) is considered. The bottom of the model is fixed in the vertical direction and the outer edges are free to expand. Material flow curves shown in Fig. 1 are used. The tool is moved incrementally downward and at each increment, the results are post-processed to determine the contact force, contact length and the asperity flattened volume. From the contact length, a contact area is determined. The contact pressure is determined as the ratio of contact force and contact area and stored as a function of percentage of flattened volume with respect to total asperity volume. Once the contact pressure database is constructed for a specific set of coating-substrate material combination and coating thickness, it can be used for modeling the flattening behavior of different surface textures and normal loads as long as the materials and mean coating thickness are unchanged. Fig. 3 shows some of the contact pressure results from the simulations as a function of percentage of flattened asperity volume for sample number 1, 2 and 4 corresponding to a coating-substrate material combination and coating thickness as per Table 2. Hardness curves derived from FE analyses are smoothed out and used in the analytical contact model (section 2.3). It is clearly seen that the contact pressure of an uncoated asperity is higher compared to the coated asperity and with increase in coating thickness the contact pressure decreases.

FE modeling and post-processing
For the coated asperities with higher aspect ratio (at base radius of 10 and 20 μm and height of 3 μm), the contact pressure initially increases however at larger indentations decreases or remains almost constant afterwards. For the coated asperities of 10 μm base radius and 3 μm height, with the increase in tool displacement the rate at which contact area increases is higher compared to the rate of increase in the contact force which results in decrease in the contact pressure. In the case of asperities with a 20 μm of base radius, 3 μm height, the rise in force is approximately equal to the increase in contact area. This leads to almost constant contact pressure. This trend is reverse for asperities with larger base radii. This can be explained by the deformation of the steel substrate in asperities with larger base radius. Fig. 4 depicts the plastic deformation of asperities with base radii of 10 and 50 μm with height of 3 μm. It is found that, for asperities of larger base radii, the plastic zone can reach to the substrate for a given tool displacement. For smaller asperities, it does not reach the substrate due to lower contact load for the same tool displacement. Hence for asperities with larger base radius, the contribution of the substrate to the contact pressure of coatingsubstrate hybrid is larger.

Multiple asperity contact model
Asperities in contact with a flat tool on rough surface at an assumed tool displacement are determined. All the contacting asperities are fitted by an ellipsoidal shape and the geometrical parameters of the ellipsoid are determined for each asperity. Material behavior of each individual contacting asperity is modeled as the contact pressure of its equivalent ellipsoid derived from the contact pressure database. It should be noted   that this is an approximation to the asperity but still the derived shape and size dependent contact pressure is expected to result in a better approximation of the asperity hardness in the flattening situation as compared to the constant relation proposed by Tabor [9]. Force balance and volume conservation equations are solved to determine the equilibrium position of the tool from which the real area of contact is determined.
A measured rough surface with a size equal to the representative area of 2 � 2 mm is used in the model to determine a real area of contact. A strategy to determine the representative area is explained in section 3.1. The rough surface is described as the measured surface points or pixels using confocal microscopy. This deterministic description of the surface is used to determine the real area of contact. An increment (e.g. 5 � 10 3 μm) in tool displacement is assumed. There are certain number of surface points or pixels that come in contact at an assumed position of the tool. The surface points which are in contact and adjacent, form a cluster referred to as contact patch as shown in Fig. 5. A contact patch is a cluster of pixels connected together by shared pixel edges or by shared pixel corner points. A procedure described by de Rooij [25] is followed to determine contact patches and the geometry of ellipsoids. The contact patches are identified by a binary image processing technique. Once the contact patches are determined, the base radius r of the fitted ellipsoid is determined from the area A of the contact patch. Knowing the base area A, the height h of the ellipsoid is determined by equating the volume of the contact patch (V) and half the ellipsoid using equations (2) and (3).
It is expected that the asperity contact pressure is influenced by noncontacting material in its vicinity. However, if the ellipsoid is fitted at a tool position as shown in Fig. 5 (right) with base area A i , the contact pressure of the asperity can not account for the adjacent material properly. The adjacency effect can be included by considering additional material adjacent to the contact patch while fitting an ellipsoid. However, considering too much adjacent material in the fitting stage can lead to the deviation in fitted geometry from the actual flattened asperity height data especially at the tip of the asperity. It is also expected that the effect of the adjacent material close to the flattened volume is most significant and reduces with increasing distance from the flattened material. A sensitivity analysis on a single asperity FE simulation with different amount of adjacent material showed that the influence of adjacent material on contact pressure corresponding to a base area � 2A i is a good approximation. Hence it is assumed that the contact pressure at a specific deformation level for the resulting contact patch (with contact area A i ) is better represented by an ellipsoid fitted on the height data with base area equal to two times the contact patch area (2A i ). The fitted ellipsoid on the asperity with base area of 2A i and contact area of A i is used to determine the contact pressure of asperity from the FE database. The schematic illustration of this strategy is shown in Fig. 5 (right). This ensures the effect of adjacent material on asperity contact pressure with little effect on asperity fitting accuracy.

Analytical contact model
To determine the fractional real area of contact α at a given nominal pressure P nom , the tool is moved at small increments on the measured rough surface such that the separation between tool and rough surface decreases. Asperities in contact are determined. Force balance and volume conservation equations are solved to determine α.

Force balance equation
The total force F int carried by the contacting asperities is determined as follows, The total force applied by the tool F ext on the rough surface is, At equilibrium position of the tool, where, H i is the contact pressure of each contacting asperity derived from the database based on asperity geometry, A i is the asperity contact area at the given tool position, N is the total number of asperities in contact and A nom is the rough surface nominal area.

Volume conservation
When the asperities deform, the displaced volume is distributed over the non-contacting area. This is taken into account in the contact pressure calculation stage using FE simulations but also required to be included in the analytical model. In the analytical model, it is assumed that, the displaced volume is uniformly distributed over the non- contacting surface [8,20] as shown in Fig. 6. Due to the rise in non-contacting surface, some of the non-contacting surface points come into contact with the tool which increases the real area of contact. The uniform rise U in non-contacting surface points for a given tool position d can be determined by equating the flattened volume (left side of equation (7)) to the rising volume (right side of equation (7)). Where φðzÞ is the normalized height distribution data of the measured surface (see section 3.1 for more details on φðzÞ).
The contact pressure value for each of the contacting asperities is extracted from the contact pressure database generated by running offline FE simulations. At each increment, the bearing force contribution from each of the contacting asperities is determined using the contact pressure and the contact area at the asperity level (equation (4)). Force balance (6) and volume conservation (7) equations are solved to determine the real area of contact.

Solution scheme
A step by step solution scheme to determine real area of contact is summarized here.
Step 1: Start with tool displacement increment d i . Determine the number of contact patches N at tool plane d i on the confocal image. Determine the base area A i and volume V i for each contact patch. Where, V i of the contact patch is determined as total volume of the surface points above the tool plane at d i and A i is the base area at tool plane d i .
Step 2: Asperity fitting and contact pressure calculation for each contact patch -Determine total fractional real area of contact α i ¼ ð -Determine the position of an imaginary plane d' i such that fractional real area of contact is equal to 2α i . d' i is determined by solving equation (9).
-Again determine the contact patches, A' i and volume V' i for each contact patch at plane d' i . Using equations (2) and (3)  Step 3: Now for assumed tool plane at d i , determine uniform rise in non-contacting surface points U by solving equation (7). Use U and update the base area A i for all contact patches identified at assumed tool plane d i . This is because the base area in some contact patches may increase.
Step 4: Determine the force carried by each contact patch F i ¼ H i A i and furthermore, the total force F int using equation (4).
Step 5: Check the relative force unbalance ratio ¼ jðF int F ext Þj= F ext . If the ratio is less than the specified tolerance of 2%, then equilibrium position of the tool is d i and the resulting fractional real area of contact α is determined using equation (8). Otherwise add next displacement increment in tool and repeat the process until the relative force unbalance ratio is less than the tolerance.

Experiments and model validation
Normal load flattening experiments are performed on coated and uncoated samples of steel with different substrate grades. The details of the samples are shown in Table 2. All the coated samples are zinc coated (GI) and electro discharged textured (EDT) with different coating thickness. Coating thickness is measured at different locations on the sample. The coating thickness measurement is done using Fischer FMP40 Dualscope coating thickness measurement probe which can measure the thickness of non-ferrous coating with an accuracy of � 1 μm. A nominal contact pressure up to 60 MPa is applied. A Sensofar confocal microscope is used to measure the surface height data with a pixel resolution of 1.3 μm. The surface height data of the sample is measured before and after loading from which the surface height distribution is derived. The height distribution function of sample 1 and 2 is bimodal with two dominant peaks (see Figs. 9 and 10) while sample 3 and 4 have a single dominant peak (see Figs. 11 and 12). This is another difference in the samples in addition to the coating thickness, substrate material and surface roughness.

Height distribution measurement
The measured data using confocal microscopy is corrected for any missing points, unrealistic sharp peaks in the measurement which is often an error introduced during the measurement. The correction is done using the same procedure as described in (1). A noise filter is used to remove sharp peaks from the measured data. A median filter is commonly used in the image or signal processing to remove noise. A median filter with kernel size of (3 � 3) is used to suppress the sharp peaks. The other type of correction required on the raw data is to remove the curvature or the tilt on the measured data. For this purpose, a polynomial function of 2nd order is fitted through the raw data. The raw data is shifted and straightened by subtracting it from the fitted function.
Furthermore, it is necessary to measure sufficient area of the rough surface which is representative of the height distribution of the whole surface. On the contrary, bigger measurement area increases overall computational time. Convergence of the height distribution data is used as a criterion to determine the size of representative area. For this purpose, height distribution density functions obtained from different sizes of the measurement areas are compared. The lowest measurement area at which the distribution smoothly converges is chosen as the representative area. Fig. 7 (left) shows the result of this study with measurement size varies from 0.5 � 0.5 mm up to 2.5 � 2.5 mm at a resolution of 1.3 μm. It can be concluded from the results that a measurement area of 2 � 2 mm can be chosen as the representative area. Therefore, all the measurements are carried out on a measurement area of 2 � 2 mm.

Stochastic surface description
In order to estimate the real area of contact, the rough surface is often modeled stochastically using a polynomial function or a normal distribution function. However, it is also observed that the height distribution function φðzÞ can be non-normal or asymmetric. For example, EDT sheet surfaces generally have a bi-modal height distribution as shown in Fig. 7. Such surfaces require flexible yet robust fitting function. For example, height distribution data can be fitted using Fourier series fit (1) or a B-spline fit [26]. The stochastic distribution facilitates an efficient translation from micro-scale contact modeling to macro-scale modeling of contact. The downside is that it eliminates the deterministic description of the surface texture and looses the adjacency information (spatial distribution of asperities) which is vital in the contact modeling. But, a stochastic description can still be used in the equations of the solution scheme where the adjacency information is not required. In this study, a B-spline function is used to fit the height distribution data which will be used in the modeling. The accuracy of the fit depends on the order and number of splines. Fig. 7 (right) shows the height distribution data and its corresponding B-spline fit using 10 cubic splines. It can be concluded that the chosen fitting function is a good choice to represent the height distribution data. The key components of the setup are a pair of polished rectangular sapphire notches (surface roughness Sa <1 nm), load cell, thrust bearing and fasteners. A sapphire notch can be considered rigid because its hardness is much higher compared to the sample material. A sheet sample of 5 mm width is placed between the pair of notches and loaded at the center. The loading notches have a total rectangular contact area of 5 � 7 mm with chamfers at the edges. The longer side (7 mm) ensures that the notches always cover the full width of the sample. A self lubricating bushing is used to set a free sliding fit between the notch holder and the load cell housing. The notch holder freely slides in the load cell housing which ensures that the whole compressive load is delivered to the load cell which is sitting at the bottom of the housing. A miniature button load cell of 2 kN is used. A thrust bearing is used for better alignment of sapphire notch against the specimen during loading. The sample is marked at the center with a square of � 2.5 � 2.5 mm to compare the asperity height data before and after loading. The bolts are tightened one by one in small increments while the applied load is continuously read during the loading until the total load reaches a magnitude corresponding to the desired nominal pressure.

Results
In this work, the real area of contact is determined using a height distribution data of the deformed rough surface [1,8]. The height distribution data of deformed and undeformed sample is overlaid to visualize the difference. When the load is applied, the asperities deform and the height distribution profile shifts towards the left giving rise to the local maximum value (see Fig. 8 right). In an ideal situation, when a flat rigid surface is deforming a soft rough surface all the taller asperities will attain the same transition height h t . However, due to elastic recovery and misalignment of any small degree the height distribution function does not show a perfect transition height (vertical line at the local maxima) but rather exhibits a steep slope. In this study, asperity height corresponding to the local maxima is used as the transition height and the real area of contact α h is determined as follows, where h t is the height of the asperity corresponding to local maximum peak and φ d ðzÞ is the normalized height distribution data of the deformed surface. As mentioned in Table 2 The 3D surface topography of the undeformed and deformed surfaces clearly show the flattening of the peak asperities with increasing applied pressure. The transition height h t at local maxima is easily recognizable which is used to determine real area of contact using equation (10). The results of the measurements for all samples are shown in Fig. 13.
The fractional real area of contact for all samples increases almost linearly with increase in nominal pressure. The real area of contact for uncoated material (sample 4) is smaller compared to the coated samples. This difference clearly demonstrates the effect of a coating on asperity deformation. The difference in real area of contact for differently coated samples at low pressure is not clearly seen but for higher pressures the difference is clearly visible. The real area of contact for a higher coating thickness is higher. This can be attributed to reduction in effective hardness of the material with increase in zinc coating thickness. Sample 2 with average coating thickness of 14 μm and sample 3 with average coating thickness of 23 μm show a similar results. This is due to the stronger substrate in sample 2 (S220) compared to sample 3 (DX56) (see Fig. 1). The key conclusion from the results is that, the fractional real area of contact is proportional to applied load and that the coating material, its thickness and underlying substrate have influence on the deformation of the asperities.

Model results and validation
The model is implemented in Matlab. Using the contact pressure database derived for each of the material combination and measured sample height data, the fractional real area of contact is determined at   Fig. 13. The results show that model is able to predict the real area of contact for the nominal contact pressure range upto 60 MPa very well. A real area of contact increases almost linearly with rise in nominal pressure. This trend is consistent with the model results by Chang et al. [13]. This is due to the fact that the rise in asperity contact pressure after initial yielding is limited for all asperity geometries and materials (see Fig. 3). Furthermore, at each tool position, asperities of different sizes (hence different contact pressure) come into contact. Hence, the resulting effective contact pressure is an area averaged quantity of all the contacting asperities and depends on the distribution of asperity geometry. For example, Fig. 14 shows the distribution of asperity base area and cut-off volume for sample 1 at equilibrium tool positions corresponding to 15 MPa and 45 MPa. The number of asperities and their size increase with nominal pressure but there is always a contribution from a large range of asperities at each pressure level. The model is able to capture the trend and the magnitude of the fractional real area of contact as determined from the experiment. Coupling the single asperity FE simulation results with the analytical contact model enables to account for the coating-substrate material combination, material hardening, neighboring not-contacting material effects and realistic asperity hardness.

Conclusions
A semi-analytical normal load flattening model is proposed which accounts for the coating and substrate material behavior, coating thickness and surface topography of the rough surface. An asperity deformation is modeled more accurately using the contact pressure under the flattening situation derived for ellipsoidal asperities of different sizes compared to previous studies using a constant relation between yield strength and hardness [8]. The model is able to predict the influence of coating thickness, surface topography coating-substrate     material combination on real area of contact. Furthermore, the model does not require any fitting parameter. The rise of non-contacting material is included using volume conservation. The model predicts a monotonic increase in real area of contact for coated and uncoated systems.
Normal load flattening experiments are performed on uncoated and coated samples with different materials and coating thicknesses. For the samples analyzed, the real area of contact increases approximately linearly relative to the applied load. The measured and calculated real area of contact for coated sheet is larger compared to the uncoated sheet. The real area of contact increases with increase in coating thickness. The predictions from the proposed model show a good agreement with the experimental results for uncoated and coated sheet samples.

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.