CFD Based Design for Ejector Cooling System Using HFOS (1234ze(E) and 1234yf)

: The ﬁeld of computational ﬂuid dynamics has been rekindled by recent researchers to unleash this powerful tool to predict the ejector design, as well as to analyse and improve its performance. In this paper, CFD simulation was conducted to model a 2-D axisymmetric supersonic ejector using NIST real gas model integrated in ANSYS Fluent to probe the physical insight and consistent with accurate solutions. HFOs (1234ze(E) and 1234yf) were used as working ﬂuids for their promising alternatives, low global warming potential (GWP), and adhering to EU Council regulations. The impact of di ﬀ erent operating conditions, performance maps, and the Pareto frontier performance approach were investigated. The expansion ratio of both refrigerants has been accomplished in linear relationship using their critical compression ratio within ± 0.30% accuracy. The results show that R1234yf achieved reasonably better overall performance than R1234ze(E). Generally, by increasing the primary ﬂow inlet saturation temperature and pressure, the entrainment ratio will be lower, and this allows for a higher critical operating back pressure. Moreover, it was found out that increasing the degree of superheat for inlet primary ﬂow by 25 K improved the entrainment ratio by almost 20.70% for R1234yf. Conversely, increasing the degree of superheat to the inlet secondary ﬂow has a relativity negative impact on the performance. The maximum overall ejector e ﬃ ciency reached was 0.372 and 0.364 for R1234yf and R1234ze(E) respectively. Comparing the results using ideal gas model, the ejector entrainment ratio was overestimated up to 50.26% for R1234yf and 25.66% for R1234ze(E) higher than using real gas model.


Introduction
Refrigeration is identified as an indispensable method that is widely used in many applications, including food storage, providing thermal comfort, and in the health care industry to preserve pharmaceuticals and medicines. The conventional vapor-compression refrigeration systems are mainly driven by electricity and are usually characterized by high energy consumption [1]. It was observed by the International Energy Agency in 2018 that the demand for air-conditioners worldwide is predicted to soar, resulting in an increase in the number of air-conditioners from 1.6 billion units today to 5.6 billion units by mid-century [2]. Therefore, there is a growing concern that the amount of electricity needed to power them will overload the electrical grids. Even though the need to reduce the cost associated with high energy consumption of conventional refrigerators has resulted in many developments in refrigeration systems such as more energy efficient compressors, there is a limit to which this can be achieved. In addition, there is an increase in global-warming emissions caused by refrigeration such as ejector, liquid pump, and heat generator. The working fluid for this system is divided into two; primary and secondary flow. The primary fluid enters a converging-diverging ejector nozzle where it expands, generates a very low pressure and subsequently accommodates a supersonic flow at the exit [13]. Due to a local pressure drop and tangential force that develop at the edge of the primary flow, the secondary fluid is sucked into the ejector mixing chamber and entrained with higher velocity. In the mixing chamber, part of the energy from primary flow in the form of kinetic energy is transmitted to secondary stream and another part converted into pressure energy. The rest of the energy from both streams is dissipated as heat due to friction and mixing. Afterwards, a normal shock wave is experienced followed by a moderate back pressure in the condenser, causing sudden pressure to increase and deceleration to subsonic flow. The compressed stream then enters the diffuser with intermediate pressure between the motive fluid featured by the primary pressure and the secondary pressure of the entrained flow with deceleration to convert its kinetic energy back into potential energy. The resulting mixed stream leaves to the condenser to reject the heat and condense to liquid. Part of the working fluid then pumped to the generator as the primary stream and the rest goes to the evaporator through the expansion valve where the heat load will be absorbed from the cooling space or products.
In the refrigeration cycles, the performance of the ejector can be analyzed by several parameters. The energy efficiency of the cycle can be quantified by assessing the COP. It is defined as the ratio between the cooling capacities at the evaporator to the gross energy input to the cycle. The ratio between the secondary flow mass flow rate and the primary flow mass flow rate has another important evaluation parameter called entrainment ratio (ω). This is because the entrainment ratio indicates the operational mode of the ejector.
For the duration of the critical mode (on design) of operation, both the primary and the secondary flow are choked upon reaching the maximum value, and it remain constant with a low range of back pressure. However, only the primary flow is choked at the sub-critical mode (off design), resulting in a change in the mass flow rate ratio with respect to back pressure. In the case of back flow mode, a reverse flow at the secondary flow is observed and the entrainment ratio declines to lower than zero. In addition, the compression ratio parameter (CR) (also called pressure lift) is a key indicator for the pumping power because it points out how effectively the pressure of the suction secondary stream can be increased and thus measure the cycle operative range. CR is defined as the ratio of the static pressure between diffuser exit flow and secondary flow. Another important dimensionless parameter is the expansion ratio ER defined as the ratio between the primary pressure to secondary pressure. Regarding the ejector efficiency itself (η ejector ), it can be defined by ASHRAE as the ratio between the actual recovered compression energy and the available theoretical energy in the primary stream [14,15].

Numerical Model Descriptions
In the analytical approach such as the 1-D model, many assumptions needed to be considered such as the friction coefficients related to the mixing losses and isentropic efficiency, hypothetical throat inside the constant area, and the location of the shock waves. CFD does not require these assumptions or any correction factors from empirical results to accurate the model.

Ejector Design
The ejector remains the critical part that influences the ERS overall performance. Therefore, the geometry has to be sensibly determined. The simulated ejector was designed based on the study proposal of 30 KW ejector refrigeration system test stand at technical university of Liberec considering the guidelines obtained from Huang et al. 1999 [16] and preceding literatures to meet the system operational conditions, heat exchangers calculations and selecting of the working fluid. In addition, the geometry is designed to encounter the optimum working condition for both HFOs (1234ze(E) and 1234yf) refrigerants. These refrigerants are promising alternatives to fluorocarbons with low GWP and generally low flammability. Most importantly, they have a positively sloped saturation vapor curve where its entropy decreases along the upper limit curve known as the dry-expansion. In other words, Energies 2020, 13, 1408 4 of 19 the expansion process of the fluid at saturated vapor will end in the superheated vapor flow region which is useful to avoid condensation within the ejector [17,18].
The general schematic drawing of the ejector simulated in this paper is shown in Figure 1. In order to reduce the complexity and time-consumption due to real gas non-linearity, an axisymmetric two-dimensional ejector model was used. In contrast, the results of axisymmetric and three-dimensional simulation are similar [19].
The ejector remains the critical part that influences the ERS overall performance. Therefore, the geometry has to be sensibly determined. The simulated ejector was designed based on the study proposal of 30KW ejector refrigeration system test stand at technical university of Liberec considering the guidelines obtained from Huang et al. 1999 [16] and preceding literatures to meet the system operational conditions, heat exchangers calculations and selecting of the working fluid. In addition, the geometry is designed to encounter the optimum working condition for both HFOs (1234ze(E) and 1234yf) refrigerants. These refrigerants are promising alternatives to fluorocarbons with low GWP and generally low flammability. Most importantly, they have a positively sloped saturation vapor curve where its entropy decreases along the upper limit curve known as the dry-expansion. In other words, the expansion process of the fluid at saturated vapor will end in the superheated vapor flow region which is useful to avoid condensation within the ejector [17,18].
The general schematic drawing of the ejector simulated in this paper is shown in Figure 1. In order to reduce the complexity and time-consumption due to real gas non-linearity, an axisymmetric two-dimensional ejector model was used. In contrast, the results of axisymmetric and three-dimensional simulation are similar [19].

CFD Implementation
In this paper, the simulation was performed with real gas modelling approach using National Institute of Standards and Technology (NIST) fluid database (REFPROP v9.1) integrated in commercial Ansys Fluent R18.1 package [20]. In this respect, the governing equations are discretized using a control volume technique. The calculated axisymmetric space domain was divided as the final structured mesh into 274,241.00 and 270,182.00 number of elements for R1234yf and R1234ze(E) respectively. The reason for that difference is to account for y+ which in all cases was less than 5 as well as to capture the relevant phenomena and probe the physical insight that occur close to the wall.
The domain was quadrilateral meshing structural as shown in Figure 2, provides higher orthogonality and good quality mesh to obtain fast and stable convergence using the real gas model. To ensure gaining accurate numerical results, mesh sensitivity analysis was performed with four-and 15-times finer mesh on both working fluids as shown in Table 1. The analysis was conducted for the entrainment ratio (ω), primary and secondary mass flow rates (ṁ p andṁ s ), and the exit temperatures T c . This indicated how the number of the mesh can affect the numerical results depending on further mesh density. It is shown from the table that the results are stable and there is no significant change with any further increase in the number of elements. Therefore, to perform the computation with less cost and time-consuming using NIST real gas non-linear model, case number 1 was selected to perform the numerical model analysis.

CFD Implementation
In this paper, the simulation was performed with real gas modelling approach using National Institute of Standards and Technology (NIST) fluid database (REFPROP v9.1) integrated in commercial Ansys Fluent R18.1 package [20]. In this respect, the governing equations are discretized using a control volume technique. The calculated axisymmetric space domain was divided as the final structured mesh into 274,241.00 and 270,182.00 number of elements for R1234yf and R1234ze(E) respectively. The reason for that difference is to account for y+ which in all cases was less than 5 as well as to capture the relevant phenomena and probe the physical insight that occur close to the wall.
The domain was quadrilateral meshing structural as shown in Figure 2, provides higher orthogonality and good quality mesh to obtain fast and stable convergence using the real gas model. To ensure gaining accurate numerical results, mesh sensitivity analysis was performed with fourand 15-times finer mesh on both working fluids as shown in Table 1. The analysis was conducted for the entrainment ratio (ω), primary and secondary mass flow rates (ṁp and ṁs), and the exit temperatures Tc. This indicated how the number of the mesh can affect the numerical results depending on further mesh density. It is shown from the table that the results are stable and there is no significant change with any further increase in the number of elements. Therefore, to perform the computation with less cost and time-consuming using NIST real gas non-linear model, case number 1 was selected to perform the numerical model analysis.
The k-ω SST turbulence model has been approved to be the most suitable for ejector calculation because of higher accuracy than another renormalization group RNG-based k-ε. In addition, it also shows better performances in term of stream mixing inside the ejector. Consequently, it has been selected as the turbulence model [1,3,[21][22][23][24]. The simulation was performed under pressure-based coupled solver and the flow was set based on steady-state. For complex flows, it may be ideal to use the second-order upwind discretization since it leads to obtain good results in case of tetrahedral and quad/hex meshes [25]. Real gas models are usually characterized by complex equations, and therefore the solution seems to converge at a much longer time than ideal gas. Hence, to ensure convergence at a shorter time, the flow courant  The k-ω SST turbulence model has been approved to be the most suitable for ejector calculation because of higher accuracy than another renormalization group RNG-based k-ε. In addition, it also shows better performances in term of stream mixing inside the ejector. Consequently, it has been selected as the turbulence model [1,3,[21][22][23][24].
The simulation was performed under pressure-based coupled solver and the flow was set based on steady-state. For complex flows, it may be ideal to use the second-order upwind discretization since it leads to obtain good results in case of tetrahedral and quad/hex meshes [25]. Real gas models are usually characterized by complex equations, and therefore the solution seems to converge at a much longer time than ideal gas. Hence, to ensure convergence at a shorter time, the flow courant number was started from 1 and the under-relaxation factors values set at 0.1 with increasing values of superheat degrees for both inlet flow streams to higher than 30 K, and then after an adequate number of iterations, these parameters can be adjusted gradually to accelerate the convergence depending on the solution stability. It is easier to converge in the single-choke condition hence it is recommended to start within non-critical condition and then step by step increase the back-flow pressure to reach the critical condition. Most of the simulation converged to the residues up to 10 −8 which is very difficult to reach if the density-based solver where used.

Boundary Conditions
The boundary conditions specify the dependent variables on the boundaries of the model to control the solution of the PDE. Basically, the simulated ejector is a passive device where the flow is confined and derived by its unique geometry. The boundary conditions were associated with two inlets linked with the motive flow from the generator and suction flow from the evaporator and one outlet-mixed flow which passed to the condenser. The walls are considered adiabatic. In this simulation, both refrigerant inlet conditions are set at saturated vapor. The primary flow stream temperatures T p were tested at 87, 80, and 70 • C while the secondary flow stream temperatures T s were at 6, 10, and 14 • C. Both flow streams were taken as the corresponding saturation pressure at their temperatures. Two pressure inlets were selected for the motive as the primary flow and the suction as the secondary flow entering the ejector and one outlet pressure for the mixed stream discharged to the condenser.

Numerical Validation
At the beginning of this work, the best available validation found was the experimental result from Huang et al. using 30 different cases at different ejector area ratios and similar operation conditions using results achieved at critical mode [16]. The reason for using Ref. 16 for validation is that the experimental data under various configurations and a wide range of operating conditions were clearly presented. The values of the entrainment ratio and critical back pressure were even specified by numbers. As well as the signature geometries were clearly mentioned and validated in the literature. The best turbulent model which can predict more precise results becomes the main challenge. Table 2 summarizes the test conditions along with the experimental and numerical entrainment ratio (ω), CFD primary and secondary mass flow rates (ṁ p andṁ s ), exit temperatures T c , the error between the CFD model and the experimental data, as well as the number of mesh cells examined. The numerical work was performed under k-ω SST turbulence model using a pressure-based solver, and the inlet flow streams are set to be at saturated vapor. The numerical simulation was observed to converge when the residues of all governing equations fall below 1 × 10 −8 . The present calculated numerical model was able to predict the experimental entrainment ration within ±14% error as seen in Figure 3, while Mohamed el at. show a similar comparison with his CFD work using density-based solver and RNG k-ε turbulence model but with result error within ±20%.     [16] and comparing result from Ref. [12].

Comparison of Real and Ideal Gas Flow Results
The entrainment ratio results of the ejector in the case of working by ideal gas have been obtained in order to evaluate the difference between ideal and the used NIST real gas model. The approach of defining the physical properties of the ideal gas model is an important consideration. For that reason, the physical properties of the ideal gas model were specified based on the primary flow [26,27]. Therefore, viscosity, specific heat, and thermal conductivity were treated as constants throughout the calculation domain and depending on the operating conditions in each case, except the density was modelled using the ideal gas law.
The results show that the performance of the ejector was overestimated using ideal gas model as illustrated in Figure 4. For example, the entrainment ratio for R1234yf in ideal gas model at T p = 70 • C and T s = 6 • C was 16.78% higher than using real gas model. This ratio increased to 50.26% when T p = 87 • C. Comparing to R1234ze(E) at the same previous working conditions, ω of the ideal gas model is about 11.40% higher than that of real gas model at T p = 70 • C, but as the inlet primary flow saturation temperature increased to 87 • C, the difference becomes 25.66%. Researchers assumed in many reviews that simulations using the ideal gas model are able to predict sufficient results but, ironically, the comparison shows how it diverges from the actual real gas model results.

Comparison of Real and Ideal Gas Flow Results
The entrainment ratio results of the ejector in the case of working by ideal gas have been obtained in order to evaluate the difference between ideal and the used NIST real gas model. The approach of defining the physical properties of the ideal gas model is an important consideration. For that reason, the physical properties of the ideal gas model were specified based on the primary flow [26,27]. Therefore, viscosity, specific heat, and thermal conductivity were treated as constants throughout the calculation domain and depending on the operating conditions in each case, except the density was modelled using the ideal gas law.
The results show that the performance of the ejector was overestimated using ideal gas model as illustrated in Figure 4. For example, the entrainment ratio for R1234yf in ideal gas model at Tp = 70 °C and Ts = 6 °C was 16.78% higher than using real gas model. This ratio increased to 50.26% when Tp = 87 °C. Comparing to R1234ze(E) at the same previous working conditions, ω of the ideal gas model is about 11.40% higher than that of real gas model at Tp = 70 °C, but as the inlet primary flow saturation temperature increased to 87 °C, the difference becomes 25.66%. Researchers assumed in many reviews that simulations using the ideal gas model are able to predict sufficient results but, ironically, the comparison shows how it diverges from the actual real gas model results.

Effect of Operation Conditions
The primary and the secondary flow inlet to the ejector were examined under their saturation conditions. Three different primary temperatures (70, 80, and 87 °C) and entrained temperatures (6, 10, and 14 °C) were selected based on the working range of the designed ejector for both refrigerants. The details of the ejector performance were considered by investigating the characteristic curves for the ejector critical and subcritical modes because the ejector performance is assessed in terms of the energy recovered by the secondary flow with respect to the energy available in the primary flow. Figure 5 represent the effect of the primary flow as the desired amount of energy to accelerate and suck the secondary flow by converting the pressure energy of the primary flow into kinetic energy while mixing. Therefore, the entrainment ratio was predicted over the working range of back pressure. The dashed lines represent that the predicted data fit the performance curves while the markers show the numerical results.

Effect of Operation Conditions
The primary and the secondary flow inlet to the ejector were examined under their saturation conditions. Three different primary temperatures (70, 80, and 87 • C) and entrained temperatures (6, 10, and 14 • C) were selected based on the working range of the designed ejector for both refrigerants. The details of the ejector performance were considered by investigating the characteristic curves for the ejector critical and subcritical modes because the ejector performance is assessed in terms of the energy recovered by the secondary flow with respect to the energy available in the primary flow. Figure 5 represent the effect of the primary flow as the desired amount of energy to accelerate and suck the secondary flow by converting the pressure energy of the primary flow into kinetic energy while mixing. Therefore, the entrainment ratio was predicted over the working range of back pressure. The dashed lines represent that the predicted data fit the performance curves while the markers show the numerical results.
It can be seen from Figure 5 that increasing the primary flow temperature caused a lower entrainment ratio and produced a higher critical back pressure. This is associated to the increase of the primary mass flow rate and serves to enlarge the expansion angle, causing a reduction of the annular effective area and increasing the resulting momentum of the mixed stream due to the higher velocity of the entrained stream attained. Therefore, the shock waves move downstream, and the ejector will operate at higher critical pressure. In addition, when the primary or secondary saturated flow temperatures is fixed, then the compression ratio will increase with increasing expansion ratio (ER), causing a lower entrainment ratio, and vice versa. = 14 °C, 10 °C, and 6 °C respectively. When Tp = 87 °C, the difference in ω between the working fluids decreases to 4.91%. These differences are not varying that much in the case of increasing Tp at the fixing Ts, as represented in Figure 5.
The operating back pressure was affected by the inlet operating temperatures. With increasing both inlet saturated flow temperatures or even by increasing one of them, the critical mode for the ejector working with two refrigerants will increase. The on design working mode range is relatively higher in case of using R1234yf. Comparatively, it was found that at fixed Ts, the value of the critical back pressure Pc-cri increased by 100 kPa when Tp increased from 70 to 80 °C or from 80 to 87 °C, while Pc-cri for R1234ze(E) increased by 75 kPa at similar conditions.
In general, R1234yf has in the region of 23.50% higher critical back pressure than R1234ze(E) and higher ejector working back pressure range. (e) (f) The maximum ejector efficiency at minimum energy cost should be considered when ejector cooling system is designed. This is linked with the energy required at the generator to increase the primary flow temperature to a certain limit where the ejection ratio reaches the maximum, nevertheless any extra energy added to the generator will have no significant effect and will be wasted. This analysis needs to be study at specific secondary and back pressure from the working R1234yf represent higher entrainment ratio than R1234ze(E). At T p = 70 • C and T s = 14 • C the ejector works at the highest entrainment ratio. R1234yf recorded ω = 0.9158 where R1234ze(E) has ω = 0.8435. Decreasing the secondary temperatures to 10 • C and 6 • C also decrease ω by (14.72%, 17.71%) for R1234yf and (16.94%, 20.28%) for R1234ze(E) which is a decrease of about 13% back pressure. In other word, decreasing T s at a fixed T p will decrease the entrainment ratio because the secondary mass flow rate,ṁ s , is decreasing whereas the primary mass flow rate,ṁ p , is fixed. Moreover, the difference in ω between the refrigerants will be decreasing in case of increasing the fixed T p . For example, at T p = 70 • C the difference decreases by 8.57%, 11.47% and15.60% in case of T s = 14 • C, 10 • C, and 6 • C respectively. When T p = 87 • C, the difference in ω between the working fluids decreases to 4.91%. These differences are not varying that much in the case of increasing T p at the fixing T s , as represented in Figure 5.
The operating back pressure was affected by the inlet operating temperatures. With increasing both inlet saturated flow temperatures or even by increasing one of them, the critical mode for the ejector working with two refrigerants will increase. The on design working mode range is relatively higher in case of using R1234yf. Comparatively, it was found that at fixed T s , the value of the critical back pressure P c-cri increased by 100 kPa when T p increased from 70 to 80 • C or from 80 to 87 • C, while P c-cri for R1234ze(E) increased by 75 kPa at similar conditions.
In general, R1234yf has in the region of 23.50% higher critical back pressure than R1234ze(E) and higher ejector working back pressure range.
The maximum ejector efficiency at minimum energy cost should be considered when ejector cooling system is designed. This is linked with the energy required at the generator to increase the primary flow temperature to a certain limit where the ejection ratio reaches the maximum, nevertheless any extra energy added to the generator will have no significant effect and will be wasted. This analysis needs to be study at specific secondary and back pressure from the working fluid conditions to avoid testing the optimum case with all parameters. Figure 6 illustrates the contours of Mach number at different primary temperatures T p inside ejector baseline for R1234ze(E) at T s = 10 • C, P c = 560 kPa and R1234yf at T s = 10 • C, P c = 730 kPa. Both refrigerants were represented at same T p = 70, 80, 87 • C from top to bottom. It can be seen that increasing T p , i.e., P p at T p saturation, will increase the oblique shock waves intensity and shift them upstream to the diffuser part. Furthermore, operating at higher P p results in double-choking condition (critical mode) because the amount of energy that drives the entrained flow is enough to accelerate it to the choking state (sonic speed) achieving the peak value of the secondary mass flow rate and lift the flow to that particular back pressure. In addition, the size of the jet core is increased with higher P p and, as a result, this has inversely proportional relation with the effective area (hypothetical throat) [28]. Consequently, it can be seen that the strong normal shock waves occur close to the diffuser at high T p which is an indicator that the ejector is working in double-choking mode of operation but away from the critical point. On the other hand, when T p is fixed and T s is increased, then higher entrained flow energy and less amount of primary flow is required to accelerate the flow. At that point theṁ s will be increasing and hence the entrainment ratio. For this reason, it is recommended to work at higher secondary flow (saturation) temperatures. The primary nozzle operates in the under-expanded regime in all cases and the jet converged toward the ejector center when increasing T p . It can be observed that the flow patterns are almost similar with slightly less intense shock waves for R1234yf. Figure 7 illustrates the contours of Mach number for both working fluids at different back pressure. The top contours of the Mach number for both refrigerants represent the results at critical region where both primary and secondary flow are choked, and the entrainment ratio remains constant. The primary flow leaves the primary nozzle at supersonic velocity and the shock train can be observed along the centerline through the fluctuation of the static pressure before the mixing. A separation layer takes place between the primary and the secondary flow because of the velocity difference in between so that the entrained flow velocity will increase until the mixing take place at the constant throat area. The velocity become subsonic and the oblique shock waves take place due to the static pressure increase. the other hand, when Tp is fixed and Ts is increased, then higher entrained flow energy and less amount of primary flow is required to accelerate the flow. At that point the ṁs will be increasing and hence the entrainment ratio. For this reason, it is recommended to work at higher secondary flow (saturation) temperatures. The primary nozzle operates in the under-expanded regime in all cases and the jet converged toward the ejector center when increasing Tp. It can be observed that the flow patterns are almost similar with slightly less intense shock waves for R1234yf.  Figure 7 illustrates the contours of Mach number for both working fluids at different back pressure. The top contours of the Mach number for both refrigerants represent the results at critical region where both primary and secondary flow are choked, and the entrainment ratio remains constant. The primary flow leaves the primary nozzle at supersonic velocity and the shock train can be observed along the centerline through the fluctuation of the static pressure before the mixing. A separation layer takes place between the primary and the secondary flow because of the velocity difference in between so that the entrained flow velocity will increase until the mixing take place at the constant throat area. The velocity become subsonic and the oblique shock waves take place due to the static pressure increase.
When the back pressure is increased to the critical value, the region of the subcritical mode starts, and that point called the critical back pressure where the Mach number contours represent at the middle contours in the figure. After that point, further increase in back pressure will cause a sharp decrease in the entrainment ratio till it reaches zero and a single choking appeared only for the primary flow as represented in the bottom contours. Afterwards, an additional increase in the back pressure will lead to the presence of reverse flow.
Moreover, it can be observed that increasing the back pressure will move the position of the shock waves towards the back of the primary nozzle direction and eliminate the existence of the shock waves when reaching the critical back pressure, as noticed in operating at back pressure of 613kPa in case of R1234ze(E) and 801 kPa for R1234yf under the inlet condition of Tp = 80 °C, Ts = 10 °C. Regardless of the higher range of critical region for R1234yf, the contours of the Mach number and the shock waves look similar. Additional increase in the back pressure after the critical point will decrease the axial velocity of the mixed flow and move the oblique shock waves to the primary nozzle exit, where it disturbs the primary jet core to the limit when the primary flow could not expand, forcing it to flow reversibly to the entrained flow entrance.  Figure 7 illustrates the contours of Mach number for both working fluids at different back pressure. The top contours of the Mach number for both refrigerants represent the results at critical region where both primary and secondary flow are choked, and the entrainment ratio remains constant. The primary flow leaves the primary nozzle at supersonic velocity and the shock train can be observed along the centerline through the fluctuation of the static pressure before the mixing. A separation layer takes place between the primary and the secondary flow because of the velocity difference in between so that the entrained flow velocity will increase until the mixing take place at the constant throat area. The velocity become subsonic and the oblique shock waves take place due to the static pressure increase.
When the back pressure is increased to the critical value, the region of the subcritical mode starts, and that point called the critical back pressure where the Mach number contours represent at the middle contours in the figure. After that point, further increase in back pressure will cause a sharp decrease in the entrainment ratio till it reaches zero and a single choking appeared only for the primary flow as represented in the bottom contours. Afterwards, an additional increase in the back pressure will lead to the presence of reverse flow.
Moreover, it can be observed that increasing the back pressure will move the position of the shock waves towards the back of the primary nozzle direction and eliminate the existence of the shock waves when reaching the critical back pressure, as noticed in operating at back pressure of 613kPa in case of R1234ze(E) and 801 kPa for R1234yf under the inlet condition of Tp = 80 °C, Ts = 10 °C. Regardless of the higher range of critical region for R1234yf, the contours of the Mach number and the shock waves look similar. Additional increase in the back pressure after the critical point will decrease the axial velocity of the mixed flow and move the oblique shock waves to the primary nozzle exit, where it disturbs the primary jet core to the limit when the primary flow could not expand, forcing it to flow reversibly to the entrained flow entrance.

Effect of Superheat.
In general, the main sources of superheat affecting the ejector operation can be divided into two parts. The first part is due to the internal superheat created by flow irreversibility as a result of the existence of normal shock waves and friction during stream mixing. Another source comes from superheating the inlet primary and/or secondary flow. The global superheat effect then will be a combination of these effects depending on the flow characteristics at inlet constant section determined by the Mach number. The input stream with superheat affects the performance of the system differently. Slightly superheating the refrigerant at the inlet may ensure that the formation of When the back pressure is increased to the critical value, the region of the subcritical mode starts, and that point called the critical back pressure where the Mach number contours represent at the middle contours in the figure. After that point, further increase in back pressure will cause a sharp Energies 2020, 13, 1408 11 of 19 decrease in the entrainment ratio till it reaches zero and a single choking appeared only for the primary flow as represented in the bottom contours. Afterwards, an additional increase in the back pressure will lead to the presence of reverse flow.
Moreover, it can be observed that increasing the back pressure will move the position of the shock waves towards the back of the primary nozzle direction and eliminate the existence of the shock waves when reaching the critical back pressure, as noticed in operating at back pressure of 613 kPa in case of R1234ze(E) and 801 kPa for R1234yf under the inlet condition of T p = 80 • C, T s = 10 • C. Regardless of the higher range of critical region for R1234yf, the contours of the Mach number and the shock waves look similar. Additional increase in the back pressure after the critical point will decrease the axial velocity of the mixed flow and move the oblique shock waves to the primary nozzle exit, where it disturbs the primary jet core to the limit when the primary flow could not expand, forcing it to flow reversibly to the entrained flow entrance.

Effect of Superheat
In general, the main sources of superheat affecting the ejector operation can be divided into two parts. The first part is due to the internal superheat created by flow irreversibility as a result of the existence of normal shock waves and friction during stream mixing. Another source comes from superheating the inlet primary and/or secondary flow. The global superheat effect then will be a combination of these effects depending on the flow characteristics at inlet constant section determined by the Mach number. The input stream with superheat affects the performance of the system differently. Slightly superheating the refrigerant at the inlet may ensure that the formation of droplets inside the ejector is prevented. Formation of droplets is caused by internal condensation after the working fluids in the primary and secondary nozzles are expanded. This should be prevented because it seriously affects the gas dynamic process and the performance of the ejector [16,29]. In addition, blockage of the effective area arises, and this may lead to bumping of the ejector wall which will cause periodically oscillating, unsteady flow, and thus the ejector operation becomes unstable [30,31].
However, the degree of superheat excess will have a negative impact on the system overall performance because of the energy being wasted. The main observations from the results obtained in case of primary flow superheated are shown in Figure 8. The results reveal that, at fixed inlet saturated pressure at T p = 80 • C and T s = 6 • C, increasing the degree of superheat in the inlet primary flow will significantly increase the entrainment ratio and thus the ejector performance. For R1234yf, ω increased by 6.31% when the primary flow temperature is superheated with 5 K. Moreover, ω extended 20.70% higher when the degree of superheating reached 25 K, and after which ω is barely improved by the excessive superheat. For example, for the last 5 K added to reach 25 K of superheat, ω increased only by 2.43% compared with 6.31% increasing of ω at the first 5 K of superheat added. R1234ze(E) has the same influence of increasing ω with increasing the degree of superheat in the primary flow but around half the values compared to R1234yf. For example, ω is increased by just 11.48% at 25 K degree of superheat. The condition of increasing superheat coincides with slight travel of critical pressure point towards higher back pressure range.
The previous result of increasing the superheat degrees in the primary flow contrasts with the secondary flow superheating which has negative impact on the performance while maintaining fixed conditions at the primary inlet flow, as illustrated in Figure 9. It can be seen that 5 K of secondary flow superheat decreased ω by 1.30% and 1.16% for R1234yf and R1234ze(E), respectively. Further, ω continued declining by around 1.17% for each 5 K degree of superheat added. This effect is very small compared to the primary flow superheat and explaining why it is negligible in most of the research papers, for example, that of Chen et al. [32].
conditions at the primary inlet flow, as illustrated in Figure 9. It can be seen that 5 K of secondary flow superheat decreased ω by 1.30% and 1.16% for R1234yf and R1234ze(E), respectively. Further, ω continued declining by around 1.17% for each 5K degree of superheat added. This effect is very small compared to the primary flow superheat and explaining why it is negligible in most of the research papers, for example, that of Chen et al. [32].

Ejector Performance Mapping
To limit the number of the numerical simulations required to analyze ejector performance that relies on different inlet and outlet operation conditions, characteristic curves are generally proposed in terms of these variables which can provide a useful and universal aid to the practical design. The analysis yields a critical system design map from which a system design can be simply studied. The map can be plotted for the ejector using data from simulations or experiments including ω and critical back pressure [13,33,34]  ω continued declining by around 1.17% for each 5K degree of superheat added. This effect is very small compared to the primary flow superheat and explaining why it is negligible in most of the research papers, for example, that of Chen et al. [32].

Ejector Performance Mapping
To limit the number of the numerical simulations required to analyze ejector performance that relies on different inlet and outlet operation conditions, characteristic curves are generally proposed in terms of these variables which can provide a useful and universal aid to the practical design. The analysis yields a critical system design map from which a system design can be simply studied. The map can be plotted for the ejector using data from simulations or experiments including ω and critical back pressure [13,33,34]

Ejector Performance Mapping
To limit the number of the numerical simulations required to analyze ejector performance that relies on different inlet and outlet operation conditions, characteristic curves are generally proposed in terms of these variables which can provide a useful and universal aid to the practical design. The analysis yields a critical system design map from which a system design can be simply studied. The map can be plotted for the ejector using data from simulations or experiments including ω and critical back pressure [13,33,34]. A complete performance map expressed by ejector entrainment ratio varying with the critical back pressure for a range of primary and secondary flow temperatures under saturation conditions for both working fluids is reported in Figure 10. varying with the critical back pressure for a range of primary and secondary flow temperatures under saturation conditions for both working fluids is reported in Figure 10.

Pareto Frontier Solution Curve
A Pareto frontier performance curve predicts the critical back pressure at which the critical working mode of the ejector ends [12]. Alternatively, it can demonstrate the boundaries of the baseline's limit, where the variation of the working design parameters is practical and can be consider as an optimum solution. Figure 11 shows the Pareto frontier performance curve for all operation cases of the ejector based on CFD results drawn on the critical points. It can be seen that R1234ze(E) performs much better, unlike as seen in a previous comparison because of the higher entrainment ratio achieved for a given compression ratio. In this nondimensional comparison, R1234ze(E) is working under higher compression ratio since the secondary flow pressure at the saturation condition is much lower relative to R1234yf. Generally, operating at high condenser temperature (high back pressure) requires high compression ratio to pressurize the mixed stream to the increased condenser pressure and thus causing low entrainment ratio. It can be concluded that by using the Pareto frontier performance curve, the critical entrainment ratio ωcri can be predicted for both refrigerants by the help of their pressure lift ratio in this ejector and specify the working mode using polynomial relations (1) and (2) with ±7.11% and ±2.26% accuracy for R1234yf and R1234ze(E), respectively. Nevertheless, at any equivalent expansion ratio under any operation condition the performance characteristics follow the same trend. For example, for R1234yf at Tp = 80 °C, Ts = 6 °C, ER = 6.54 and Tp = 87 °C, Ts =10°C, ER = 6.63, the difference in critical value of compression ratio is ±1%. Moreover, R1234ze(E) has ±0.77% difference at Tp = 80 °C, Ts = 6 °C, ER = 7.47, and Tp = 87 °C, Ts = 10 °C, ER = 7.54.
The expansion ratio (ER) of both refrigerants has been accomplished in linear relationship using their critical compression ratio within ±0.30% accuracy. These correlations are concluded in equations (3) and (4). By using these relations, the ejector expansion ratio can be easily achieved with the help of the critical compression ratio and simplified the techniques to calculate the entrainment ratio at any operation condition. This approach could save computational time and cost for predicting the ejector performance. R1234ze(E) R1234yf T p =80°C Figure 10. Performance map of the ejector based on CFD data.

Pareto Frontier Solution Curve
A Pareto frontier performance curve predicts the critical back pressure at which the critical working mode of the ejector ends [12]. Alternatively, it can demonstrate the boundaries of the baseline's limit, where the variation of the working design parameters is practical and can be consider as an optimum solution. Figure 11 shows the Pareto frontier performance curve for all operation cases of the ejector based on CFD results drawn on the critical points. It can be seen that R1234ze(E) performs much better, unlike as seen in a previous comparison because of the higher entrainment ratio achieved for a given compression ratio. In this nondimensional comparison, R1234ze(E) is working under higher compression ratio since the secondary flow pressure at the saturation condition is much lower relative to R1234yf. Generally, operating at high condenser temperature (high back pressure) requires high compression ratio to pressurize the mixed stream to the increased condenser pressure and thus causing low entrainment ratio. It can be concluded that by using the Pareto frontier performance curve, the critical entrainment ratio ω cri can be predicted for both refrigerants by the help of their pressure lift ratio in this ejector and specify the working mode using polynomial relations (1) and (2) with ±7.11% and ±2.26% accuracy for R1234yf and R1234ze(E), respectively. Nevertheless, at any equivalent expansion ratio under any operation condition the performance characteristics follow the same trend. For example, for R1234yf at T p = 80 • C, T s = 6 • C, ER = 6.54 and T p = 87 • C, T s =10 • C, ER = 6.63, the difference in critical value of compression ratio is ±1%. Moreover, R1234ze(E) has ±0.77% difference at T p = 80 • C, T s = 6 • C, ER = 7.47, and T p = 87 • C, T s = 10 • C, ER = 7.54.
The expansion ratio (ER) of both refrigerants has been accomplished in linear relationship using their critical compression ratio within ±0.30% accuracy. These correlations are concluded in Equations (3) and (4). By using these relations, the ejector expansion ratio can be easily achieved with the help of the critical compression ratio and simplified the techniques to calculate the entrainment ratio at any operation condition. This approach could save computational time and cost for predicting the ejector performance.
ER  Figure 11. Pareto Frontier performance curve for all operation cases of the ejector based on CFD data.

Ejector Efficiency
The ejector efficiency was first introduced by Köhler et al. [35]. The great advantage of this represented efficiency is that it can be calculated by using only the external parameters which are easy to measure. The same expression for the total ejector efficiency was presented by Elbel et al. deriving by different approach as shown in Equation (5) [36]. This formula was used in many other literatures, as in [37][38][39]. It was defined as the ratio of the amount of the recovered ejector expansion work rate with maximum possible expansion work rate recovery potential. It can also describe the effectiveness of the ejector as the ratio of actual work recovered to the maximum available work recovery.
where hs,is and hp,is are referred to specific enthalpy for an assumed isentropic state from the inlet primary and secondary conditions to the ejector outlet pressure. hs and hp are the specific enthalpy for the inlet primary and secondary conditions. Moreover, Dvorak et al. proposed the ejector efficiency based on the relation between the specific compression work acquired by the secondary stream and the specific work exerted by the primary flow as Equation (6) [40]. This formula has been used in many researches [41,42]. The only problem in Dvorak formula is the adiabatic exponent which needs to be constant and represents the ideal gas flow. Therefore, Equation (5) is used in this work to analyze the efficiency for the ejector model.
In Figure 12 the ejector efficiency based on the numerical results were presented for different working conditions. The markers on the figure show the CFD obtained results while as the lines predict the efficiency that fit these results. The same trends were revealed for both refrigerants' efficiencies. In addition, the characteristic curve of the ejector efficiency ends at which the critical flow occurs. A further reduction in the back pressure at critical mode will not lead to an increase in the entrainment ratio through the ejector and will only decline the efficiency until zero. It can be seen that R1234ze(E) has higher efficiency compared to R1234yf under the same entrainment ratio while R1234yf could provide higher ejector efficiency at critical point. The maximum efficiency reached at the critical point were 0.372 and 0.364 for R1234yf and R1234ze(E) respectively. Furthermore, at fixed Ts, the ejector has high efficiency when it works under low Tp. On the other hand, at subcritical mode condition where the back pressure increases, the efficiency decreases which is as a result of the decreasing of the mass entrainment ratio at these operation conditions. Comparing the results in the case of using Dvorak's proposal as an ideal efficiency of the gas flow, the efficiency of the ejector will   Figure 11. Pareto Frontier performance curve for all operation cases of the ejector based on CFD data.

Ejector Efficiency
The ejector efficiency was first introduced by Köhler et al. [35]. The great advantage of this represented efficiency is that it can be calculated by using only the external parameters which are easy to measure. The same expression for the total ejector efficiency was presented by Elbel et al. deriving by different approach as shown in Equation (5) [36]. This formula was used in many other literatures, as in [37][38][39]. It was defined as the ratio of the amount of the recovered ejector expansion work rate with maximum possible expansion work rate recovery potential. It can also describe the effectiveness of the ejector as the ratio of actual work recovered to the maximum available work recovery.
where h s,is and h p,is are referred to specific enthalpy for an assumed isentropic state from the inlet primary and secondary conditions to the ejector outlet pressure. h s and h p are the specific enthalpy for the inlet primary and secondary conditions. Moreover, Dvorak et al. proposed the ejector efficiency based on the relation between the specific compression work acquired by the secondary stream and the specific work exerted by the primary flow as Equation (6) [40]. This formula has been used in many researches [41,42]. The only problem in Dvorak formula is the adiabatic exponent which needs to be constant and represents the ideal gas flow. Therefore, Equation (5) is used in this work to analyze the efficiency for the ejector model.
In Figure 12 the ejector efficiency based on the numerical results were presented for different working conditions. The markers on the figure show the CFD obtained results while as the lines predict the efficiency that fit these results. The same trends were revealed for both refrigerants' efficiencies. In addition, the characteristic curve of the ejector efficiency ends at which the critical flow occurs. A further reduction in the back pressure at critical mode will not lead to an increase in the entrainment ratio through the ejector and will only decline the efficiency until zero. It can be seen that R1234ze(E) has higher efficiency compared to R1234yf under the same entrainment ratio while R1234yf could provide higher ejector efficiency at critical point. The maximum efficiency reached at the critical point were 0.372 and 0.364 for R1234yf and R1234ze(E) respectively. Furthermore, at fixed T s , the ejector has high efficiency when it works under low T p . On the other hand, at subcritical mode condition where the back pressure increases, the efficiency decreases which is as a result of the decreasing of the mass entrainment ratio at these operation conditions. Comparing the results in the case of using Dvorak's proposal as an ideal efficiency of the gas flow, the efficiency of the ejector will have the same change trends, with lower values up to 16% in the case of R1234yf and 5.36% for R1234ze(E).
Energies 2020, 13, x FOR PEER REVIEW 12 of 20 condition where the back pressure increases, the efficiency decreases which is as a result of the decreasing of the mass entrainment ratio at these operation conditions. Comparing the results in the case of using Dvorak's proposal as an ideal efficiency of the gas flow, the efficiency of the ejector will have the same change trends, with lower values up to 16% in the case of R1234yf and 5.36% for R1234ze(E).

Conclusions
This paper is mainly devoted to computational fluid dynamics simulation conducted for a 2D axisymmetric ejector to present the performance over the entire operation conditions using NIST real gas model integrated in commercial ANSYS Fluent. The experimental data were used to validate Energies 2020, 13, 1408 16 of 19 how effective the model is. The comparison results have shown that the model could accurately be able to predict the experimental results with error less than ±14%. The numerical analysis study was applied using k-ω SST turbulence model and pressure-based coupled solver. Most of the simulations converged to the residues fall below 1×10 −8 . HFOs (R1234ze(E) and R1234yf) were selected as the working fluids because of their promising alternatives to fluorocarbons with low global warming potential and flammability as well as adhering to EU Council regulations. Mesh sensitivity was performed to ensure the result stability via different number of elements.
The inlet primary and the secondary flow were examined under their saturation conditions. Given the energy available in the primary flow, the ejector performance is assessed based on how much of this energy the secondary flow can recover. This was done by investigating the characteristic curves for the ejector critical and subcritical modes. It was shown that increasing the primary flow temperature caused lower entrainment ratio and produced higher critical back pressure. When the primary or secondary saturated flow temperatures are fixed then the compression ratio will increase with increasing the expansion ratio and causes lower entrainment ratio indeed. As a result, R1234yf represented higher entrainment ratio than R1234ze(E) by 4.91% to 15.60%. The contours of Mach number for both working fluids at different inlet primary pressure and different back pressure were represented. Several phenomena in the ejector were discussed, like the shock waves, effective mixing, and operational modes. In order to prevent the formation of droplet inside the ejector, adding degree of superheat in the inlet flow is important but the excess of value will have a negative impact on the system overall performance because of the energy being wasted. Therefore, the amount of superheat to be added should be carefully studied. From the performance evaluation point of view, the analysis demonstrated that the entrainment ratio can be increased up to 20.70% with 25 K superheat for R1234yf, which is almost double the percentage compared to R1234ze(E). Secondary flow superheating has a negative impact on the performance, which slightly decreases ω within 1.17% for both refrigerants for each 5 K degree of superheat added.
Ejector performance map was proposed using the critical back pressure to obtain the characteristic curve and show the large critical mode which R1234yf can work in comparing to R1234ze(E). A Pareto frontier performance curve was performed. This approach could save computational time and cost for predicting the ejector performance by calculating the critical entrainment ratio using ejector pressure lift ratio. Moreover, the expansion ratio (ER) of both refrigerants was accomplished in linear correlation using their critical compression ratio. It was also shown that regardless of any operation condition, the refrigerant working at equivalent expansion ratio have similar performance characteristics curve trend. The effectiveness of the modelled ejector was determined to compare both working fluids efficiency. The maximum efficiency reached at the critical point was 0.372 and 0.364 for R1234yf and R1234ze(E) respectively. In other words, the efficiency was lower when using ideal gas flow formula. The ideal gas model has been simulated and compared with the real gas model. It was concluded that the performance of the ejector will be overestimated in the case of using the ideal gas model with maxima of 50.26% and 25.66% for R1234yf and R1234ze(E), respectively. This explains the necessity of using a real gas model in CFD and why the ideal gas model could fail to predict the performance of the ejector. In a nutshell, R1234ze(E) can achieve relatively better overall performance than R1234yf.
This study was conducted to design a high-performance ejector that has been manufactured already and currently under test. The main essential parameters for energy-efficient HFOs ejector cooling systems have been investigated to estimate the suitability of the modeled ejector. However, more results and greater accuracy, for example, the effect of the NXP, pressure distribution along the ejector, and the entropy generation will be considered and validated with experiments in future work.