Numerical Modeling of Ejector and Development of Improved Methods for the Design of Ejector-Assisted Refrigeration System

Abstract

An ejector is a simple mechanical device that can be integrated with power generation or the refrigeration cycle to enhance their performance. Owing to the complex flow behavior in the ejector, the performance prediction of the ejector is done by numerical simulations. However, to evaluate the performance of an ejector integrated power cycle or refrigeration cycle, the need for simpler and more reliable thermodynamic models to estimate the performance of the ejector persists. This research, therefore, aims at developing a single mathematical correlation that can predict the ejector performance with reasonable accuracy. The proposed correlation relates the entrainment ratio and the pressure rise across the ejector to the area ratio and the mass flow rate of the primary flow. R141b is selected as the ejector refrigerant, and the results obtained through the proposed correlation are validated through numerical solutions. The comparison between the analytical and numerical with experimental results provided an error of less than 8.4% and 4.29%, respectively.

1. Introduction

An ejector is a simple mechanical device that uses the low pressure created by the expansion of high-pressure motive fluid (primary) and generates a vacuum, which is then used to entrain and subsequently compress a secondary fluid [1]. The widespread application of an ejector includes: the removal of ash from flue gases [2]; in desalination plants [3]; solid oxide fuel cells [4]; refrigeration machines [5]; fueling of hydrogen vehicles [6]; and power generation cycle [7]. Ejectors could be integrated into various energy conversion processes to improve the energy efficiency, and have been actively researched since the 1950s [8]. However, the ejector-assisted refrigeration systems, or the ejector expansion refrigeration cycle (EERC), have gained the most interest, due to their superior performance, operational flexibility, and ease of control [9]. Recently, the interest in ejector-assisted refrigeration machines has been greatly renewed, since they can be operated using low-grade heat [10,11].

In an EERC, the ejector is either applied to replace the mechanical compressor, as in a steam jet cooling system, or reduce the irreversibility during the throttling process by harnessing the kinetic energy released during the process [9]. In the latter case, the kinetic energy is utilized to raise the suction pressure higher than the evaporation pressure, thereby reducing the compression power. A schematic of a typical ejector expansion refrigeration cycle (EERC) and the corresponding pressure enthalpy (P-h) diagram is shown in Figure 1 [12,13]. A 17% increase in coefficient of performance (COP) has been reported using EERCs [13]. In contrast to the conventional refrigeration machine, in an EERC, the expansion valve is replaced by an ejector. Instead of throttling the high-pressure refrigerant (motive fluid) from the condenser, it is expanded in the convergent divergent (CD) primary nozzle of the ejector. The expansion of the refrigerant in the CD generates a vacuum and entrains the refrigerant from the evaporator (secondary fluid). The two streams subsequently mix in the mixing chamber, undergo a shock wave, and are fed into a diffuser to increase the pressure of the mixed stream. The ejector pre-compresses the working fluid to back pressure in the EERC, and, consequently, reduces the compression load [9].

The ejector is at the core of the EERC, and an integral step in the design and performance estimation of an EERC is the modeling of the ejector. The performance of the ejector is quantified in terms of two main global parameters; (i) entrainment ratio, and (ii) compression pressure ratio. The entrainment ratio is the ratio of secondary fluid flow rate to the primary (motive) fluid flow rate, while the compression pressure ratio reflects the pressure gain of the secondary fluid [14]. Many thermodynamic models were developed to predict the ejector performance. Huang et al. [15] pioneered the development of a 1D model for the prediction of ejector performance. The model was based on ideal gas assumption, isentropic flow relations with constant isentropic efficiency, and the constant value for specific heat, and showed 15% deviation from the experimental results. Over time, many improved models to estimate the ejector performance have been developed and proposed [9]. In the Huang model, the choice of isentropic efficiencies is the major source of error, therefore, Haghparast et al. [16] replaced the isentropic efficiencies with polytropic efficiency to improve the model accuracy. Selvaraju et al. [17] proposed another analysis code that induces the effect of friction in the ejector, which is normally neglected by its predecessors. The accuracy of the 1D analytical model has been improved; nonetheless, they are based on assumptions of ideal gas flow and normal shock, which does not reflect the actual flow behavior in the ejector [9].

The fluid flow through the ejector is compressible and supersonic, involving a series of complex flow interactions, such as shock wave–boundary layer interaction and flow mixing [18]. The 1D analytical models do not incorporate the complex flow interactions occurring in the ejector; however, the complex flow structure plays a significant role in estimating the performance of the ejector [18]. Therefore, various researchers [12,18] employed computational fluid dynamics (CFD) models to better visualize the fluid flow and provide a detailed description of complex fluid flow inside the ejector. These models essentially solve the Reynolds-averaged Navier-Stokes (RANS) equation by using turbulence models and numerical methods [19]. Desevaux et al. [20] used a standard k-ε turbulence model in FLUENT to study the ejector flow. However, their study was based on ideal gas assumption, and underpredicted the back pressure by about 20%. Del Valle et al. [21] presented a CFD model that considered real gas properties, and the results showed better agreement with the experimental results. Han et al. [22] utilized the CFD model to comprehensively investigate the mechanism of boundary layer separation in the ejector and its effect on ejector performance. Their study concluded that either a too small or too big mixing chamber diameter can induce boundary layer separation, which is the main cause of back flow in the ejector. Elbarghthi et al. [23], as well, employed CFD to explore the performance of an EERC with alternative refrigerants, such as R-1234ze(E) and R-1234yf. An insightful review of system and component level numerical modeling of ejector was done by Little et al. [24], in which they highlighted that the pathway to improving the ejector analysis is the advance modeling of turbulence effects and the phase changes.

The 1D, as well as CFD, analysis of ejectors have been continuously evolving. Although the 1D models pose certain limitations, however, they are readily adaptable and compatible for the system level analysis of EERCs. The constant isentropic efficiencies are a major source of discrepancies in the 1D analytical models. Therefore, to improve the choice of isentropic efficiencies for the sizing and performance prediction of the ejector, Haghparast et al. [16] used CFD tools to approximate the isentropic efficiencies. In their study, the appropriate efficiencies are extracted from the CFD model, and subsequently used in the 1D model to estimate the ejector performance. Riaz et al. [12] extended the same integration of CFD tools and 1D models to design and optimize the low-grade waste heat driven ejector refrigeration system. Rogié et al. [25] utilized a similar technique in the investigation of the ejector as a potential replacement of expansion valve in the hydrogen fueling station to minimize the energy expenditure and the fueling time. The CFD analysis provides more credibility in the analysis results; however, they lack the ability to be readily integrated for a system level analysis. The varying operating conditions leading to high computation and time requirements make CFD models unfit for real time implementation in the practical systems [18,19].

The CFD models considerably improve the quality of ejector analysis; however, they are focused on the component level study and, hence, are not compatible for the performance analysis of a complete system. On the other hand, the accuracy of 1D models for ejectors is still questionable. Therefore, in this study, empirical correlations, based on experimental data of Huang at al. [15], have been developed to predict ejector performance parameters. The study describes how the performance of an ejector is mainly governed by two parameters: (i) area ratio of the initial CD nozzle and the mixing chamber, and (ii) the ratio of motive fluid inlet pressure to the secondary inlet pressure. Therefore, a multivariate polynomial regression technique is applied to predict the entrainment ratio and compression pressure ratio of the ejector as a function of major performance influencing parameters. Furthermore, a CFD model is developed to investigate the ejector performance. Zhu et al. [26] explored the ejector flow by using experimental and numerical methods. Their study investigated four turbulence models and found out the k-ε model agrees best with measurements for predictions of mass flow rates, which is adopted in this study. The CFD model is validated against the experimental results. The validated CFD model is subsequently used to evaluate the accuracy of the results calculated using the empirical correlations. The empirical correlations have the advantage of easy integration with the overall system model, and facilitate the real time implementation of ejectors into practical refrigeration systems.

2. Ejector Description and Modeling

In the following section, the 1D model of ejector analysis, with the necessary assumptions, is given. Subsequently, the development and implementation of the numerical (CFD) model is also discussed in detail, with the focus on geometry construction, meshing, selection of turbulence model, and the related FLUENT settings.

2.1. Mathematical Solution of Ejectors

The thermodynamic models are the early models developed to understand and estimate the ejector performance [27]. Basically, these 1D and 2D mathematical models analytically solve the gas dynamic relations for compressible flow to model the flow field inside the ejector [28]. Although these mathematical models provide a certain level of understanding of fluid flow through the ejector, these models are still far behind in fully explaining the ejector flow field, because of involvement of certain flow discrepancies, such as normal shock waves, primary and secondary flow mixing, and boundary layer separation, as well as use of isentropic coefficients [29]. These flow features obstruct the development of detailed models that could fully analyze the ejector flow, as then the complexity of governing equations go beyond the capability of analytical procedures. The existing analytical models simulate the flow by making various assumptions such that the governing system of equations is simplified [15]. All these models are used to evaluate the ejector performance in terms of global flow parameters, such as entrainment ratio (ω) and compression ratio (P_c) [12]. The geometry of an ejector is shown in Figure 2. The important boundary conditions include throat area of primary nozzle (A_t), area of mixing chamber (A₃), and the inlet pressures of motive fluid (P_g) and secondary fluid (P_e).

In this study, the 1D thermodynamic model of Huang et al. [15] is utilized for describing the performance of the ejector at various values of P_g corresponding to saturated vapor temperature from 351 K to 368 K, which lies within the range of ejector refrigeration machinery. The various geometrical specifications of ejector that were used in this study are given in

Table 1. Geometrical specifications of ejector.

Ejector Geometry
Primary Nozzle			Mixing (Constant Area) Section
			Serial No.	Diameter (D3)
Serial No.	Throat Diameter (Dt)	Exit Diameter (D4)	A	6.70 mm
A	2.64 mm	4.50 mm	B	6.98 mm
			C	7.60 mm
			D	8.10 mm
E	2.82 mm	5.10 mm	E	8.54 mm
			G	7.34 mm
			H	9.20 mm

.

The system of equations used to solve the analytic ejector model shown in Figure 2 are summarized in

Table 2. Conventional calculation procedure.

Step	Inputs	Equations	Output	Comments
1	$P_g,T_g,$$A_t$	${\dot{m}}_p=\frac{P_g{A}_t}{\sqrt{T_e}}\sqrt{\frac{\gamma }{R}{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma +1}{\gamma -1}}}\sqrt{{\eta }_p}$	${\dot{m}}_p$	At choking condition, the mass flow rate through the primary nozzle follows a gas dynamic relation.
2	$A_{p4}$	$\frac{A_{p4}}{A_t}=\frac{1}{M_{p4}^2}{\left[\frac{2}{\gamma +1}\left(1+\frac{\gamma -1}{2}{M}_{p4}^2\right)\right]}^{\frac{\gamma +1}{\gamma -1}}$ $\frac{P_g}{P_{p4}}={\left[1+\frac{\gamma -1}{2}{M}_{p4}^2\right]}^{\frac{\gamma }{\gamma -1}}$	$M_{p4},$$P_{p4}$	The Mach no. M4 is calculated by using the Newton–Raphson method.
3	$P_e$	$P_{sH}=\frac{P_e}{{\left(1+\frac{\gamma -1}{2}{M}_{sH}^2\right)}^{\frac{\gamma }{\gamma -1}}}$	$P_{sH}$	Referring to the assumptions made, the Mach number of secondary flow at hypothetical throat is MeH = 1.
4	$P_{sH}=$$P_{pH}$	$\frac{P_{sH}}{P_{p4}}=\frac{{\left[1+\frac{\gamma -1}{2}{M}_{p4}^2\right]}^{\frac{\gamma }{\gamma -1}}}{{\left[1+\frac{\gamma -1}{2}{M}_{pH}^2\right]}^{\frac{\gamma }{\gamma -1}}}$ $\frac{A_{pH}}{A_{p4}}=\frac{\left(\frac{{\phi }_p}{M_{pH}}\right){\left[\frac{2}{\gamma +1}\left(1+\frac{\gamma -1}{2}{M}_{pH}^2\right)\right]}^{\frac{\gamma +1}{2\left(\gamma -1\right)}}}{\left(\frac{1}{M_{p4}}\right){\left[\frac{2}{\gamma +1}\left(1+\frac{\gamma -1}{2}{M}_{p4}^2\right)\right]}^{\frac{\gamma +1}{2\left(\gamma -1\right)}}}$	$M_{pH},$$A_{pH}$	${\phi }_p$: an isentropic coefficient that represents flow losses as primary fluid flow from section 4-4 to section H-H.
5	$A_{pH}$	$A_{sH}=A_3+A_{pH}$	$A_{sH}$	If AsH < 0, calculate A3 by using A3 = ApH + ΔA3, otherwise return to step 4 to recalculate ApH, and again the condition is checked.
6	$A_{sH},P_e,$ $T_e,\gamma ,$ $R, {\eta }_s$	${\dot{m}}_s=\frac{P_e{A}_{sH}}{\sqrt{T_e}}\sqrt{\frac{\gamma }{R}{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma +1}{\gamma -1}}}\sqrt{{\eta }_s}$	${\dot{m}}_s$	${\eta }_s$: Isentropic efficiency of entrained flow.
7	$T_g,T_e,$ $M_{pH},M_{sH}$	$\frac{T_g}{T_{pH}}=1+\frac{\gamma -1}{2}{M}_{pH}^2$ $\frac{T_e}{T_{sH}}=1+\frac{\gamma -1}{2}{M}_{sH}^2$	$T_{pH,}{T}_{sH}$	Value of Tg and Te can be taken from step 1 and 3, respectively.
8	$T_{pH,}{T}_{sH},$ $\gamma ,R,$ $M_{pH},$ $M_{sH},$ ${\dot{m}}_p,{\dot{m}}_s$	${\phi }_m\left[{\dot{m}}_p{V}_{pH}+{\dot{m}}_s{V}_{sH}\right]=\left({\dot{m}}_p+{\dot{m}}_s\right)V_m$ $V_{pH}=M_{pH}{a}_{pH},V_{sH}=M_{sH}{a}_{sH}$ $a_{pH}=\sqrt{\gamma R{T}_{pH}},a_{sH}=\sqrt{\gamma R{T}_{sH}}$	$V_M,$ $V_{pH},$ $V_{sH}$	${\phi }_M$: mixed flow friction coefficient, PM = PpH = PsH, MsH = 1, and MpH can be taken from step 4.
9	$V_M,\gamma ,R,$ $C_p,T_{pH},$ $T_{sH},V_{pH},V_{sH},$ ${\dot{m}}_p,{\dot{m}}_s$	${\dot{m}}_p\left(C_p{T}_{pH}+\frac{V_{pH}^2}{2}\right)+{\dot{m}}_s\left(C_p{T}_{sH}+\frac{V_{sH}^2}{2}\right)=\left({\dot{m}}_p+{\dot{m}}_s\right)\left(C_p{T}_m+\frac{V_m^2}{2}\right)$ $a_M=\sqrt{\gamma R{T}_M},M_M=\frac{V_M}{a_M}$	$T_M,M_M$	The first equation gives TM, which is then used to find value of αM and MM.
10	$P_M,M_M,\gamma$	$P_3=P_M\left(1+\frac{2\gamma }{\gamma +1}\left(M_M^2-1\right)\right)$ $M_3=\sqrt{\frac{1+\frac{\gamma -1}{2}{M}_M^2}{\gamma M_M^2-\frac{\gamma -1}{2}}}$	$P_3,M_3$	Flow is solved after the shock wave, and value of PM can be taken from step 8.
11	$P_3,M_3$	$P_{Cn}=P_3{\left(1+\frac{\gamma -1}{2}{M}_3^2\right)}^{\frac{\gamma }{\gamma -1}}$	$P_{Cn}$	Flow pressure at diffuser exit is calculated.
12	$P_{Cn},P_{Cn}^{*}$	$if P_{Cn}\ge P_{Cn}^{*} then A_3=A_3-\Delta A_3$ $if P_{Cn}<P_{Cn}^{*} then A_3=A_3+\Delta A_3$	$A_3$	$P_{Cn}^{*}$: critical condenser pressure, and A3 must be equal to A3 in step 5, otherwise procedure starts again from step 5.
13	${\dot{m}}_p,{\dot{m}}_s$	$\omega =\frac{{\dot{m}}_s}{{\dot{m}}_p}$	$\omega$	Entrainment ratio is calculated by using ṁs and ṁp from step 6 and step 1, respectively.

, with the flowing assumptions.

• The working fluid acts as an ideal gas, having constant specific heat (C_p) and specific heat ratio (γ).
• Steady, adiabatic, and 1D flow.
• Negligible kinetic energy at secondary flow inlet, primary nozzle inlet, and diffuser exit.
• Use of isentropic relations and constant mixing chamber efficiency.
• The primary and secondary fluid flow mixes at hypothetical throat located within constant area section (Section 2–3).
• Constant pressure mixing (P_pH = P_sH).
• Choking of entrained flow at hypothetical throat and Mach number of M_sH = 1 is assumed.
• Adiabatic ejector walls.

The limitations of abovementioned analytical model include (a) constant values of isentropic coefficients and isentropic efficiencies of various components, and (b) the procedure for analytically solving the governing equations, mentioned above, is an iterative process.

2.2. Mathematical Solution of Ejectors

In the present study, a systematic approach was proposed to develop a simple method to calculate the optimal values of ω and P_c without the need of iterations. A regression analysis was conducted to identify the likely relation between ejector performance parameters and ejector operating conditions and geometrical specifications. To correlate the inputs (A₃, A_t, P_g, P_e) with outputs (ω and C_P), regression equations were developed using the least square method. As ω and P_c are two dependent variables, two regression equations were proposed, which required both area ratio (A₃/A_t) and pressure ratio (P_g/P_e) as an input (independent variables). The relationship between dependent variable ω with A₃/A_t and P_g/P_e is shown Figure 3.

The figure clearly shows a nonlinear relation between the ω with P_g/P_e and A₃/A_t, and there exists a negative nonlinear association between ω and P_g/P_e, that is, decreasing the value of P_g/P_e results in increased entrainment ratio. Similarly, the association between the other dependent variable (P_cn) and independent variables is also shown by Figure 4. Each graphical representation shows that the ejector performance parameters are highly influenced under the variations in inlet flow pressures (P_g, P_e) and ejector geometry (A₃, A_t).

The correlations were developed using the experimental data available in the study conducted by Huang et al. [15], and is provided in Appendix A. This data includes experimentally determined values of ω, compression ratio (P_cn/P_e), A₃/A_t, and P_g/P_e. An empirical correlation of second-degree polynomial (that predicts ejector performance in terms of ω as a function of A₃/A_t and P_g/P_e), was developed using polynomial regression techniques because of the existence of the nonlinear behavior among the predictor and response variables, as shown in Figure 3 and Figure 4. Using a forward selection method, a second-degree polynomial, in terms of area ratio and pressure ratio, was selected. The general equation with second-degree polynomial [30] in the present case becomes:

$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_{11}{x}_1^2+{\beta }_{22}{x}_2^2+{\beta }_{12}{x}_1{x}_2+\mathsf{\epsilon }$$

(1)

where y: dependent variable; x: independent variable; β₀: intercept; β₁ and β₂: linear effect parameters; β₁₁, β₂₂: quadratic effect parameters; β₁₂: interaction parameter; and ε: error. The above equation can be represented in terms of specific independent variables, as given below:

$$\begin{gathered} E\left(y\right)=\omega ={\beta }_0+{\beta }_1\left(A_3/A_t\right)+{\beta }_2\left(P_g/P_e\right)+{\beta }_{11}{(A_3/A_t)}^2+{\beta }_{22}{(P_g/P_e)}^2 \\ +{\beta }_{12}\left(A_3/A_t\right)\left(P_g/P_e\right) \end{gathered}$$

(2)

The above expression is also termed as response surface [31] where β₀, β₁, β₂, β₁₁, β₂, and β₁₂ are coefficient parameters of the equation. These equation parameters were found by applying the least square approach on the matrix form of the above equations, which generated the following matrix relation:

$$\left[\beta \right]={\left({\left[X\right]}^{*}\left[X^t\right]\right)}^{-1}\times \left[X^t\right]\left[Y\right]$$

(3)

where the matrix [X] contained n values of all the independent variables (A₃/A_t and P_g/P_e), and matrix [Y] contained all n values of dependent variable (ω) associated with specific values of independent variables. The above matrix expression was then solved in MATLAB, and the following empirical correlation was obtained:

$$\begin{gathered} \omega =0.1705+0.1479\left(\frac{A_3}{A_t}\right)-0.07002\left(\frac{P_g}{P_e}\right)+0.0014{\left(\frac{A_3}{A_t}\right)}^2+0.0035{\left(\frac{P_g}{P_e}\right)}^2 \\ -0.0074\left(\frac{A_3}{A_t}\times \frac{P_g}{P_e}\right) \end{gathered}$$

(4)

Similarly, the second empirical correlation for predicting the compression ratio was obtained following the same procedure, except that, due to the linear relation between the response and predictor variables, the general regression equation is shown below, and the values of coefficients were calculated as explained above and the developed empirical equation is also shown below [30]. In the later section, both empirical equations are validated against the CFD model developed, to predict ejector performance at varying ejector specifications.

$$\mathrm{y}={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\mathsf{\epsilon }\frac{P_{cn}}{P_e}=2.3204-0.12314(\frac{A_3}{A_t})+0.1713(\frac{P_g}{P_e})$$

(5)

Figure 5 shows a comparison between performance evaluation of ejector done through the conventional 1D thermodynamics model and the empirical correlations. The figure clearly demonstrates how the complex solution procedure of the 1D model can be replaced by a single-step empirical correlation.

CFD is a robust tool that helps in better visualization of ejector flow phenomena, such as supersonic flow, flow mixing, shock trains, and boundary layers, as discussed above. In this study, a CFD model has been developed, utilizing the commercially available simulation package ANSYS FLUENT 2019. The software package uses finite volume techniques, that are based on discretization of governing equations, by dividing the physical geometry into smaller elements, forming a control volume mesh. Originally, the fluid flow through the ejector is essentially unsteady in a 3D space, mainly due to the turbulent flow nature. However, by utilizing the RANS equations, which are used for determining averaged values of flow quantities by time averaging over long intervals [32], the problem is assumed to be steady and axisymmetric, providing acceptable level of accuracy in both global and local flow phenomena [33]. Therefore, in this study, the ejector is modeled utilizing axisymmetric geometry, as literature studies have shown that a 2D axisymmetric model will produce identical results to a 3D flow model, but with less computation effort [34]. By utilizing axis symmetric condition, 2D governing equations defined in terms of radial components are solved using a finite volume discretization method, with second-order upwind scheme [12]. These axisymmetric equations for simulating ejector flow include standard conservation equations of energy, mass, and momentum, as used in studies by Zhang et al. [35] and Bartosiewicz et al. [36].

Discretization is an integral part of a numerical study. The accuracy and quality of results of any CFD model is highly influenced by the mesh density. A well-refined mesh is a necessity to have negligible discretization error. A mesh independence study was performed by refining the mesh size, and results showing variation of entrainment ratio with number of cells is given in Figure 6. It can be seen that when the number of cells are increased further from 500,000 to 600,000 for a structured mesh, the entrainment ratio rises only by 0.1%, and, therefore, mesh density of 500,000 cells is shortlisted in this study as a tradeoff between accuracy of results and corresponding computation effort (time). Furthermore, use of structured mesh allows much better control over mesh quality, and quadrilateral elements are used for meshing the fluid domain. The visualization of flow domain is shown in Figure 7. The mesh quality at the onset of mixing of primary and secondary fluid is refined to capture the complex flow interaction and oblique shocks. After the mixing, the mixed flow undergoes normal shock, which reduces the flow velocity from supersonic to subsonic [37]. The magnified view of mixing region and normal shock is also given in Figure 7.

Generally, a density-based solver is recommended widely in the literature for simulating compressible fluid flows. However, the recent studies show that a pressure-based solver is also equally capable of solving compressible flows in less time than a density-based solver [28]. Hence, a pressure-based solver with coupled algorithm for pressure–velocity coupling is selected. The stagnation pressure boundary condition was applied at both the primary and secondary inlet. The CFD study was performed using experimental data of Huang et al. [15], and for R141b working fluid, the thermodynamic properties were obtained from the REFPROP v9.1.

The solution of RANS equations involves various turbulent fluid models, and proper selection of these models highly influences the flow behavior inside the ejector. Two turbulence models are shortlisted through the literature: k-ε (k epsilon) and k-ω (k omega) [28]. Among the two models, the k-ω model is widely recommended in the literature, as compared to the k-ε model, because the former provides a more accurate global and local flow feature estimate than the latter [38]. In addition, the k-ε model requires higher number of iterations, and provides poor performance, as compared to the k-ω-sst (k omega sheer stress transport) model in terms of convergence, thermal, and flow field prediction [39]. Also, according to various other CFD studies [26,40], the k-ω-sst model provides better prediction for shock wave location and intensity. Hence, in this paper, a CFD study was conducted using a k-ω-sst turbulence model for all simulations, and summary of numerical setup adopted, boundary conditions detail, and mesh employed is provided in

Table 3. CFD model setup summary.

Mesh
Mesh Type	Structured
Number of elements	500,000
Element type	Quadratic quadrilateral
Boundary Conditions
Primary flow inlet	Pressure inlet
Secondary flow inlet	Pressure inlet
Discharge flow outlet	Pressure outlet
Numerical Model Setup
Solver	Pressure based
Turbulence model	k-ω-sst
Method of initialization	Hybrid
Fluid density	Ideal gas
Working fluid	R141b
Discretization scheme	Second order upwind
Convergence criteria	Residuals <10−6

. The convergence criteria were satisfied typically around 1500 iterations when all the residuals fell below 1 × 10⁻⁶.

2.3. CFD Model Validation

For validating the CFD model, the geometry (AA, AB, AC, AD) from already published experimental work [15] is utilized, where Figure 8 shows fixed geometrical configuration [41] and values of variable parameters, such as A_t, A₄, and A₃, are taken from Table 1.

Under these geometrical specifications and experimental operating conditions (P_g: 0.4–0.604 MPa, P_e: 0.04–0.047 MPa and P_C), the developed CFD model is employed to validate the results obtained through empirical correlation of entrainment ratio and compression ratio. In addition, these comparisons provide the information about the reliability of the developed correlations.

3. Results

For validation, five different geometries were simulated from Huang et al.’s paper, with the same boundary conditions. Simulations were performed, and the Mach number contours and pressure contours are shown in Figure 9 and Figure 10, respectively.

From Mach number contours, it is evident that flow in ejectors is a very complex phenomenon that includes shock waves, oblique shock waves, and turbulent effects. The figure shows the varying flow conditions along the different sections of ejector geometry. Initially, flow enters the convergent part of the primary converging diverging (CD) nozzle, where the flow decelerates until it achieves the choking condition (Ma = 1) at the contraction (throat) of the CD nozzle. The motive flow, while passing through the divergent section of the nozzle, accelerates and continues accelerating, even after exiting the primary nozzle, which is as per the assumption made for the analytical solution of the ejector. After the two flows (primary and secondary) mix together in the mixing section, the velocities of these flows sum up into a single averaged velocity (of mixed flow) along the ejector length. The mixed flow undergoes a normal compression shock wave, which transforms the supersonic flow (Ma > 1) to subsonic (Ma < 1). This mixed flow is then fed into a diffuser, where it is further compressed until it attains condenser pressure.

The secondary flow is accelerated and is then entrained by the motive flow, due to the formation of shear layer between primary and secondary flows, under the influence of high-velocity gradient in between the motive and entrained flow at the CD nozzle exit, as also observed in Figure 9, where the velocity transforms to dark blue, followed by a light blue, and finally to yellow. After the motive flow exits the primary nozzle, it slows down, due to its interaction with the entrained flow, starting from the outer boundary of flow, and, consequently, the entrained flow gets accelerated. As the flow passes through the constant area mixing section, more uniformity in color is observed, which points to good flow mixing.

The variation in flow pressure along the ejector centerline is shown in Figure 10. The motive flow, as it passes through the CD nozzle, undergoes a sudden drop in pressure while it remains steady inside the constant area mixing section. Inside the mixing section, the frictional effects slightly reduce the flow velocity, which, in turn, increases the pressure of flow. At the end, when this fluid (supersonic) is fed into a diffuser, the pressure drops due to the supersonic flow condition, followed by an abrupt and rapid rise of pressure due to the occurrence of normal shock wave, which brings the supersonic flow to subsonic state. The deviation between values of entrainment ratio obtained through the CFD model, empirical correlations, and experimental data are shown in Figure 11.

The comparison of empirical correlation and the experimental data is shown in Figure 12, where the maximum deviation from experimental results is 8.40%.

3.1. Case Study: Simple Refrigeration Machine

A simple refrigeration machine, as shown in Figure 13, was designed and parametrically investigated for a case study. The system is designed for a generator temperature (T_Pri) and evaporator temperature (T_Eva) of 70–100 °C and 10–20 °C, respectively. The ambient sink is considered for condenser, so the condensing temperature (T_Cond) is fixed to 40 °C. The working fluid is taken as R141b, thus the developed correlations are valid. The conditions are summarized in

Table 4. Specifications of refrigeration system.

Parameters	Values
Refrigeration capacity (QCool)	300 W
Required condensing saturation temperature (TCond,Req)	40 °C
Generator saturation temperature (TPri)	70–100 °C
Evaporator saturation temperature (TEva)	10–20 °C

. For the case study, the refrigeration capacity (Q_Cool) is fixed to 300 W. Using the known Q_Cool and T_Eva, the msec can be calculated. Subsequently, the correlations are applied to calculate m_Pri and P_Cond. The calculation procedure is summarized in

Table 5. Summary of calculation procedure.

Component	Input	Output	Equations	Comment
Evaporator	QCool, TEva	msec	$m\_\sec =Q_{Cool}/\Delta h_{fg}$	Since the saturation temperature is known, Δhfg can be calculated.
Ejector	msec, A3/AT, PPri/PEva	ω, PCond, mPri		Using the developed co-relation, the ω and PCond can be calculated. If $P_{Cond}\ne P_{Cond,Req}$ then iterate A3/At.
Generator	TPri, mPri	QAdd	$m_{Pri}=Q_{Add}/\Delta h_{fg}$	Since the saturation temperature is known, Δhfg can be calculated.

. After solving the system, the coefficient of performance (COP) can be calculated. The pumping work was too small, and is, therefore, neglected in calculating the COP.

$$CO{P}_{System}= \frac{Q_{Add}}{Q_{Cool}}$$

(6)

3.2. COP Variation with Generation Temperature and the Evaporation Temperature

The T_Eva and T_Pri have a significant impact on the COP of the system, shown in Figure 13. Both can vary, depending on the heat source and the target temperature. For the system shown in Figure 13, the Q_Add can be supplied using any low-grade heat source, such as waste heat or the solar energy [39]. Furthermore, P_Eva is governed by the required target temperature and can also vary, depending on the requirement. Therefore, a parametric study with respect to T_Eva and T_Pri was conducted, and the results are shown in Figure 14. As can be seen from the figures, higher values of T_Pri and T_Eva both favor the performance of the proposed design. The results are rather obvious, since the increase in evaporation or generation temperature both reduce the load on a refrigeration machine. However, the investigation demonstrates the ease of parametric investigation using empirical correlations for the ejector.

4. Conclusions

The interest in ejector-assisted refrigeration systems has been rekindled, since these compact devices can be operated using low-grade heat. The ejector is one of the crucial components in the ejector-assisted machines, and the need for a reliable model for estimating ejector performance persists. However, the ejector flow is comprehensible, and involves complex flow features like normal shock wave and boundary layer interaction, and, so far, all the 1D thermodynamic models are rather complex and lack the adequate accuracy. Therefore, this study proposed an empirical correlation to predict the performance of an ejector operated with R141b in the range of 0.4–0.604 MPa primary pressure and 0.040–0.047 MPa secondary pressure. The correlation can predict both global flow parameters of ejectors, i.e., entrainment ratio and compression ratio, with an accuracy of 8.4%, and 6.3%, respectively. Furthermore, due to the complex flow behavior in the ejector, a k-ω CFD model was developed. The model was validated against the experimental results. Using the CFD model, Mach number and pressure contours were plotted. The results show that the maximum error in the CFD was 4.29%. Moreover, the validity of correlation was also studied, with the validated CFD results also shown. In the end, a case study was shown that successfully demonstrated the ease of using the correlations for the design of an ejector-assisted refrigeration machine. A parametric study was conducted that concluded that both high primary fluid temperature and evaporation temperature increase the performance of the machine.