Application of modelling approaches of twin-screw compressors: thermodynamic investigation and reduced-order model identification

Refrigeration is an essential part of the food chain. It is used in all stages of the chain, from industrial food processing to final consumption at home. In these processes, mechanical refrigeration technologies are employed, where compressors increase gas pressure from evaporation to condensation. In industrial refrigeration systems, twin-screw compressors represent the most widely used technology. A detailed mathematical model of a twin-screw compressor has been developed in Simulink® using differential equations for energy and mass balances to simulate the compression cycle that includes suction, compression and discharge phases. Gas pressure and enthalpy can be calculated as time functions during the cycle. However, the computational times obtained limit the possibility to extend the use of the model in the development of control strategies for the whole refrigeration plant in its real operating conditions. Therefore, the detailed model has been used to train a simplified model developed in Matlab®: the simulated mass flow rate, shaft power and the fluid discharge temperature have been employed to identify several geometrical and thermodynamic parameters of the simplified model. The latter relies on nonlinear algebraic equations and, thus, requires a very short computational time. A limited performance dataset has been used to train the model, and a different dataset to test it: the results of the models have been compared, and small errors in mass flow rate, shaft power and fluid discharge temperature have been observed.


Introduction
Nowadays, refrigeration has become an important part of people's daily lives because it plays a key role in countless sectors, e.g. in the food industry, air conditioning and healthcare. In particular, refrigeration is vital for the food chain in preserving food products, and in the healthcare sector in preserving pharmaceutical products and medicines, especially vaccines. Hence, from an energy and economic point of view, the importance of the refrigeration sector is paramount. The International Institute of Refrigeration estimates that the number of refrigeration systems in operation worldwide is about 5 billion, which account for approximately 20 % of the overall electricity consumed worldwide [1]. Moreover, the global electricity demand for refrigeration is expected to double by 2050 [1]. In this context, significant energy saving opportunities can derive from the optimal management and control of refrigeration devices.
System modelling and simulation are widely used to predict performance and to optimise the system, since they allow money and time to be saved in the development of prototype units and their continuous testing. When modelling the overall refrigeration system, all components have to be modelled. One of the main components is the compressor. It represents the heart of refrigeration systems because it increases the pressure of the working fluid from evaporation to condensation. Therefore, over the years, engineers and researchers have made great efforts to study and model all types of compressor. In industrial refrigeration systems, twin-screw compressors represent the most widely used technology. Many types of screw compressor models appear in the literature. Stosic and Hanjalic developed a geometrical and thermodynamic model for the simulation of screw machines [2,3]. They presented an algorithm for generating rotor profiles [2] which are used to simulate the thermodynamic compression process and to optimise design parameters [3]. Spille-Kohoff et al. presented a 3-D CFD simulation of a dry screw compressor to investigate leakage flows and rotor heating [4], while Ding and Young presented a 3-D transient CFD model of a screw compressor with oil injection to analyse its effects on compressor performance [5]. These CFD models require a large computational time, firstly for creating machine geometry or grid to mesh, and then for solving differential equations and, therefore, they are not suitable for whole-system simulations. Chamoun et al. presented a mathematical model of the compression cycle for twin-screw compressors using spatial discretisation according to the rotational angle of the male rotor to study the dynamic behaviour of the compressor [6]. Giuffrida described a semi-empirical model of an open drive twin-screw compressor by dividing the compression process into a number of steps [7]. The latter is a very simplified model that requires a short computational time. However, this model needs the performance map provided by compressor manufacturers.
In order to create a suitable twin-screw compressor model for the simulation of entire refrigeration systems, a new procedure that links both modelling approaches has been presented in this work. In particular, a mathematical model of a twin-screw compressor has been implemented to simulate the compression cycle in detail. Simulation results of the detailed model have been used to create the compressor performance map on which a new simplified model is identified. This cross-identification procedure allows the calculation of several geometrical and thermodynamic parameters of the simplified model. Finally, the procedure has been tested and validated by comparing the results of both models using several performance datasets.

Modelling approaches
The twin-screw compressor represents the most widely used technology in industrial refrigeration systems. It is a rotary positive displacement machine and, as such, it has distinct working phases, namely suction, compression and discharge, which take place simultaneously. Compression is achieved by the intermeshing of the male and female rotors. When the male rotor starts to move out of mesh with the female rotor, the inlet port is uncovered and gas flows into the compressor until the entire interlobe space is filled (suction). As the rotor continues to turn, the volume of the gas is progressively reduced and its pressure increases (compression). Finally, further rotation uncovers the discharge port and the compressed gas starts to flow out of the compressor (discharge). The pressure increase is determined by the Built-in Volume Ratio (BVR). It is a geometrical feature of a specific compressor and it is given by the ratio between the gas volume when the loading section closes and the final volume at discharge. Hence, it determines the internal compression ratio and leads to under-or over-compression phenomena when the designed internal pressure is lower or higher than the desired discharge pressure. These phenomena affect compressor performance, increasing the required compression power. Moreover, twin-screw compressor efficiency is influenced by internal fluid leakages between two chambers which are at different stages of the compression process. In the following sections, a detailed dynamic model of a twin-screw compressor, based on differential equations for energy and mass balances, and a simplified thermodynamic model of this component are developed and described.

Detailed dynamic model
As shown in Fig. 1a, this modelling strategy divides the compression process into three phases that take place simultaneously: suction, compression and discharge. The compression process is modelled through a set of differential equations for the conservation of mass and energy to describe how the thermodynamic and flow properties vary during the compression cycle. Fig. 1b shows the control volume, which is defined in analogy with a reciprocating compressor, with its mass and energy flows. A variation of the thermodynamic properties with the rotational angle of the male rotor θ is considered. However, if ω is the compressor angular speed, mass and energy balances can be written as time functions. Therefore, the mass balance for the control volume shown in Fig. 1b can be expressed as follows: where ṁi and ṁo represent the inlet and outlet mass flow rates, respectively. The first law of thermodynamics for the same control volume can be written as follows: where U, dW /dt and Q represent the total internal energy, the mechanical power and the thermal power, respectively. Considering that: Equations (1) and (2) can be rewritten as: The effect of the heat transfer between the gas and the compressor body ̇ is very low compared with the gas internal energy. Therefore, it can be neglected as described by Ignatiev et al. [8]. Equations (7) and (8) have been particularised for each compression phase through mass and energy flows that characterise a specific process (e.g. the suction mass flow rate for the suction process). When rearranging Equations (7) and (8), it is possible to calculate gas pressure and specific enthalpy as time functions during the cycle: The model has been implemented in Simulink ® by approximating partial derivatives with forward difference quotients and calculating properties through CoolProp [9]. The time derivative of the volume is assumed to have constant values in the aspiration phase and in the compression and discharge phases.
An artifice has been used to model the three phases that occur simultaneously and have different durations, as suggested by Kauder et al. [10]. In particular, a unified angle phase δ has been used. It represents the angle for which a new chamber is created. For twin-screw compressors, it is possible to calculate this angle as follows: where Nm represents the number of lobes of the male rotor. Moreover, another particular angle β has been considered. It represents the angle for which the chamber completes its cycle and disappears. The number of chambers that have to be implemented in the Simulink ® scheme is given by the ratio between β and δ. If β is not an integer multiple of δ, the number of chambers is the closest multiple. Each chamber communicates with the next chamber transferring all thermodynamic data (pressure, specific enthalpy and volume) that then become the initial condition of the integrator of the latter when a new suction chamber is created. It is possible to model the fluid leakages between the two chambers which are at a different stage of the process by calculating the leakage mass flow rate ṁleak between the two chambers as described in [6]: where Aleak represents the leakage area. Each chamber receives as input the leakage mass flow rate calculated by the next chamber and provides as output the leakage mass flow rate to the preceding chamber. The suction mass flow rate ṁsuc, the compression power Pc and the discharge temperature Tdis are calculated by means of Equations (13) to (15): The input parameters of the model are the suction temperature and pressure, the discharge pressure, the rotational speed and some geometrical features of a specific compressor such as the volume curve, the BVR, the leakage area, and the intake and exhaust port areas. This dynamic model simulates the compression process in a very detailed way and predicts compressor performance. Nevertheless it requires a long computational time (⁓ 30 min for simulating 0.12 s) that limits its application to whole system simulation and to the development of plant control strategies. A simplified thermodynamic model of a twin-screw compressor has been developed and described in the next section in order to simulate the overall performance of the compressor in terms of ṁsuc and Pc when it is included in entire refrigeration systems.

Simplified thermodynamic model
The compression process can be modelled by dividing it into a number of steps, as described by Giuffrida for an open-drive twin-screw compressor [7] and Winandy et al. for a hermetic scroll refrigeration compressor [11]. The following modelling assumptions are made: the kinetic energy of the working fluid, the pressure drop at the inlet and discharge ports and the presence of oil in the working fluid are neglected. Thus, the working fluid only consists of the refrigerant. Moreover, fluid leakages are considered adiabatic. Fluid experiences some thermodynamic transformations flowing through the compressor, as shown in Fig. 2. The first step, from section 1 to section 2, considers the mixing of the fluid entering the compressor at the inlet port and the fluid leakages. This process leads to an increase in fluid specific enthalpy, expressed through Equation (16): where ṁtot represents the total mass flow rate: ̇t ot =̇s uc +̇l eak Fluid leakages are treated by considering that all leakage paths gather into a unique fictitious path which links the exhaust and suction lines. In order to calculate the leakage mass flow rate, the equation of an isentropic flow through a convergent nozzle is used as described in [7]: where ρleak and hleak are the leakage density and enthalpy which are calculated considering the isentropic condition (sleak = sdis) and the highest value between the suction pressure psuc and the critical pressure pcrit,leak: The second step, from 2 to 3, considers a heat transfer between the fluid and the compressor body, which is assumed to be an isothermal envelope at temperature Tw. The thermal power Q suc can be written as: where AUsuc,nom is the suction heat transfer coefficient for the nominal mass flow rate ṁnom.
The refrigerant mass flow rate that is compressed can be calculated as follows: As schematised in Fig. 2, the compression process, from 3 to 5, is divided into two steps. The first step, from 3 to 4, is related to an isentropic compression that reduces the refrigerant volume according to the BVR: Therefore, the refrigerant experiences a pressure increase up to pressure p4. The second step (i.e. 4 to 5) considers an adiabatic process at constant volume to simulate the under-or overcompression phenomena that occur when pressure p4 is lower or higher than discharge pressure pdis. These phenomena lead to some losses and to an increase in the mechanical power input Pin, expressed by the second term in the square brackets of Equation (24) as described in [7]: The last step (i.e. 5 to 6) refers to heat transfer Q dis between the refrigerant and the compressor body. If Tw is higher than T5 the refrigerant is heated; otherwise it is cooled. To model this step, Equations (20) and (21) where τvis is the torque loss due to the viscous friction, μ is the dynamic viscosity of lubricant oil and atl,1 and atl,2 are two dimensionless constants that are used to introduce a relation with the internal compression power Pin and the compressor displacement Vsw [7]. Hence, the overall compression power Pc can be calculated as follows: c = in + loss,1 + loss,2 Another loss occurs in the heat transfer between the compressor body and the external ambient: The temperature Tw can be calculated through the energy balance considering the compressor body as the control volume: loss,1 + loss,2 −̇s uc +̇d is −̇a mb = 0 Finally, the overall energy balance can be written as follows: The inputs of the model are suction pressure psuc, inlet temperature Tsuc, rotational speed n, discharge pressure pdis and ambient temperature Tamb. This simplified model is implemented in Matlab ® and, being based on a set of non-linear algebraic equations, it requires a very short computational time. Hence, it is suitable for system-wide simulation.

Identification of the simplified model parameters
In order to use the simplified model for simulation purposes, the geometrical and thermodynamic parameters Aleak, AUsuc,nom, ṁnom, Vsw, BVR, atl,1, atl,2, AUdis,nom, AUamb have to be identified. The identification of the model parameters presupposes the availability of the compressor performance maps and it is obtained through the minimisation of the following error function: in which the simplified model results, in terms of simulated mass flow rate ṁsuc,sim, compression power Pc,sim, and discharge temperature Tdis,sim, are compared to the compressor map data (ṁsuc,cat, Pc,cat, Tdis,cat). It is a constrained nonlinear multivariable function where the variables are the model parameters. In order to solve this optimization problem, the Matlab ® function fmincon is used. It finds the minimum of constrained nonlinear multivariable functions by using the interior-point algorithm, which is a large-scale algorithm that satisfies bounds at all iterations.
Therefore, by choosing a number of datasets Nd (nine at least) with different external operating conditions (suction and discharge pressures) on the compressor maps, the model parameters are identified.
In this section, the results of the detailed model and the procedure for the identification of the simplified model parameters are presented. Finally, a comparison between the results of the two models is made in order to validate the cross-identification procedure.

Simulation results
The detailed dynamic model calculates gas pressure and enthalpy as time functions through Equations (9) and (10). In Table 1, the input data for the detailed model are listed. The leakage and port area are assumed to be constant and the fluid used for these analyses is ammonia.
The evolutions in time of volume, pressure and temperature in one of the compression chambers that form the Simulink ® model are shown in Fig. 3. The gas pressure increase is due to the decrease in volume, according to the operating principle of twin-screw compressors. The discontinuity of the volume occurs when a new chamber is created and it is due to the eulerian approach used for the detailed model implementation (i.e. the process continues in the next chamber). Moreover, Fig. 3 highlights that the pressure increase leads to an increase in the compression chamber temperature. The most interesting results of the model are suction mass flow rate, mechanical compression power and fluid discharge temperature, as represented in Fig. 4. Indeed, in addition to defining the performance of the compressor, these results represent the parameters of the simplified model in the crossidentification procedure described in the next section.

Cross-identification procedure
In order to overcome the limits of both models discussed in the previous section, a crossidentification procedure has been fine-tuned. In the absence of the compressor performance maps, the detailed model has been used to create them, by varying the external operating conditions of the compressor. Thus, fourteen datasets have been generated to train the simplified model and are shown in Table 2. The results of the detailed model, in terms of suction mass flow rate, mechanical compression power and fluid discharge temperature, have been used to identify the thermodynamic and geometrical parameters of the simplified model through the minimisation of the error function in Equation (33), as described in Section 2.3. Time-averaged values have been calculated with the detailed model in stationary operating condition and have been considered in the analyses. Note that the geometrical parameter BVR has been removed from the identification procedure because it is shared by the two models.
The results of the identification procedure are listed in Table 3.  Finally, fourteen datasets with completely different boundary conditions have been used to test the model and to validate the cross-identification procedure. Fig. 5 and Fig. 6 show the comparison between the results of the two models. The light grey dash dot lines represent the relative errors of  1 % and  4 % for mass flow rate (Fig. 5) and power (Fig. 6), respectively. It is possible to note that all the values are within these boundaries. Therefore, the procedure can be considered validated.

Conclusions
Two modelling approaches for simulating twin-screw compressor performance have been implemented. In order to overcome the limits of both detailed and simplified models, a crossidentification procedure has been fine-tuned. The dynamic model has been developed in Simulink ® and calculates the suction mass flow rate, the mechanical compression power and the fluid discharge temperature in great detail for varying operating conditions. The results of the detailed model have been used to generate compressor performance maps which are necessary for the identification of the simplified model developed in Matlab ® . In order to validate the procedure, a comparison between the results of both models has been made: the suction mass flow rate and the mechanical compression power are predicted with an absolute error less than 1 % and 4 %, respectively. The simplified model is, therefore, reliable and suitable to be used in the simulation and optimisation of entire refrigeration systems.