Three-dimensional finite element analysis of shallow indentation of rough strain-hardening surface

: Three-dimensional finite element modeling of the contact between a rigid spherical indenter and a rough surface is presented when considering both the loading and unloading phases. The relationships among the indentation load, displacement, contact area, and mean contact pressure for both loading and unloading are established through a curve fitting using sigmoid logistic and power law functions. The contact load is proportional to the contact area, and the mean contact pressure is related to the characteristic stress, which is dependent on the material properties. The residual displacement is proportional to the maximum indentation displacement. A proportional relationship also exists for plastically dissipated energy and work conducted during loading. The surface roughness results in an effective elastic modulus calculated from an initial unloading stiffness several times larger than the true value of elastic modulus. Nonetheless, the calculated modulus under a shallow spherical indentation can still be applied for a relative comparison.


Introduction
There has been a plethora of research on the surface roughness and its dependence on the machining parameters [1,2] owing to its influence on the contact heat transfer [3], wear [4], adhesion, stiction, electrical conductivity of the interface, and surface functions (e.g., coating performance, frictional behavior, and fluid load capacity [5]), and because it plays a key role in the physical properties of micro/nano structures. Adhesion between particles and surfaces, which is closely associated with various technological processes (e.g., pharmaceutical processes and powder painting), depends on the surface roughness because it reduces the true contact area leading to reduced body interaction [6]. Surface roughness also plays a significant role in the dispersion forces, which can generate severe problems such as spontaneous stiction or permanent adhesion between separate elements (e.g., suspended structures with a very small gap distance) in micro-electromechanical systems [7]. Lyon et al. [8] theoretically studied the effects of roughness at the atomic scale on the surface plasmon excitation, and found visible effects of the surface roughness on the image potential and stopping power. Through a molecular dynamics simulation, Liu et al. [9] found that the surface roughness plays a significant role in the plasticity initiation of silicon nanowires. Nunez and Polycarpou [10] experimentally investigated the effects of the counterpart surface roughness on the formation of a transfer layer from the polymer films for application to dry or solid lubrication, and determined the dependence of the friction coefficient and wear rate on the surface roughness. Zhu et al. [11] studied the influence of the surface roughness on wheel-rail adhesion and wear through the use of a rollingsliding tribometer. They found that increasing the surface roughness will increase the wear, and that the surface roughness influences the adhesion recovery. Dawood et al. [12] studied the effects of the surface roughness during the friction stir welding of 6061 aluminum alloy workpieces. They found perceptible influences of the surface roughness of a workpiece on the quality of the weld surface, as well as the heat flux, grain refinement, tensile properties, micro-hardness, and fracture of the joints. Curry et al. [13] found that the structure, failure, and thermal conductivity of a columnar coating generated by suspension plasma spaying are strongly influenced by the topography of the surface onto which the coating is deposited.
The important effects of the surface roughness have motivated extensive studies on the contact between the rough surfaces of fractal or sinusoidal features [14−19] (e.g., frictionless normal contact between rough surfaces [20], elastic frictionless non-adhesive contact between rough surfaces [21−23], plastic deformation of rough rolling contact [24], contact stiffnesses of rough surfaces [25], and plastic contact with/without strain hardening [26−28]). Greenwood and Williamson [29] assumed a Gaussian distribution in the heights of spherical asperities through the application of Hertz theory [30] to each asperity. In addition, variations in radius, height, and ellipticity were considered by Bush et al. [31]. A scaling approach was used by Persson [32] as a way to consider the roughness of successive length scales owing to the self-affine fractal characteristic of the surfaces. Because the contact pressure is often sufficiently high to induce plastic deformation, statistical models considering elastic-plastic asperities have also been developed [33−35]. These theoretical models ignore the interaction between neighboring asperities as well as the bulk deformation of the base material, however, and suffer from many shortcomings (e.g., a discontinuous elastic-plastic transition [34]). Because the actual surface topography, interactions among asperities, and the plasticity can be considered, the finite element method is a deterministic approach to investigating the elastic-plastic contact of rough surfaces. Pei et al. [36] adopted a fully three-dimensional (3D) finite element (FE) modeling for contact of an elastic-plastic rough surface and a rigid flat surface, and found that the results for elastic-plastic solids were in contrast to previous studies on elastic solids [37]. Poulios and Kilt [38] studied the frictionless contact of nominally flat rough surfaces by applying a 3D FE model with a real surface topography and elasticplastic material behavior. Dong and Cao [39] simulated the deformation behavior of the elastic-linear plastic asperity of a sinusoidal profile using a 3D FE model, and correlated the subsurface stress conditions with a fatigue layer. Liu et al. [40] studied the normal contact of elastic asperities with an elliptical contact area, and investigated the effects of the geometrical and material parameters on contact stiffness. This work is devoted to addressing the contact mechanism of a rough surface through the finite element modeling of a spherical indentation of an experimentally measured rough surface. The contact variables (e.g., normal load and contact area) are quantified, and explicit relations are presented based on the simulation results.

Finite element modeling of rough surface contact
The material used herein is gold, and the surface topography is extracted from atomic force microscopy (AFM). Figure 1 shows an AFM image of a gold surface with a scan size 5 μm × 5 μm, data array of 512 × 512, and surface resolution Δ = 9.8 nm, which is sufficient for modeling the contact behavior because a high resolution of about tens of nanometers is required to properly represent the surface topography [41]. The data are visualized and analyzed using Gwyddion 2.41: the arithmetical mean deviation R a ≈ 6.71 nm, the root-mean-square (rms) roughness R q ≈ 8.21 nm, the skewness R sk = -0.0285, the kurtosis R ku = -0.436; and the rms slope      2 | / | 0.94 z x (the brackets indicate the average of all values, where z is the surface height). The average peak-to-valley and curvature radius of the asperities are 35 nm and 240 nm, respectively [42].
The one-dimensional height-height correlation function in this study is defined as where N = M = 512 (the same as a data array of 512 × 512 for the AFM scan), and  x is the distance from point (n + m, l) to point (n, l). As expected,  x (H is the Hurst or roughness exponent) when  x is smaller than the lateral length (≈ 10R q for the case shown in Fig. 2   with the reasonable range of 1 < D < 2 [46]. Figure 3 shows a 3D FE mesh with half a million nodes for the sample to be indented. A fine mesh is used on the top surface, and the mesh becomes increasingly coarser away from the surface. The 3D FE simulation is carried out using the finite element package Z-set, applying its parallel solver [47−49]. The FE mesh is constructed using the AFM data shown in Fig. 1. The in-plane dimensions are 2.5 μm × 2.5 μm because only the central region of the scanned image is used for generating a rough surface in finite element modeling. The thickness of the indented sample is about 2.5 μm, which is sufficiently high to represent a bulk material deformation under the displacements considered (≤ 20 nm). The smallest element at the surface of the mesh is 9.8 nm, which is the same as the surface resolution in an AFM scan. Linear elements are used because they were found to be more accurate owing to a greater consistency in terms of the lumped mass matrix [50], although Mesarovic and Fleck [51] found that second-order elements show better convergence and accuracy than linear elements. The shapes of the nano-scaled asperities cannot be determined because the radii of the AFM tip and peaks are comparable. A 3D finite element mesh with a rough surface is generated from a 3D finite element mesh with a flat surface. Each node on the top surface corresponds to one measurement in the AFM map, and the surface roughness is generated by displacing each node on the surface in the direction normal to the top surface based on the height from this map. The finite element model consists of approximately half a million nodes. The contact area predicted by the finite element method is closely related to the number of nodes in contact, and thus depends strongly on the finite element mesh [48]. A reasonable mesh convergence requires 32 elements to describe the shape of one single asperity [48]. However, the difference between various meshes is  | https://mc03.manuscriptcentral.com/friction more sensitive when the displacements or forces are small [48]. The force is overestimated by a few percentages for the present mesh density compared with a denser mesh with an element size half of the present size. The present mesh is used for computational efficiency. To ensure good precision, the increment of the contact interference is set to 1 nm. A rigid spherical indenter with a radius of 50 μm is fixed above the center of the indented surface. The sample is brought into contact with the indenter by applying a uniform displacement on the bottom surface of the sample. Zero in-plane displacements are applied to the bottom surface. The symmetric boundary conditions are applied to the lateral faces. In the initial state, the highest point of the rough indented surface and the lowest point of the spherical indenter both have the same coordinate of z = 0, and no contact occurs between the indenter and sample.
A frictionless (i.e., perfectly slip) contact condition, which corresponds to a fully lubricated contact condition in a real experiment, is typically adopted in a contact analysis [51−58] because no discernable differences between frictionless and frictional contacts have been found [59] in spite of the difference in the stress field near the contact surface [60] between a frictionless contact and a fully stick contact [61]. The classical non-penetration/non-adhesion contact condition with a rigid surface is used in the following form [62]: where      n n n is the normal stress on the contact surface, n g is the normal gap between the deformable solid and the rigid plane, and  t is the tangential stress.
The contact algorithm enforces an impenetrability constraint on the contact surfaces. Strain hardening is considered because very few materials exhibit an elastic-perfectly plastic behavior [63]. Isotropic rate-independent J 2 plasticity is used for the constitutive equation, where ij s is the deviatoric part of the stress tensor,  is the effective Mises stress, and  is an internal hardening variable. Because a generalized solution applicable to all types of materials has yet to be developed [64], and the contact responses differ from one hardening law to another [51], three different strain hardening laws are considered: linear, exponential, and power hardening: where  p and  p are respectively the effective plastic strain and effective plastic strain rate where  ij is the plastic strain rate tensor.
The material properties are listed in Table 1. The material used for the linear strain-hardening law is gold (Au) owing to its linear relationship between load and reduction area, and the similarity between the load-reduction area and stress-strain curves. Linear strain hardening is expressed in terms of the tangent modulus, which is the slope of the uni-axial stressstrain curve beyond the elastic limit. The poisson's ratio of Au is from [41] and [65]. Two different sets of material parameters are used for linear strain hardening: Au-1 with a small strain hardening, and Au-2 with a large strain hardening, which is beneficial in reducing the friction and wear [66] from a larger stiffness. The tangent modulus K = 1.07 GPa for Au-1 is less than 0.02E, which is the upper limit of many practical materials [67,68], and Au-1 can be considered an elastic-perfectly plastic material because isotropic linear hardening renders the same results as elasticperfectly plastic behavior [69], and a tangential modulus of up to 10% of the elastic modulus has only a small effect on the frictionless and non-adhesive contact [70]. The material used for the exponential hardening law is Table 1 Mechanical properties of elastic-plastic materials.  [48]. The elastic modulus and poisson's ratio for Fe are from [71]. The material used for the power hardening material is copper, the properties of which are from [58].

Results and discussion
Uncertainty exists in the first contact point, and there is difficulty in capturing the first contact. The first appearance of a detectable contact area is regarded as corresponding to zero indentation displacement (h = 0). The indentation displacement is the displacement of the bottom surface of the bulk sample minus this initial displacement corresponding to the smallest detectable reaction force. Only the loading process has been previously studied [48, 72−74]. Both loading and unloading are considered in the present work, similar to the analysis of the indentation cycle in [26]. contact with each other is quite small), and the area is strongly dependent on the specific realization of the random surface [36,37] as well as the specific material properties. A power function can be used to express the relation between the contact area and indentation displacement, namely,  2 1 D A D d , with parameters D 1 and D 2 being dependent on the material properties. A difference in the A-d curves for different materials occurs because the plasticity can affect the distribution of the contact area [36].

Relation between contact area and indentation displacement
The unloading A versus d curves for different materials can be expressed using a sigmoid logistic type function (see Figs. 4(b)−4(e)): where the subscript "max" indicates the value at the final loading (i.e., the maximum indentation displacement, d max ), and D 3 , D 4 , and D 5 are dependent on both the material properties and d max . At the initial stage of unloading (d > 0.98d max ), the contact area remains unchanged with a decrease in the indentation displacement. The presence of residual indentation displacement at the final unloading is due to the occurrence of plastic deformation. Other types of sigmoid logistic functions are also applicable for the unloading curve fitting (see the example in Fig. 4(f)).   where  r is the characteristic strain, and   r 0.1 for Au-1,   r 0.05 for Au-2, and   r 0.17 for Cu. In addition,  r for Au-2 is higher than that for Au-1

Relation between mean contact pressure and contact area
because Au-2 is stiffer owing to a larger tangent plasticity modulus, which is consistent with a larger contact pressure for a stiffer material [66]. The criterion used for choosing the characteristic strains for different materials is that, under a large contact area, the normalized stresses for different materials become almost the same, as shown in Fig. 5(a). A material with exponential strain hardening is used as a reference because it has no characteristic strain. Changing the strain-hardening rate is an effective way to modify the deformation resistance capability of a material under elastic-plastic deformation, and the characteristic stress is an imprint of the material stiffness. As expected, the mean contact pressure is sufficiently large to produce plastic deformation [36]. If asperities constituting a rough surface are approximated using sinusoidal shapes, the aspect ratio (i.e., the ratio of height to width) is about 0.05 for most asperities [48]. Under such a small aspect ratio (i.e., short-wide type), the plastic strain reaches a stable platform for contact of an isolated asperity [39]. The plateau of the mean contact pressure for contact of a rough surface is in line with the plateau of the contact stress and the plastic strain on the peak of each individual asperity [39]. For a material with small strain hardening (Au-1, Fe, and Cu), the contact pressure reaches approximately 3 times the characteristic stress very quickly after initial contact occurs, and a plateau of the mean contact pressure can be approximated during loading. For a material with large strain hardening (Au-2), p m increases gradually during the initial stage of contact, and approaches the material hardness under sufficient contact, although strictly speaking, there exists a small increasing tendency of p m owing to interactions of neighboring asperities, in line with the findings in Refs. [38,41,52,72]. Many models of plastic contact of a rough surface assume that p m corresponds to the material hardness, H, which is proportional to the yield stress [75]. Contact mechanics predicts   y / 3 m p for a simple isolated asperity [76]. However, neighboring asperities can decrease the deviatoric components of the stress tensor, and tend to increase  y / m p ; thus, the interactions between neighboring asperities allow   m y 3 p [36]. Our results suggest that, under well developed plastic deformation, the hardness or mean contact pressure of a strain hardening material is better for associating with a characteristic stress, which is a function of the yield stress and a characteristic strain, and thus    m r

. H p
The unloading A versus p m curves for different materials can also be expressed through asigmoid logistic type function (see Figs. 5(b)−5(e)): where 2 a and 3 a are dependent on the material properties, but independent of max 0 / A A . The value of 1 a is directly related to the contact area, and  1 a 0 100( / ) A A (see Fig. 5(f)); in addition, note that the unit of measurement for a dimensionless contact area 0 / A A is the percentage (%). At the initial stage of unloading, the contact area remains unvaried, although the mean contact pressure continues decreasing while unloading.

Relation between contact area and load
Both the area and geometry of the contact regions affect the interfacial stiffness and area-dependent properties (e.g., adhesion, contact stiffness, and electrical and thermal conductivity [36,75,77]). The contact area topology at different loading steps is displayed in Fig. 6, which shows the expansion of the equivalent plastic strain on the contact surface of Au-1 with increasing indentation displacement. Plastic deformation is an important factor in predicting the generation of debris caused by residual stress and micro-cracking [66]. As expected, the contact occurs at a finite number of discrete and isolated spots, leading to a complex contact morphology [36] because the entire contact region is composed of many tiny discrete regions [73]; the asperities experience plastic deformation as soon as contact occurs [41], demonstrating the invalidity of an elastic model in a real application because plastic deformation occurs even at an indentation displacement of as small as 1 nm owing to the fine roughness. Experiments showed that the interference under a fully plastic contact state is much smaller than the characteristic size (e.g., about 0.1% of the radius of the ball for contact between a rigid ball and a deformable space) [78], which is very difficult, if not impossible, to detect. As expected, because the rms slope      2 0.5 | / | 0.94 z x is much lar- is the reduced elastic modulus), plastic deformation is induced at very small loads as soon as peaks of the asperities touch the indenter [36]. Neither the first contact nor the maximum contact pressure occurs at the surface center of the indented sample owing to the surface roughness. Figure 7(a) shows a variation in the true contact area A with indentation load F during loading. Here, F is normalized by . The contact area is proportional to the contact load, which is in line with previous studies [31, 32, 36−38, 72, 73, 79−81], because the area increases with the load to maintain a nearly constant contact pressure. For a rigid-perfectly plastic material, the real contact area is proportional to the load because the mean contact pressure is constant and equals the flow stress [72]. Similarly, for a strain hardening material, the mean contact pressure is also constant, and is equal to the characteristic stress. The ratio of the load to the projected area is 0.92 GPa for gold thin films under a spherical indentation [42], which is close to H (≈ 1.3 GPa) for Au-1. It can be approximated that for the loading, and the effect of the plastic constitutive law on the relation between the contact area and normal load during loading is based on the characteristic stress  r .
The real contact area is a highly significant quantity because it dominates the creep curve of a rough surface [82], controls the electrical contact resistance [41], and provides insight into the adhesion-induced pull-off force [83]. The friction from adhesion is proportional to the real contact area A, and the formation probability of a wear particle increases with an increase in the real contact area [84]. It is well known that the frictional force, T, is proportional to a normal load, F, when sliding between a pair of contacting bodies occurs. The proportionality between F and A for the contact of rough surfaces is consistent with  .
T F If the contact load is proportional to the contact area, a calibration of the area function of the indenter prior to the indentation test is unnecessary, provided that the proportionality, H = F/A, is known because the load can be known directly during the experiment. In addition, the contact area during loading can be calculated using A = F/H.
The unloading A versus F curves for different materials can be fitted using a power function (see Figs. 7 where D 6 and D 7 are dependent on the material properties and the maximum indentation variables. The non-linearity of the unloading F-A curves manifests  Figure 7(f) compares A versus F curves under different conditions: a flat surface, rough surface, and purely elastic material (the reference material is Au-1). Under a shallow indentation, a proportional relationship between the load and contact area also exists for a purely elastic solid. Under the same load, a purely elastic material underestimates the contact area, whereas a perfectly flat surface overestimates the contact area, because the contact stiffness for an elastic material is the largest, and for a flat surface is the smallest. A purely elastic material can sustain larger contact pressures than an elastic-plastic material, and therefore the same load requires a smaller contact area for an elastic material. A reduced contact area and higher mean contact pressure for a purely elastic solid was also found in [38], which describes an effective way to modify the surface roughness to adjust the sur-face stiffness. The mean contact pressure is dependent on both the material and roughness because   is not a constant but is dependent on the material properties and geometrical parameters [85]. For a spherical indentation of a flat surface, when a plastic deformation is initiated, Hertz theory predicts that   m y / 1.07 p [76], which is close to the value of for a spherical indentation of a linear hardening material, and indicates the role of characteristic stress in the contact of a strain hardening material. function can be used for the F-d relation during loading for both purely elastic and elastic-plastic solids:

Indentation load-displacement curves
For elastic contact of two nominally flat rough surfaces, the contact stiffness (i.e., the slope of the F-d curve during loading) increases proportionally with the normal load F [46], suggesting a power-law relationship [86], although this has yet to be confirmed [87]. A power-law relation exists between the contact stiffness and F for contact between a nominally flat surface and a sphere during loading, demonstrating the critical role of the contact geometry.
The unloading indentation load-displacement curves can be fitted using a plateau curve (see where parameters D 10 and D 11 depend on both the material properties and the maximum loading variables. Another power-law function [88] can also be used for the unloading curves (see Fig. 8(f)): where subscript "f" denotes the value upon final unloading. The exponent D 13 is distinctly larger than 1, implying a nonlinear unloading behavior different from that of a flat punch. It was found that  f m a x / d d 0.75 is in line with the contact deformation of an isolated asperity [85], and exponent D 13 is within the range of 1.6 and 2.3, which is larger than the normal value [89].

Calculated variables
The residual displacement was found to be proportional to the maximum displacement, and  f m a x / d d 0.75 (see Fig. 9(a)). The plastic energy dissipated (i.e., / d d and p t / U U , from indentation tests can be used to calibrate the surface topography, as proposed by [36].
For a smooth indented surface, the slope S of F versus the d curve at the initial unloading is related to the effective modulus, E eff , as in Refs. [81,89,91]: where S can be calculated using the Oliver and Pharr power-law function as As expected, a proportional relationship exists between contact stiffness S and the square root of the contact area (see Fig. 9(c)). Although, the calculated effective elastic modulus E eff is larger than the input value E r used in the simulation, under a shallow indentation (i.e., a small A), the calculated E eff from Eq. (12) is proportional to the input values E r (  eff E r 3.3E ); under a large contact area, the relation between E eff and E r is dependent on the specific strain-hardening rule. If a circular contact region is assumed for contact between a smooth sphere and a rough flat surface [92], the contact area can be approximated as [89] where R is the radius of the spherical indenter, and Figure 9(d) shows similar results as Fig. 9(c), except that the calculated A cal is used instead of the measured A. Under a shallow indentation, the relation between normal contact stiffness S and the calculated contact area A cal is independent of the material, whereas the relation is affected by the strain hardening law under a large indentation. Although the magnitude of the effective elastic modulus calculated from a shallow indentation test is larger than the real value, the results from a shallow indentation can be used for a relative comparison. Under a large indentation, the strain hardening law has a prominent effect on the calculated elastic modulus. Figure 9(e) compares the measured A and calculated A cal . The real contact area for a rough surface is smaller than that calculated when assuming a circular contact. The assumptions (e.g., a circular contact region, and the maximum contact pressure at the center of the surface for contact between a smooth sphere and a rough flat surface) used in theoretical models of the contact of a rough surface [92] should be improved for consistency between the calculated and real contact areas.

Conclusion
Contact between a rigid sphere and a nominally flat rough surface was studied using the finite element method while considering different strain-hardening rules. The surface roughness, measured experimentally using AFM, was applied to generate a 3D FE mesh. Contact variables such as the indentation load, displacement, contact area, and mean contact pressure were obtained for both loading and unloading, and their relations were quantified using a curve fitting of the power law and sigmoid logistic types. The characteristic stress, which depends on both the material properties and the contact geometries, was found to be associated with the mean contact pressure. The deterministic finite element modeling shows a proportional relation between the contact area and load. The surface roughness makes the effective elastic modulus calculated from the initial unloading loaddisplacement curve several times larger than the material value for a shallow indentation. The effects of the wavelengths on the surface profiles with various degrees of surface roughness and surface measurement parameters (e.g., AFM scan resolution and range) will be studied in a future work. The calculated results under a small contact area can only serve as a relative comparison. The normalized contact responses (e.g., contact area and load) are material dependent. The ratios of residual displacement over the maximum displacement, and the plastic energy dissipated after the final unloading over the work applied after loading, were found to be independent of the material properties. The dependence of these ratios on the surface roughness will be studied as a future work.