Numerical Study on Thin Film Lubrication Performance with Imperfect Coating under the Effect of the Electrical Double Layer in Point Contact

Abstract

Revealing the interface effect on lubrication is great significance for designing high performance water-based lubricating friction pairs. In this study, a new thin film lubrication numerical model considering the electrical double layer (EDL) effect with an imperfect coating in point contact is proposed, and the validity and effectiveness of the model is verified. The numerical results show that increasing the zeta potential increases the film thickness, which indicates that the EDL effect can improve the firm lubricated film formation. However, the effect of the zeta potential on the von Mises stress is very small because the liquid film pressure remains constant at different surface potentials. On the other hand, with an increase in coating thickness, the effects of the coating on the film thickness and pressure gradually become evident, but the von Mises stress is affected by the imperfectly bonded interface and coating thickness. Finally, the effect of surface roughness on thin film lubrication was analyzed. The liquid film thickness slightly decreased and the occurrence of stress concentration on the surface coating was evident. The proposed model was expected to provide a calculation basis for solving the thin film lubrication problem of coatings affected by the EDL.

1. Introduction

In nanolubrication, the solid–liquid interface effects play a fundamental role, especially in water-based liquid lubrication. Among the interface effects, the electrical double layer (EDL) effect is a typical interfacial effect in confined gaps that has received considerable attention [1,2,3]. In the EDL, the charge on the wall surface absorbs the nearby counter-ions and repels the same ions in the flowing fluid, forming an adsorption layer near the surface during the movement of liquid molecules. Far away from the surface, positive and negative ions diffuse along the flow direction to form a diffusion layer. The EDL effect exists in the lubrication system, as ceramic lubrication [4], electrotunable friction [5] and liquid superlubricity, for instance [6]. Therefore, an EDL structure can significantly affect the lubrication performance.

Relevant studies [7,8,9] have shown that when the thickness of the lubricating film is less than 100 nm, the EDL effect cannot be ignored. Therefore, analyzing the effect of the EDL on lubrication is crucial in thin film lubrication (TFL). Wang et al. [10] studied in detail the action mechanism of ceramics in water lubrication and indicated that the EDL produced on a ceramic surface is an important factor affecting the lubrication performance. Funari et al. [2] used a quartz crystal microbalance to investigate the effect of the EDL structure on the rheological properties. Based on the EDL effect and friction-induced electric field on the electroosmosis behavior, Luan et al. [3] studied the electro-osmosis effect of a ceramic ball surface under a water-based lubricant and analyzed the EDL effect on the lubrication performance. In their study, free ions in the EDL diffusion layer drove liquid molecules to form the electro-osmotic flow (EOF) and surfactants mainly affected lubrication performance through EOF. With in-depth research, researchers have obtained excellent lubrication performance in polytetrafluoroethylene (PTFE)/steel friction pairs by using ionic lubricants. Their experiments show that thick and dense electrical double layers and the formation of friction chemical films can achieve superlubricity [11]. The above studies show that the EDL is an important factor in TFL, which can promote the development of high-performance lubrication.

Theoretical modeling is one approach for studying the effects of EDL on TFL. Chen et al. [12] considered the influence of asymmetric electrical potential and established an overlapping EDL model to study the effect of κh (κ is the Debye length and h is the film thickness) and surface potential on the electrical viscosity. Jing et al. [13] proposed a numerical model for discussing the effect of the zeta potential on fluid lubrication in sliding bearing. Bai et al. [8] deduced a mathematical model of electrical viscosity based on the Poisson–Boltzmann equation and proposed a modified Reynolds model, considering the effect of the EDL. In addition, Zuo et al. [1,9] considered the influence of different material surfaces and studied the effect of an asymmetric EDL and temperature on film lubrication by numerical simulation. Fang et al. [14] established a numerical analysis method for investigating the TFL performance of rough surfaces, considering the EDL effect and the conditions for achieving ultralow friction in water-based lubrication. These methods can be used to analyze the effect of the EDL on lubrication performance.

On the other hand, a single material cannot satisfy the requirements of complex working conditions. Coatings [15,16,17,18] are often used to improve the performance of construction machinery. The study of coating lubrication has received much attention. Wang et al. [19] analyzed the effect of shear stress on the coating stress under elastohydrodynamic lubrication (EHL) by establishing a mathematical model derived from the Papkovich–Neuber potentials and fast Fourier transform (FFT). Based on a finite element method, Habchi [20] studied the influence of the thermal and mechanical properties of coatings on EHL. Alakhramsing et al. [21] modified the coating lubrication model by Habchi [20] and analyzed the EHL of coatings under finite line contact. Wu et al. [22] established a mixed friction load model for the coated angular contact ball bearings. The aforementioned studies considered the coating layer to be perfect; however, coatings often have interface flaws due to technical defects, materials, and other factors. In practical applications, the coating layer structure can lead to variations of the stress [23], resulting in delamination and lubrication failure. The existing imperfect coating interface models can be divided into spring-like [24], frictionless [25], dislocation-like, and force-like [26,27] interfaces. Based on the interface model proposed by Wang et al. [27], He et al. [28] considered the defect interface established an imperfect coating–substrate interface model under EHL. Researchers used models to study the effect of imperfect coatings on lubrication performance under EHL. Furthermore, in coatings layer system, such as ceramic–metal coating structure, it can form a more stable charge distribution in the lubricated region, thereby affecting the formation and characteristics of the EDL, which play an important role in lubrication. However, under thin film lubrication, there is still a lack of corresponding research on the impact of the electrical double layer effect on the load-bearing capacity of the lubricating film with coating. At the same time, in regard to micro–nano motion devices, the stress distribution under different conditions must be considered, especially in terms of the influence of rough surfaces.

Considering the above research gap, a point-contact TFL model based on an imperfect coating–substrate interface, which considers an asymmetric EDL structure, is proposed. Section 2 introduces the theoretical method, and Section 3 presents the numerical method and its quantitative validation. Section 4 discusses the effects of the EDL, coating properties, imperfect coating–substrate interface, and surface roughness on the lubrication behavior of TFL. Finally, Section 5 summarizes the conclusions of this study.

2. Theoretical Method

Figure 1 shows a point-contact lubrication model between a sliding ball and coating surface. The EDL effect can appear on the solid–liquid interface. The coating–substrate interface is used the interface jumping coefficient t_i (i = 1, …,6) to express the discontinuities of the displacement and stress. Here, $t_i=0$ represents complete detachment, whereas $t_i=1$ represents complete adhesion. The coating thickness is h₁. The average velocity along the x-direction in this model is $U=(u_s+u_b)/2$, where $u_b$ and $u_s$ are the ball and substrate velocity, respectively.

2.1. Lubrication Model

In the thin film lubrication, considering the EDL effect, the isothermal Reynolds equation can be modified to obtain the liquid pressure p as follows:

$$\frac{\partial }{\partial x}\left(\frac{\rho h^3}{{\eta }_a}\frac{\partial p}{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{\rho h^3}{{\eta }_a}\frac{\partial p}{\partial y}\right)=12U\frac{\partial \rho h}{\partial y},$$

(1)

where h and ρ represent the film thickness and density, respectively, and ${\eta }_a$ denotes apparent viscosity of the lubricant considering the effect of the EDL. The apparent viscosity ${\eta }_a$ consists of the lubricant viscosity and electrical viscosity, that is, ${\eta }_a=\eta +{\eta }_e$. The pressure p must satisfy the following boundary conditions:

$$\left\{\begin{array}{l}p(x_0,y)=0\\ p(x_e,y)=0,\frac{\partial p(x_e,y)}{\partial x}=0\end{array}\right..$$

(2)

The relationship between the lubricant viscosity and pressure is based on the Roelands model, as follows:

$$\eta ={\eta }_0\exp \{[\ln {\eta }_0+9.67][-1+{(1+5.1\times {10}^{-9}p)}^{0.68}]\}.$$

(3)

The asymmetric potential EDL model is used for the electrical viscosity, as this model is more consistent with the physical reality [1,9]. When the upper and lower surface potentials are ${\xi }_1$ and ${\xi }_2$, respectively, the total surface potential can be determined as ${\xi }_a={\xi }_1+{\xi }_2$. Subsequently, the electrical viscosity can be calculated as:

$${\eta }_e=\frac{3{\epsilon }_e^2{\xi }_a{}^2}{4{\pi }^2{\kappa }^2{h}^4\lambda }\left(\frac{\cosh \kappa h-1}{\sinh \kappa h}-\frac{\kappa h}{2}\right).$$

(4)

The relationship between the lubricant density and pressure is provided by the Dowson–Higginson model [29]:

$$\rho ={\rho }_0(1+\frac{0.6\times {10}^{-9}p}{1+1.7\times {10}^{-9}p}).$$

(5)

The applied load $w$ is balanced with the lubrication film pressure in the contact area:

$$w={\displaystyle {\iint }_{\Omega }pdxdy}.$$

(6)

The film thickness equation between two surfaces is expressed as:

$$h\left(x,y\right)=h_0+h_s\left(x,y\right)+u_3\left(x,y\right)+\delta \left(x,y\right),$$

(7)

where $h_0$ is the gap between the upper and lower surfaces, $h_s$ describes the geometry of the two contact surfaces, $\delta$ is the roughness amplitude of the upper and lower surfaces, and $u_3$ indicates the surface deformations of the upper and lower surfaces under pressure. While the influence coefficient between the deformation and normal pressure is not directly available from the analytic solution, the response function is derived in the frequency domain in [27]. This study ignores the influence of shear stress on the surface displacement. Therefore, the deformations of the sphere and coating surface are calculated as follows:

$$u_3=IFFT({\widehat{\widehat{C}}}_p^{u_3^{}}\ast \widehat{\widehat{p}}),$$

(8)

where ${\widehat{\widehat{C}}}_p^{u_3^{}}$ denotes the response function of the pressure displacement in the frequency domain and IFFT denotes the inverse FFT. This calculation deformation process is described in Section 2.2.

2.2. Surface Deformation and Subsurface Stresses

The Papkovich–Neuber potentials are used to express the displacement $u_i^{\left(k\right)}$ and stress ${\sigma }_{i,j}^{\left(k\right)}$ in the coating–substrate system (k = 1 denotes the coating, whereas k = 2 indicates the substrate). The detailed procedure can be found in the reference [27]. This study provides a simplified introduction to the calculation process. The displacement and stress are calculated as follows:

$$\begin{array}{l}{u}_i^{\left(k\right)}=\frac{1}{2G_k}[{\phi }_{,i}^{\left(k\right)}+x{\psi }_{1,i}^{\left(k\right)}+z_k{\psi }_{3,i}^{\left(k\right)}-(3-4v_k){\psi }_i^{\left(k\right)}]\\ {\sigma }_{i,j}^{\left(k\right)}={\phi }_{i,j}^{\left(k\right)}-2v_k({\psi }_{1,1}^{\left(k\right)}+{\psi }_{3,3}^{\left(k\right)}){\delta }_{ij}-(1-2v_k)({\psi }_{i,j}^{\left(k\right)}+{\psi }_{j,i}^{\left(k\right)})+x{\psi }_{1,ij}^{\left(k\right)}+z_k{\psi }_{3,ij}^{\left(k\right)}\end{array},$$

(9)

where G_k represents the shear modulus of the coating or substrate, and ${\delta }_{\mathrm{i}\mathrm{j}}$ is the Kronecker delta. The Papkovich–Neuber potentials $\varphi$ and $\psi ({\psi }_1,{\psi }_2,{\psi }_3)$ are transformed using the FFT as follows:

$$\begin{array}{l}\tilde{\tilde{\varphi }}=A^{(k)}{e}^{-\alpha z_k}+{\overline{A}}^{(k)}{e}^{\alpha z_k}\\ {\tilde{\tilde{\psi }}}_1=B^{(k)}{e}^{-\alpha z_k}+{\overline{B}}^{(k)}{e}^{\alpha z_k}\\ {\tilde{\tilde{\psi }}}_3=C^{(k)}{e}^{-\alpha z_k}+{\overline{C}}^{(k)}{e}^{\alpha z_k}\end{array},$$

(10)

where $a=\sqrt{m^2+n^2}$; m and n correspond to the frequency variables in the x- and y-directions, respectively; and $A^{\left(k\right)}$,${\overline{A}}^{\left(k\right)}$,$B^{\left(k\right)}$,${\overline{B}}^{(k)}$,$C^{(k)}$, and ${\overline{C}}^{(k)}$ are unknown coefficients. According to the force transmissibility, the boundary conditions for the coating surface stress are given by Equation (11), which ignores the effect of shear stress in the y-direction.

$$\begin{array}{l}{\tilde{\tilde{\sigma }}}_{33}^{(1)}(m,n,0)=-\tilde{\tilde{p}}(m,n)\\ {\tilde{\tilde{\sigma }}}_{31}^{\left(1\right)}\left(m,n,0\right)=-{\tilde{\tilde{q}}}_x\left(m,n\right)\\ {\tilde{\tilde{\sigma }}}_{32}^{(1)}(m,n,0)=0\end{array}$$

(11)

At the coating–substrate interface, the boundary conditions are:

$$\begin{array}{l}{u}_i^{\left(2\right)}\left(x, y, 0\right)= t_{ij}^u{u}_j^{\left(1\right)}\left(x, y, h_1\right)\\ {\sigma }_{3i}^{\left(2\right)}\left(x, y, 0\right)=t_{ij}^{\sigma }{\sigma }_{3j}^{\left(1\right)} \left(x, y, h_1\right)\\ (t_{ij}^u,t_{ij}^{\sigma }\in [0,1];i,j=1,2,3)\end{array},$$

(12)

where $t_{ij}^u$ and $t_{ij}^{\sigma }$ are jumping coefficients defined to quantify the displacement and traction discontinuities across the coating–substrate interface, respectively. Due to the certain bonding strength between the coating and the substrate, the value is usually greater than 0, which can be transformed into a diagonal matrix form as follows:

$$\left[\begin{array}{l}{u}_1^{\left(2\right)}\left(x,y,0\right)\\ u_2^{\left(2\right)}\left(x,y,0\right)\\ u_3^{\left(2\right)}\left(x,y,0\right)\\ {\sigma }_{31}^{\left(2\right)}\left(x,y,0\right)\\ {\sigma }_{32}^{\left(2\right)}\left(x,y,0\right)\\ {\sigma }_{33}^{\left(2\right)}\left(x,y,0\right)\end{array}\right]==\left[\begin{array}{cc}\begin{array}{ccc}{t}_1& & \\ & t_2& \\ & & t_3\end{array}& 0\\ 0& \begin{array}{ccc}{t}_4& & \\ & t_5& \\ & & t_6\end{array}\end{array}\right]\left[\begin{array}{l}{u}_1^{\left(1\right)}\left(x,y,h_1\right)\\ u_2^{\left(1\right)}\left(x,y,h_1\right)\\ u_3^{\left(1\right)}\left(x,y,h_1\right)\\ {\sigma }_{31}^{\left(1\right)}\left(x,y,h_1\right)\\ {\sigma }_{32}^{\left(1\right)}\left(x,y,h_1\right)\\ {\sigma }_{33}^{\left(1\right)}\left(x,y,h_1\right)\end{array}\right],$$

(13)

where t₁, t₂, and t₃ are the jumping coefficients in the x-, y-, and z-directions of the dislocation-like interface, respectively, whereas t₄, t₅, and t₆ are the jumping coefficients in the x-, y-, and z-directions of the force-like interface, respectively. In this study, t_i = 0 indicates that the coating and substrate are independent of each other, whereas t_i = 1 indicates that the coating and substrate are perfectly bonded.

The frequency response function and other formulas are based on the reference [27]. The normal displacement frequency response function is as follows:

$${\tilde{\tilde{u}}}_3^{(k)}=\frac{1}{2G_k}\left\{\begin{array}{l}-\alpha [A^{(k)}{e}^{-\alpha z_k}-{\overline{A}}^{(k)}{e}^{\alpha z_k}]+im{z}_k[B^{(k)}{e}^{-\alpha z_k}+{\overline{B}}^{(k)}{e}^{\alpha z_k}]-\\ im{\alpha }^{-1}[B^{(k)}{e}^{-\alpha z_k}-{\overline{B}}^{(k)}{e}^{\alpha z_k}]-(3-4{\nu }_k)[C^{(k)}{e}^{-\alpha z_k}+{\overline{C}}^{(k)}{e}^{\alpha z_k}]-\\ \alpha z_k[C^{(k)}{e}^{-\alpha z_k}-{\overline{C}}^{(k)}{e}^{\alpha z_k}]-i\alpha [B_{,m}^{(k)}{e}^{-\alpha z_k}-{\overline{B}}_{,m}^{(k)}{e}^{-\alpha z_k}]\end{array}\right\}.$$

(14)

To obtain the normal displacement, the rectangular element of the grid node is calculated, and the discrete Fourier transform of the shape function is obtained:

$$\tilde{\tilde{Y}}(m,n)=\frac{4\sin (m{\Delta }_x / 2)\sin (n{\Delta }_y / 2)}{mn},$$

(15)

where ${\Delta }_x$ and ${\Delta }_y$ are the space grid sizes along the x- and y-directions, respectively. The shape and frequency response functions are convolved to obtain $C_{sub}^{sup}$, where sub represents the excitation and sup denotes the response.

$${\tilde{\tilde{C}}}_p^{u_i^{\left(k\right)}}={\tilde{\tilde{u}}}_i^{\left(k\right)}(m,n)\tilde{\tilde{Y}}(m,n)$$

(16)

The normal displacement frequency response function is aliased to calculate the influence coefficients (ICs) of the displacement. The ICs are converted from a continuous form to a discrete form, and the methods of discrete convolution and FFT are used as follows:

$$\begin{array}{l}{\widehat{\widehat{C}}}_p^{u_3^{\left(k\right)}}=\frac{1}{{\Delta }_x{\Delta }_y}{\displaystyle \sum _{rx=-AL}^{rx=AL}{\displaystyle \sum _{ry=-AL}^{ry=AL}{\tilde{\tilde{C}}}_p^{u_3^{(k)}}}}\left(\frac{2\pi }{M_e{\Delta }_x}i-\frac{2\pi }{{\Delta }_x}{r}_x,\frac{2\pi }{N_e{\Delta }_y}i-\frac{2\pi }{{\Delta }_y}{r}_y\right)\\ \left(-M_e / 2<i\le M_e / 2,-N_e / 2<j\le N_e / 2\right)\end{array},$$

(17)

where $M_e=2^{\gamma }M$ and $N_e=2^{\gamma }N$ are used to refine the original mesh; M and N are number of grid nodes in x- and y-directions, respectively; and $\gamma$ is used to control the aliasing level. The calculation method for the stress influence coefficient is the same as that for the normal displacement. Therefore, the formulas for deformation calculation and stress calculation are expressed as follows:

$$\begin{array}{l}{u}_3^{\left(k\right)}=C_p^{u_3^{\left(k\right)}}*p=IFFT[{\widehat{\widehat{C}}}_p^{u_3^{\left(k\right)}}\cdot \widehat{\widehat{p}}]\\ {\sigma }_{ij}^{(k)}=C_p^{{\sigma }_{ij}^{(k)}}*p+C_{q_x}^{{\sigma }_{ij}^{(k)}}*q_x+C_{q_y}^{{\sigma }_{ij}^{(k)}}*q_y=IFFT[{\widehat{\widehat{C}}}_p^{{\sigma }_{ij}^{(k)}}\cdot \widehat{\widehat{p}}+{\widehat{\widehat{C}}}_{q_x}^{{\sigma }_{ij}^{(k)}}\cdot {\widehat{\widehat{q}}}_x+{\widehat{\widehat{C}}}_{q_y}^{{\sigma }_{ij}^{(k)}}\cdot {\widehat{\widehat{q}}}_y]\end{array},$$

(18)

In the THL, since shear stress is relatively small compared with normal pressure, the influence of shear stress on the von Mises stress is ignored. The stress calculation can be simplified as follows:

$${\sigma }_{ij}^{(k)}=C_p^{{\sigma }_{ij}^{(k)}}*p=IFFT[{\widehat{\widehat{C}}}_p^{{\sigma }_{ij}^{(k)}}\cdot \widehat{\widehat{p}}].$$

(19)

2.3. Gaussian Surface Roughness

In nanoscale film lubrication, surface roughness is also an important factor affecting the lubrication performance [22,30,31,32,33]. Therefore, this study considers the influence of surface roughness on TFL. The surface roughness can be obtained by experimentally measuring the surface morphology or by digitally generating a rough surface geometry. In this study, a Gaussian-distributed rough surface is used as the morphology to determine the applicability. Using two-dimensional digital filtering techniques and computer simulations, rough surfaces with specific autocorrelation functions were generated. The autocorrelation function is given by:

$$\delta \left({{\displaystyle \tau }}_x,{{\displaystyle \tau }}_y\right)=Rm{s}^{2}exp\left(-2.3{\left({\left(\frac{{{\displaystyle \tau }}_x}{{\beta }_x}\right)}^2+{\left(\frac{{{\displaystyle \tau }}_y}{{\beta }_y}\right)}^2\right)}^{\frac{1}{2}}\right),$$

(20)

where Rms is the root mean square of the surface roughness, ${{\displaystyle \tau }}_x$ and ${{\displaystyle \tau }}_y$ are decay length in the x- and y-directions, respectively; ${\beta }_x$ and ${\beta }_y$ correspond to the relevant lengths in the x- and y-directions, respectively. ${\beta }_x={\beta }_y$ means that the rough surface is isotropic; otherwise, it is anisotropic. The specific rough surfaces can then be generated quickly using Fourier transform. The detailed process of the surface morphology generation can be found in [34,35]. To study the effect of surface roughness on the lubrication performance, the values of the correlation lengths were set as ${\beta }_x$ = ${\beta }_y$ = 2, 4, 6 µm, and the surface roughness morphology is depicted in Figure 2. The effect of the rough surface of the sphere is ignored in this study.

3. Methods and Verification

3.1. Numerical Method

In this study, to avoid convergence and accuracy issues in the numerical solution of the lubrication control equations, a nondimensionalization process on the lubrication control equation is first performed. The dimensionless parameters are selected as follows:

$$X=\frac{x}{a}, Y=\frac{y}{a}, H=\frac{h{R}_x}{a^2}, \epsilon =\frac{{\rho }^{\ast }{H}^3}{{\eta }^{\ast }\lambda }, \lambda =\frac{12{\eta }_{0}U{R}_x^2}{a^3{P}_H}, {\rho }^{\ast }=\frac{\rho }{{\rho }_0}, {\eta }^{\ast }=\frac{\eta }{{\eta }_0},P=\frac{p}{P_H},$$

where a is the Hertzian contact radius of the circular contact model, R is the synthetic radius of curvature in the x-direction, and P_H is the maximum Hertzian contact pressure.

The dimensionless lubrication equations can be obtained by applying the dimensionless parameters. The lubrication equations are discretized using the finite difference method. The three-dimensional calculation domain is $256\times 256\times 200$, whose calculation range is specified as $-2\le x/a\le 2 ,-2\le y/a\le 2$, and $0\le z/a\le 2$. The corresponding element dimensions are ${\Delta }_x/a=0.0156, {\Delta }_y/a=0.0156, {\Delta }_z/a=0.01$, where ${\Delta }_z$ is the grid size in the z-directions. In the computation program, the basic computation parameters are inputted, and the influence coefficients, which are used to calculate the surface deformation and subsurface stress, are then obtained using the Papkovich–Neuber potentials and FFT. Next, the lubricating film thickness, equivalent viscosity, and density parameters are calculated, and the Reynolds equation is then solved to obtain the lubricant film pressure distribution. Subsequently, the pressure convergence is checked, and if it converges, the load calculation and load balance judgment are performed. If the load is not balanced, the rigid contact clearance is adjusted, and the next calculation cycle is initiated. When both the load and pressure reach the computational accuracy, the cycle ends, and the subsurface stress calculation results are calculated. To satisfy the pressure and load convergence requirements, the pressure and load convergence accuracies are set to ${\epsilon }_P=1.0\times {10}^{-4}$ and ${\epsilon }_w=1.0\times {10}^{-4}$, respectively. Among them, $\overline{P}$ is the dimensionless pressure obtained from the previous calculation, w_new is the calculated load.

In the actual calculation, the convergence accuracy is often higher than the value shown, and load convergence is often reached when the full set of pressures converges; thus, the calculation accuracy can be guaranteed. The calculation process is displayed in Figure 3.

3.2. Model Validation

In order to verify the validity and correctness of the numerical calculation method in this study, we compared the calculation results with those in the literatures [27,28]. It is set that the radius of ball is 19.05 mm, elastic modulus of ball material is 210 GPa, and Poisson’s ratio is 0.3; that the elastic modulus of the substrate is 210 GPa and Poisson’s ratio is 0.3; and that the thickness of the coating is 0.5a, the elastic modulus of the coating is 52.5 GPa, and Poisson’s ratio is 0.3. The applied load was 50 N and the average velocity was 0.6 m/s. The Hertz contact radius and Hertz contact pressure were 0.1836 mm and 708 MPa, respectively. The viscosity of the liquid lubricant was 0.0262 Pa∙s and the viscosity coefficient was $1.25\times {10}^{-9}{ \mathrm{Pa}}^{-1}$. The Barus viscosity model [36] was applied in [27,28], so this viscosity model is replaced the Roelands viscosity model for verification of this study proposes the calculation.

To verify the stress change in the imperfect interface, we set the average speed to U = 0.01 m/s, friction coefficient to ${\mu }_f=0.5$, the elastic modulus of the coating is E₁/E₂ = 2.0, and the jumping coefficients are t₁ = t₂ = t₃ = t₄ = t₅ = t₆ = 0.5. Based on the above conditions, the obtained pressure distribution was similar to that under dry contact. According to the stress change results in [27], the validation results are presented in Figure 4a. Considering the bond between the coating and substrate, the dislocation-like interface [28] was also used to verify the model. The average velocity of the upper and lower surfaces was 0.6 m/s. The elastic modulus of the coating is E₁/E₂ = 0.25, and the jumping coefficients are t₁ = t₂ = 0.5, t₃ = t₄ = t₅ = t₆ = 1.0. The pressure and film thickness were plotted, as shown in Figure 4b. Here, the dotted line represents the data from published paper [28] and the solid line denotes the results of the current model.

As depicted in Figure 4a, the stress at the substrate interface is equivalent to half of the coating interface stress. Moreover, the change trend of the stress tensor matches with the result in [27]. Furthermore, the predicted numerical results of the pressure and film thickness under the interface condition considering dislocation are also consistent with the published results, as shown in Figure 4b. The film thickness and pressure distribution in [28] are highly coincident with the current model calculation results. The comparison between the calculation method presented in this study and the two published cases is able to sufficiently demonstrate that the proposed calculation method is effective and correct, so that it can be applied to investigate the effect of the EDL in point contact with imperfect coating on thin film lubrication performance.

4. Results and Discussion

In this study, a smooth sphere with a radius of 6.35 mm was used, where the elastic modulus and Poisson’s ratio of the sphere were 320 GPa and 0.26, respectively. The Poisson’s ratio of the substrate was 0.3 and the elastic modulus was 210 GPa. The Poisson’s ratio of the coating was 0.3 and the applied load was 1 N. The maximum Hertz contact pressure P_H = 450.41 MPa was calculated using the elastic modulus of the substrate, and the Hertz contact radius was a = 0.0326 mm. The applied load and substrate remained unchanged, and the coating thickness was dimensionless.

The average velocity of the upper and lower surfaces along the x-direction was set to 1 m/s. The viscosity of the lubricant was 0.002 Pa∙s, and the pressure–viscosity coefficient was $1.2\times {10}^{-8} {\mathrm{Pa}}^{-1}$.

Table 1. Model parameter.

Parameter	Value	Parameter	Value
Radius of sphere (Rb)	6.35 mm	Maximum pressure of Hertz contact (PH)	450.41 MPa
Poisson’s ratio of sphere $({\nu }_b)$	0.26	Radius of Hertz contact (a)	0.0326 mm
Elastic modulus of sphere (Eb)	320 GPa	Viscosity $({\eta }_0)$	0.002 Pa∙s
Poisson’s ratio of the substrate and coating $({\nu }_1,{\nu }_2)$	0.3	Pressure–viscosity coefficient	$12.00\times {10}^{-9} {\mathrm{Pa}}^{-1}$
Elastic modulus of the substrate (E2)	210 GPa

lists the specific parameters. The absolute electrical constant of the fluid was $\epsilon =7.08\times {10}^{-10} \mathrm{F}/\mathrm{m}$, the bulk electrical conductivity was $\lambda =1.9\times {10}^{-4} \mathrm{S}/\mathrm{m}$, and the Debye length parameter was ${\kappa }^{-1}=10 \mathrm{nm}$. Furthermore, the mesh refinement parameters, $\gamma =6$ and $\gamma =4$, were used to calculate the surface displacement and material stress, respectively. The von Mises stress was normalized to the maximum Hertz contact pressure P_H.

4.1. Effect of Zeta Potential

This study mainly considers the effect of the zeta potential on the pressure and film thickness, ignoring the discontinuities of the displacement and stress at the coating–substrate interface. The jumping coefficients were set to t₁ = t₂ = t₃ = t₄ = t₅ = t₆ = 1.0 and the coating thickness to h₁ = 1.0a. The zeta potential was varied from 0 to 300 mV. The pressure and film thickness were calculated and plotted for different coatings, as depicted in Figure 5. The variation of the central film thickness with increasing zeta potential is displayed in Figure 6. Figure 7 shows the von Mises stress when the coating elastic moduli are E₁/E₂ = 0.25, 1.0, and 4.0 under surface zeta potential ${\xi }_a$, ranging from 0 to 300 mV.

Figure 5 shows the pressure and film thickness distribution under different zeta potentials and elastic modulus of coatings. There is a significant secondary pressure peak in the rigid coating (E₁ = 4E₂), as shown in Figure 5a. With the elastic modulus of the coating increases, it can be seen that pressures are increased. Figure 5b shows the variation of film thickness under the zeta potential of EDL. The generation of the electrical double layer can increase the film thickness. Based on the central film thickness variation with coating the elastic modulus and electrical double layer, as shown in Figure 6, it can obviously be seen that the increase in zeta potential significantly increases the lubricant film thickness. While the elastic modulus of the coating increases, the film thickness decreases.

As shown in Figure 7, owing to the difference in the elastic modulus of the coating, there is a significant difference in the von Mises stress along the z-direction, and the variation of stress in the rigid coating is more pronounced. Based on Figure 5a, the maximum pressure under the flexible coating is lower than that of the rigid coating, so the maximum stress under the flexible coatings is smaller than that under rigid coatings. Meanwhile, the pressure distribution of lubricant film changes are extremely small with the increase in zeta potential, so the effect of electrical double layer on subsurface stress is minimal.

The above analyses show the effects of the zeta potential and coatings on TFL. Of these, the pressure under the rigid coating is much greater than that of the flexible coating, but the increase in the zeta potential significantly improves the load-bearing capacity of the lubricant film, indicating that the influence of the electrical double layer structure on the coating is more important.

4.2. Effect of Imperfectly Bonded Interface

This section discusses the effects of the coupled imperfectly bonded interface and coating thickness and analyzes the von Mises stress under thin coating. The transfer of the normal displacement of the coating–substrate interface should be continuous during the application, so the jumping coefficient t₃ = 1.0, whereas the other coefficients varied from 0.0001 to 1.0. The effect of the EDL was considered, and the zeta potential was set as 200 mV.

4.2.1. Effect of Coating Thickness and Interface Continuity on Pressure

The pressure and film thickness curves at different jumping coefficients are compared in Figure 8. There is an obvious difference in the pressure between the jumping coefficient t₆ and other jumping coefficients in Figure 8a,c. In terms of film thickness, the trend of film thickness curve changes in Figure 8b is different from that in Figure 8d. These results indicate that the pressure is less affected by the coating–substrate interface slip and tangential stress transmission (t₄, t₅) but mainly by the continuity of the normal stress transmission (t₆).

Figure 9a shows the maximum pressure curves for coating thicknesses, ranging from 0.01a to 2.0a at jumping coefficients of t₆ = 1.0, 0.5 for different coatings. When the coating thickness is small, the maximum pressure at t₆ = 0.5 is significantly greater than the maximum pressure at the perfectly bonded interface. With the increase in coating thickness, the pressure curve with t₆ = 0.5 gradually approaches that of perfectly bonded interface and these two curves tend to be consistent. When the coating thickness is 2.0a, the maximum pressures under E₁/E₂ = 0.25, 1.0, 4.0 becomes close to 0.5P_H, 1.0P_H, and 1.5P_H, respectively.

The central film thickness curves in Figure 9b have similar trends to the maximum pressure curves in Figure 9a. However, under flexible coatings, there is a significant change in the film thickness, and with an increase in the elastic modulus of the coating, the change rate of the film thickness gradually decreases. Meanwhile, the change in the film thickness of the rigid coating in this model is more stable, which can be attributed to the influence of the electrical double layer effect on film thickness.

This indicates that as the thickness of the coating increases, the influence of the coating elastic modulus on lubrication performance becomes increasingly important. Under the rigid coating, a smaller lubrication area is generated, which increases surface pressure and reduces the film thickness. Under flexible coatings, the increase in film thickness is more pronounced.

4.2.2. Effect of Coupling Interface Continuity and Coating Thickness on Subsurface Stress

This section discusses the effect of an imperfect interface on the von Mises stress when the coating thickness varies from 0.1a to 1.0a. Figure 10 and Figure 11, respectively, display the contours on the cross-section of xoz at y = 0 for the von Mises stress. The data in Figure 10 were calculated for a perfect interface under a flexible coating, whereas those in Figure 11 were calculated for a perfect interface under a rigid coating. The coating thicknesses were set to 0.1a, 0.5a, and 1.0a.

There is an obvious stress concentration in the flexible coating with a decrease in E₁ when the coating thickness is 0.1a, as shown in Figure 10. When the coating thickness increases, the von Mises stress of the subsurface increases, owing to the increase in the elastic modulus of the flexible coating, whereas the von Mises stress in the rigid coating is mainly concentrated in the coating, as depicted in Figure 11. However, the maximum stress decreases with the increase in coating thickness in rigid coatings, which is confirmed by results in other papers [27,28,37].

The non-uniformity of material properties can lead to changes in stress distribution. When the characteristics of the same material are consistent, the change in stress is continuous. However, the differences between the coating and the substrate, especially the differences in elastic modulus, will result in significant stress changes at the junction of these two parts. Furthermore, under the condition of the perfectly bonded interface, when the coating reaches a certain thickness, the stress concentration area is mainly concentrated in the interior of the coating. At the same time, in the rigid coating, this phenomenon will cause great stress at the interface.

The effect of the dislocation-like interface on the stress can be observed clearly in Figure 12. The decrease in the jumping coefficients t₁ and t₂ leads to a stress concentration in the coating when the coating thicknesses are 0.1a and 1.0a. As the coating thickness coincides with the stress concentration, there is a significant difference in the stress between the coating and substrate. Compared with the effect of the dislocation-like interface on the stress, that of the force-like interface is more evident, as displayed in Figure 13. When the coating thickness is small, there is no significant difference in the stress between the rigid and flexible coatings, which is different from that at t₁ = t₂ = 0.001. When the coating thickness exceeds a certain range, the discontinuity of the tangential force transfer leads to a large decrease in the stress in the substrate. This phenomenon becomes evident as t₄ and t₅ decrease. Figure 14 illustrates the effect of the discontinuity in the transmission of the normal force on the von Mises stress. With a decrease in t₆, the coating and substrate become considerably independent. Therefore, the effect of lubrication on the coating is much greater than that on the substrate at discontinuous interfaces, especially under the effect of normal stress transmission. In most cases, the coating exhibits a significant stress concentration, and the maximum stress decreases with increasing coating thickness.

Figure 15 shows the effect of interface slip on the coupled force-like and dislocation-like coating–substrate interfaces. It can be seen that the slip of the coating–substrate interface reduces the von Mises stresses of the coating and substrate to varying degrees compared with those in Figure 10 and Figure 11. The decrease in stress in the flexible coating is greater than that in the rigid coating, owing to the difference in their elastic moduli. When the coating thickness is small, the von Mises stress in the rigid coating shows significant stress concentration.

4.3. Effect of Surface Roughness

This section focuses on the effects of different rough surfaces on film lubrication and ignores the discontinuities of displacement and stress at the coating–substrate interface. Therefore, the jumping coefficients were set to t₁ = t₂ = t₃ = t₄ = t₅ = t₆ = 1.0 and the coating thickness to h₁ = 1.0a. The parameters of the EDL were the same as those in the upper section. The pressure and film thickness at surface roughness values of 0, 3, 6, and 9 nm were calculated and plotted in Figure 16. Figure 17 shows the effect of the density of rough peaks on the film thickness. Figure 18 shows the von Mises stress in smooth and rough surface with different elastic modulus of the coating under the zeta potential ${\zeta }_a$ = 200 mV of TFL.

With an increase in Rms, the pressure fluctuates considerably, and the maximum pressure increases as depicted in Figure 16. At high pressures, the film thickness fluctuation is small owing to the effect of the EDL, which ensures the liquid film existence to realizes TFL state. The changes in the density and size of the rough peak significantly affect the film thickness, as displayed in Figure 17, which is consistent with the viewpoint of Prajapati [35]. As the elastic modulus of the coating increases, the pressure and film thickness fluctuations of the rigid coating are often greater than those of the flexible coating. From these results, it is deduced that the lubrication state is destroyed with an increase in the roughness peak and elastic modulus of the coating, which promotes the transition from lubrication to boundary lubrication or dry friction and leads to lubrication failure. Analyzing the von Mises stress in Figure 18, there is a significant stress concentration and deviation in the roughness surface, which is more evident in the rigid coating. Thus, a conclusion can be drawn that the formation of rough peaks leads to pressure fluctuations that increase with the Rms, the density of the rough peaks, and the elastic modulus of the coating, further affecting the von Mises stress in subsurface.

5. Conclusions

The currently developed model can quantitatively analyze the impact of electrical double layer and imperfect coatings on lubrication. By analyzing the composite effect of imperfect coating and electrical double layer under thin film lubrication, the changes in pressure, film thickness, and stress distribution are studied, and the lubrication design of precision parts can be guided based on this study. The main conclusions are as follows:

(1)

The increase in zeta potential significantly affects the film thickness of the thin film lubrication. This phenomenon forms substantially thicker films and reduces the liquid shear stress.

(2)

A decrease in the jumping coefficient t₆ promotes an increase in the lubrication pressure. With an increase in the coating thickness, the influence of the interface discontinuities on the liquid film pressure decreases gradually. When the coating thickness exceeds 2.0a, the pressure is primarily affected by the coating material properties.

(3)

The von Mises stress of the rigid coating is more concentrated on the coating–substrate interface than that of the flexible coating, which can easily lead to coating failure. When the coating is near the stress concentration area, the stresses under different coatings are affected by the interface continuity and exhibit considerable differences. As the coating thickness increases, the stress is mainly affected by the elastic modulus of the coating.

(4)

With the increase in the surface roughness, it can lead to an increase in the pressure peak in a small thickness range, and stress concentration distributes in the shallow coating layer.

The research results provide a theoretical basis for coating lubrication design under the EDL effect. Experimental research on lubrication performance with low friction, considering the frictional chemical reaction film, coating structure, and electrical double layer effect will be conducted in future work.