## 3 Numerical Procedure  
The hybrid numerical model consists of the suction-compression-discharge chamber thermodynamic sub-model and the pipe gas dynamic sub-model, which achieves the mutual interaction of the two systems. The former one simulates the thermodynamic processes based on an integral approach. The resolution of the latter one follows two different numerical techniques depending on points considered (internal points or end-points of pipe). The flow properties of internal points are updated by solving the fundamental conservation equations using the LW2 numerical scheme [33; 34]. The flow properties at end-points are computed by means of the Trapezoidal Method of Characteristics (TMOC) developed by authors [33]. The boundary conditions applied follow the [PERSON]'s quasi-steady physical models [42] for specific boundaries. For the real gas model, it is noteworthy here that the TMOC formulation for ideal gas is applied but gas properties are calculated by calling the NIST REFPROP database. An improvement of evaluating the influence of friction and heat transfer along the characteristic lines and the path-line is thus obtained due to the non-ideal formulation of the TMOC without substantially increasing the complexity and computation effort [43]. In the developed numerical model, the coupling between pipelines and cylinder is realized through estimating the temperature and pressure in the suction and discharge chambers. For each time step, the independent points of the pipelines are solved by using the LW2 numerical scheme, but the boundary condition connected to chambers must be solved simultaneously with the thermodynamic sub-model in an iterative fashion. Firstly, the mass flow rate and energy flow through the ports contiguous to chambers are used as inputs to the thermodynamic integral solution. Differential Equations (9), (10) or (25) and (20) correspond in sequence to the first law of thermodynamics, the continuity equation and valve dynamic are calculated simultaneously by using the standard 4 th order Runge-Kutta method. Specific internal energy e and gas mass m inside control volumes can be obtained and one can compute the instantaneous density by:  
\[\rho_{cv}=\frac{m_{cv}}{V_{cv}} \tag{30}\]  
The two independent thermodynamic properties (specific internal energy and density) are enough to select the rest of the properties in the control volumes handily. Subsequently, the new solutions for the temperature and pressure in the suction and discharge chambers are applied as boundary conditions for updating the flow properties of pipe end-points, which are transferred as new inputs to the thermodynamic integral solution. This re-evaluation process could be repeated as many times as required until further significant improvements are not reached. It is found, however, that application of the correcting step only once is accurate enough since the computational time grid for thermodynamic solution is usually short [18]. Through this methodology, the pressure pulsation propagation in the pipeline system is coupled into the compressor thermodynamic processes. The time-step of the numerical simulation must be properly set. It is common sense that for pipe calculation, the spatial mesh size \(\Delta x\) is usually defined by the user to achieve a suitable trade-off between accuracy and computational cost. However, the time interval \(\Delta t\) will be determined by the stability criterion of Courant-Friedrichs-Lewy (CFL) applied to each pipe in the system [41]:  
\[\Delta t\leq\frac{\Delta z}{a+|u|} \tag{31}\]where \(a\) is the speed of sound. It is general that the time-marching procedure will be most efficient when Equation (31) has an equal sign. In addition, the working processes of the reciprocating compressor are carried out in a cyclic pattern. At each cycle, the simulation must start just at the beginning of the period and finish at the end of the period.