Numerical simulation model of a reciprocating seal based on the mixed-lubrication theory

: Based on the mixed-lubrication theory, a numerical simulation model of a reciprocating seal is established. The model includes finite element, fluid mechanics, contact mechanics and deformation mechanics analyses. The hydrodynamic analysis uses elastohydrodynamic lubrication theory and inverse hydrodynamic lubrication to obtain the relationship between fluid pressure and film thickness, while the contact mechanics analysis obtains the relationship between the asperity contact pressure and film thickness. The influence of the coefficient matrix and static contact pressure is calculated by the finite element analysis. The coupling relationship among the fluid pressure, asperity contact pressure and static contact pressure is provided by the influence coefficient method of the deformation mechanics analysis. The calculation results of the numerical model are verified by experimental results, which are measured by the bench test developed by the authors.


Introduction
A reciprocating seal is a key component in modern industry and plays an important role in aerospace and automobile fields. Sealing occupies very little weight and volume in the whole mechanical system, but it is one of the most important parts to ensure the normal and safe operation of the system. Sealing failure increases the maintenance cost and many accidents are related to it, which will cause great loss. The research on sealing technology is required to improve the reliability and life of sealing.
The study on reciprocating seals dates back to the late 19th century. At that time, only a few simple rubber O-rings were studied [1]. In 1935, Gronau published his doctoral thesis on reciprocating seals [2], and more and more studies have focused on the mechanism of reciprocating seals, but there is no major progress for a long time. The first breakthrough occurred between 1944 and 1947, when Denny of Imperial College and his mentor White conducted extensive experimental and theoretical work on reciprocating seals [3,4]. In 1969, Dowson and Swales designed a special reciprocating test bench that can be used to measure the contact pressure and film thickness during the reciprocating test [5]. After that, along with the development of the computer simulation technology and numerical calculation method, a growing number of researchers have begun to study the numerical simulation of a reciprocating seal [6].
The core of the numerical calculation model is solving the Reynolds equation of fluid mechanics in the sealing region. The leakage and friction, two most important parameters for evaluating sealing performance, can be calculated from the film thickness and fluid pressure. There are two theories to solve the Reynolds equation, elastohydrodynamic lubrication (EHL) theory and inverse hydrodynamic lubrication (IHL) theory.
The EHL theory is a method in which the oil film thickness is known, and then, the oil pressure can be calculated by solving the Reynolds equation in a forward direction. A representative study using the EHL method employs a numerical model constructed by Yang and Salant, his mentor [7][8][9]. The numerical models are coupled with fluid mechanics, contact mechanics, deformation mechanics and thermal analyses. The leakage rate, the friction of the seal against the piston rod and the distribution of oil film thickness, fluid pressure and contact pressure are calculated.
The IHL theory is a method for resolving the thickness of the oil film by using the known oil film pressure. Kanters used the IHL method to analyse the rectangular sealing ring [10]. First, the finite element simulation analysis was used to calculate the static contact pressure. Since the thickness of the liquid film is on the order of micrometres, the deformation of the seal ring is small; it can be assumed that the contact stress in the absence of the liquid film state is equal to the fluid pressure, and then, the distribution of film thickness is calculated by solving the Reynolds equation in the reverse direction.
The numerical simulation model in this paper uses both EHL and IHL methods, and couples finite element, fluid mechanics, contact mechanics and deformation mechanics analyses together. In the process of model building, some factors that have little impact need to be ignored, and the following reasonable simplifications and assumptions are proposed: i. The stroke length is much larger than the contact width of the sealing lip, ignoring the transient effect and the viscoelastic effect of the sealing material, only the sealing performance under steady-state conditions is studied. ii. The piston rod is considered smooth, but the seal has a certain roughness. iii. Assume that there is sufficient oil on the air side, so the flow of the inner stroke can always be calculated. iv. The sealed geometry model is considered to be axisymmetric. v. The friction between the rough peak and the piston rod is calculated using the empirical friction coefficient [8].

Mixed-lubrication numerical model
According to the mixed lubrication theory, the asperity contact exists in the sealing lip contact area. So, the model needs to consider the effect of the asperity contact pressure; the sum of the fluid pressure and asperity contact pressure is equal to static contact pressure.

Finite element analysis
First of all, the finite element modelling of the sealing ring is required to obtain the static contact pressure and the influence of the coefficient matrix. Fig. 2 shows the finite element model of a Trelleborg sealing ring.
In order to reduce the calculation amount and to be able to see the stress state in the whole sealing ring, the axisymmetric model was chosen. As the deformation of the piston rod and groove is very small, these can be set as a rigid body. The O-ring material is Nitrile-Butadiene Rubber; its stress-strain curve was measured by the uniaxial tensile compression test, which is imported into the hyperelastic model in the ABAQUS software. The Mooney-Rivlin model was selected for modelling the strain potential energy relationship. The V-ring material is Polytetrafluoroethylene; its stress-strain curve is also measured by the uniaxial test, which is imported into the plastic model in the ABAQUS software.
After the model is established, the overfill assembly is simulated by moving the rods and grooves, and the pressure permeation method in ABAQUS is used to pressurise the sealing ring to simulate the effect of high-pressure oil.
Then, the stress-strain state of the sealing system and the contact stress of the sealing lip can be obtained by the model. Finally, the influence of the coefficient matrix is estimated by the finite element analysis; it is used to describe the relationship between the deformation of the seal contact area and the force on the joints of the seal contact area.
Before extracting the matrix, the piston rod needs to be removed so as not to affect the deformation of the sealing ring. Then, a unit force is applied to each node in the sealing contact area. After each unit force is applied, the radial displacement of each node in the sealing contact area is recorded in a vector form. The deformation vectors caused by each unit force are combined to form the influence coefficient matrix. The specific process is shown in Fig. 5.

Fluid mechanics analysis
Since the sealing model is axisymmetric and the film thickness is far less than the contact area width, the oil film is analysed by quasi-one-dimensional Reynolds equation [11] where h 0 is the thickness of the membrane at the maximum pressure point.

Forward solution of the quasi-one-dimensional Reynolds equation (EHL):
The quasi-one-dimensional Reynolds equation is reduced to The integral of both sides is Here, h 0 was calculated according to the known pressure boundary conditions, and then, the fluid pressure distribution was obtained.
Reverse solving quasi-one-dimensional Reynolds equation (IHL): Fig. 6 shows a sketch of the fluid film flow.
Differentiating (1) yields Let the point where the second derivative of the pressure is 0 be the point a; then, since the dh/dx of the point is not 0, 6ηU − 3h a 2 dp dx a = 0 Let w a = dp dx a h a = 2ηU w a Substituting (6) into (1) gives Substituting (6) and (7) into (1) gives A = (dp/dx)(h 0 2 /6ηU), H(x) = h(x)/h 0 ; the cubic equation is used to find the root formula to get H(x).

Contact mechanics analysis
According to the mixed lubrication theory, when the oil film thickness is less than or equal to three times the RMS roughness value, the distance between two surfaces becomes too small, and the asperity contact stress generated by the asperity contact must be considered [12]. Based on the classic Greenwood-Williamson (G-W) model, the asperity contact pressure of the sealed contact zone is analysed. Under the joint action of the asperity contact pressure and oil film pressure, the interface of the seal ring reaches the balance of film thickness and pressure.
The G-W model simplifies the contact of two rough surfaces into a purely elastic contact between an absolutely smooth rigid surface and a rough surface, and considers the asperity of the rough surface as a regular hemisphere, as shown in Fig. 7, as a sketch of a single-asperity contact.
The load generated by a single-asperity contact is E′ on behalf of the two contact surfaces of equivalent stiffness is calculated as where E and υ are the elastic modulus and the Poisson ratio of the sealing ring materials. If the asperity height distribution is in line with the Gaussian distribution, the pressure generated by the asperity contact can be expressed as where n represents the asperity density on the contact surface. In addition, φ(z) is the Gaussian distribution probability density function of the wave peak height, and it is expressed as

Deformation mechanics analysis
The purpose of the deformation mechanics analysis is to couple the finite element simulation of solid mechanics with the analysis of fluid mechanics and contact mechanics so that the final film thickness and stress can be balanced. It can calculate the deformation caused by the unbalanced stress value to the contact area of the sealing ring, so it is called the deformation mechanics analysis.
In this paper, the influence coefficient method is used to analyse the deformation mechanics of the seal ring. The influence coefficient method assumes that the deformation of the seal ring can be approximated as linear when the received force is small. An additional unit force is then applied to each of the sealed contact zones, and the normal deformation of the sealed contact zone when the unit force is applied is recorded to form an influence coefficient matrix I ii . Thus, if the calculation does not reach equilibrium, the normal deformation of the sealed contact zone can be calculated from the imbalance force as δ i is the deformation vector, and p f + p con − p sc i is the unbalanced force vector.

Calculation process
According to the above-mentioned analytical process, a mixed lubrication numerical calculation model of a reciprocating seal coupled with solid mechanics, fluid mechanics, contact mechanics and deformation mechanics analyses is obtained. Friction is calculated by the following formula: The calculation process is as follows: i. Establish a finite element simulation model to calculate the contact pressure distribution p sc of the sealing lip under a given working condition and extract the influence coefficient matrix I ii . ii. Assuming that the static contact pressure distribution p sc is the same as the oil film pressure distribution p f , then the film thickness distribution h is obtained by the IHL method. iii. Calculate the value of p f + p con − p sc . If the value is zero, it means that the balance between p f + p con − p sc is reached, and the friction and the amount of leakage are calculated from the obtained results. If there is any difference, further calculations are needed. iv. Using the influence coefficient matrix I ii and the unbalanced stress p f + p con − p sc to calculate the deformation δh of the lip contact area, the film thickness h = h + δh is corrected, and the third step is recalculated.

Experimental verification on the test rig
In order to verify the simulation model, some experiments were carried out on the high-pressure reciprocating seal test rig built by the authors. Fig. 8 shows the system of the test rig. The test rig consists of three parts: a test cylinder and a driving cylinder, a high-pressure oil supply system and a hydraulic system for the driving cylinder. Fig. 9 shows the structure of the test cylinder. The piston rod is fixed to the frame by a force sensor on the left side, and the test cylinder reciprocates. The purpose of this design is to remove the effects of inertial forces generated when the piston rod moves. Fig. 10 shows the test result of the reciprocating seal test.
The ordinate value is the result of the force sensor, which represents the sum of the internal stroke friction and the external stroke friction.
Then, the reciprocating seal numerical simulation model established above is used to simulate the test process under the same working conditions. The calculated results are compared with the test results, as shown in Fig. 11.

Conclusion
A numerical simulation model of a reciprocating seal under high pressure has been established, which couples solid mechanics, fluid mechanics, contact mechanics, and deformation mechanics. Solid mechanics is analysed by finite element software, and any other type of seal can be analysed by changing the model in the software. The quasi-one-dimensional Reynolds equation is solved, which is based on EHL and IHL theories, for fluid mechanics analysis. The contact pressure of asperities on the sealing lips, which is calculated by the classic G-W model, is also considered. Then, unbalanced stress is eliminated by the influence coefficient method. The model is suitable for all types of reciprocating seals and can be calculated at a high pressure of 35 MPa.