Numerical algorithm for fluid–solid coupling in reciprocating rod seals

Journal Pre-proof Numerical algorithm for fluid–solid coupling in reciprocating rod seals Chong Xiang, Fei Guo, Xiang Liu, Yijie Chen, Xiaohong Jia, Yuming Wang PII:

S0301-679X(19)30592-4

DOI:

https://doi.org/10.1016/j.triboint.2019.106078

Reference:

JTRI 106078

To appear in:

Tribology International

Received Date: 26 July 2019 Revised Date:

5 November 2019

Accepted Date: 15 November 2019

Please cite this article as: Xiang C, Guo F, Liu X, Chen Y, Jia X, Wang Y, Numerical algorithm for fluid– solid coupling in reciprocating rod seals, Tribology International (2019), doi: https://doi.org/10.1016/ j.triboint.2019.106078. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Numerical

algorithm

for

fluid–solid

coupling

in

reciprocating rod seals Chong Xianga,b, Fei Guoa,b*, Xiang Liuc, Yijie Chend, Xiaohong Jiaa,b, Yuming Wanga,b a

State Key Laboratory of Tribology, Tsinghua University, Beijing 100084, China

b

Joint Research Center for Rubber and Plastic Seals, Tsinghua University, Beijing 100084, China

c

Nanyang Cijian Auto Shock Absorber Co., LTD., Nanyang, 473000, China

d

Science and Technology on Vehicle Transmission Laboratory, China North Vehicle Research

Institute, Beijing 100072, China *Corresponding author. E-mail address: [email protected]

ABSTRACT By analyzing the principle of pressure balance in the sealing area, a more stable and faster convergence fluid–solid coupling numerical algorithm is proposed. The overall flexibility matrix of fluid film and asperities is extracted from both fluid and contact mechanics analyses. Then, multiplying the flexibility matrix by the current unbalanced pressure yields the required film thickness correction. Through some example analysis, it can be found the present model can be set with a larger relaxation factor (0.1) or even do not need any relaxation, while conventional model usually need to sacrifice efficiency and choose a smaller relaxation factor (0.001) to achieve stable convergence, meaning that the stability of the present model can be ensured at a higher convergence rate. Keywords: Flexibility matrix, Film thickness correction, Fluid–solid coupling, Reciprocating rod seals 1. Introduction In hydraulic systems, reciprocating rod seals play a critical role in ensuring reliability in function and reducing pollution by preventing hydraulic fluid escaping into the surroundings. Nowadays, with working conditions of reciprocating rod seals becoming increasingly stringent, the sealing systems designed by qualitative forecasting methods have been difficult to meet the actual requirements. Numerical simulation technologies thus have been more widely used in the design of sealing systems. However, establishing an accurate simulation model of reciprocating rod seals is 1

not easy because of the complexity in the sealing mechanism. From the experimental results of Karaszkiewicz [1], when the rod is in motion a thin fluid film would adheres to the rod surface from hydrodynamic effects with the seal ring floating on the fluid film. When the rod moves toward the air side, called the outstroke, hydraulic oil leaks out of the cylinder; when the rod returns toward the fluid side, called the instroke, oil is pumped back into the cylinder. Therefore, the average leakage is actually equal to the flow during the outstroke minus that during the instroke. Hence, to calculate the average leakage, the fluid film thickness of the two strokes must be calculated. For this purpose, the stress–strain state of the seal system, the flow state of the fluid film, and the contact characteristics of the sealing interface all should be analyzed. With the benefit from previous research in solid mechanics, fluid mechanics, and contact mechanics, these physical processes mentioned above all can be modeled, but coupling them is very difficult because these processes are highly nonlinear. Fortunately, during the last decade or so, many sophisticated numerical simulation models has been developed [2] and two coupling approaches have emerged: the inverse hydrodynamic lubrication (IHL) theory proposed by Kanters, Nau, and Nikas et al. [3–6], and the elastohydrodynamic lubrication (EHL) theory by Salant et al. [7–10]. In IHL theory, it is assumed that the pressure applied directly by the sealing ring to the rod when the fluid film is not present is the same as the pressure applied to the fluid film when the seal ring is floating on it because the fluid film is very thin. Through solid mechanics analysis, the static contact pressure can be obtained, which means the fluid film pressure also can be determined. Then, the film thickness distribution is obtained by solving the Reynolds equation in reverse. This coupling method reasonably ignores the change in stress of the sealing ring so that the film thickness distribution can be calculated directly, promising a high calculation efficiency. Nevertheless, this method is of low accuracy, because the contact between surface asperities, which was shown to be present by Yang [11] and Crudu [12], is ignored. Different from IHL theory, EHL theory is based on mixed lubrication. It first realizes the forward solution of fluid film pressure and asperities contact pressure and then the thickness distribution of the film when the fluid pressure plus asperities contact pressure is equal to the static contact pressure need to be calculated. 2

Nonetheless, the film thickness at pressure equilibrium is very difficult to calculate because the changes in all the three pressures with film thickness are highly nonlinear. With regard to this problem, Yang [11] proposed the influence coefficient method to correct the film thickness during each calculation cycle. This correction is obtained by multiplying the unbalanced pressure by the flexibility matrix of the sealing ring. In addition, Huang [13] applied the pressure balance relationship to calculate the asperities contact pressure, after the fluid film pressure and static contact pressure are obtained. Then he inversely solved the film thickness distribution given the asperities contact pressure. Finally, he took the obtained fluid film thickness as the input condition for the next iteration, and repeated the above process until the film thickness converged. However, both of these methods have drawbacks. When the film thickness changes, the static contact pressure, fluid film pressure, and asperities contact pressure also changes, so the film thickness correction should be equal to the unbalanced stress divided by the sum of the stiffness of the sealing ring, fluid film, and asperities. However, Yang and Huang only used the stiffness of the sealing ring and rough peak, respectively. This would result in a larger correction in the film thickness, leading to oscillating or even divergent results. Beyond that, the distribution of stiffness of the three parts may not be consistent, making convergence of the computation described above difficult. In the present study, a high stability fluid–solid coupling numerical algorithm for reciprocating rod seals is established. Through the analyses of solid mechanics, fluid mechanics, and contact mechanics, the static contact pressure, fluid film pressure, and asperities contact pressure are determined. In addition, the total flexibility matrix of the seal system is extracted. According to the current flexibility matrix and unbalanced pressure, the film thickness correction can be obtained. The corrected film thickness is carried into the next step of calculation, and the iteration continues until convergence of the film thickness is achieved. 2. Decoupling analysis of sealing problem For this study, the O-ring (Fig. 1), the most commonly used reciprocating seal, was chosen on which to perform numerical simulations. In the model, the sealing performance parameters are obtained by analyzing the lubrication state of the sealing interface. However, as the lubrication process is a complicated fluid-solid coupling 3

problem, it is very difficult to solve the problem directly. The current solution is to analyze the single physical process involved in the problem first and then consider the coupling between each physical field. The sealing interface lubrication problem involves the deformation of the sealing ring, the flow of the lubrication medium and the contact of the asperities, and these physical processes would produce the static contact pressure, fluid film pressure, and asperities contact pressure in the sealing zone, which can be obtained from analyses of the solid mechanics, fluid mechanics, and contact mechanics, respectively.

Fig. 1. Schematic of the O-ring seal and sealing zone.

C10 , C01 , d

Mooney-Rivlin coefficients

dR

Rod diameter

p pa pcon psc

dc

Height of asperities

p

Er

Young’s modulus of the rod

q$

Dimensionless flow rate

Es

Young’s modulus of the seal ring

R

Average radius of asperities

U

Rod velocity Dimensionless axial coordinate,

Nomenclature

E' F

Equivalent modulus Cavitation index

h

Fluid film thickness

hA h0

Young’s

Dimensionless

h$ T

Dimensionless film thickness

I

Influence coefficient matrix

Kcon

Stiffness of asperities

thickness, h /σ

fluid

film

truncated

Asperities contact pressure Static contact pressure

p / pa

x$

x/L

φ xx The film thickness at the φs .c. x inflection point of pressure φ f , φfss , φfpp The film thickness at Φ

h$

Ambient pressure

Dimensionless fluid film pressure,

z

maximum pressure

Fluid film pressure

The height of asperities Pressure flow factor Shear flow factor Shear stress factor

µ

Average density/pressure function Lubricant viscosity RMS roughness

σ τ$ υr υs η

σ R1/3η 2/3

σ

4

Dimensionless fluid shear stress Poisson’s ratio of the rod Poisson’s ratio of the seal ring Density of asperities

Ksc Kf

Stiffness of seal ring

L

Length of the sealing zone

Stiffness of fluid film

R1/3η 2/3 EL / pa

ς ξ

Dimensionless rod velocity,

( µ0UL) / [( pa )σ 2 ]

2.1. Analysis of solid mechanics

To extract the static contact pressure distribution psc , a finite-element (FE) model of the O-ring sealing system (Fig. 2) was established. Because the structure and loading of the sealing system are axisymmetric, a 2D axisymmetric model is configured. This reduces both the amount of computation and more intuitively displays the stress–strain state of the sealing system. The O-ring is made of nitrile rubber (NBR), and to simulate its hyper-elastic mechanical properties, a three-parameter Mooney–Rivlin model is applied. The rod and groove are assumed rigid because their stiffness is much higher than the O-ring. By fitting the measured uniaxial nominal stress–strain curve, model parameter settings are obtained:

C10 = 2.51399769MPa , C01 = −0.269075684MPa , and d = 0.008968787398 . The preloaded state is simulated by applying a radial displacement to the rod and groove, and then the medium pressure on the oil side is applied through the pressure permeation method. Finally, the static contact pressure distribution under a given working condition is obtained.

Fig. 2. Finite-element model of the O-ring sealing system

2.2. Analysis of fluid mechanics

With a high calculation efficiency but low accuracy, IHL theory provides an initial distribution of the film thickness. As the static contact pressure is known, the 5

initial fluid film thickness is calculated by solving the Reynolds equation in reverse [14]

h3

dp f

− 6µU (h − h0 ) = 0 ,

dx

(1)

where h0 is the film thickness at maximum pressure. Differentiating yields,

h3

d2 pf dx

2

+

dp dh (3h 2 f − 6µU ) = 0 . dx dx

(2)

At the inflection point A where the second derivative of the pressure is zero, (d 2 p f / dx 2 ) A = 0 and (dh / dx) A ≠ 0 , the film thickness hA is then hA =

2 µU , wA

(3)

where wA = (dp f / dx) A . At inflection, substituting Eq. (3) into Eq. (1), h0 =

2 8µU . hA = 3 9 wA

(4)

Substituting Eq. (4) into Eq. (1), the original equation simplifies to

BH 3 ( x) − H ( x) + 1 = 0 ,

(5)

where B = (dp f / dx)(h02 / 6 µU ) and H ( x) = h( x) / h0 . Finally, the initial film thickness distribution is found by solving the above cubic equation. The next step is to calculate the fluid film pressure distribution from the film thickness distribution. Following Salant et al. [7], a dimensionless-one-dimensional Reynolds equation in which the effect of cavitation and surface roughness are considered is applied to describe the relationship between film thickness and pressure, 3 ∂ ( F Φ) dφ  ∂ ∂ (φxx h$ ) = 6ξ  1 + (1 − F ) Φ  h$ T + F s.c. x  , ∂ x$ ∂ x$ ∂ x$   ∂ x$

{

}

(6)

where F is the cavitation index, Φ average density/pressure function, and hT the truncated film thickness,  Φ ≥ 0, F = 1 and p = Φ in the liquid zone   Φ < 0, F = 0 and p = 0, ρ = 1 + Φ in the cavitation zone

(7a)

h$ h$ h$ 1 − h$ 2 / 2 . h$ T = + erf ( ) + e 2 2 2 2π

(7b)

In applying the finite volume method, Eq. (6) is discretized and linear equations 6

for the pressure of three adjacent nodes are obtained. Then the fluid film pressure distribution is obtained by solving the linear equations using the tri-diagonal matrix algorithm procedure. Finally, the dimensionless flow rate and dimensionless fluid shear stress are calculated from the pressure and thickness of the fluid film, 3 ∂F Φ q$ = −φ xx h$ + 6ξ [1 + (1 − F )Φ ] ( h$ T + φs .c . x ) , ∂ x$

τ$ =

−σξ σ h$ ∂ ( F Φ ) . (φ f − φ fss ) − φ fpp 2ς ∂ x$ ς h$

(8)

(9)

The flow factors φ xx , φs .c. x , φ f , φ fss , and φ fpp are obtained in accordance with the results of [15, 16]. 2.3. Analysis of contact mechanics Although there is a fluid film between the sealing ring and the rod, the pressure arising from the contact of the asperities cannot be ignored when the film thickness is less than three times the equivalent root-mean-square (RMS) roughness of the two contact surfaces. In this section, the asperities contact pressure is calculated using the classical Greenwood–Williamson (G-W) contact model [17], widely used in numerical calculations of seals [7–10, 18–20]. In the model, the contact pressure of each single contact between the asperities is given by W =

4 E ' R1/2 ( d c − h ) 2/3 , 3

(10)

where R and d c is the average radius and height of the asperities, respectively. Generally, this height obeys a Gaussian distribution, hence the contact pressure is ∞ 2 2 4 1 pcon = η E ' R1/ 2 ∫ e − z / 2σ ( z − h ) 3/ 2 dz , h 3 2πσ

(11)

where η , z and σ are the density, height and RMS roughness of the asperities, respectively. E ' the equivalent Young modulus,

1 1 − υs2 1 − υr2 = + , E' Es Er

(12)

where Es and υs are Young’s modulus and Poisson’s ratio of the seal ring; Er and

υr are those of the rod. Because the Young modulus of the seal ring is much smaller 7

than that of the rod, E ' is given by

E'=

Es . 1 − υs2

(13)

3. Fluid-solid coupling calculation According to the above analysis, it is found that only one unknown quantity determines the above physical processes, that is, the fluid film thickness. Thus, the way to couple these physical fields is to find the fluid film thickness when the static contact pressure, fluid film pressure, and asperities contact pressure are in equilibrium. 3.1.Correcting the fluid film thickness According to mixed lubrication theory, static contact pressure, fluid film pressure, and asperities contact pressure are all determined by the film thickness; the relationship between the three pressures and film thicknesses can be simplified as three springs (Fig. 3). When the film thickness changes, the three pressures also change, and therefore the stiffness of the three springs, denoted Kcon, Kf, and K sc , must be considered in the correction of film thickness, which is calculated from

∆h = ( p f + pcon − psc ) / ( K f + K con + K sc ).

(14)

Fig. 3. Simplified model of pressure in sealing zone

In Yang’s model [11], only Ksc is considered, and in Huang’s model [13], only

Kcon is considered; hence, the film thickness corrections calculated by their models are

8

∆h = ( p f + pcon − psc ) / K sc , in Yang's model; ∆h = ( p f + pcon − psc ) / K con , in Huang's model.

(15)

Therefore, the film thickness correction obtained by them are too large, thereby creating the following problems: (1) When

K sc <

1 1 ( K f + K con + K sc ) and K con < ( K f + K con + K sc ) , the calculated 2 2

film thickness correction is more than twice the actual required correction, which increases the deviation of film thickness from the equilibrium film thickness. (2) When the current pressure state is further from the equilibrium pressure state, the error in the thickness correction is larger, and hence when the current fluid film is much thicker than the equilibrium fluid film, the modified film thickness may be negative, which would disrupt the calculation. (3) For Huang's model, if fluid film pressure is higher than static contact pressure,

pcon is negative, then the calculation breaks down. In fact, since the dynamic pressure effect is inversely proportional to the cubic of film thickness, at twice the thickness of the equilibrium film, the dynamic pressure effect is already very weak, and when film thickness is three times roughness, asperity contact pressure is negligible. However, the purpose of fluid-solid coupling is find a thickness profile of the film when the fluid pressure plus asperity contact pressure is equal to the static contact pressure. What’s more, the equilibrium film thickness and roughness are all microns, so Pf + Pcon would drop from Psc to very low pressure in a few microns. Hence, near the equilibrium film thickness, K f + K con is on the same order of magnitude as

Psc . Moreover, because the seal ring is through the 1µ m

millimeter deformation to form a static contact pressure, K sc is on the same order of magnitude as

K sc

Psc , so in the range of film thickness that we need to find, 1mm

K f + Kcon . What’s more, whether the stiffness of asperities is higher than that

of the fluid film depends on the seal ring structure and actual operating conditions. Therefore, without the relaxation factor, both Yang’s and Huang’s models would experience convergence problems. Moreover, consider the second and third points, the initial film thickness must be chosen with great care to avoid negative fluid film 9

thickness and negative pressure. With both K f and Kcon considered, the film thickness correction of our model is

∆h = ( p f + pcon − psc ) / ( K f + K con ). Because K sc

(16)

K f + Kcon , the calculated correction is very close to the theoretical

deviation. Hence, theoretically speaking, our model is more reasonable and calculates the fluid film thickness under balanced pressure more quickly and stably. Based on the analyses of fluid mechanics and contact mechanics, the fluid film pressure and asperities contact pressure at any fluid film thickness can be calculated. When the fluid film thickness changes, both fluid film and asperities contact pressures change (the change in static contact pressure is ignored), and the amount of change in pressure per unit change in film thickness represents the stiffness of the sealing system. Therefore, the total stiffness matrix can be obtained by the following steps: first, apply a unit increment to the film thickness of each node successively; second, calculate the sum of the fluid film pressure and asperities contact pressure after the film thickness changes; finally, extract the pressure increment after each change to obtain the total stiffness matrix K . Finally, the total flexibility matrix I is obtained by inverting the total stiffness matrix; the fluid film thickness correction ∆h can be then be calculated (Fig. 4).

10

Fig. 4. Flow chart of the calculation for the fluid film thickness correction. Boldface lowercase letters represent column vectors, boldface uppercase letters represent matrices; other letters represent scalars.

3.2. Computational procedure In outline (see Fig. 5), the computational procedure is as follows: (1) Establish the finite-element simulation model of the reciprocating sealing system and obtain the static contact pressure psc . (2) Input the static contact pressure into the EHL theory to get the initial film thickness h . (3) Input the film thickness h into the IHL theory and the G-W contact model to obtain the fluid film pressure p and asperities contact pressure pcon . (4) Calculate the total flexibility matrix I . (5) If psc − p − pcon < ε1 , output the calculation result, otherwise, calculate the film thickness correction ∆h , and repeat Step (3) with the new film thickness.

11

Fig. 5. Computational procedure of the algorithm

4. Results and discussion Table 1 lists the basic parameters and settings for the simulated sealing system. To reflect the practicality and advantages of this algorithm presented in this paper, the calculation principles and results of different models are compared. Table 1 Basic parameters and settings Seal material

NBR

Elastic modulus of NBR, Es

43 MPa

Poisson’s ratio of NBR, υs

0.49

Rod diameter, dR

46 mm

Hydraulic oil viscosity, µ

0.0387 Pa s

Dry friction coefficient,

fc

0.1

Asperity radius of V-ring, R

4 µm

Asperity density of V-ring, η

1 × 1011

Atmospheric pressure,

pa

0.1 MPa 12

Sealed pressure,

psealed

0.1 MPa(outstroke),1.0 ~ 7.0 MPa(instroke)

0.1 ~ 0.7 m / s

Rod velocity, v

4.1. Calculation without relaxation factor With K sc

K f + Kcon , Yang's model as was noted earlier does not converge

without a relaxation factor. Furthermore, the final convergence results of the different models are identical as the only difference between the models is the coupling method. Therefore, only the simulation processes of the present model and Huang's model are compared. Simulations under different medium pressures and rod velocities were performed; the convergence condition for the calculation is that the unbalanced pressure is below 1 × 10 −5 MPa. Fig. 6 shows the results of the static contact pressure, fluid film pressure, and asperities contact pressure under a medium pressure of 5 MPa (instroke: 5 MPa, outstroke: 0.1 MPa) and a reciprocating rod velocity of 0.1 m/s. Cavitation is observed at the end of the outstroke, agreeing with earlier results [7].

(a)

(b)

Fig. 6. Simulation results of a 0.1 m/s stroke: (a) 5 MPa, instroke; (b) 0.1 MPa, outstroke.

4.1.1. Calculation process of the present model without relaxation factor Fig. 7 shows the unbalanced pressure at different stages in the correction during the convergence process. Evidently, even the present model does not eliminate unbalanced pressures after the first correction. This is because the model only uses the stiffness in the current state to correct the film thickness. However, the sealing system is highly nonlinear, and the stiffness of both the fluid film and asperities is not 13

constant, but varies with film thickness. Nevertheless, with the initial film thickness obtained by IHL method, the initial film thickness is close to the equilibrium film thickness, which means the initial unbalanced pressure is low, and therefore the film thickness correction is small. When the film thickness changes little, the stiffness changes little, and then the calculation error is small, so the present model converges rapidly within ten cycles.

(a)

(b)

Fig. 7. Unbalanced pressure at different stages in the correction, 0.1 m/s: (a) 5 MPa, instroke; (b) 0.1 MPa, outstroke.

In addition, during convergence, as the unbalanced pressure falls lower and lower, the calculation error becomes smaller and smaller, and then the rate of convergence becomes quickens. Fig. 8 shows maximum unbalanced pressure under successive stages in the correction, and to show the convergence speed more intuitively, the curve is plotted on a logarithmic scale. It is obvious that the maximum unbalanced pressure decreases faster with successive corrections.

14

(a)

(b)

(c) Fig. 8. Maximum unbalanced pressure calculated by the present model without relaxation factor: (a) 0.1 m/s, 1–7 MPa, instroke; (b) 0.1–0.7 m/s, 5 MPa, instroke; (c) 0.1–0.7 m/s, 0.1 MPa, outstroke.

4.1.2. Calculation process of Huang's model without relaxation factor For Huang's model, if the film thickness calculated by EHL method is used as the initial film thickness, the problem of negative film thickness or negative pressure (Section 2.4) occurs readily, causing the computation to fail. Therefore, to investigate the convergence characteristics of this model, the film thickness calculated (Section 3.2.1) is chosen as the initial film thickness in simulations. The maximum unbalanced pressure calculated by Huang’s model (Fig. 9) is seen to rise faster indicating that the calculation results diverge. Presumably this is due to the low stiffness of the rubber material, resulting in K con <

1 ( K f + K con + K sc ) ; the film thickness correction then 2

becomes so large that the film thickness deviates more from the equilibrium film thickness.

15

(a)

(b)

(c) Fig. 9. Maximum unbalanced pressure calculated by Huang’s model without relaxation factor: (a) 0.1 m/s, 1–7 MPa, instroke; (b) 0.1–0.7 m/s, 5 MPa, instroke; (c) 0.1–0.7 m/s, 0.1 MPa, outstroke.

4.2. Calculation with relaxation factor From this above analysis (Section 3.1), without the relaxation factor, the present model oscillates because of the high nonlinearity of the system, and Huang’s model diverges even though the initial film thickness is very close to the equilibrium film thickness. Therefore, to ensure stability of the calculation, sometimes a relaxation factor must be added to the calculation. 4.2.1. Calculation process of the present model with relaxation factor The relaxation factor ensures the stability of the calculation by reducing the thickness correction of each step, but a value too small affects the calculation efficiency strongly, so a larger relaxation factor should be selected on the premise of ensuring the stability of the calculation. 16

For the present model, when the relaxation coefficient is 0.1, the calculation process is very stable, The maximum unbalanced pressure at different stages of correction (Fig. 10) shows that, on the logarithmic coordinate axis, the maximum unbalanced pressure is approximately linear with successive corrective iterations. Moreover, after each correction, the maximum unbalanced pressure is approximately 0.9 times that of the previous correction. These phenomena suggest that, in this instance, the calculation is fairly stable and the rate of convergence is independent of the operating conditions and only depends on the selection of initial fluid film thickness.

(a)

(b)

(c) Fig. 10. Maximum unbalanced pressure calculated using the present model when the relaxation factor is 0.1: (a) 0.1 m/s, 1–7 MPa, instroke; (b) 0.1–0.7 m/s, 5 MPa, instroke; (c) 0.1–0.7 m/s, 0.1 MPa, outstroke.

In addition, to compare with the Huang model, a calculation was performed with relaxation factor of 0.001. The corresponding calculation results (Fig. 11) show that 17

the change in the maximum unbalanced pressure is very similar to those when the relaxation factor is 0.1, but now the pressure ratio before and after each correction is approximately 0.999.

(a)

(b)

(c) Fig. 11. Maximum unbalanced pressure calculated by the present model when the relaxation factor is 0.001: (a) 0.1 m/s, 1–7 MPa, instroke; (b) 0.1–0.7 m/s, 5 MPa, instroke; (c) 0.1–0.7 m/s, 0.1 MPa, outstroke.

4.2.2. Calculation process of Huang's model with relaxation factor For Huang’s model, the calculation converges only if the relaxation factor is about 0.001. For this instance, the maximum unbalanced pressure (Fig. 12) oscillates at the beginning of the calculation, but later ceases; the trend in the rate of change is also different. Because of the oscillation in the early stage of the calculation, the choice of initial film thickness affects its convergence. In addition, the difference in the slope of the curve at the later stage indicates different convergence characteristics under different operating conditions. The above analysis suggests that even if the 18

relaxation factor is small, Huang's model is not stable enough. It is necessary to select an appropriate initial film thickness and further reduce the relaxation factor in some extreme conditions to achieve convergence.

(a)

(b)

(c) Fig. 12. Maximum unbalanced pressure calculated by Huang’s model when the relaxation factor is 0.001: (a) 0.1 m/s, 1–7 MPa, instroke; (b) 0.1–0.7 m/s, 5 MPa, instroke; (c) 0.1–0.7 m/s, 0.1 (MPa), outstroke.

4.3.Total computation time The results of the above examples show that the present model has better convergence and fewer iterations than traditional algorithm, but whether the calculation loss is lower also depends on the calculation amount of a single cycle. Huang's model includes two steps in the single cycle: solving the Reynolds equation and reverse solving G-W model. Reynolds equation is solved by TDMA method, so the calculation cost of first step is proportional to the number of sealing nodes. The process of reverse solving G-W model can be equivalent to finding the root of 19

higher-order equations N times, its calculation cost is also proportional to the number of sealing nodes. However, the present model needs to obtain the stiffness matrix and invert it. To obtain the stiffness matrix, the Reynolds equation needs to be solved N times, its calculation cost is proportional to the square of the number of nodes. What’s more, since the stiffness matrix is a dense matrix, the amount of computation required to invert it is proportional to the third power of the number of nodes. Therefore, the calculation amount of a single cycle in the present model is larger than that in Huang's model, and the more the number of nodes, the greater the difference. Fig. 13 shows the CPU time costed by the two models under different pressures during instrokes, and the results under other working conditions are not listed because they are very similar to those shown in Fig. 13. It can be found that the total CPU time costed by the present model is much less than that costed by Huang’s model but the CPU time costed in a single cycle is the opposite. In addition, it can be found from Fig 13(a) that in the same model, the number of cycles determines the calculation amount. According to Fig 13(b), it can be found that the total CPU time increases as the pressure increases in the case of almost the same number of cycles, this is because the greater the pressure, the higher the pressure, the larger the contact area and the more nodes.

(a)

(b)

20

(c) Fig. 13. CPU time and correction times during instrokes, .1m/s, 1–7 MPa: (a) the present model without relaxation; (b) the present model with relaxation factor is 0.1; (c) Huang’s model with relaxation factor is 0.001.

5. Conclusions A new fluid–solid coupling numerical algorithm for reciprocating rod seals was established to improve the stability and rate of convergence of the numerical simulation. In this model, the static contact pressure is obtained by solid mechanics analysis, and through fluid mechanics analysis and contact mechanics analysis, the fluid film pressure and asperities contact pressure at any film thickness is calculated. In consequence, the overall flexibility matrix of the sealing system is obtainable from which the required film thickness correction is calculated for the current state of unbalanced pressure. The film thickness is modified iteratively until the pressure in the sealing zone is balanced. The convergence process of the present model and the Huang’s model were compared for various conditions. Without the relaxation factor, the unbalanced pressure of the present model oscillates when the initial film thickness is obtained by the IHL method, but does converge in a short time (less than ten cycles). However, for Huang’s model, the calculation results are divergent no matter how small the initial unbalanced pressure is. With the relaxation factor, the present model converges very smoothly when the relaxation factor is not very small (0.1), indicating that the stability of the algorithm is ensured with a high rate of convergence. By contrast, Huang's model provide no advantages in either respect. It converges only when the relaxation coefficient is very 21

small (0.001), and even then, the calculation process may oscillate. The biggest defect of the current model is that the calculation time to invert the stiffness matrix for a single cycle is proportional to the third power of the number of sealing nodes, which would make the computational cost increases rapidly as the number of nodes increases. However, for reciprocating rod seals, an axisymmetric model can be used to reduce the number of nodes. Moreover, because of its high convergence rate, its total computation time is still far less than that of Huang's model.

Acknowledgements The work described in this paper has been supported by the National Key R&D Program of China (Grant No. 2018YFB2001001) and the National Natural Science Foundation of China (Grant No. 51735006).

References [1] Karaszkiewicz A. Hydrodynamic lubrication of rubber seals for reciprocating motion; leakage of seals with an O-ring. Tribology international 1988, 21(6): 361-367. https://doi.org/10.1016/0301-679x(88)90113-2. [2] Nikas G. Eighty years of research on hydraulic reciprocating seals: review of tribological studies and related topics since the 1930s. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 2010; 224(1): 1-23. https://doi.org/10.1243/13506501jet607. [3] Kanters A F C, Verest J F M, Visscher M. On reciprocating elastomeric seals: calculation of film thicknesses using the inverse hydrodynamic lubrication theory. Tribology transactions 1990; 33(3): 301-306. https://doi.org/10.1080/10402009008981959. [4] Nau B S. An historical review of studies of polymeric seals in reciprocating hydraulic systems. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 1999, 213: 215-226. https://doi.org/10.3403/00374960. [5] Nikas G K, Sayles R S. Study of leakage and friction of flexible seals for steady motion via a numerical approximation method. Tribology International 2006, 39: 921-936. https://doi.org/10.1016/j.triboint.2005.09.003. 22

[6] Nikas G K, Sayles R S. Modelling and optimization of composite rectangular reciprocating seals. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 2006, 220(4): 395-412. https://doi.org/10.3390/lubricants6030077. [7] Yang B, Salant R F. Elastohydrodynamic lubrication simulation of O-ring and U-cup hydraulic seals. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 2011, 225(7): 603-610. https://doi.org/10.1016/j.triboint.2018.08.044. [8] Salant R F, Maser N, Yang B. Numerical model of a reciprocating hydraulic rod seal. Journal of Tribology 2007, 129(1): 91-97. https://doi.org/10.1115/ijtc2007-44066. [9] Salant R F. Reciprocating Lip Seal Analysis. Encyclopedia of Tribology 2013: 2748-2752. https://doi.org/10.1007/978-0-387-92897-5_142. [10] Salant R F, Yang B, Thatte A. Simulation of hydraulic seals. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 2010, 224(9): 865-876. https://doi.org/10.1243/13506501jet709. [11] Yang B. Elastohydrodynamic model of reciprocating hydraulic rod seals. Georgia Institute of Technology; 2010. [12] Crudu M, Fatu A, Cananau S, et al. A numerical and experimental friction analysis of reciprocating hydraulic ‘U’ rod seals. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 2012, 226(9): 785-794. https://doi.org/10.1177/1350650112445973. [13] Huang Y, Salant R F. Numerical analysis of a hydraulic rod seal: flooded vs. starved conditions. Tribology international 2015, 92: 577-584. https://doi.org/10.1016/j.triboint.2015.07.037 [14] Müller HK, Nau BS. Fluid sealing technology principles and applications. M. Dekker; 1998. [15] Patir N, Cheng H S. An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication. Journal of lubrication Technology 1978, 100(1): 12-17. https://doi.org/10.1115/1.3453103. [16] Patir N, Cheng H S. Application of average flow model to lubrication between rough sliding surfaces. J. Lubrication Technology, Trans. ASME, 1979, 101: 220-230. https://doi.org/10.1115/1.3453329. 23

[17] Greenwood J A, Williamson J B P. Contact of nominally flat surfaces. Proceedings of the royal society of London. Series A. Mathematical and physical sciences 1966, 295(1442): 300-319. https://doi.org/10.1098/rspa.1966.0242. [18] Guo F, Jia X, Suo S, et al. A mixed lubrication model of a rotary lip seal using flow factors. Tribology International 2013, 57: 195-201. https://doi.org/10.1016/j.triboint.2012.08.008. [19] Guo F, Jia X, Lv M, et al. The effect of aging in oil on the performance of a radial lip seal. Tribology International 2014, 78: 187-194. https://doi.org/10.1016/j.triboint.2014.05.017. [20] Guo F, Jia X, Huang L, et al. The effect of aging during storage on the performance of a radial lip seal. Polymer degradation and stability 2013, 98(11) :2193-2200. https://doi.org/10.4271/841145.

24

Highlights 1. The total flexibility matrix of the system is needed to modify the film thickness. 2. The present model do not need any relaxation. 3. The calculation process is more stable than traditional method. 4. The calculation time cost by the present model is much less.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: