A statistics view of contact pressure distribution for normal contact of fractal surfaces

Abstract

Due to inherent nonlinear stiffness and damping characteristics, dry friction interfaces have a significant impact on the dynamics of jointed structures. When two surfaces are brought into purely normal contact, contact pressure distribution is of high concern. This work focused on the statistics of the contact pressure distribution between a rough surface and a smooth rigid plane, which provides new insight into the interface contact behaviour. First, the fractal rough surface is generated using the Weierstrass–Mandelbrot function with measured roughness parameters. With meshed rough surfaces, an elastic–plastic finite element contact analysis is performed to determine the contact pressure distribution. Then, the results of contact pressure are statistically analysed. The effects of roughness and contact load on contact area, contact stiffness and mean contact pressure are thoroughly investigated. The probability distribution of contact pressure is determined by fitting a continuous function using a twofold Weibull mixture model. The proposed probability distribution function is found to be capable of describing contact pressure. The contact pressure distribution is affected by the surface fractal characteristics and evolves with the contact load.

Keywords

Dry friction rough surface fractal description contact pressure statistical distribution

Introduction

Dry frictional joints, such as bolted joints and riveted joints, are widely applied in mechanical and aerospace engineering.¹ Subjected to vibration excitation, contact and slip between frictional surfaces result in the nonlinear dynamic response of the jointed structure. It is very important to understand the interface contact and friction characteristics for accurately predicting the dynamics of structures with frictional joints.

The most recent investigations have revealed that the nonlinear dynamic response is highly sensitive to the contact pressure distribution pattern at the contact interface of the friction joint.^2–7 Schwingshackl et al.⁴ divided the frictional interface of an under-platform damper (UPD) into five regions with uniform pressure distribution and assigned five different constant pressure values at each region. Although the contact pressure distributions in their studies were randomly generated and do not represent real contact status, the results emphasized the effect of contact pressure distributions on vibration amplitude. Following this, Schwingshackl et al.⁵ conducted a close visual inspection of a flange joint and identified that a slightly convex shape presented in the interfaces was due to manufacturing imperfections. Thus, they modified the pressure value of individual nodes located along the top and bottom edges of the interface to get a more realistic contact pressure distribution. The modified contact pressure distribution shows higher accuracy in predicting damping than that of not considering manufacturing imperfections. Furthermore, an analytical pressure distribution, which is the sum of two exponentials (one with a positive exponent and the other with a negative exponent), is proposed to describe the concave shape and convex shape of the frictional interface due to manufacturing imperfections.⁶

These investigations highlight the effect of contact pressure distribution patterns on the mechanical response of jointed structures. However, these works only explore the contact pressure distribution resulting from macroscale shape deviation or mesoscale waviness. It is well-known that the real contact surface is covered with microscopic asperities. Therefore, it is vital to explore the real contact pressure distribution caused by these microscale rough characteristics.

Contact pressure distribution due to the macroscopic interface can be easily measured using the pressure-sensitive film.⁸ However, the lateral resolution of the pressure-sensitive film is restricted by the size of the ink microcapsules. Thus, the local distribution of the contact pressure due to the microscopic feature of the contact interface can hardly be observed. The optical method^9,10 is an alternative approach for in-situ measurement of the contact pressure distribution. One type of these methods relies on pre-calibration to accurately set up the correlation between the light intensity and the contact pressure. The other methods measure the three-dimensional (3D) displacement field and then the contact pressure distribution is inversely calculated based on the Boussinesq solutions. Due to the requirement that one of the contact bodies should be transparent, the optical approach is not feasible for mechanical joints.

There are many investigations on elastoplastic rough surface contact^11,12 and loading–unloading contact between rough surfaces.^13–15 However, only Persson's theory can directly predict the probability density function (PDF) of the contact pressure distribution.^16,17 The constant hardness assumption is adopted to approximate the elastic-perfectly-plastic material model. This assumption states that the contact pressure cannot exceed $3 σ_{Y}$ , where $σ_{Y}$ is the yield strength of the weakest contact body. Persson's theory tends to underestimate the PDF of the non-zero contact pressure less than the hardness of the contact body following the J2 plasticity model with strain hardening.^18,19

The finite element (FE) method and boundary element method are two numerical approaches for solving rough surface elastoplastic contact. Boundary element method^20–22 is less popular since the development of the elastoplastic contact model is less straightforward than its linear elastic counterpart, especially when a more realistic J2 plasticity model is implemented. Therefore, the FE method has mainly been regarded as the dominant approach for analysing the contact behaviour of rough surfaces.

The first step of FE rough contact analysis is to generate the rough surface data for establishing the FE model. Generally, there are mainly two kinds of methods applied to obtain the rough surface data for the FE simulation. One is collecting the measured surface data directly, and the other is utilizing geometry description. Thompson et al.^23,24 used the white light interferometer to measure surface data for generating rough surface and found that it is important to consider the sub-surface mesh density for obtaining convergence and reliable results in a rough FE model. On the other hand, the multiscale fractal nature of rough surfaces has long been observed.²⁵

One of the dominant approaches is to generate the fractal rough surfaces using the Weierstrass–Mandelbrot (WM) function.²⁶ Sahoo et al.^27,28 generated rough surfaces using a modified two-dimensional WM function and performed the FE analysis of an elastic–plastic contact with a rigid flat. Then parametric studies, including the effect of variation of fractal parameters on contact behaviour involving contact area, load and displacement, were conducted to describe the general relations between contact properties and surface parameters. Hyun et al.²⁹ generated fractal rough surfaces with a range of roughness amplitudes and roughness exponents using the successive random midpoint algorithm, which was also used by Pei et al.³⁰ to conduct the FE analysis of elastic–plastic contact between rough surfaces. Wang et al.³¹ proposed a spectral interpolation method to construct rough surfaces in conjunction with the measured points and systematically studied elastic–plastic friction contact between rough surfaces and a rigid flat under the combined normal and tangential loading. The results of this study showed that the reduced real contact area due to roughness also governed the friction between the surfaces.

In summary, surface roughness is recognized as having a significant effect on contact friction behaviours and with the aid of FE simulation the contact behaviour, such as the relationship among normal load, normal displacement and real contact area are well investigated. However, the works studying the contact pressure distribution of rough surfaces with microscale asperities, which have an important influence on the nonlinear friction and dynamic response of the assembled structures, are still not sufficient. At present, there is no explicit method to describe the PDF of contact pressure. It is also noted that for interface nonlinear tangential friction modelling, several popular tangential models such as the original Iwan model,^32,33 the four-parameter Iwan model,³⁴ the six-parameter Iwan model^35,36 and so on, are derived by artificially predefining a probability density function (PDF) of critical slipping force or local normal contact pressure. Thus, this article focuses on establishing a mathematical description model with continuous function to determine the form of PDF explicitly.

To address these issues, the FE method is used to simulate the contact characteristics between a rough surface and a rigid plane, and the fractal method is used to describe the geometry of the rough surface. Considering the statistical characteristics and the microcontact mechanism of rough surfaces, the contact pressure distribution is analysed from a new statistical point of view to provide insight into contact pressure distribution and its mathematical description.

The article is organized as follows. In the ‘FE simulation methods’ section, a fractal rough surface is generated using the WM function with measured roughness parameters and the FE modelling method is introduced. In the ‘Simulation cases’ section, three FE models with different roughness are established. In the ‘Results and discussion’ section, the effects of roughness and the normal load on the contact area, contact stiffness and mean contact pressure are investigated. The probability density distribution of contact pressure is determined by mathematical statistics, and the effect of the roughness on the contact pressure distribution is discussed. In the ‘Conclusions’ section, the main conclusions are presented.

FE simulation methods

Fractal description of rough surface

It is well-known that all engineering surfaces are not perfectly smooth and different surface topography characteristics have a significant impact on the contact friction behaviour of connection surfaces and the dynamic response of assembled structures. In order to effectively predict the normal contact and tangential friction behaviour, it is important to accurately obtain the rough surface data. One of the dominant approaches to simulate the rough surface topography is based on the real surface microscopic geometric characteristics.

Microscopic observations have shown that engineering surfaces exhibit a self-affine structural property revealing similar geometrical features at various scales.^37,38 Fractal geometry is characterized by continuity, non-differentiability, and self-affinity. These characteristics are considered to be related to the intrinsic physical laws that govern the generation of surfaces.³⁹ Thus, fractal theory can be used as an effective way to characterize multi-scale surface topography.

The WM function is an appropriate function for fractal rough surface representation satisfying all the previous mathematical properties. Mathematically, the WM function can be seen as sine waves of different wavelengths and amplitudes superimposed on each other with random phase angles. A practical engineering surface can be generated by the modified two-variable WM function, which can be written as⁴⁰:

\begin{aligned} z (x, y) = & L {(\frac{G}{L})}^{(D - 2)} {(\frac{\ln γ}{M})}^{1 / 2} \sum_{m = 1}^{M} \sum_{n = 0}^{n_{m a x}} γ^{(D - 3) n} \\ \times {\cos ϕ_{m, n} - \cos [\frac{2 π γ^{n} {(x^{2} + y^{2})}^{\frac{1}{2}}}{L} \times \cos (\tan^{- 1} (\frac{y}{x}) - \frac{π m}{M}) + ϕ_{m, n}]} \end{aligned}

(1)

where x and y are the Cartesian coordinates of a surface point with height

z (x, y)

; L denotes the length of the square sampling area; D denotes the fractal dimension (

2 < D < 3

); G denotes the fractal roughness parameter;

γ (γ > 1)

is a frequency density control parameter (considering surface flatness and frequency distribution density,

γ = 1.5

⁴⁰; M denotes the number of superposed ridges used to construct the surfaces. n denotes the frequency index (the upper limit of n is given by

n_{m a x} = int [\log (L / L_{s}) / \log γ]

, where

L_{s}

is the cut-off length and

int []

represents the integer part of the number in the bracket; and

ϕ_{m, n}

is a random phase distributed within

[0, 2 π]

The fractal roughness G is the height scaling parameter, which is independent of frequency within the scale range. Physically, the value of fractal roughness G reflects the surface flatness and higher G values represent the rougher surface topography. The fractal dimension D represents the contribution of the high-frequency and low-frequency components in the surface function $z (x, y)$ . Higher D values indicate that high-frequency components are more dominant than low-frequency components in the surface topography profile. Physically, the fractal dimension D is the complexity and irregularity of the fractal rough surface, which is a measurement parameter describing the similarity of surface profiles. The magnitude of the fractal dimension D represents the difference of the surface profiles at different scales. The fractal dimension D is significantly affected by the machining methods and process parameters. Generally, the fractal dimension D of the grinding surface is higher than that of the turning surface.

According to Nayak,⁴¹ any surface represents a real random process with spatial frequency variation. The power spectrum density (PSD) method,^37,42 which decomposes a rough surface into contributions from different spatial frequencies (wavelength), is an effective method to identify fractal parameters. After determining the fractal parameters, MATLAB is used to generate the surface height $z (x, y)$ from equation (1) as per the supplied x and y coordinates. Figure 1 shows an example of a rough surface generated by WM function with m. Here, the length of the sample was specified to be $L = 1 \times 10^{- 3} m$ with $M = 10$ . Figure 2 shows the height probability density (HPD), which indicated that the height of asperity is approximately a Gaussian distribution with mean of 0 and standard deviation of 8.173 µm. From a microscopic point of view, contact between rough surfaces is the interaction between asperities, and asperities of real surfaces are usually postulated to not be sharp peaks but more rounded or curved. To ensure that the simulated surface contains frequency content of the asperity scale, an equivalent asperity radius of $R \approx 40 μ m$ measured by Eriten et al.⁴³ was chosen as a reference. The relationship between cut-off length $L_{s}$ and asperity radius R is $L_{s} = 4 R$ . To ensure that the simulation results remain physically sound, the cut-off length $L_{s}$ is set to be a smaller value $L_{s} = 2 R = 80 μ m$ . Since the surface size is $1 mm \times 1 mm$ , the sampling period in the x and y direction is $10 μ m$ , which is sufficiently lower than the cut-off length.

Figure 1.

Fractal surface generated by Weierstrass–Mandelbrot (WM) function.

Figure 2.

Height probability density (HPD) of the generated fractal surface.

FE modelling

As shown in Figure 1, a 3D fractal rough surface has been generated by equation (1). Then, the coordinates of the fractal rough surface $(x, y, z)$ are imported to Abaqus as key points and those key points are connected in an ordered manner to create the rough surface. Finally, the surface is made solid, which can be used as a 3D model for contact analysis. FE simulations are performed for the 3D elastic–plastic contact between a rough surface and a smooth rigid flat. The smooth rigid flat surface is created to just touch the contact surface from the top. The rigid top surface is set as the master surface and the upper surface of a rough solid is set as the slave surface. As shown in Figure 3, a meshed model for the contact between a rigid plane and a fractal rough surface is created.

Figure 3.

Finite element meshing model for contact between a rigid plane and a rough surface.

The parameters for the rough surface are the same as shown in Figure 1. The mesh resolution of the rough surface is $100 \times 100$ and the rough deformable solid body is meshed with 110000 8-node elements (C3D8). The interaction between contact surfaces is defined as normal ‘Hard’ contact and the tangential behaviour is set as ‘Frictionless’. The movement of the rigid plane is controlled by the normal load P, which is applied from small to large to investigate the evolution of the contact pressure distribution and contact patch distribution. The bottom of a rough solid is fixed. To guarantee the geometric periodicity of the rough surface, the computational domain is expanded in the x and y directions by two additional lines. In this way, two newly generated lines have the same height as its opposite edges. All degrees of freedom of the facing nodes on opposite side surfaces are coupled to achieve a strict periodic boundary condition. To be specific, all the nodes on the side surfaces are put together to create a set, and constraint equations are established between the sets using Python script. It is worth noting that the order of the nodes in the set must satisfy the one-to-one relationship of two facing nodes. Although the copying line operation may probably induce surface discontinuity at the local region, it has little influence on the contact behaviour and the contact pressure distribution.

Simulation cases

According to Zong et al.,⁴⁴ three rough surfaces with different roughness levels were machined by the fine grinding process, and the surface topographies were measured by a 3D white light interferometric profilometer. According to the surface topography measurement results by Zong et al.,⁴⁴ the surface roughness $R_{a}$ is directly calculated, and the fractal dimension D and fractal roughness parameter G are obtained by curve-fitting the power spectral density, as listed in Table 1. By comparing the arithmetical mean deviation $R_{a}$ of the experimentally measured value by Zong et al.⁴⁴ with the numerically simulated fractal rough surface in this article, it is observed that the simulated surface for specimens 2 and 3 are in agreement with the experimentally measured results by Zong et al.,⁴⁴ while $R_{a}$ of simulated and measured surfaces for Specimen 1 have a relative difference of 40%. Overall, the simulated surfaces are related to the surface by Zong et al.⁴⁴

Table 1.

Fractal parameters of three rough surfaces.

Specimen	Fractal dimension $D$	Fractal roughness parameter $G$ ( $m$ )	Arithmetical mean deviation $R_{a}$ ( $μ m$ )
Specimen	Fractal dimension $D$	Fractal roughness parameter $G$ ( $m$ )	Experimental measured	Simulated surface
Specimen 1	2.5428	$3.9933 \times 10^{- 9}$	0.499	0.701
Specimen 2	2.4696	$4.4354 \times 10^{- 8}$	4.625	5.334
Specimen 3	2.4305	$4.7087 \times 10^{- 8}$	8.951	8.616

As shown in Figure 4, three fractal rough surfaces were generated using the fractal parameters listed in Table 1. The two-dimensional (2D) rough surface profile at $x = 0.5 mm$ was also plotted in Figure 4. Then, three FE models with different roughness levels are obtained according to the ‘FE modelling’ section. For these three models, the elastic modulus is $E = 200 GPa$ and Poisson's ratio is $ν = 0.3$ . The composite Young's modulus $E^{*}$ is evaluated as $E / (1 - ν^{2})$ , which is equal to $219.78 GPa$ . The FE simulation is based on elastic-perfectly plastic material with yield stress $σ_{Y} = 400 MPa$ .

Figure 4.

Three fractal rough surfaces with different roughness levels: (a) $R_{a} =$ 0.499 µm, (b) $R_{a} =$ 4.625 µm and (c) $R_{a} =$ 8.951 µm.

Results and discussion

Elastoplastic contact behaviour of rough surfaces

Figure 5 shows the simulated contact pressure of the rough surface with $D = 2.4696, G = 4.4354 \times 10^{- 8} m$ at different normal loads. It is clearly shown that initial contact occurs exclusively around the highest asperities. As normal load increases, plastic deformation becomes increasingly important and the contacting nodes are clustered into large connected regions around these peaks of rough surfaces.

Figure 5.

The contact pressure of the rough surface at different normal loads: (a) $P = 100 N$ , (b) $P = 200 N$ , (c) $P = 400 N$ and (d) $P = 600 N$ .

The profile height of the rough surface has a significant impact on the contact pressure, which determines the elastoplastic contact status of the surface. The Mises yield criterion⁴⁵ is applied to assess whether the surface is elastic or plastic contact. Figure 6 presents the Mises stress at different normal loads, where the blue and red areas are elastically and plastically deformed, respectively. Comparing Figure 6 with Figure 5 it is seen that the projected plastic region is larger than that of the contact region. It is well known that when the asperity is in contact with the rigid plane, a plastically deformed region is first nucleated below the contact region.⁴⁶ As the load increases, the plastic region expands to the surface to form an annulus plastic region until the entire contact area is plastically deformed.^47,48 When a rough surface is in contact with a rigid plane, the contact between one asperity and the rigid plane first occurs. As load increases, plastic deformation of the asperity continues to develop and the other asperity participates in contact. The neighbouring asperities deform in a similar manner to the former asperity, and the two asperities interact with each other. Consequently, the deformation mechanism of the asperity and the interaction of asperities together result in a larger projected plastic region than the contact area. It is noted that in Figure 6, the Mises stress of some nodes exceeds the yield stress, which is caused by the extrapolation of the integral points to the nodes in the FE model.

Figure 6.

Elastic–plastic contact status of the rough surface at different normal loads: (a) $P = 100 N$ , (b) $P = 200 N$ , (c) $P = 400 N$ and (d) $P = 600 N$ .

Based on the simulation results, the dimensionless contact load ( $P^{*}$ ), dimensionless contact area ( $A^{*}$ ) and dimensionless displacement $(δ^{*})$ are normalized as follows:

P^{*} = \frac{P}{E^{*} A_{0}}, A^{*} = \frac{A}{A_{0}}, δ^{*} = \frac{δ}{L}

(2)

where P is the contact load,

A_{0} = L^{2}

is the nominal contact area (

L

being the length of the sample), A is the real contact area, and

δ

is the displacement of the rigid plane.

Figure 7 shows the variation of the dimensionless contact area $A^{*}$ with dimensionless contact load $P^{*}$ . To better analyse the effect of roughness, an additional FE model of the contact between a perfectly smooth surface and a rigid plane is established. As expected, the real contact area of the perfectly smooth surface is the entire nominal contact area under any compressive normal load. For three rough contact surface models, it is seen that the smoother surface has a larger contact area at the same load. At small loads, the relationship between the contact area and the normal load is approximately linear. But as the contact load increases the curve increasingly deviates from this linear proportionality. A linear relationship between contact area and the load was predicted by most statistical contact models,^49,50 which assumed that there is statistically a significant number of asperities in contact and only the tops of the asperities are in contact. As the contact load increases, the latter assumption breaks down and asperities undergo elastoplastic deformation, which contribute to the deviation from linearity at larger loads. This behaviour is also observed in purely elastic contact.²⁹

Figure 7.

Dimensionless contact area versus dimensionless load for different surfaces.

Figure 8 shows the plots of dimensionless load as a function of dimensionless displacement of the rigid plane and the slope of the load–displacement curve represents the contact stiffness. For a perfectly smooth surface, the contact stiffness exhibits a bilinear behaviour. Before the turning point, only elastic deformation occurs. Following the turning point, the material experiences plastic deformation, resulting in a decrease in stiffness. Due to the contact status involving just linear elastic contact of several asperities, the contact stiffness for rough surfaces grows gradually at the initial loading stage. Some asperities undergo plastic deformation as the load increases, whereas others retain elasticity. During this stage, the contact stiffness is approximately constant with slight fluctuations caused by the randomness of the rough surface. When the load is high enough, the deformable solid body progressively becomes fully plastic, resulting in a further drop in contact stiffness. Finally, the asperity undergoes complete plastic deformation, and thus the contact stiffness of the rough surface approaches that of the smooth surface.

Figure 8.

Dimensionless load versus dimensionless displacement for different surfaces.

The variation of the mean contact pressure with dimensionless contact load is shown in Figure 9. The mean contact pressure is the ratio of the normal load to the real contact area. Therefore, the contact pressure of a perfectly smooth surface increases linearly with the load because the contact area is the nominal contact area under any compressive normal load. However, the mean contact pressure of rough surfaces is greater than that of perfectly smooth surfaces because of the stress concentration at the asperity contact spots. For surface with $R_{a} = 8.951 μ m$ and $R_{a} = 4.624 μ m$ , the mean contact pressure first remains constant with increasing normal load. However, for surfaces with $R_{a} = 0.499 μ m$ , the mean contact pressure keeps increasing as normal load increases. This may be probably due to the different elastic–plastic deformation of the asperity caused by surface roughness. Under a small normal load, the asperities of a smoother surface predominantly undergo elastic deformation, while the asperities of a rougher surface mainly occur plastic deformation. When the normal force is sufficiently large, the mean contact pressure converges to the contact case of the perfectly smooth. The convergence trend is mostly attributable to the fact that under large normal loads almost all asperities are involved in the contact and the rough surface is nearly flattened completely.

Figure 9.

Mean contact pressure versus dimensionless load for different roughness levels.

In the original GW model,⁴⁹ plasticity index $Ψ = (\frac{E^{*}}{H}) \sqrt{\frac{σ}{R}}$ was defined to determine the mode of deformation (elastic or plastic), where $E^{*}$ is the composite Young's modulus, H is the hardness, $σ$ is the standard deviation of the height distribution and R is the average asperity radius of curvature. The 9-point peak criterion (9 PP criterions) was used to identify an asperity and the asperity radius was determined using equation (1) by Zhou et al.⁵¹ The asperity radius of each rough surface was determined by averaging asperity radius of all asperity peak. For three rough surfaces ( $R_{a} =$ 0.499, 4.625 and 8.951 µm) generated in this article, the asperity radii are 202, 121 and 65 µm, respectively. Therefore, the plasticity indices are 6.22, 43.31 and 71.91, respectively. The statistics-based rough contact model¹⁴ showed that for the same normal load the contact area gradually decreases and the dimensionless displacement increases as the plasticity index increases, which is consistent with simulation results in Figures 7 to 9. In general, the surface roughness reduces the real contact area, resulting in an entirely different relationship between contact pressure and contact load compared to smooth surfaces. The contact occurs in a more localized region due to the roughness, which leads to a reduction in the overall contact stiffness.^52,53

Description of contact pressure by statistical method

Statistical modelling of contact pressure

The CPRSS of the field output in software Abaqus describes the contact pressure at each node. The contact pressure data of all contacting nodes are extracted to conduct the statistical analysis. The PDF of contact pressure distribution is determined by the contact pressure data within the contact area, not the entire contact interface. Taking Specimen 2 with $R_{a} = 4.625 μ m$ as an example, the histogram in Figure 10 depicts the probability density of contact pressure under 400 N preload. It is apparently seen that the contact pressure distribution of the real surface is neither uniform distribution nor exponential distribution, but a bimodal distribution. The two peaks in the histogram are mainly caused by the contact pressure within the elastically deformed and plastically deformed regions, respectively. However, the contact pressure distribution, which essentially determines the probability distribution of slider critical slip force in the Iwan model, was artificially assumed to follow a uniform distribution or an exponential distribution without a clear physical explanation.^54,55 To facilitate the derivation of the continuous Iwan density function, it is necessary to know the continuous distribution function of the contact pressure.

Figure 10.

Contact pressure distribution under 400 N preload for Specimen 2.

Three popular distributions (Gaussian, Weibull and twofold Weibull mixture) are employed to model the probability distribution function (PDF) and cumulative distribution function (CDF) for the contact pressure data. The Gaussian distribution (normal distribution) and Weibull distribution are extensively used in continuous statistical functions. The formula for the Weibull distribution function with three parameters is written as

f (x) = \frac{β}{α} {(\frac{x - γ}{α})}^{β - 1} e^{- {(\frac{x - γ}{α})}^{β}}

(3)

where

f (x) \geq 0

represents the probability density, x is the independent variable under observation,

α

denotes the scale parameter,

β

denotes the shape parameter and

γ

denotes the location parameter.

The corresponding estimated parameters of Gaussian distribution and Weibull distribution are calculated with the use of the software Matlab by the toolbox function ‘fitdist’. As shown in Figure 11, both Gaussian and Weibull distributions do not fit the contact pressure histogram very well. Noticing the bimodal nature of the histogram, a twofold Weibull mixture model with each representing one type of distribution is proposed to fit the distribution of contact pressure. Mathematically, the bimodal distribution is the superposition of two unimodal distributions. A twofold Weibull mixture model is written as

g (x) = \sum_{i = 1}^{2} w_{i} \frac{β_{i}}{α_{i}} {(\frac{x - γ_{i}}{α_{i}})}^{β_{i} - 1} e^{- {(\frac{x - γ_{i}}{α_{i}})}^{β_{i}}}

(4)

where

α_{i}

β_{i}

and

γ_{i} (i = 1, 2)

are the parameters of ith Weibull probability density function,

w_{i} (i = 1, 2)

is the weight of ith distribution with

w_{1} + w_{2} = 1

Figure 11.

The probability density function of contact pressure.

For simple Weibull distribution, the nonlinear relationship between the sample data x and the corresponding function value $f (x)$ can be converted to a linear relationship by simple mathematical transformation. Then, the linear regression method can be applied to estimate the parameters. However, this approach is incapable of handling a mixture model. Parameter estimation of the twofold Weibull mixture model is more complicated compared to simple Weibull distribution since the twofold Weibull mixture model contains more independent parameters. For the twofold Weibull mixture model, the parameter estimation problem is transformed into a data-fitting problem in this article. A nonlinear least squares regression technique was applied to the probability density data of the contact pressure to establish an optimized objective function and the unknown parameters were solved by means of the genetic algorithm. Discrete samples of contact pressure probability density are denoted as $(p_{j}, P_{j}) (j = 1, 2, \dots, m)$ , where m represents the number of samples, $p_{j}$ denotes the discrete contact pressure and the corresponding probability density is denoted as $P_{j}$ . Then the objective function can be written as

min F (θ) = \sum_{j = 1}^{m} [P_{j} - g (p_{j})]^{2}

(5)

where

θ = (w_{1}, α_{1}, β_{1}, γ_{1}, w_{2}, α_{2}, β_{2}, γ_{2})

denotes the vector of unknown parameters. A set of optimal parameters can be found by minimizing the value of the objective function using the genetic algorithm. Considering that the value of the contact pressure is positive, the constraint

γ_{i} \leq 0

is imposed in optimization to ensure that the function remains defined. Figure 11 reveals that, among three distributions, the twofold Weibull distribution appears to predict very well the contact pressure distribution of real rough surfaces: a slight decrease when the contact pressure is small, followed by an increase and then a decrease when the contact pressure is large.

Checking the validity of the twofold Weibull mixture model

The CDFs of the Gaussian distribution and Weibull distribution can be found by Murthy et al.⁵⁶ and Ahsanullah et al.⁵⁷ and the CDF of the twofold Weibull mixture model (equation (6)) can be derived using its probability density function (equation (4)).

G (x) = 1 - \sum_{i = 1}^{2} w_{i} e^{- {(\frac{x - γ_{i}}{α_{i}})}^{β_{i}}}

(6)

By presenting the simulated contact pressure in ascending order, the cumulative distribution of the contact pressure can be calculated. Figure 12 presents the cumulative distribution corresponding to the probability density function shown in Figure 11.

Figure 12.

Cumulative distribution function of contact pressure.

$D_{n}$ is defined to examine how well a hypothesized distribution function $G (x)$ fits an empirical distribution function $F_{n} (x)$ , where $G (x)$ denotes the CDF of three hypothesized distributions in this article and $F_{n} (x)$ is the empirical distribution function for the observations.

D_{n} = sup_{x} | G (x) - F_{n} (x) |

(7)

where sup denotes supremum.

The values of $D_{n}$ for the Gaussian distribution, Weibull distribution and the twofold Weibull distribution are 0.3198, 0.3854 and 0.0540, respectively. Therefore, based on finite element simulation results, the probability density function of the contact pressure distribution might be approximated as a twofold Weibull distribution. Following this, the goodness-of-fit of contact pressure distribution varying with normal loads and roughness is investigated based on the twofold Weibull distribution. To illustrate the goodness-of-fit of the twofold Weibull distribution for different normal loads and different roughness, Table 2 lists the value of $D_{n}$ for nine sets of contact pressure data (three specimens under three different normal loads). The maximum value of all nine sets is < 0.1, which suggests that the twofold Weibull distribution is an adequate choice to describe the PDF of the contact pressure.

Table 2.

The value of $D_{n}$ for all three specimens under three normal loads.

Normal load	200 N	400 N	600 N
Specimen 1	0.0823	0.0823	0.0997
Specimen 2	0.0963	0.0540	0.0418
Specimen 3	0.0798	0.0352	0.0338

Contact pressure distribution characteristics

Effect of load and roughness on contact pressure distribution

As shown in Figures 13 to 15, the simulated discrete values of contact pressure at varying applied loads for different surfaces are extracted, and statistical analysis is performed to determine the discrete values of probability density. The horizontal axis represents the range of contact pressure. The vertical axis represents the probability density corresponding to the contact pressure. A good agreement is observed between simulated results and the twofold Weibull mixture model for each of the applied normal loads. It is also noted that some fitted curves do not represent well the distribution around 0, which might be related to the optimization process that minimizes the global error without ensuring the minimization of the local error. The contact pressure distribution of a perfectly smooth surface is exactly a uniform distribution. However, as shown in Figures 13 to 15, it is seen that the contact pressure distribution of the rough surface is significantly different from that of the smooth surface. Compared with the Dirac function of the contact pressure PDF for the perfectly smooth surface, the PDF of the contact pressure of the rough surface is broadened with a larger variance towards the low contact pressure. With the increase of roughness, the variability of contact pressure increases and the proportion of high contact pressure increases significantly. The presence of the small contact pressure is mostly attributable to the elastic contact pressure at the margin of the real contact area.

Figure 13.

Contact pressure distribution at varying applied loads for $R_{a} =$ 0.499 µm: (a) $P = 200 N$ , (b) $P = 400 N$ and (c) $P = 600 N$ .

Figure 14.

Contact pressure distribution at varying applied loads for $R_{a} =$ 4.625 µm: (a) $P = 200 N$ , (b) $P = 400 N$ and (c) $P = 600 N$ .

Figure 15.

Contact pressure distribution at varying applied loads for $R_{a} =$ 8.951 µm: (a) $P = 200 N$ , (b) $P = 400 N$ and (c) $P = 600 N$ .

Figure 16 shows the fitting results of twofold Weibull mixture model for different surfaces at varying applied loads. For the same normal contact load, the peak value of the contact pressure probability density moves towards higher contact pressures as the roughness increases. The probability density function changes significantly when the roughness increases from $R_{a} =$ 0.499 to 4.625 µm, while the difference of probability density function between rough surfaces with $R_{a} =$ 4.625 and 8.951 µm is smaller. One of the possible reasons is that surfaces with $R_{a} =$ 0.499 µm quickly reach nearly full contact, while the process of approaching full contact is more similar for surfaces with $R_{a} =$ 8.951 and 4.625 µm.

Figure 16.

Fitting results of twofold Weibull mixture model for different surfaces: (a) $P = 200 N$ , (b) $P = 400 N$ and (c) $P = 600 N$ .

Effect of yield stress on contact pressure distribution

To analyse the effect of yield stress $σ_{Y}$ on contact pressure redistribution, three models with different yield stresses ( $σ_{Y} = 200, 400 and 600 MPa$ ) are established. As shown in Figure 17 under three different normal loads, the peak of contact pressure probability density shifts to the right with the increase of yield stress, and the peak value of probability density decreases. Additionally, the contact pressure corresponding to peak probability density is about $3 σ_{Y}$ , which is approximately the value of hardness.^11,46,47

Figure 17.

Variation of contact pressure distribution with yield stress: (a) $P = 200 N$ , (b) $P = 400 N$ and (c) $P = 600 N$ .

Conclusions

This work simulated the normal contact behaviour of fractal rough surfaces by employing the FE method. The aim of this article is to obtain a more realistic or physically-based contact pressure distribution form and investigate the effect of roughness and normal load on contact pressure distribution. The fractal rough surface was generated using the modified two-variable WM function. The contact pressure distribution between a fractal rough surface and a rigid plane was investigated by combining the FE method and statistical analysis.

Surface roughness reduces the real contact area, resulting in a completely different relationship between contact pressure and contact load compared to perfectly smooth surfaces, as determined by numerical simulations. Due to the roughness, the contact occurs in a more localized region, which reduces the overall contact stiffness. The mean contact pressure of a perfectly smooth surface increases with the load increasing. In contrast, for rough surfaces, the mean contact pressure initially remains constant as the normal load increases, and when the normal force is sufficiently large, the mean contact pressure converges to the contact case of the perfectly smooth. The probability density distribution of contact pressure is determined using the statistical method and might be approximated as a twofold Weibull distribution.

This article proposed a mathematical description model with a continuous function to explicitly determine the PDF of the contact pressure distribution. In fact, the contact pressure distribution is largely influenced by the geometry of contact bodies as well as the waviness and roughness of surfaces. The influence of geometry deviation and waviness on contact pressure distribution will be considered in the future.

Footnotes

Acknowledgements

The authors would like to acknowledge support from the National Natural Science Foundation of China (Nos. 12072268 and 52105179) and the China Science Challenge project (No. TZ2018007).

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (grant no. 12072268 and 52105179).

ORCID iDs

Huiyi Yang

Chao Xu

Appendix A

Figure 18(a) shows the FE model for Specimen 2 ( $R_{a} = 4.625 μ m$ ) used in this article. As shown in Figure 18(b), a model with fine mesh along the vertical direction was constructed to investigate the mesh convergence along the vertical direction. The element number along the vertical direction of the fine mesh model is about 1.45 times that of the original model. Figure 19 shows the relationship between the dimensionless contact area and dimensionless normal load for two models. It is clearly seen that the variation of contact area with normal load is almost the same for two different vertical mesh resolutions. The probability density functions of the contact pressure distribution are also compared in Figure 20 and there is no significant difference of the simulation results between the two models. Therefore, the contact pressure results are converged for 11 vertical element number used in this article.

To analyse the influence of height the height of the bulk, two additional models with height = 0.75 and 1 mm are established for comparative analysis, as shown in Figure 18(c) and (d). Figures 19 and 20 show the comparative results among three FE models with height = 0.5, 0.75 and 1 mm. It is seen that the curve between the contact area and normal load and the contact pressure distribution for the three models are in good agreement. Therefore, it is feasible to investigate the contact pressure distribution using the FE model with height = 0.5 mm.

Since the computation cost is less for surfaces with smaller roughness, three FE models ( $with different mesh resolutions 100 \times 100$ , $150 \times 150$ and $200 \times 200$ ) for Specimen 1 ( $R_{a} = 0.499 μ m$ ) were created to investigate the mesh convergence along the x and y directions. The contact pressure of three models with different mesh resolutions is compared in Figure 21. There is no significant difference of the simulation results among the three models. Therefore, it is feasible to simulate the contact behaviour with 100 × 100 mesh resolution considering both the convergence and computation time.

References

Gaul

Nitsche

. The role of friction in mechanical joints. Appl Mech Rev 2001; 54: 93–106.

Schwingshackl

Petrov

. Modeling of flange joints for the nonlinear dynamic analysis of gas turbine engine casings. J Eng Gas Turbines Power Trans Asme 2012; 134: 122504. DOI: 10.1115/1.4007342.

Pesaresi

Salles

Jones

, et al. Modelling the nonlinear behaviour of an underplatform damper test rig for turbine applications. Mech Syst Signal Process 2017; 85: 662–679.

Schwingshackl

Petrov

Ewins

. Effects of contact interface parameters on vibration of turbine bladed disks with underplatform dampers. J Eng Gas Turbines Power-Trans Asme 2012; 134: 032507. DOI: 10.1115/1.4004721.

Schwingshackl

Di Maio

Sever

, et al. Modeling and validation of the nonlinear dynamic behavior of bolted flange joints. J Eng Gas Turbines Power Trans Asme 2013; 135: 122504. DOI: 10.1115/1.4025076.

Schwingshackl

. Identification reassembly uncertainties for a basic lap joint. Dynamics of coupled structures. 4. Switzerland: Springer, 2017, pp. 53–61.

Schwingshackl

. Measurement of friction contact parameters for nonlinear dynamic analysis. Topics in modal analysis. 5. Switzerland: Springer, 2012, pp. 167–177.

Ise

Kobayashi

Itoh

, et al. Pressure sensitive paint visualization of pressure distribution of externally pressurized circular thrust gas bearing with a single gas supply hole: experimental validation of variation of pressure distribution by supply pressure. J Tribol Trans Asme 2017; 139: 054501. DOI: 10.1115/1.4035156.

Sharp

Poole

Kleiman

. Optical measurement of contact forces using frustrated total internal reflection. Phys Rev Appl 2018; 10: 034051. DOI: 10.1103/PhysRevApplied.10.034051.

10.

Bai

Cao

, et al. Imaging dynamic three-dimensional traction stresses. Sci Adv 2022; 8: eabm0984. DOI: 10.1126/sciadv.abm0984.

11.

Tiwari

Almqvist

Persson

BNJ

. Plastic deformation of rough metallic surfaces. Arxiv 2020; 68: 129. DOI: arXiv:2006.11084.

12.

Chen

. Fractal modeling of elastic-plastic contact between three-dimensional rough surfaces. Ind Lubr Tribol 2018; 70: 290–300.

13.

Yuan

Zhao

. The loading-unloading model of contact between fractal rough surfaces. Int J Prec Eng Manuf 2020; 21: 1047–1063.

14.

Kadin

Kligerman

Etsion

. Unloading an elastic-plastic contact of rough surfaces. J Mech Phys Solids 2006; 54: 2652–2674.

15.

Wang

Jia

Wang

, et al. Research on the loading-unloading fractal contact model between two three-dimensional spherical rough surfaces with regard to friction. Acta Mech 2020; 231: 4397–4413.

16.

Persson

BNJ

. Elastoplastic contact between randomly rough surfaces. Phys Rev Lett 2001; 87: 116101. DOI: 10.1103/PhysRevLett.87.116101.

17.

Mulvihill

. Revisiting the Persson theory of elastoplastic contact: a simpler closed-form solution and a rigorous proof of boundary conditions. Tribol Lett 2022; 70: 91. DOI: 10.1007/s11249-022-01633-z.

18.

Frerot

Anciaux

Molinari

J-F

. Crack nucleation in the adhesive wear of an elastic-plastic half-space. J Mech Phys Solids 2020; 145: 104100. DOI: 10.1016/j.jmps.2020.104100.

19.

Venugopalan

Irani

Nicola

. Plastic contact of self-affine surfaces: Persson's theory versus discrete dislocation plasticity. J Mech Phys Solids 2019; 132: 103676. DOI: 10.1016/j.jmps.2019.07.019.

20.

Almqvist

Sahlin

Larsson

, et al. On the dry elasto-plastic contact of nominally flat surfaces. Tribol Int 2007; 40: 574–579.

21.

Weber

Suhina

Junge

, et al. Molecular probes reveal deviations from Amontons’ law in multi-asperity frictional contacts. Nat Commun 2018; 9: 888. DOI: 10.1038/s41467-018-02981-y.

22.

Frerot

Bonnet

Molinari

J-F

, et al. A Fourier-accelerated volume integral method for elastoplastic contact. Comput Methods Appl Mech Eng 2019; 351: 951–976.

23.

Thompson

. A comparison of methods to evaluate the behavior of finite element models with rough surfaces. Scanning 2011; 33: 353–369.

24.

Thompson

. Considerations for the incorporation of measured surfaces in finite element models. Scanning 2010; 32: 183–198.

25.

Mandelbrot

. The fractal geometry of nature. New York: WH Freeman, 1982.

26.

Majumdar

Tien

. Fractal characterization and simulation of rough surfaces. Wear 1990; 136: 313–327.

27.

Sahoo

Ghosh

. Finite element contact analysis of fractal surfaces. J Phys D-Appl Phys 2007; 40: 4245–4252.

28.

Sahoo

Adhikary

Saha

. Finite element based elastic-plastic contact of fractal surfaces considering strain hardening. J Tribol Surf Eng 2010; 1: 39–56.

29.

Hyun

Pei

Molinari

, et al. Finite-element analysis of contact between elastic self-affine surfaces. Phys Rev E 2004; 70: 026117. DOI: 10.1103/PhysRevE.70.026117.

30.

Pei

Hyun

Molinari

, et al. Finite element modeling of elasto-plastic contact between rough surfaces. J Mech Phys Solids 2005; 53: 2385–2409.

31.

Wang

, et al. The effect of resolution on the deterministic finite element elastic–plastic rough surface contact under combined normal and tangential loading. Tribol Int 2020; 144: 106141. DOI: 10.1016/j.triboint.2019.106141.

32.

Iwan

. On a class of models for the yielding behavior of continuous and composite systems. 1967.

33.

Iwan

. A distributed-element model for hysteresis and its steady-state dynamic response. J Appl Mech 1966; 33: 93.

34.

Segalman

. Modelling joint friction in structural dynamics. Struct Control Health Monit 2006; 13: 430–453.

35.

Hao

Feng

, et al. Investigation into discretization methods of the six-parameter Iwan model. Mech Syst Signal Process 2017; 85: 98–110.

36.

Hao

. A six-parameter Iwan model and its application. Mech Syst Signal Process 2016; 68-69: 354–365.

37.

Jacobs

TDB

Junge

Pastewka

. Quantitative characterization of surface topography using spectral analysis. Surf Topogr Metrol Prop 2017; 5: 013001. DOI: 10.1088/2051-672X/aa51f8.

38.

Gujrati

Sanner

Khanal

, et al. Comprehensive topography characterization of polycrystalline diamond coatings. Surf Topogr Metrol Prop 2021; 9: 014003. DOI: 10.1088/2051-672X/abe71f.

39.

Komvopoulos

Yan

. A fractal analysis of stiction in microelectromechanical systems. J Tribol Trans Asme 1997; 119: 391–400.

40.

Yan

Komvopoulos

. Contact analysis of elastic–plastic fractal surfaces. J Appl Phys 1998; 84: 3617–3624.

41.

Nayak

. Random process model of rough surfaces. J Lubr Technol 1971; 93: 98.

42.

Persson

BNJ

. On the fractal dimension of rough surfaces. Tribol Lett 2014; 54: 99–106.

43.

Eriten

Polycarpou

Bergman

. Physics-based modeling for fretting behavior of nominally flat rough surfaces. Int J Solids Struct 2011; 48: 1436–1450.

44.

Zong

Qin

Chu

. Modeling of frictional stick-slip of contact interfaces considering normal fractal contact. J Appl Mech Trans Asme 2022; 89: 031003. DOI: 10.1115/1.4052882.

45.

Lubliner

. Plasticity theory. London: Pearson Education, 2008.

46.

Johnson

. Contact mechanics. Cambridge: Cambridge University Press, 1987.

47.

Kogut

Etsion

. Elastic–plastic contact analysis of a sphere and a rigid flat. J Appl Mech Trans Asme 2002; 69: 657–662.

48.

Jackson

Green

. A finite element study of elasto–plastic hemispherical contact against a rigid flat. J Tribol Trans Asme 2005; 127: 343–354.

49.

Greenwood

Williamson

. Contact of nominally flat surfaces. Proc R Soc A 1966; 295: 300–319.

50.

Yuan

Long

Ding

, et al. Statistical contact model of rough surfaces: the role of surface tension. Int J Solids Struct 2018; 138: 217–223.

51.

Zhou

Wang

, et al. Three-dimensional asperity model of rough surfaces based on valley-peak ratio of the maximum peak. Meccanica 2021; 56: 711–730.

52.

Kartal

Mulvihill

Nowell

, et al. Measurements of pressure and area dependent tangential contact stiffness between rough surfaces using digital image correlation. Tribol Int 2011; 44: 1188–1198.

53.

Akarapu

Sharp

Robbins

. Stiffness of contacts between rough surfaces. Phys Rev Lett 2011; 106: 204301. DOI: 10.1103/PhysRevLett.106.204301.

54.

Song

Hartwigsen

McFarland

, et al. Simulation of dynamics of beam structures with bolted joints using adjusted Iwan beam elements. J Sound Vib 2004; 273: 249–276.

55.

Segalman

. A four-parameter Iwan model for lap-type joints. J Appl Mech-Trans Asme 2005; 72: 752–760.

56.

Murthy

Xie

Jiang

. Weibull models. New Jersey: John Wiley & Sons, 2004.

57.

Ahsanullah

Kibria

Shakil

. Normal distribution. In: Normal and student st distributions and their applications. Switzerland: Springer, 2014, pp. 7–50.