A Novel Approach of Studying the Fluid–Structure–Thermal Interaction of the Piston–Cylinder Interface of Axial Piston Pumps

Abstract

The friction in the swash plate type axial piston pumps is mainly influenced by the fluid film in the friction interface. The piston–cylinder interface is one of the key friction interfaces in the pumps. The film geometry is determined by the gap between the piston and the cylinder. The dimensions of the parts determine the gap geometry, and the deformation of the structure also influences the gap geometry. The fluid viscosity is strongly influenced by temperature. Thus, a novel approach of studying the fluid film, the structure, and temperature interaction is provided in this paper. A full and quick fluid–structure–thermal interaction simulation is realized. Then, a dynamic model of the piston–cylinder interface, which integrated the fluid–structure–thermal interacting effects, has been developed. Finally, an approach for calculating the extra friction force between the piston and the cylinder is provided. Compared with the measurement data, the simulation results of the axial friction force achieve a good fit. The present work allows a fast prediction and detailed support for designing the piston–cylinder interfaces.

1. Introduction

Swash plate type axial piston pumps are prevalent in fluid power systems due to its high pressure, high volumetric efficiency, and compact structure [1]. There are three key friction interfaces in the axial piston pump: the slipper–swash plate interface, piston–cylinder interface, and valve plate–cylinder interface. The wear of the friction interfaces is the primary failure mode of the axial piston pumps. In particular, high pressure and high sliding velocity are usually exerted for the piston–cylinder interface. In addition, the piston supported by the swash plate is always inclined in the cylinder block, and the lubrication condition is even worse. The fluid film in the gap between the piston and the cylinder greatly influences the lubrication performance. The main challenge of studying the interface is the complex multi-domain physical phenomena—fluid–structure–thermal interaction. The sealing and bearing gap design is generally completed without considering the fluid–structure–thermal effects. However, the elastic deformation of the structure and the variation of temperature significantly influence the fluid film. Therefore, it is desirable to develop a model which can realize the fluid–structure–thermal interaction simulation to study the influence of various factors on friction. Based on this model, the friction force can be predicted, and the design of the gap geometry can be optimized.

From the literature review, some other simulation models have been developed. Pelosi et al. [2,3] developed a fully coupled multi-body dynamic model. The model can capture the complex fluid–structure interaction phenomena by solving the Reynolds and the energy equations on a finite volume discretized domain. The simulation is relatively difficult and time-consuming, as the averaged heat fluxes are used to determine the solid parts’ temperature distribution. It is challenging to achieve a convergent temperature field in this way. The contact force between the cylinder and the piston was not mentioned in this research, however the friction force in the simulation results was much greater than the flow friction force [4]. Wang et al. [5] provided a two-dimensional adiabatic method to calculate the temperature field. This method neglected the elastic deformation of the structure and greatly simplified the calculation process. A dynamic model [6] was also developed, but the film squeeze term in the Reynolds equation was neglected. This term is considered to be one of the main sources of film’s load-bearing capacity, and it will have an effect on the dynamic analysis of the piston. Hu et al. [7,8] developed an elastohydrodynamic lubrication model of piston–cylinder interface. In this model, the deformation matrix of the bushing was obtained through ANSYS, and then the bushing’s elastic deformation was calculated by the linear superposition of the deformation matrix. However, the influence of temperature and the piston’s elastic deformation was not considered. The research also discovered that the conical bushing could reduce the friction force. This proved that the proper film geometry can improve the lubrication performance. Wang [9] developed a fluid–structure interaction model for the piston–cylinder interface. Similar to the deformation matrix, the influence matrix was proposed to calculate the elastic deformation of the structures. In this paper, it seems that only the power loss caused by the fluid’s viscous friction has been considered. In fact, the contact friction also has a great influence on power loss. In addition, Jiang et al. [10,11] investigated the chamfer and groove’s influence on the film characteristics. Song et al. [12] presented a coupled multi-disciplinary model which included the temperature field, and analyzed the influence of the thermal effect on the contact time between the piston and cylinder. The research placed extra emphasis on the influence of temperature, and proposed that there was a critical temperature beyond which the contact time will increase rapidly. Zhang et al. [13] established a mathematical model of the film, and analyzed the dynamic characteristics of the film under 20 MPa and 70 MPa. In particular, the film pressure characteristics under ultra-high pressure conditions were obtained by experimentation. Li et al. [14] established a nonlinear mathematical model of the Reynolds equation and energy equation, and obtained the film thickness and pressure field. The research discussed the influence of hydraulic system variables and structural parameters on the leakage flow of the piston–cylinder interface. Jiang et al. [15] developed a leakage model of the piston–cylinder interface and used the genetic algorithm and the finite volume method to obtain the pressure field. Jiang et al. also analyzed the influence of multiple parameters on the leakage.

In this paper, the pressure and temperature field of the oil gap are obtained by solving the Reynolds equation and energy equation simultaneously. The finite volume method realized by the MATLAB scripts is used to solve the equations. The numerical tool COMSOL is introduced to calculate the elastic deformation of the structure. In this way, the interaction simulation can be realized more quickly and easily. A dynamic model of the piston–cylinder interface integrated the fluid film, the structure, and the temperature is developed. The Newton iteration method is used to compute the shifting velocity of the piston radial movement. Then, a numerical integration algorithm is used to calculate the position of the piston. Based on the approach, a full understanding and quick prediction of piston–cylinder interface is realized for axial piston machines.

2. Modeling

2.1. Film Geometry

In an axial piston pump, there is a gap between the piston and the cylinder. The fluid film is determined by the gap geometry, which plays the role of lubrication and sealing during the piston movement. The piston supported by the swash plate is inclined in the cylinder bushing, as shown in Figure 1.

Due to the inclined piston, the film geometry in the gap is complex. As shown in Figure 1, the eccentricity of the piston cross-sections 1 and 2 are [$e_1,e_2$] and [$e_3,e_4$], respectively. Assuming that the eccentricity of the section m changes linearly with $L_m$, the eccentricity of the section m can be calculated with:

$$\{\begin{array}{l}{e}_{xm}=e_1+\frac{\left(e_3-e_1\right)}{L}{L}_m\\ e_{ym}=e_2+\frac{\left(e_4-e_2\right)}{L}{L}_m\end{array}$$

(1)

As shown in Figure 2, c denotes the average height of the gap when the piston’s axis and the cylinder coincide. When the piston is inclined, the film height at the point [m, n] can be derived:

$$h(L_m,{\theta }_n)=c-\sqrt{e_{xm}{}^2+e_{ym}{}^2}\cos \left({\theta }_n-{\psi }_m\right)+\Delta h$$

(2)

$${\psi }_m=\arctan \left(e_{ym}/e_{xm}\right)$$

(3)

where $\Delta h$ denotes the elastic deformation of the piston and the cylinder.

Similarly, the film squeeze rate can be calculated with:

$$\dot{h}\left(L_m,{\theta }_n\right)={\dot{e}}_{xm}\cos {\theta }_n+{\dot{e}}_{ym}\sin {\theta }_n+\Delta \dot{h}$$

(4)

$$\{\begin{array}{l}{\dot{e}}_{xm}={\dot{e}}_1+\frac{\left({\dot{e}}_3-{\dot{e}}_1\right)}{L}{L}_m\\ {\dot{e}}_{ym}={\dot{e}}_2+\frac{\left({\dot{e}}_4-{\dot{e}}_2\right)}{L}{L}_m\end{array}$$

(5)

where the film squeeze rate $\Delta \dot{h}$ caused by the elastic deformation is neglected.

The film height is small compared with the dimension of the piston, so the curvature of the film can be neglected. The film can be transformed into a rectangular coordinate system, as shown in Figure 3. The gray surface denotes the piston surface, and the yellow one denotes the bushing surface.

In the rectangular coordinate system, the coordinates of the point can be expressed as:

$$\widehat{x}=R{\theta }_n,\widehat{y}=L_m,\widehat{z}=h$$

(6)

where R denotes the radius of the piston.

2.2. The Reynolds Equation and Energy Equation

The hydraulic oil used in the piston pumps is assumed incompressible. Therefore, the Navier–Stokes equation of incompressible fluid can be obtained:

$$\partial \rho u/\partial t+\nabla \left(\rho uu\right)=-\nabla p+\nabla \left(\mu \nabla u\right)$$

(7)

The following assumptions can be made: (1) The flow is under the steady-state; (2) compared with the viscous force, the inertial force can be neglected; (3) as the film height is small, the pressure does not change in the direction of the film height; and (4) the velocity gradient of fluid in the direction of film height is much greater than that in other directions.

As shown in Figure 3, in the $\widehat{z}$-axis direction, as the upper boundary of the film contacts with the piston surface, the velocity of the upper boundary is the same as that of the piston. The velocity includes the reciprocating velocity in the $\widehat{y}$-axis direction and the rotation velocity in the $\widehat{x}$-axis direction. On the other hand, the lower boundary of the film contacts the surface of the bushing. As the bushing is fixed, the velocity of the lower boundary should be 0.

The Reynolds equation can be obtained based on the above conditions:

$$\frac{\partial }{\partial \widehat{x}}\left(\frac{\partial p}{\partial \widehat{x}}\frac{h_p{}^3}{\mu }\right)+\frac{\partial }{\partial \widehat{y}}\left(\frac{\partial p}{\partial \widehat{y}}\frac{h_p{}^3}{\mu }\right)=6\left(u_p\frac{\partial h_p}{\partial \widehat{x}}+v_p\frac{\partial h_p}{\partial \widehat{y}}+2\frac{d{h}_p}{dt}\right)$$

(8)

where $h_p$ denotes the film height; $u_p$ is the rotation velocity (approximate to the revolution velocity of the piston); and $v_p$ denotes the reciprocating velocity.

According to Newton’s law of viscosity, the shearing stress per unit area of the piston/bushing surface can be defined by:

$$\{\begin{array}{l}{\tau }_{pa}={\left(-\mu \frac{\partial v}{\partial \widehat{z}}\right)}_{z=h_p}\\ {\tau }_{ba}={\left(-\mu \frac{\partial v}{\partial \widehat{z}}\right)}_{z=0}\end{array}$$

(9)

where ${\tau }_{pa}$ is the axial shearing stress exerted on the piston; ${\tau }_{ba}$ denotes the axial shearing stress exerted on the bushing.

The viscous friction force exerted on the piston/bushing can be calculated with:

$$F_{pa}\left(\phi \right)={\displaystyle \iint {\tau }_{pa}\left(\phi \right)}dA$$

(10)

$$F_{ba}\left(\phi \right)={\displaystyle \iint {\tau }_{ba}\left(\phi \right)}dA$$

(11)

where $\phi$ represents the position of the piston.

After obtaining the film geometry and the pressure distribution, the temperature distribution can be obtained by solving the energy equation. According to the law of conservation of energy, suppose the effects of volume force and thermal radiation are not considered. In that case, the energy equation for the viscous fluid flow can be expressed as [5]:

$$\rho \frac{D\left(c_{p}T\right)}{Dt}=\nabla \left(k_f\nabla T\right)-\frac{T}{\rho }\frac{\partial T}{\partial \rho }\frac{Dp}{Dt}+\varphi$$

(12)

where T denotes the temperature; $c_p$ denotes the specific heat at constant pressure; $k_f$ is the thermal conductivity coefficient of hydraulic fluid; and $\varphi$ denotes the energy dissipation term.

The convection and heat conduction in the direction of the film height can be neglected. The energy equation of the film can be simplified to:

$$c_p\left(v_p\frac{\partial T}{\partial \widehat{y}}+u_p\frac{\partial T}{\partial \widehat{x}}\right)=k_f\left(\frac{{\partial }^{2}T}{\partial {\widehat{y}}^2}+\frac{{\partial }^{2}T}{\partial {\widehat{x}}^2}\right)+\varphi$$

(13)

where the energy dissipation term $\varphi$ can be calculated with [16]:

$$\varphi =\mu \left({\left(\frac{\partial v}{\partial \widehat{z}}\right)}^2+{\left(\frac{\partial u}{\partial \widehat{z}}\right)}^2\right)$$

(14)

Both Equations (8) and (13) belong to the partial differential equation, and the expression of film height is complex. Hence, it is challenging to obtain the analytical solution of these two equations. In this paper, the numerical method is used to solve the equations. Firstly, the partial differential equations are discretized by the finite volume method, and the linear system of equations is obtained. Then, the equations are solved by the relaxation method. Finally, the numerical method is realized with the MATLAB program.

2.3. Finite Volume Method

The finite volume method is used to discrete the Reynolds equation. Figure 4 represents a part of the film grid. Integrating Equation (8) within the shaded area in Figure 4, Equation (15) becomes:

$$\Delta \widehat{y}\frac{\left(p_E-p_C\right)}{\delta \widehat{x}}{\left(\frac{h_p{}^3}{\mu }\right)}_e-\Delta \widehat{y}\frac{\left(p_C-p_W\right)}{\delta \widehat{x}}{\left(\frac{h_p{}^3}{\mu }\right)}_w+\Delta \widehat{x}\frac{\left(p_N-p_C\right)}{\delta \widehat{y}}{\left(\frac{h_p{}^3}{\mu }\right)}_n-\Delta \widehat{x}\frac{\left(p_C-p_S\right)}{\delta \widehat{y}}{\left(\frac{h_p{}^3}{\mu }\right)}_s=6S\left(h_p\right)\Delta \widehat{x}\Delta \widehat{y}$$

(15)

Equation (15) can be written as:

$$a_W{P}_W+a_N{P}_N+a_E{P}_E+a_S{P}_S=a_C{P}_C+6S\left(h_p\right)\Delta \widehat{x}\Delta \widehat{y}$$

(16)

$$a_W=\frac{\Delta \widehat{y}}{\delta \widehat{x}}{\left(\frac{h_p{}^3}{\mu }\right)}_W,a_N=\frac{\Delta \widehat{x}}{\delta \widehat{y}}{\left(\frac{h_p{}^3}{\mu }\right)}_N,a_E=\frac{\Delta \widehat{y}}{\delta \widehat{x}}{\left(\frac{h_p{}^3}{\mu }\right)}_E,a_S=\frac{\Delta \widehat{x}}{\delta \widehat{y}}{\left(\frac{h_p{}^3}{\mu }\right)}_S$$

$$a_C=a_W+a_N+a_E+a_S$$

Such a linear equation can be written for any grid point. These linear equations form a linear equation system, which is solved by the relaxation method under given boundary conditions (the outlet pressure and inlet pressure in the direction of $\widehat{y}$, and the continuity condition in the direction of $\widehat{x}$).

In the same way, Equation (13) can be discretized and solved. However, to avoid negative coefficient in the discretized equation, and ensure the convergence of the iterative solution scheme, the power-law scheme should be used [2].

2.4. Fluid–Structure–Thermal Interaction

Due to the eccentricity of the piston, the high peak will appear in the pressure field, as shown in Figure 5. The pressure field Figure 5a is calculated by the fluid model, and the pressure field Figure 5b is calculated by the fluid–structure–thermal interaction model. The maximum pressure is significantly reduced. It indicates that the influence of elastic deformation and temperature on the pressure field of the film is very significant. The high pressure exerted on the piston and the cylinder will lead to elastic deformation $\Delta h_p$, as shown in Figure 5c,d. The elastic deformation cannot be neglected. In the meantime, the change of temperature can also lead to elastic deformation, which is called thermoelastic deformation $\Delta h_T$. In addition, the temperature has a significant influence on the fluid viscosity.

The pressure can lead to the elastic deformation of the structure. On the other hand, the elastic deformation of the structure changes the film height and the pressure distribution until the pressure distribution and the elastic deformation do not change. The elastic deformation due to the pressure field is instantaneous, but it is an iterative process to obtain this equilibrium state. When the pressure distribution changes, the elastic deformation changes immediately. However, the change of temperature due to the heat transfer takes time. Therefore, before achieving a thermal equilibrium state, the temperature of the structure is constantly changing. If we couple the process that the structure temperature changes from the initial temperature to steady-state temperature into the dynamic model, the calculation is quite time-consuming. Thus, assume that the temperature distribution of the piston and the cylinder has been in a steady-state. This paper uses the Roelands equation to determine the viscosity [5]:

$$\mu ={\mu }_0\exp \left\{\left(\ln {\mu }_0+9.67\right)\left[-1+\left({\left(1+5.1\times {10}^{-9}p\right)}^Z{\left(\frac{T-138}{T0-138}\right)}^{-W}\right)\right]\right\}$$

(17)

The viscosity–pressure relationship and the viscosity–temperature relationship of the fluid are shown in Figure 6. It indicates that the viscosity of the fluid is more sensitive to the temperature than pressure. As shown in Figure 6c, at T = 293.15 K, the viscosity of the fluid is 0.070 Pa·s; at T = 323.15 K, the viscosity of the fluid is 0.019 Pa·s. The viscosity dropped by 73%. On the contrary, as the pressure increases, the viscosity of the fluid increases, as shown in Figure 6b. The proper viscosity can improve the lubrication performance and reduce the leakage [9].

In this paper, MATLAB is used to solve the pressure field and temperature field, and COMSOL is used to calculate the elastic deformation. As shown in Figure 7, the thermoelastic deformation will be calculated before the fluid–structure interaction simulation.

2.5. Dynamic Model

2.5.1. The Force Exerted on the Piston

It is desirable to analyze the motion of the piston to develop a dynamic model of the piston–cylinder interface. First, the piston revolves around the cylinder axis. At the same time, the piston makes a reciprocating motion in the cylinder along the axis. The velocity of reciprocating is related to the rotation velocity of the cylinder block. As shown in Figure 8, the z-axis coordinates (this coordinate system is different from that of the film) of one point on the piston can be calculated with:

$$z=-R_p\cos \phi \tan \beta +const$$

(18)

The axial velocity and acceleration of the piston can be given as:

$$\{\begin{array}{l}\dot{z}=\omega R_p\sin \phi \tan \beta \\ \ddot{z}={\omega }^2{R}_p\cos \phi \tan \beta \end{array}$$

(19)

where $\omega$ is the angular velocity of the piston revolution; $R_p$ denotes the distribution radius of the piston; $\beta$ is the inclination angle of the swash plate; and const denotes a constant.

It is also desirable to analyze the forces exerted on the piston.

The piston is pressurized by the hydraulic fluid in the displacement chamber. The force can be calculated with:

$$F_d=\pi R^2\left(p_s-p_0\right)$$

(20)

At the same time, the piston reciprocates along the z-axis. The inertia force is derived as:

$$F_a=\left(m_p+m_s\right)\ddot{z}=\left(m_p+m_s\right){\omega }^2{R}_p\cos \phi \tan \beta$$

(21)

where $p_s$ denotes the pressure of the displacement chamber; $p_0$ denotes the pressure of hydraulic fluid in the piston pump shell; $m_p$ denotes the mass of the piston; and $m_s$ denotes the mass of the slipper.

As shown in Figure 8, a viscous friction force $F_{pa}$ applied on the piston due to hydraulic fluid flow. Thus, the forces in the direction of the piston axis can be summarized as:

$$F_n\cos \beta +F_a+F_d+F_{pa}=0$$

(22)

In Figure 8, the supporting forces of the film exerted on the piston are not shown. The supporting force of the film in the x-axis direction is defined as $F_{px}$; the force in the y-axis direction is defined as $F_{py}$. Relative to the coordinate system of the piston cross-section 1, the torque of the film exerted on the piston can be defined as $M_{px}$ and $M_{py}$. $F_{st}$ denotes the friction force exerted on the slipper. $F_{\omega }$ denotes the centrifugal force exerted on the piston, which can be calculated with:

$$F_{\omega }=\left(m_p+m_s\right){\omega }^2{R}_p$$

(23)

As shown in Figure 8, in the x–y plane, two force balance equations can be written as:

$$\{\begin{array}{l}-F_{st}\cos \phi +F_{\omega }\sin \phi +F_{px}=0\\ F_{st}\sin \phi +F_{\omega }\cos \phi +F_n\sin \beta +F_{py}=0\end{array}$$

(24)

In addition, two torque balance equations can be written as:

$$\{\begin{array}{l}-F_{st}\sin \phi L_z-F_{\omega }\cos \phi \left(L_z-L_c\right)-F_n\sin \beta L_z+M_{px}=0\\ -F_{st}\cos \phi L_z+F_{\omega }\sin \phi \left(L_z-L_c\right)+M_{py}=0\end{array}$$

(25)

where $L_z$ denotes the distance from the center of the piston ball head to the piston section 1 and $L_c$ denotes the distance from the center of the piston ball head to the piston center of mass.

The forces of the film can be obtained by integrating the pressure field. According to the Reynolds equation, the film pressure field is determined by the film height, squeeze rate, and piston position. The film height and squeeze rate can be determined by the eccentricity and the shifting velocity of the piston radial movement. Therefore, the forces of the film exerted on the piston can be assumed a function of the eccentricity e, the shifting velocity $\dot{e}$, and the piston position $\phi$. Except for the forces of the film, other forces are only related to the piston position $\phi$. Thus, Equations (24) and (25) can be written as:

$$\{\begin{array}{l}{f}_1\left(\dot{e}\right)=F_x\left(\phi \right)+F_{px}\left(\phi ,e,\dot{e}\right)=0\\ f_2\left(\dot{e}\right)=F_y\left(\phi \right)+F_{py}\left(\phi ,e,\dot{e}\right)=0\\ f_3\left(\dot{e}\right)=M_x\left(\phi \right)+M_{px}\left(\phi ,e,\dot{e}\right)=0\\ f_4\left(\dot{e}\right)=M_y\left(\phi \right)+M_{py}\left(\phi ,e,\dot{e}\right)=0\end{array}$$

(26)

When the eccentricity e and the piston position $\phi$ are known, the Newton iteration method can solve Equation (26).

As the pressure distribution of the film is obtained by the numerical method, there is no specific expression for the forces of the film. However, the Jacobian matrix in the Newton iteration method can be calculated by the difference quotient method.

After obtaining the shifting velocity at the current time, the eccentricity at the next time can be determined by:

$$e_{n+1}=e_n+\Delta t\cdot {\dot{e}}_n$$

(27)

where $e_n$ denotes the eccentricity at the current time; ${\dot{e}}_n$ denotes the shifting velocity at the current time; and $e_{n+1}$ denotes the eccentricity of the next time. $\Delta t$ denotes the time required for the piston to rotate a slight angle around the cylinder axis, which should be small.

2.5.2. The Extra Friction Force

We did not consider the contact force of the cylinder to the piston in the above analysis. This is as we believe that there will always be a film between the piston and the cylinder. However, if there is only flow friction, the calculated result of friction force does not meet the actual situation. Therefore, we propose a calculation method for extra friction force.

If the piston is always in a state of fluid lubrication during movement, the friction force exerted on the piston is mainly the shear stress of the hydraulic fluid. When the film height is small, extra friction force will act on the piston due to the influence of the surface roughness of the piston and the bushing, as shown in Figure 9. In areas where the film height is small, there will be very high-pressure peaks. This paper argues that the extra friction force is related to the pressure in the area and can be calculated with:

$$f_{ext}={\mu }_{ext}{\displaystyle \iint pdA\left(h\le h_{min}\right)}$$

(28)

where ${\mu }_{ext}$ is the friction coefficient when there is lubricating fluid between the piston and the bushing.

In this case, the total friction force acting on the piston can be expressed as:

$$f_{pa}=F_{pa}+dir\cdot f_{ext}$$

(29)

where dir represents the direction of movement of a point on the piston.

$$dir=v_p/\sqrt{v_p{}^2+u_p{}^2}$$

In this case, Equation (22) needs to be modified as:

$$F_n\cos \beta +F_a+F_d+F_{pa}+dir\cdot f_{ext}=0$$

(30)

The complete calculation process is shown in Figure 10.

3. Results and Discussion

3.1. The Influence of Temperature

The critical simulation parameters are shown in

Table 1. Simulation parameters.

Parameters	Symbol	Value
Piston outer diameter (mm) [2]	$D_p$	20.700
Bushing inner diameter (mm) [2]	$D_b$	20.724
Piston height (mm)	$L_p$	54
Swash plate inclination (°)	$\beta$	17
Hydraulic fluid dynamic viscosity at 40 °C (cst)	$\mu$	32
Friction coefficient	${\mu }_{ext}$	0.02
Minimum film height (μm)	$h_{min}$	2

.

When the temperature of the structure changes by small amounts, COMSOL can calculate the thermoelastic deformation of the piston and the bushing. As the bushing thickness is slight, assume that the bushing temperature is the same as the film. Energy dissipation occurs due to friction of the viscous fluid in the gap, making the film temperature higher than the hydraulic oil in the high/low-pressure chamber. The temperature distribution can be obtained by solving the energy equation, and the average temperature can be calculated.

The next step is to calculate the thermoelastic deformation of the piston and the cylinder block based on the temperature boundary conditions, as shown in

Table 2. Thermal analysis results.

Operating Conditions	Condition 1	Condition 2	Condition 3
The average temperature of the film (°C)	53.0	54.9	55.7
Change of the gap height (μm)	0.95	1.01	1.20
Average dynamic viscosity (Pa⋅s)	0.020	0.019	0.019

. The pressure of the condition 1 is 8.1 MPa, the pressure of the condition 2 is 12.1 MPa, the pressure of condition 3 is 15.1 MPa, and the angular speed of the three conditions is 1000 rpm [2]. It indicates that, as the temperature increases, the gap between the piston and the bushing does not change significantly, as the piston and the bushing expand simultaneously. However, the dynamic viscosity of the oil changes significantly compared with the viscosity at 40 °C.

3.2. Axial Friction Force

The results of the axial friction force refer to the measurement data of the Maha Fluid Power research center of Purdue University [2]. The simulation parameters are the same as the experiment. The comparison between the simulation results and the measurement results of the axial friction force exerted on the bushing is shown in Figure 11.

The change of temperature will influence the performance of the piston pump. There are two main effects of temperature. One is the influence on thermoelastic deformation of the structure, and the other is the influence on the dynamic viscosity of the hydraulic fluid. According to the simulation results, the average height of the gap slightly increases with the increase in temperature. However, the dynamic viscosity of the hydraulic fluid decreases rapidly with the increase in temperature. This leads to a rapid decrease in the load-bearing capacity of the film. Therefore, it is difficult to maintain fluid lubrication at a high temperature.

The simulation results are obviously inconsistent with the measurement results if only the flow friction force is calculated.

The results of the fluid–structure–thermal interaction model show a good fit compared with the measurement results. It shows that the calculation method of the extra friction force is valid.

As shown in Figure 11, the maximum friction force tends to occur at 60°. The distribution of the film pressure field at different positions in the process of one revolution of the piston is shown in Figure 12. As expected, it indicates that the maximum pressure peak in the film pressure field appears between 45° and 90°.

3.3. Leakage of the Piston–Cylinder Interface

There is a gap between the piston and the bushing, through which there will be leakage. There are two main causes for the leakage. One is the pressure difference between the chamber and the shell, which is called pressure flow. The other is the fluid’s flow, which is called velocity flow. The velocity flow is related to the axial reciprocating velocity of the piston. As shown in Figure 13, in the three operating conditions, the angular speed is the same, so the velocity flow does not change much. On the contrary, the pressure flow changes a lot. The positive leakage represents that the hydraulic fluid flows from the chamber to the shell, and the negative leakage represents that the fluid flows from the shell to the chamber. In the high-pressure chamber, as the pressure flow and the velocity flow have opposite effects, the leakage is small. On the contrary, the velocity flow plays a leading role in the low-pressure chamber, and the leakage is relatively large.

3.4. Discussion

As the above results show, the model can accurately predict the friction force of the piston–cylinder interface. As shown in Figure 5 and Figure 6, the structure’s elastic deformation and temperature have a great influence on the film pressure field and the fluid’s dynamic viscosity. Therefore, to accurately calculate the film characteristics, it is necessary to consider the influence of the elastic deformation and temperature. Compared with the fluid model, the fluid–structure model, and the fluid–thermal model, the fluid–structure–thermal model is more comprehensive. We used COMSOL’s finite element solver to reduce the difficulty of calculating the elastic deformation. We set the average film temperature as the boundary condition to calculate the thermoelastic deformation of the structure, which greatly simplified the calculation process. In addition, we also gave a novel method to calculate the extra friction force.

However, compared with other models, it is still more complicated. The process of the structures’ temperature rise is ignored, and it is difficult to calculate the characteristics of the film in the process.

The calculation process of the fluid–structure–thermal model is time-consuming, which is related to the size of the grid, the time step, and the performance of the computer. As mentioned above, we set the average temperature of the film as the boundary condition to calculate the thermoelastic deformation, skipping the process of the temperature rise. It is also easy to calculate elastic deformation with COMSOL. In other fluid–structure–thermal models, such as the Pelosi’s model, the averaged heat fluxes are used to determine the solid parts temperature distribution. This process will cost a lot of time. Therefore, the calculation of our model should be quick theoretically, under the same conditions.

4. Conclusions

In this paper, a novel approach of studying the interaction of the fluid film, structure, and temperature is provided. Compared to previous models, the present work allows a fast simulation by solving the deformation using a software module. Additionally, the film height effect is particularly considered to improve the prediction of the friction force. A complete and quick fluid–structure–thermal interaction simulation is realized. The finite volume method is used to solve the Reynolds equation and energy equation. COMSOL is introduced to calculate the elastic deformation of the structure. Then, a dynamic model of the piston–cylinder interface integrated the temperature, and elastic deformation effects are developed. The numerical methods are used to the shifting velocity and the position of the piston. Finally, a novel approach for calculating the extra friction force between the piston and the cylinder is provided. Based on the above results, the following conclusions can be obtained:

(1)

The temperature greatly influences the lubrication performance. The dynamic viscosity will drop by 73% when the temperature rises from 293.15 K to 323.25 K at the pressure of 0.1 MPa. This will lead to the decrease in the film’s load-bearing capacity and the increase in the leakage;

(2)

Insufficient film height brings extra friction force, which is related to the pressure in the area;

(3)

A dynamic model of the piston–cylinder interface integrated the temperature and elastic deformation effects can accurately predict the friction force.

The purpose of the paper is to accurately predict the friction force and the leakage of the piston–cylinder interface with a comprehensive and easy-to-implement model. The model can be used to optimize the structural parameters and shape of the piston in the future.