Numerical simulation of gas leakage through stuffing box packing of reciprocating compressors

1. Introduction

Compressors are the heart of various industrial plants (e.g., natural gas processing plants). The duty of the compressors is to transmit (from a low-pressure point to a high-pressure point) and compress gasses (Wu et al., 2022). Compressors are divided into two categories, namely dynamic and positive displacement (Wahren, 1997). The latter is split into rotary and reciprocating. The reciprocating compressors are widely used in natural gas processing plants, chemical plants, offshore platforms, liquidification, gas pipeline systems, petrochemical industries, and CNG stations (Farzaneh-Gord et al., 2015; Elhaj et al., 2008). In a typical reciprocating compressor, the piston moves forward in the cylinder and pushes gas to a smaller space. Therefore, the pressure of the working gas is increased (Bloch & Hoefner, 1996). It is critical to seal the cylinder chamber to avoid leakage of the pressurized gas (Simons et al., 2019). One of the stationary components of any reciprocating compressor for performing this duty is the stuffing box (Zahorulko, 2015; Zainal et al., 2019). The schematic of the stuffing box is depicted in Fig. 1. Each row of packing contains radial and tangential rings. The rod has a reciprocating motion, and the rings have to stick on the cup walls to avoid the leakage of high-pressure gas (Mackay, 2004).

Figure 1. Schematic of a typical stuffing box.

Many researchers utilized computational fluid dynamics (CFD) to study fluid flow in different parts of compressors. Habing (2005) simulated the flow through the valve of a typical reciprocating compressor with an open-source code (i.e., EFD-Flow). He performed an FSI analysis considering the moving valve plate and its stationary seat. He neglected viscous forces and captured the vortexes at the plate's edge. He honestly admitted his modeling was not converged. However, the obtained results were interesting, and they were in good agreement with the experimental results. Kinjo and his co-workers (2010) investigated the compression process in a reciprocating compressor utilizing their own CFD code. They applied the virtual flux method and obtained promising results. Balduzzi et al. (2013) performed a three-dimensional steady-state simulation for predicting pressure loss from the valves of a reciprocating compressor. They utilized the finite volume method through ANSYS CFX software and obtained a correlation for the flow coefficient of the valves. Balduzzi et al. (2018) exploited FLUENT in another research to simulate flow in a complete reciprocating compressor. They reported mass flow inlet and outlet during one cycle. However, they do not consider leaking. Zhao et al. (2018) performed a 3D CFD simulation. They investigated the interaction between valve motion and pressure pulsation. Mohammadi-Amin and his co-workers (2020) numerically investigated a two-dimensional chamber of a reciprocating compressor by COMSOL-Multiphysics. Also, they reported pressure distribution in the compress chamber for different frequencies of the piston. Coskun et al. (2020) numerically simulate fluid flow in a reciprocating compressor of a refrigerator with its suction muffler by FLUENT. They utilized the SIMPLE method to solve unsteady two-dimensional governing equations. Zhao et al. (2021) utilized FLUENT for 3D simulation of the pipeline and cylinder of a reciprocating compressor by the PISO algorithm. They verified the acquired results with their experiments. Wu et al. (2023) assessed the valve FSI problem with ADINA. They found promising results in their numerical studies. Pashak et al. (2023) numerically analyzed the impacts of valve lift and piston position on the discharge mass flow rate (for a typical reciprocating compressor). They exploited a three-dimensional model through ANSYS/Fluent. They showed the double-variable paradigm was the most reliable strategy for discharge loss prediction. Bacak and his co-authors (2023) utilized CONVERGE 3.0 for 3D FSI simulation. They found good agreement with the experimental results. Most recently, Braga and Deschamps (2023) numerically simulated leakage flow between the outer surface of a piston and the inner surface of its cylinder by FLUENT. They found that the leak rate highly depends on the chamber's pressure. They reported that piston speed did not affect the leakage magnitude. Through this brief literature survey, it could be concluded the stuffing box has not been considered in previous research. Therefore, such a study is motivated in this paper.

Due to its importance, the designation and manufacturing of robust and efficient stuffing boxes is aimed at Hanyag. In doing so, numerical simulating of leaking flow through stuffing box packing is studied in this study in detail. In this section, the importance of the subject is clarified. The second section provides the mathematical formulation of the problem. The third section is devoted to illustrating the details of the numerical approach. The results are presented in the fourth section. Finally, the concluding remarks are stated.

2. Governing Equations

The main goal of the present work is to study the leaking flow through the stuffing box. Two challenges are encountered in this way: the leakage passage is very narrow (Sun et al., 2020; Kovacevic & Rane, 2021), and the working fluid is compressible. According to the existing geometric dimensions, the characteristic length of the leakage path, taking into account the smoothness of the surfaces of the ring and ring cap, is about 0.02 mm, and the average free path length, for example, for air molecules under standard conditions in the range is 0.00001mm. It means the Knudsen number is smaller than 0.001, and the flow is in the continuum regime (Karniadakis et al., 2001). Therefore, the equations governing the present problem are the familiar equations of the conservation of mass, momentum, and energy, which are solved with the equation of state of the working gas. Considering that the gas inside the cylinder piston of the reciprocating compressor has a compressible behavior, the governing equations should also be considered in the compressible form. Although the real problem is unsteady, the steady simulation is performed with the maximum pressure as the worst case to be on the safe side (i.e., the designation of the proposed stuffing box will have a good reliability factor). The schematic of the problem is displayed in Fig. 2.

Fig. 2. The schematic of the problem.

The continuity (conservation of mass) for a compressible gas is as follows (Roudgar Saffari et al., 2021),

The conservations of momentum (i.e., Navier-Stokes equations) for a Newtonian fluid in three directions are written in the following form. Note that the leakage velocity is small, and the flow is laminar (White, 2000).

Noting the problem is solved in the steady state condition; therefore, the first left terms of the above equations are equal to zero. Also, the conservation of energy is stated as (Bergman et al., 2011),

where e, k, and  are specific energy, thermal conductivity, and rate of energy generation in control volume, respectively. Note that in the present study . Assuming the gas is calorically perfect, the conservation of energy could be simplified,

where c_v, and μФ are resembles specific heat and viscous dissipation, respectively. Again, the most left term is set to zero because the flow is considered to be steady. Viscous dissipation is formulated as follows,

where λ is the second viscosity coefficient which is considered to be  according to Stokes’ law (its application for polyatomic gas is almost correct). The working fluid is assumed to be an ideal gas, and the famous equation of state is utilized,

Now, a system of six coupled equations and six unknowns is obtained. The system is solved by utilizing the no-slip and no-penetration boundary conditions on the walls. 

3. Numerical Method

To solve the prescribed equations, the well-known finite volume method (FV) is applied by utilizing ANSYS/FLUENT 18.0. In the finite volume method, the domain is fully covered by a number of non-overlapping computational cells (Patankar, 1980). Then, the conservation laws are applied to all of these cells, and the obtained integrals are converted to algebraic equations by approximating the integrant to be constant over each of them (i.e., equal to its magnitude in the center of the cell). Two strategies are available for the finite volume method: density-based and pressure-based. According to initial tests, it is demonstrated that the Mach number is under 0.8 in the majority of the domain. Therefore, the pressure-based strategy is selected to have faster convergence. Also, the staggered grid is utilized to remedy the check-board problem. Indeed, the coupled scheme is applied for splitting pressure-velocity coupling. Furthermore, the second-order scheme is exploited for spatial discretization.

Two distinguished geometries are considered in this study, as shown in Fig. 3. The first geometry is utilized as a model for one-row packing rings. The second one is exploited for analyzing the effects of the second row of rings. The width of the leakage path is dependent on the surface's roughness. However, herein, a severe condition is considered that the width of the passage everywhere is set to the maximum allowable tolerance (i.e., 0.02mm). All walls are assumed to be stationary, and applying no-slip conditions is suitable. The inlet boundary (i.e., the left line in the drawing) is set to pressure-inlet in FLUENT, and the maximum pressure of the chamber is applied to it. On the other hand, the right line is set to pressure-outlet condition with atmosphere magnitude (usually, the pressure in the cross-head chamber is atmospheric air). It is assumed that the working fluid properties are similar to air, and obey the ideal gas rule. Also, the symmetric two-dimensional modeling is performed.

Figure 3. The created geometry in Design Modeler of ANSYS software. a) First geometry, b) Second geometry.

4. Results and Discussion

4.1. Verification

Before applying the numerical method, it should be verified. The problem introduced by Braga and Deschamps (2023) is selected for verification purposes. In this test case, the leakage flow between the outer surface of the piston and the inner surface of the cylinder of a small reciprocating compressor is simulated three-dimensionally. Accordingly, the pressure in the inlet and outlet sides are set to 7.6 bar and 0.629 bar, respectively. Also, the piston's radius, length, and velocity are 10.5 mm, 18.1 mm, and 3m/s, respectively. The radial clearance between the cylinder and piston is 0.005 mm. The mass flow rate of the leakage is obtained below 0.15 g/s, which is in good agreement with Braga and Deschamps (2023). Therefore, the applied approach is validated for the original problem. For a better understanding, the pressure contours are displayed in Fig. 4. Unfortunately, Braga and Deschamps (2023) only report the mass flow leakage. In other words, there is no available data for comparison of the pressure contours. However, the pressure contours are presented for future comparison by other researchers.

Figure 4. Contours of pressure between the piston and cylinder of a reciprocating compressor with p_in-p_out=7bar.

4.2. Results for one-row and two-row packing

Now, it is time to simulate the leakage flow through the stuffing box. Firstly, the geometry of the stuffing box with only one row of packing is analyzed. The pressure difference is changed from 7 bar to 70 bar, and the mass flow rate of leakage is evaluated. The acquired results are presented in Fig. 5. As seen, the leakage from the compression chamber is increased from 3.056 to 49.6 g/s by increasing the pressure difference. The cause of the flow is the pressure difference. Therefore, the obtained trend is logical. Note that utilizing the second row of packing is required to reduce the leakage.

Figure 5. The leakage rate versus various pressure differences for one-row packing.

For better understanding, the contours of the pressure, Mach number, and stream function are obtained for various simulations (i.e., various pressure differences). For instance, these contours for Δp = 70 bar are shown in Fig. 6. It is worth mentioning the prominent part of the pressure reduction occurred in the narrow passage between the ring and cap. Almost in the whole domain, the Mach number is below 0.7 (see Fig. 7b). Hence, applying the pressure-based strategy is validated. Also, the flow pattern, as understood by stream function contours, shows the creation of vortexes before and after the narrow passage (see Fig. 7c). The backflow before the passage means that the packing effectively prevents the gas leakage to the cross-head chamber. It is worth mentioning that the presented results are obtained by 439000 cells. A mesh study is done by utilizing five different resolutions (i.e., 109750, 219500, 439000, and 878000), and it is concluded using more computational cells does not improve the results. So, the results are independent from the grid (see Fig. 8).

Figure 6. Flow field contours for one-row packing with Δp = 70 bar. a) pressure distribution, b) Mach number contours, and c) stream function.

Figure 7. The mesh study for the 1^st geometry. Mass leakage versus cell numbers for Δp = 70 bar.

Secondly, the two-row packing is simulated. Again, the pressure difference is altered between 7 bar to 70 bar. The leakage rate for each case is plotted in Fig. 8. As expected, the leaking is resonated for higher pressure differences. The mass flow leakage is 1.2 g/s and 41 g/s for pressure differences of 7 bar and 70 bar, respectively. 

The distributions of various properties are acquired for different pressure differences. Among them, the flow field contours (including pressure, velocity magnitude, and Mach number distribution) for Δp = 70 bar are displayed in Fig. 9. The presence of vortexes before and after each leaking passage occurred again. The backflow before each passage verified that the packing performed its duty effectively and prevented the gas from flowing to the cross-head chamber. The results are acquired by utilizing 878000 computational cells.

Figure 8. The leakage rate versus various pressure differences for two-row packing.

Figure 9. Flow field contours for two-row packing with Δp = 70 bar. a) pressure contours, b) velocity magnitude, and c) stream function.

For comparison purposes, the leakage rates for one-row and two-row packing are given in Table 1. Furthermore, the reduction percentage caused by the second row of the rings is presented in the upper right column of this table. The results show that utilizing the second row for low-pressure differences is more efficient. For high-pressure differences, another row of packing rings is required to keep the leakage below the admissible rate (i.e., 25 g/s).

Table 1. The effect of second-row packing rings on leakage.

5. Conclusions

In this study, the three-dimensional steady flow of an ideal gas through the very thin passage of the stuffing box is numerically simulated by applying the finite volume method. Two different geometries are considered, and pressure differences are altered from 7 to 70 bar. The simulations are performed with ANSYS/Fluent commercial software. The main findings are summarized as follows,

• The leakage flow is in the range of the continuum regime, and application of the usual form of continuity, Navier-Stokes, and energy equation is admissible.
• For one row packing, the leakage mass flow rate is increased from 3.056 to 49.6 g/s by increasing the pressure difference from 7 to 70 bar. 
• The addition of the second row leads to a decrease in the leakage rate, especially for lower pressure differences. 
• Although the first row has a crucial role in leakage control, the presence of the second row is also mandatory to keep the leakage rate in allowable condition.

Experimental work is suggested for future research.

Acknowledgments

The first and second authors wish to thank the Hanyag Corporation for their constructive data and crucial support.