A detailed review on CO2 two-phase ejector flow modeling

Ejector-equipped vapor-compression systems for refrigeration and cooling, relying solely on CO2 (R744) as a natural working fluid, are perceived to be an eco-friendly and highly efficient solution for many applications. However, the complexity of two-phase ejector flows makes it very challenging to find realiable and efficient ejector designs. Improved design methods are necessary in order to achieve higher performance in R744 units compared to the traditional compressor-based systems with refrigerants that put a high strain on the environment. Consequently, the development of advanced models and tools for an accurate design of the R744 ejectors has been a highly prioritized research topic. To the best of the authors’ knowledge, the current status of R744 ejector models and their limitations has not been thoroughly evaluated yet. To summarise the current state of the art and knowledge gaps, this work presents an exhaustive overview of the available numerical models applied to R744 two-phase ejectors, i.e. multiphase flow modeling, turbulence aspects, numerical solution methods, applications of models, to further encourage the adoption of R744 vapor-compression solutions. Finally, a thorough discussion of different focus points for future research as well as the main challenges in the field is presented.


Introduction
Growing concern for human impact on the environment has brought about a major shift in the field of Heating, Ventilation, Air Conditioning, and Refrigeration (HVAC&R). The transition from high global warming potential (GWP) working fluids, e.g. hydro-fluorocarbons (HFCs), to environmentally friendly refrigerants is a critical milestone to achieve carbon neutrality. Legislative agreements aiming at HFC phase-down, such as the Kigali Amendment to the Montreal Protocol and the EU F-gas Regulation 517/2014 [1] (European Commission, 2014), have been in force since 2015 and apply for both developed and developing countries. Particular emphasis has been placed on developing solutions using natural refrigerants as these are cheap, available, and have no unknown bi-effects on the atmosphere.
Vapor compression refrigeration systems consist of a circulating refrigerant (working fluid) that absorb latent and sensible heat at a low temperature and releases it at a higher temperature. This is possible by compressing and decompressing the refrigerant to reach appropriate temperature levels. A compressor supplies mechanical work, increasing the refrigerant temperature and pressure. At the higher temperature, the heat is released to the ambient air through a condenser or gas-cooler. The refrigerant is then reduced in pressure and temperature in an expansion device and collects heat through an evaporator. In the simplest cycles, the expansion device is a simple expansion valve. However, these devices do not recover any of the supplied mechanical energy and therefore impose a throttling loss. Alternatively, work-recovery devices, such as expanders or two-phase ejectors are used. This is especially important for refrigerants that operate at high pressures (such as CO 2 ). Of the natural and environmentally friendly refrigerants, CO 2 (R744) stands out as an efficient, long term solution. R744 is characterized by negligible GWP, non-flammability, non-toxicity, low cost, favorable thermophysical properties [2]. Additionally, the R744 cycle can use smaller and more compact components due to its thermodynamic properties. These characteristics make R744 refrigeration systems a hot research topic. Significant work is currently being carried out, highlighting the benefits of R744 for many applications, such as supermarkets [3,4], heat pump units [5][6][7], vehicles [8][9][10], light commercial refrigeration [11,12], tumble dryers [13,14], chillers [15]. R744 HVAC&R systems are rapidly becoming widely accepted for many purposes worldwide. The wide recognition of this ascending technology is furthermore a consequence of its high performance in any operating https://doi.org/10. 1016 T mode, such as at high ambient temperatures (in warm climates). The high efficiency of R744 HVAC&R systems is highly dependent on support by a two-phase ejector. It has been shown that the recovery of part of the available expansion work allows energy savings up to 25% compared to HFC-based systems in supermarkets [16].
Compared to the expanders, two-phase ejectors are characterized by low cost, absence of moving parts (i.e. great reliability), and the ability to handle two-phase flows without risks of damage. Consequently, Lawrence and Elbel claimed that the adoption of an ejector should be favored for HVAC&R units [17].
The conventional trans-critical R744 cycle with ejector and its pressure-enthalpy (P-h) diagram are shown in Fig. 1a and b, respectively. Furthermore, the refrigerant flow and a qualitative pressure and velocity profile are presented in Fig. 1 (a). In this solution, CO 2 exiting the gas cooler as vapor (thermodynamic state 3, identifying the high pressure) is referred to as the primary or motive flow.
The high pressure at the motive inlet accelerates the flow in the motive nozzle converging section to sonic conditions (Ma = 1) at the throat and further accelerates it to super-sonic flow (Ma >1) in the motive nozzle diverging section. The acceleration is coupled to a pressure reduction, which initiates a phase change process of the liquid CO 2 in the nozzle, called flashing [19]. The motive flow fans out from the divergent part of the nozzle creating a low-pressure region (thermodynamic state 4) which drives the CO 2 from the secondary inlet (thermodynamic state 9) into the suction chamber (thermodynamic state 5). The entrained CO 2 stream is generally referred to as the secondary or suction flow. The primary and secondary streams are then mixed in the mixing chamber (thermodynamic state 6) due to the exchange of mass, momentum, and energy (heat). As the flow subsequently decelerates in the diffuser, a part of the residual kinetic energy is converted into a pressure increase (thermodynamic state 7, identifying the intermediate pressure). Thus, the energy lost due to the expansion of the motive flow is recovered and used to increase the pressure of the suction flow to produce the wanted pressure lift of the suction flow.
An ejector-equipped vapor compression system presents two main advantages when compared to the conventional R744 cycle: (I) Reduction in energy consumption thanks to pre-compression of R744 through the ejector from the low evaporator pressure to the intermediate one; (II) Increase in the refrigerating capacity, as the R744 enters into the evaporator with a lower vapor quality and enthalpy. The performance of a two-phase ejector is generally evaluated using the following indicators: the mass entrainment ratio ( ), the suction pressure ratio ( ), the pressure lift (P lift ) and the ejector isentropic efficiency ( ejector ). The mass entrainment ratio ( ) is defined as the ratio of the entrained mass flow rate (m) to the motive mass flow rate and describes the ability of the ejector to entrain the low-pressure secondary flow: = m m 9 3 (1) where the subscript indicates the thermodynamic state referred to in Fig. 1. The suction pressure ratio ( ) is used to evaluate the compression ratio between the ejector outlet pressure and the ejector suction pressure. The pressure lift (P lift ) is used to evaluate the difference between the ejector outlet pressure and the ejector suction pressure.
Lastly, the ejector efficiency ( ejector ) is the ratio between the actual amount of work recovered by the ejector (W r ) and the total work recovery potential for an isentropic process (W r,max ) [ Here, h is the specific enthalpy. The importance of the two-phase ejector for R744 systems was highlighted by Elbel and Lawrence [17], who showed that the adoption of a two-phase ejector is much more beneficial to R744 systems compared to HFC-based units. The R744 ejector efficiency generally ranges between 0.2 and 0.4 [20], whereas the efficiency that of both R410A and R134a ejectors is normally less than 0.2 [17]. Therefore, the characteristics of R744 favorable for expansion recovery, combined with the high energy efficiency offered by the adoption of a two-phase ejector, are further promoting the application of R744 refrigeration system at an industrial scale [21][22][23][24][25][26]. However, due to the significant complexity of ejector flow, connected to the interactions between the motive and entrained flow, the performance of a two-phase ejector is highly dependent on its mechanical design. Furthermore, R744 ejector designs are significantly different from ejectors with other working fluids, due to different thermodynamic properties and component size.
Experimental work is the most reliable approach to identify the optimal design of the ejector and has so far been the most used design methodology. However, it is characterized by some limitations. Primarily, the large time-and resource costs of experimental design are  [18].
K.E. Ringstad, et al. Thermal Science and Engineering Progress 20 (2020) 100647 prohibitive for large scale test campaigns since the performance of each ejector is defined for at least 5 degrees of freedom: motive nozzle pressure and density, suction nozzle inlet pressure, and density, and outlet pressure. Extensive experimental work has been carried out to investigate efficient ejector designs, looking into the geometry of the mixing section [27,28], motive nozzle [29,18], diffuser [28,18], adjustable needles [30,31]. However, as stated by Elbel and Lawrence [17] the specific geometry obtained in each study is only optimal for that specific case and operation. Furthermore, Elbel and Lawrence [17] point out the need for studies that cover generalized ejector geometries. However, it is costly to produce results with large variations in operating conditions, ejector design, and system layout.
Fortunately, many of these limitations can be bridged using experimentally validated numerical modeling. Additionally, such models can provide a better understanding of the flow patterns. Such a solution can allow identification of optimal operation and ejector design based on optimization algorithms.
The first interest in two-phase ejector modeling started in the 1990s using a one-dimensional approach [32] and has since then been extended with additional experimental data and newly developed models to achieve better prediction of two-phase ejector flow. The current state of the art within fluid flow modeling involves solving the full equations of fluid motion in three dimensions, i.e. the Navier-Stokes equations, using computational fluid dynamics (CFD). Computational fluid dynamics is a powerful numerical tool that can capture the local flow behavior realistically and can yield more accurate predictions and more physically realistic solutions than simpler models. Furthermore, CFD is not dependent on extensively tuned parameters based on experiments, which improves the applicability of the models.
For literature on the modeling of other flows with relevance for R744 ejectors, the authors refer to research investigating atomizers [33], diesel injector [34], and accidental leaks in devices such as nuclear reactors, engines or hazardous gas containers [19].
Even though two-phase ejector technologies has been thoroughly reviewed in the past, covering such topics as cycles and applications [16,17,[35][36][37][38], thermodynamic modeling [17,39,40,37,36,41] geometric design [37,38] and historic developments [38,41,42], to the best of the authors' knowledge limited attention was given to the determination of appropriate modeling approaches for two-phase ejectors [43]. R744 two-phase ejectors models have been previously reviewed by Nowak et al. [43], however, several modeling aspects were not covered by this overview. This work attempts to further detail modeling aspects previously left untouched, such as in-depth multiphase modeling, turbulence approaches, model applications, new experimental validation data, numerical methods, and model accuracy comparisons. This work intends to fill this knowledge gap by providing a complete state of knowledge, summarizing the recent progress within R744 ejector modeling.
In this review, first, the state of the art of multiphase flow models are presented and compared. A thorough discussion about the available turbulence models such as different two-equation models and more advanced multiphase models for R744 ejector models is given. Different numerical methods and their significance for accuracy and convergence are examined, and currently available experimental data for model validation is exhaustively reviewed. Lastly, the available modeling strategies and suggestions for further work are summarised.

Flow characteristics
The fluid flow in an R744 two-phase ejector is characterized by multiple two-phase flow regimes in the different sections of the ejector. These flow regimes are dependent on operating pressures and ejector design. Due to this complexity, there is a lack of knowledge of these flow regimes, which poses a modeling challenge.
The flow characteristics in the motive nozzle depend on the nozzle inlet conditions. Two types of expansion paths can be identified: expansion from a supercritical state that intersects the phase envelope close to the critical point, and an expansion that intersects the phase envelope far away from the critical point. The former will be referred to as "near-critical" expansion and the latter as "off-critical" expansion. At near-critical conditions bubble nucleation in the motive nozzle is almost instantaneous, while the phase change is delayed at off-critical conditions due to non-equilibrium effects, discussed in Section 2.2. The bubbles quickly reach the velocity of the surrounding liquid flow. As the bubbly flow leaves the motive nozzle, the flow is supersonic and can go through a series of two-phase shocks. At some point, the flow becomes predominantly vapor and the flow regime inverts from a bubbly flow to a droplet flow. An illustration of such a jet flow is presented in Fig. 2 (Jiang et al. [44]).
Dimensionless numbers are commonly used to evaluate flow characteristics in fluid mechanics. These numbers help in quantifying the scales of different flow phenomena in the flow. A limited number of the most relevant dimensionless numbers will be presented below.
The Reynolds number, Eq. (5), is defined as the ratio between inertial and viscous forces, and is commonly used to quantify turbulent flows. The Weber number, Eq. (6), is defined as the ratio between inertia and surface tension, and is used for discussing the breakup and coalescence of bubbles and droplets. The Ohnesorge number, Eq. (7), is a combined dimensionless number that describes the effects of surface tension, inertia, and viscous forces. The Prandtl number, Eq. (8), is defined by the ratio between momentum and heat diffusion. The gas and liquid Prandtl numbers of R744 vary from Pr 1.5 6 in the temperature range 0-30°C. The Biot number, Eq. (9), is a number compare the heat transfer within and at the surface of a bubble or droplet. For low Biot numbers, Bi 1 the temperature inside a bubble or droplet can be assumed homogeneous. The Nusselt number, Eq. (10), is similarly defined, however, compare convective and conductive heat transfer at the surface.
= Pr c µ k , K.E. Ringstad, et al. Thermal Science and Engineering Progress 20 (2020) 100647 Here, U is the characteristic velocity, l is a characteristic length, is the surface tension, c p is the specific heat, µ is the fluid viscosity, h c is the convective heat transfer coefficient, and is the density. k is the thermal conductivity. Subscripts are used for numbers where two bodies are in contact, i.e. bubbles or droplets, i and e denote internal and external heat transfer for that body, respectively.
Based on numerical values for an R744 ejector by Smolka et al. [45], the dimensionless numbers were calculated at two locations; the centerline after the motive nozzle exit and the end of the mixing chamber/ start of the diffuser. The presented quality, pressure, and velocity plots in Smolka et al. [45] were used to estimate the dimensionless numbers. The REFPROP library [46] Table 2.
From these calculations the estimated Reynolds numbers are very high due to the high velocities in the flow, indicating that the flow is highly turbulent, both in the mixing chamber and in the motive nozzle. Furthermore, the very high Weber numbers indicate that bubble breakup will happen very rapidly.

A discussion on non-equilibrium
Knowledge of non-equilibrium effects is critical to properly understand two-phase ejectors. In a two-phase flow, there exists non-equilibrium states can be divided into thermodynamic and transport nonequilibriums. Thermodynamic non-equilibrium relates to the superheated or subcooled state of a liquid or gas, respectively. Transport nonequilibrium relates to differing temperatures, pressures, or velocities of the phases.
During rapid depressurization, as in the motive of an R744 ejector, the saturation temperature will drop below the liquid temperature, superheating the liquid. The liquid starts to evaporate until the equilibrium state is reached. The degree of superheat is therefore limited by the phase change mechanism. Beyond the homogeneous nucleation line any perturbation will instantly force the phase change, imposing the upper limit on superheating. This is illustrated in Fig. 3, where the saturation and homogeneous nucleation lines of R744 are illustrated. Three cases of ejector motive conditions are plotted with an illustrated isenthalpic expansion process, corresponding to the expansion process in the motive nozzle (points 3-4 in Fig. 1). Case 1 is a near-critical expansion, and cases 2 and 3 are off-critical expansions, in terms of the previously established notions. This figure demonstrates that as the degree of sub-cooling at the motive nozzle inlet increases (lower enthalpy) and moves towards off-critical expansion, the theoretically achievable superheat increases, and therefore the relevance of chemical potential non-equilibrium. As mentioned in Section 2.1, it is clear that for near-critical expansion the possible thermodynamic non-equilibrium is very small, and therefore the phase change will occur very rapidly. The figure was generated based on the framework established by Wilhelmsen et al. [47] and Aursand et al. [48].
Transport (or mechanical) non-equilibrium relates to flow characteristics, as well as fluid properties. The commonly considered disequilibria are velocity, temperature, and pressure. Thermal non-equilibrium refers to a state where the phases exist at differing temperatures. In this case, heat transfer between the bubbles and the surrounding liquid must be considered. These states can also affect the thermodynamic non-equilibrium as two phases can co-exist at thermodynamic equilibrium, however at different temperatures. Furthermore, within a bubble or droplet the temperature is not homogeneously distributed. A droplet or a bubble may, for example, contain regions where the fluid is at equilibrium and other regions where it is at non-equilibrium. For low Biot numbers, the heat transfer inside a bubble or droplet is much faster than heat transfer to the surrounding fluid. In this case, the temperature within the particulate can be assumed homogeneous. Still, limited knowledge is available on the Biot number for R744 ejector flows, and further research is needed. This was investigated by Bartosiewicz and Seynhaeve [49], where it is shown that a non-equilibrium temperature front may in some cases only penetrate a short distance into a droplet. This is illustrated in Fig. 4, where different control volume resolutions and the corresponding resolved temperature distributions are presented. The effects of pressure non-equilibrium are often neglected due to the very short time scale of these disequilibria.
A summary of the currently available experimental literature discussing the importance of non-equilibrium modeling in R744 ejectors is given in this paper. Firstly, measurements of wall pressure in a converging-diverging nozzle with R744 were conducted by Nakagawa et al. [50]. It was found that the observed pressures were close to equilibrium conditions. However, for low operating temperatures, the experiments suggested the occurrence of non-equilibrium phase change as pressures lower than saturation pressures were observed in the nozzle. Secondly, shock waves in the nozzle were further discussed by Berana et al. [51]. They put forward that equilibrium flow is an ideal limiting case, and that for some conditions dispersed and pseudo shocks were observed. These shocks are weaker and thicker than equilibrium shocks. Berana et al. [51] suggested that velocity non-equilibrium could significantly impact shock solutions due to its effect on the speed of sound. Thirdly, Li et al. [52] used the wall pressure and a visualization technique to investigate the phase change point for different operating conditions in an R744 ejector. Also here, the measurements showed conditions approximately at equilibrium conditions. The aforementioned studies were conducted at supercritical motive pressures. These investigations agree that the expansion in a nozzle is in general at homogeneous equilibrium. However, non-equilibrium effects have been suggested by the aforementioned authors to be important to capture transcritical flow at certain operating conditions.

Multiphase flow models for R744 ejectors
Accurate multiphase flow modeling is generally considered as the most challenging aspect of R744 ejector modeling and the key factor for successful modeling of R744 ejectors. Compared to single-phase flows, multiphase flows introduce complex interactions at the interface between the liquid and vapor phases. In addition, R744 ejector flow involves sudden and complex phase change through flashing. The reader is referred to Polanco et al. [53] and Liao and Lucas [19] for additional information on flashing flow modeling. As previously mentioned, R744 ejector modeling is not yet fully explored, thus many simplifying assumptions are used in the available models. Assumptions such as equilibrium of velocity i.e. homogeneous flow (assumed in most currently available models), thermal equilibrium (assumed in mixture and equilibrium models), thermodynamic/chemical equilibrium (assumed in Equilibrium models) and pressure equilibrium (assumed in all currently available models) are considered to achieve well-posed and consistent models.
Ejector models are predominantly based on inputs of operating conditions, i.e. temperature and pressure, at the motive, suction, and outlet of the ejector. The two primary objectives of an ejector are to provide pressure lift and entrainment of the secondary fluid. Therefore, most models attempt to predict the entrainment ratio ( ) from the pressures and temperatures at the boundaries. These variables can easily be obtained experimentally, and are the variables used for control of most refrigeration systems.
The accuracy of current R744 ejector models varies significantly based on operating conditions, ejector geometry, and model complexity. As an example, thermodynamic relation models (0-D) can achieve accuracy within a 10-15% error in the motive mass flow rate [54]. However, such models have typically a narrow range of validity in terms of varying ejector geometry and operating conditions, and are based on tuning from previous experimental data. On the other hand, more complex models (CFD) can achieve higher accuracies with a much wider range of validity, at the expense of much higher computational costs.
The understanding of two-phase ejector models is further complicated by their high sensitivity to the model inputs. Smolka et al. [45] discussed the effects of small changes in the numerical boundary conditions. Changing the boundary conditions within the experimental uncertainty had a large impact on the model's predicted mass flow rate. This result is summarised in Table 3, and highlights the importance of improved models and high accuracy experimental data. The reason for the motive mass flow rate being easier to predict is due to the flow becoming choked when reaching supersonic speeds. A chocked flow is  insensitive to the downstream conditions and is, therefore, easier to predict. On the other hand, the prediction of the suction mass flow rate is highly sensitive to downstream conditions. This is because the downstream conditions, such as mixing chamber pressure, is one of the main driving mechanism of the suction flow. It is therefore hard to say whether the low model accuracy for suction flow is due to wrong modeling of the mixing process, i.e. turbulence, or the wrong modeling of the motive flow downstream of being chocked. This will be further discussed in the following sections. An illustrative overview of models found in R744 ejector literature is presented in Fig. 5. The models are classified according to the multiphase aspects considered. Firstly, the models are classified according to the number of fluid flow equation sets solved; in the two-fluid approach, one set of equations is solved for each phase, while in the pseudo-fluid approach the equations of both phases are averaged into one set. Thus far, all currently available ejector models fall under this latter category. Secondly, the models are classified according to whether or not the pseudo fluid transport properties are evaluated assuming equilibrium or not. Thirdly, the models are classified according to their treatment of velocity non-equilibrium. Only one study has thus far evaluated velocity slip [55], and most models consider the flow as homogeneous.

Homogeneous flow models
In multiphase flow modeling, a commonly used assumption is that of homogeneous flow. This approach assumes the phases to be at mechanical equilibrium, i.e. both phases are described with a single velocity-(u ) and pressure-field (P) in order to reduce the complexity of the model, as it neglects the modeling of the slip velocity. Essentially, the two phases can then be treated as a single pseudo-fluid with transport properties derived according to an averaging procedure. This pseudo-fluid will be governed by the equations of fluid motion, Eqs. (11)-(12) as well as an energy equation, Eq. (13).
Here the Einstein notation is used with subscript-indexes i and j, and the subscript m indicates the pseudo-fluid mixture properties. u p E h q , , , , , refer to the density, velocity, pressure, total energy, enthalpy and heat flux, respectively. The effective stress tensor ij,eff is the laminar (Newtonian) and turbulent stress tensors combined, The homogeneous flow approach is prevalent in R744 two-phase ejector models and will, therefore, be a focal point of this review. However, the transport properties of a mixture of two phases have to be defined as they are derived from the averaging of the equations. Therefore, estimates of these properties must be carefully considered. For example, the mixture viscosity is typically defined as a volumeweighted average of the two phases µ l , and µ v . In this case, the mixture viscosity, µ m , is a function of the liquid and vapor viscosities and the vapor fraction: To evaluate this property (Eq. 14), three factors must be considered: (I-Phase-fraction) The evaluation of phase fraction of the mixture, (II-Properties) The fluid properties of the phases, µ l and µ v in Eq. (14). (III-Averaging) Mixture averaging procedure i.e. the function f in Eq. (14). As an example is a volume-weighted average presented below (Eq. (14) These three factors will be discussed in terms of each model. Factors (II-Properties) and (III-Averaging) are discussed further in Section 3.3.

Homogeneous equilibrium
One simple, yet reasonable solution to evaluate the phase fraction (I-Phase-fraction) and the fluid properties (II-Properties) is to assume thermodynamic equilibrium of the phases at all points in the flow. This is the main concept adopted in the Homogeneous Equilibrium Model (HEM). This model has been extensively used in the literature to model R744 ejector flow structure and characteristics [31,45,56,57].
Bulinski et al. [58] conducted an early exploratory investigation in 2010 to investigate two different multiphase CFD models for an R744 ejector. A heterogeneous model, similar to the mixture model, and a temperature-based homogeneous model was implemented into the ANSYS Fluent software and tested. Unfortunately, the models performed poorly and the model results could not reproduce the experimental data. This led to investigating a homogeneous equilibrium model formulation based on total enthalpy. One preferable property of equilibrium flow is that the pressure and enthalpy uniquely define the thermodynamic state, and thus the properties in the two-phase dome. Properties are typically divided into thermophysical-(Eq. 16) and transport (Eq. 17) properties. , where µ k c , , , , p are the pseudo-fluid density, vapor volume fraction, kinematic viscosity, thermal conductivity and heat capacity, respectively. In general, to achieve the most accurate thermodynamic relations, it is preferred to use density-energy formulations as they always conserve mass and energy, however these are less practical for comparisons with experimental results and require more computations.
This property was exploited in the work by Smolka et al. [45]. Furthering their previous work [58], a full 3D steady-state CFD model for R744 ejectors was implemented such that all properties were defined using the pressure and total enthalpy. To obtain the total enthalpy, an alternative formulation of the energy equation, Eq. (18), was implemented using user-defined functions in the ANSYS Fluent software [59].
In this equation, h is the specific enthalpy, u is the velocity vector, eff the effective diffusion coefficient. The source terms S h1,2,3 describe the mechanical energy, the irreversible dissipation of the kinetic energy variations, and the dissipation of the turbulent kinetic energy, respectively [45]. The two-phase mixture properties were evaluated in Smolka et al. [45] using the REFPROP library [46]. The results were validated with experimental data obtained from a lab-test campaign carried out by the SINTEF Energy Research laboratory in Trondheim (Norway) [60][61][62]. The experimental results were obtained using an ejector implemented in an R744 transcritical heat pump system under varied conditions. It was found that the non-symmetries were significant and that 3D flow modeling should be performed, in contrast to previous investigations with 2D flow [63]. The disagreeing results indicate that the significance of 3D effects depend on the ejector design and should be verified on a case to case basis. This is further discussed in Section 5.2. The authors found that global parameters (e.g. mass flow rate) were well approximated by the HEM approach for operation near the critical point, i.e. for near-critical expansion. The average deviation from experimental data was 5.6% and 10.1% for the motive and suction mass flow rates (MFR), respectively.
Later, Lucas et al. [57] presented a similar CFD model using the total enthalpy HEM approach. This model was implemented in the OpenFOAM framework [64] with the TEMO-media library. The investigations compared the results of the pressure recovery with those carried out by Lucas et al. [65] in earlier experimental work. When operating without suction flow up to 10% deviation in pressure distribution between model and experimental results were observed. However, this error increased to 20% when the suction flow was included. It was concluded that the larger deviation when simulating with the suction flow was due to inaccurate modeling of the pressure losses associated with flow mixing.
Furthermore, a HEM was implemented in Giacomelli et al. [66,67 using look-up tables implemented into ANSYS Fluent from the RE-FPROP library. They conclude that the HEM is an efficient tool for achieving reasonable results. Still, they state that further model development, in terms of meta-stable properties, is required to properly describe the flow physics.
Recently, Fang et al. [68] implemented a HEM based on a thermodynamic look-up table [69] with an inviscid density-based solver in the OpenFOAM framework, (discussed in Section 5), to model converging-diverging (CD) nozzle with R744. The results were compared to the wall pressure profile [27] with reasonable accuracy for the larger expansion angle. Interestingly, the model performed similarly for the sub-cooled as for the super-critical motive conditions. They state that work is ongoing for an R744 ejector simulation.
Following their previous work [45], Palacz et al. [56] investigated the application range of the HEM approach within typical supermarket refrigeration conditions. The model was validated with experimental results for 24 running modes. The results showed that for operating conditions close to or above the critical point of R744, the HEM approach can accurately reproduce the experimental results (within 5% error of mass flow rate). However, as the temperature and pressure are reduced below the critical point, to off-critical expansion, the HEM approach accuracy is reduced to approximately 50% deviation in motive mass flow rate. The significant errors were attributed to the nonequilibrium effects that are not considered by the HEM approach, which had been argued for in previous ejector models [70].

Homogeneous relaxation
Palacz et al. [56] concluded that non-equilibrium modeling is necessary to model the whole range of common supermarket R744 applications. At equilibrium the pressure and enthalpy define the vapor quality of the flow. However, at non-equilibrium an additional transport equation for the vapour fraction, Eq. (19), with a phase change model is needed to evaluate (I-Phase-fraction).
where is a phase source term modeling phase change. One potential extension to the HEM for non-equilibrium conditions is the homogeneous relaxation model (HRM) introduced by Bilicki and Kestin [71].
Similarly to the HEM, the HRM assumes homogeneous multiphase flow. The HRM treats the phase change as a relaxation process toward the equilibrium vapor-quality: where µ is a relaxation factor, is the mixture density, is the instantaneous vapor mass fraction and¯e q the time-averaged equilibrium vapor mass fraction. The discussed relaxation introduces a delay to the onset of phase-change with a time-scale referred to as the relaxation time, µ . In relation to the previously established notation, this yields a more accurate description of the phase fraction (I-Phase-fraction). This improvement is, however, based on an appropriate estimate of the relaxation time scale.
Downar-Zapolski et al. [72] defined the relaxation time as described in presented in Eq. (21), validated for the case of flashing water.
where 0 is the initial relaxation, a and b are empirical coefficients, and is the mixture void fraction parameter defined as: = , sl sl sv (22) where sl , sv are the saturated liquid and vapor density. This formulation was extended for R744 by Angielczyk et al. [70], who investigated a formulation of , the pressure parameter in Eq.
where p c is the critical pressure at the given conditions and the subscript sat indicates saturated conditions. The mentioned authors found that the appropriate empirical coefficients were, 1.76, 2.14 07 0 . Such an HRM was implemented by Colarossi et al. [63] into a 2D CFD model using the OpenFOAM framework [64] for R744 ejectors. The study aimed to investigate the presence of non-equilibrium effects. As an initial approach, the authors used the relaxation time for flashing water [72]. A comparison with the experimental results obtained for the pressure recovery performed by Nakagawa et al. [27] revealed an average error of 18.6% and a maximum error of 50%. It was also stated that the expansion follows a path of equilibrium states, and concluded that non-equilibrium effects were negligible. The study was evaluated at supercritical conditions so these results are not contradicting the results of Palacz et al. [56].
In a later study, Palacz et al. [73] compared an HRM and a HEM approach. The constant relaxation time parameters suggested by Angielczyk et al. [70] was employed. The investigation showed that the HRM outperformed the HEM for operating regimes distant from the critical point. However, the HEM was more accurate than the HRM approach at supercritical conditions. Overall up to 5% improved accuracy was observed for the HRM only in the sub-critical range. This improvement was found unsatisfactory, and it was concluded that further investigation into the relaxation time parameter should be conducted to achieve higher accuracies.
Recently, Colombo et al. [74] presented an HRM approach based on the relaxation parameter of Angielczyk et al. [70] and compared with the experimental results by Palacz et al. [75]. Relatively low errors were observed for the three operating points (2-14% error in motive K.E. Ringstad, et al. Thermal Science and Engineering Progress 20 (2020) 100647 MFR, 2-16% error in suction MFR), however, investigations with at wider operating ranges are needed for further validation. The investigations carried out by Palacz et al. [73] revealed that the assumption of a constant relaxation time parameters decreased the model accuracy for some operating conditions. To generalize the HRM for a wider range of operating conditions investigation into a variable relaxation time parameters was conducted by Haida et al. [61]. The parameters in Eq. (21), ie.
a , 0 and b, were thoroughly investigated to understand their impact on the flow. Furthermore, an optimization algorithm was used to construct an optimal set of relaxation time parameters for different operating ranges, to better replicate experimental results and achieve higher accuracy. The motive nozzle conditions were divided into a high, a medium, and a low-pressure range. In each pressure range, the parameters were optimized to minimize the discrepancy in the results. The empirical coefficients that best reproduced the experimental results were identified as follows: The optimization was based on a set of operating conditions and validated with experimental results outside this set. The modified HRM was able to produce similar accuracy to the HEM approach for transcritical operating conditions while retaining the benefits of the HRM for sub-critical conditions. However, further work is needed to validate the transferability of this model to other geometries.

Mixture models
A modeling approach that has recently garnered attention from the research community is the mixture model. These 4-equation models aim to model the phase change mechanisms in the flow. Here, the transport equation for the volume fraction (Eq. 19) is solved with the mass transfer terms for evaporation and condensation explicitly included: where e and c are the mass transfer sources due to evaporation and condensation, respectively. Such an improved estimation of the phase composition, i.e. (I-Phase-fraction), allows for a more accurate estimate of the mixture properties. The properties are calculated based on a mass weighted average (Eq. 25) for the thermodynamic variables, e.g. enthalpy or total energy, or volume weighted average (Eq. 26) for the variables mixture density, molecular viscosity or thermal conductivity: where , and are generic properties, and are the volume and mass fractions respectively. The subscripts l, v, and m correspond to the liquid, vapor and mixture properties, respectively. Different approaches have been proposed to model the phase change mechanisms. Yazdani et al. [55] implemented a mixture model into the ANSYS Fluent framework to model the R744 ejector. Here, twophase transfer mechanisms were considered (I-Phase-fraction), ie. in the terms e and c of Eq. (24); a cavitation model adopted from the work by Singhal et al. [76], and a boiling model based on kinetic theory [77,78]. However, as noted by Giacomelli et al. [62] the superposition of these two phase-change mechanisms is not justified during flashing, as cavitation and boiling cannot be considered as two independent mechanisms.
Yazdani et al. [55] implemented different methods to estimate the thermodynamic fluid properties that enter the interpolation schemes, Eq. (25)- (26). The thermal conductivity and viscosity were set as a constant mean value for each phase, the liquid and vapor specific heats were set at saturated conditions, and the density of the liquid and vapor phase were estimated based on the REFPROP-database [46] and a Peng-Robinson equation of state. Additionally, this paper introduced a drift flux model to account for velocity-slip between the phases. This will be further discussed in Section 3.4.
Recently, a more extensive mixture model was presented by Giacomelli et al. [67] to model the nozzle flow in the R744 ejector ( [50]). In this work, the full set of thermodynamic fluid properties for both phases were found through look-up tables based on the NIST REFPROP library. This was implemented by introducing two non-interacting "species", where one represents the R744 vapor phase and one represents the R744 liquid phase. This allows for each phase to be defined with a User Defined Real Gas Model (UDRGM). Giacomelli et al.
The two parameters e and c require tuning, and can be analogously compared to the relaxation time. This model is similar to the boiling model employed by Yazdani et al. [55], through the Clausius-Clapeyron equation.
The phase change models used Yazdani et al. [55] and Giacomelli et al. [67] indirectly assume an interface area, A i , given the form: ie. assuming spherical bubbles or droplets of uniform diameter d B modified by the phase fraction. This assumption is likely to become invalid, especially in the motive nozzle. Giacomelli et al. [62] state that further work will be devoted to including interface density tracking to overcome the limitation of this assumption. Similar mixture models have also been implemented for flashing steam in converging-diverging nozzles [79], two-phase R600a [80] and LNG-BOG [81] ejectors. The discussed mixture model was recently applied to simulate an R744 ejector by Giacomelli et al. [62]. Only a 2% error was observed for the motive mass flow rate showing the high accuracy of the model, however, significant deviations are still observed for the suction mass flow rate (10-17%). The parameter e was found to have a significant impact on the model accuracy in terms of mass flow rate. The parameter c was found to have a negligible impact on the mass flow rate, however, up to 11% higher values of outlet volume fraction were observed as this value was varied from = 0.1 c to = 10000 c . Giacomelli et al. [62] compared their results with those obtained using a HEM and found that the mixture model performed significantly better (19% error for the mixture compared to 48% error for the HEM in terms of entrainment ratio). The difference is flow pattern prediction is depicted in Fig. 6, where the difference between the density field predicted by the HEM and mixture is highlighted. This illustrates the effect of relaxing the phase transition, where the HEM shows sharp discontinuities in density, the mixture model produces a smoother transition. However, the model by Giacomelli et al. [62] displayed problems with numerical stability and slow convergence. Giacomelli et al. [62] reported convergence times up to 8-10 days on a 12-core workstation, which was close to ten times higher than a similar HEM. This slow convergence rate is prohibitively high for applications such as optimization.

Delayed equilibrium
The mixture model allows for the existence of meta-stable and saturated states of liquid and vapor phases, whereas the HEM allows only for saturated conditions of both phases. An alternative approach is the Delayed Equilibrium Model (DEM), where the liquid phase is treated as a combination of saturated and meta-stable liquid. This coexistence is supported by the assumption is that only a fraction of the liquid is held where the subscripts s m l , , , and g indicate saturated condition, metastable condition, liquid, and gas, respectively. This approach has been used in 1D modeling in an R744 converging-diverging nozzle to predict the critical mass flow rate by Angielczyk et al. [82] and to model an R744 ejector by Banasiak and Hafner [83]. Angielczyk et al. [82] compared the pressure distribution of the HEM and DEM in a converging-diverging nozzle. Their results showed that the non-equilibrium introduced in the DEM reduced the model accuracy, however, the limiting case of equilibrium (the HEM approach) still revealed significant discrepancies.
Later, Banasiak and Hafner [84] combined the DEM with both homogeneous and heterogeneous nucleation theory to model an R744 ejector. The authors divided the flow path of the expanding refrigerant into three sections, based on the fluid's thermodynamic state and degree of thermal equilibrium reached: (i) the thermodynamically stable single-phase section (containing either pure liquid or pure vapor), (ii) the metastable two-phase section (containing a mixture of the metastable fluid, saturated liquid, and saturated gas), and (iii) the equilibrium two-phase section (containing only the equilibrium liquid-gas mixture). The fluid bulk temperature in the different regions was defined and calculated according to the fluid properties except in the twophase metastable zone, where two temperatures should be distinguished: the metastable liquid temperature and the equilibrium temperature. The authors replaced the isentropic expansion of the liquid as formulated in the work by Banasiak and Hafner [85], with a more general approach, independent of any transition trails, allowing for a more realistic assumption for the non-isentropic expansion profiles.

Modeling non-equilibrium
A full compressible two-phase model introduces seven degrees of freedom, independent of the turbulence modeling [86]. Such a model allows for disequilibrium of pressure (subscript p), velocity (subscript u), temperature (subscript T) and chemical potential (subscript µ) [87]. The time scale to reach equilibrium for each of these variables is often denoted . The relaxation time discussed for the HRM and mixture model is the chemical potential relaxation time, µ . An analysis of the magnitudes of these time scales [88,89] for flashing water found that the thermal and chemical relaxation times were dominant in comparison with velocity and pressure relaxation times.
The complexity of the model and the required sub-models increase as additional non-equilibrium are introduced. The models are typically arranged into a hierarchy [90], describing the different combinations of relaxation models available. Model choice also changes the predicted speed sound [90,87]. The introduction of a velocity disequilibrium is discussed in Section 3.4. Pressure non-equilibrium is often neglected in most multiphase flow models, due to the rapid response of pressure waves. However, it is considered a challenge to show that problems with pressure non-equilibrium can be well-posed [91]. Thermal nonequilibrium is treated in the two-fluid models, where one energy equation is solved for each of the phases, and coupled by sub-models for heat transfer. The effects of thermal non-equilibrium are dependent on the mesh resolution (i.e. the control volume size). If the mesh is refined up to a resolution where each cell contains primarily one phase, see Fig. 4, thermal disequilibrium can be neglected. Potentially, it may be that both liquid and vapor phase exist in chemical potential equilibrium, however at differing temperatures (i.e. thermal non-equilibrium). Such a flow exists at equilibrium however will be poorly described by any current R744 model. Table 4 presents an overview of current R744 models and which non-equilibrium states are considered in each model.

Two-phase mixture properties
Different models and data sets for R744 properties are available, for a thorough overview see Banasiak and Hafner [92]. When approximating the fluid properties of R744 (II-Properties), the property library REFPROP [46] has been extensively used. REFPROP uses the equation of state (EOS) by [93] based on the Helmholtz equation for the determination of R744 properties. This EOS is widely regarded for its high accuracy [92], however, it is highly CPU intensive. Thus, the published literature has used look-up tables for more efficient simulations [69,68,45,67]. The accuracy of these tables must be considered, which is a trade-off between storage size and accuracy. As an example, Giacomelli et al. [67] found errors in property estimate up to 1.6% in their coarse look-up table. This look-up table was, however, preferred to reduce computation time. Alternatively, look-up tables based on  variable step size is preferred, such as implemented by Banasiak and Hafner [69]. This is especially important near the critical point of R744, where large variations in properties occur. Since the work by Banasiak and Hafner [93], further work has improved upon the modeling of fluid properties of R744. Especially on the modeling of liquid and near-critical viscosity has been improved [94][95][96], however uncertainties up to 3-4% in liquid viscosity are still observed in newer correlations. The effect of bulk viscosity (volume viscosity) can be significant for compressible flow with poly-atomic gases. Fang et al. [68] studied the effect of bulk viscosity on a converging-diverging nozzle with R744 using CFD, concluding that the bulk viscosity had no noticeable effect on the flow. When mixture properties are considered, it is important to consider appropriate averaging laws for the mixture properties [97], which corresponds to (III-Averaging) in the previously established notion. Due to the high velocity and often supersonic flow in the ejector, the evaluation of the speed of sound can significantly impact model results and/or convergence. A recent discussion on critical two-phase speed of sound for different models in pipe-flow was presented by De Lorenzo et al. [89], Fang et al. [69], and Lund [90]. In density-based models [68], the speed of sound directly modifies the propagation speed of acoustic waves and is critical for proper model results. On the other hand, for pressure-based solvers, such as the mixture model [67] and HRM [61], the speed of sound does not explicitly enter the solution algorithm. However, improper speed of sound models can still cause slow convergence or instabilities [62]. In the mixture model presented by Giacomelli et al. [67], the speed of sound could not be set as an independent variable due to limitations in the ANSYS Fluent software [59]. In this software, the Wallis speed of sound model is built into the UDRGM approach [98]. Giacomelli et al. [67] compared the Wallis model with a model by Brennen [99] by calculating the solution speed of sound in post-processing. Their investigation suggests that the Brennen model may be more accurate than the Wallis model. Zhu and Elbel [24] state that the presence of fluids other than R744 in an ejector can dramatically affect the flow physics. As an example, oil flow through the ejector could impact ejector performance. Numerical simulations can be useful to model these effects and can help develop novel solutions to the treatment of oil in R744 ejector systems.

Velocity slip models
In Section 3 the homogeneous flow assumption, where all phases move with the same velocity, was considered. However, velocity slip (velocity non-equilibrium) can have a significant impact on the speed of sound [100,90], two-phase turbulence [101], and shock-wave pressure distribution [55].
To incorporate velocity slip, the models can be extended by including additional terms in the void-fraction (Eq. 19) and momentum equations (Eq. 12) to model the phase velocity slip. Velocity slip is considered in terms of the drift velocity, ie. the relative velocity between a phase (p) and the mixture (m). This is formulated as: or formulated in terms of relative velocities, where subscript k is the phase index 1 to n, subscript m referres to the mixture velocity. Several relations have been presented to estimate the velocity slip in the literature on flashing flows, see the review by Liao and Lucas [19]. Typically the velocity slip models account for drag related forces, neglecting other two-phase interactions. Still, research considering the phase velocity slip condition is very limited. Yazdani et al. [55] introduced a simple velocity slip model based on the formulation for mixture models proposed by Manninen et al. [102]. The drag formulation used was based on the relations of Schiller and Naumann [103], where drag is considered as a modified Stokes drag coefficient for low particle Reynolds number < Re 1000 p , and a constant drag coefficient at > Re 1000 p . The particle Reynolds number is the Reynolds number, Eq. (5), based on the slip-velocity and particle diameter. Yazdani et al. [55] demonstrated that the slip model produced pressure waves at the nozzle exit, which were smoothed by the homogeneous flow model. However, the slip model was found to only have a minor effect on global ejector performance.
The relative importance of velocity slip depends on the velocity relaxation time, v . An estimate of this time scale can be found by considering Stokes drag on a spherical particle. The relevant time scale for velocity slip is the drag relaxation time, ie. the time to accelerate a particulate to the surrounding fluid velocity. Considering the Newtonian (laminar) flow regime at a constant drag coefficient, C D , yields [102]: Hence, the ratio of the particulate to continuum density is critical for velocity slip. For bubbly flows, where the particulate is at a lower density than the surrounding fluid, this ratio is low and bubbles will quickly be accelerated to the velocity of the surrounding liquid. As the flow inverts from a bubbly flow regime with increasing vapor fraction, the flow regime turns into a droplet flow. Compared to bubbles, droplets, on the other hand, will be largely independent of the surrounding gas velocity.

Advanced two-phase modeling
The models presented so far (HEM, HRM, mixture, DEM) have considered the two-phase problem by treating the phase change and phase slip indirectly using supplementary models, coupled to the pseudo-fluid solver. Alternatively, these issues can be treated more directly through the use of a two-fluid model (TFM) [91]. The two-fluid model treats each of the two phases as a separate fluid. This involves the use of separate equations for mass, momentum, and energy transfer for both fluids. In 3D CFD models, this yields a set of 10 equations; 2 for mass, 6 for momentum (one for each direction and phase), and 2 for energy. Such an approach has the benefit of being able to capture nonequilibria (such as temperature and velocity) between the phases directly. However, additional modeling is required to capture the interactions between the phases. The TFM was tried by Menegay [104] in 1998 for R134a ejectors, however, the CFD model was eventually significantly simplified due to its complexity.
Many additional effects of two-phase flow have not yet been thoroughly investigated for R744 ejectors, such as the effects of bubble and droplet collisions, interphase drag, and jet break up (atomization). As an example, the breakup process after a nozzle displays a complex flow pattern, illustrated in Fig. 2. Furthermore, these effects are interdependent with multiphase turbulent effects, which increases problem complexity. Additionally, the close relationships for multiphase flows are strongly dependent on the flow regime, as discussed in Section 2.1. Unfortunately, more advanced investigations will typically require more accurate experimental data for validation. Such experimental data is not yet available, see Section 6.

Comparison of multiphase models
This paper has mainly discussed three approaches for advanced modeling of R744 ejectors: HEM, HRM, and the mixture approach. A comparison of the errors is presented in Fig. 7 and summarized in Appendix A. These data are compiled from the studies carried out by Giacomelli et al. [62], Palacz et al. [73], and Haida et al. [61]. To the best of the authors' knowledge, no such comparison has been done in the literature. Ejector geometry has not been considered in this K.E. Ringstad, et al. Thermal Science and Engineering Progress 20 (2020) 100647 comparison and could potentially have a significant impact on model accuracy.
In Fig. 7 the modeling error, i.e. the discrepancy between model and experimental data, is presented for the different studies. The points are color-coded according to the percentage model error in the motive MFR. Fig. 7 illustrates that current R744 ejector models are able to reasonably capture near-critical motive flow conditions. The recently developed mixture model [62] shows a very low motive MFR discrepancy (error 0.1-2%). However, as can be observed in Fig. 7 this model has not yet been tested for pressures below critical. Moreover, the method shows promise for future development. The motive MFR errors increase for lower motive pressures, as the expansion moves from near-critical to off-critical. As previously discussed in this review, this is likely due to the increased effect of non-equilibrium conditions during off-critical operation.
The average errors of the mentioned models are shown in Table 5  ). The lowest average MFR error was observed for the mixture model both for suction and motive flow rate at supercritical conditions. Nevertheless, this approach should be further tested for sub-critical pressure conditions for a full comparison. Currently, the modified HRM formulation has demonstrated the best accuracy in the sub-critical pressures. However, this approach still significantly deviates from experimental measurements, especially for suction MFR at low pressures.
From Appendix A and Table 5 it is clear that the suction MFR is still not well described by any current model. The relative error in suction MFR varies from 20-100%. The challenge of correctly capturing the suction mass-flow rate stems from its complicated governing mechanisms. Firstly, the suction mass flow rate is induced by the motive flow, which means that errors in the motive flow compound for the suction flow. The motive flow rate can be predicted by a wide range of models, as it is primarily affected by the upstream and downstream pressures. In comparison, the suction flow is more sensitive to local flow phenomena, such as the local velocity, vapor quality, turbulence, and pressure. This means that even if the global mass flow rate of the motive flow is well predicted, if the local solution of the motive flow is not correctly reproduced, it will give cause errors for the suction flow. Secondly, the physics that produces the motive and suction flow is very different. The motive flow is decided primarily by thermodynamics, phase-change, and the motive pressure. While the suction flow is entrained by the motive flow, which is governed primarily by turbulence, the local velocity field, and multiphase momentum interaction.
A summary of the main advantages and challenges of each of these approaches is presented in Table 6. The more recently implemented methods (Modified HRM and mixture model) have been less extensively tested in the literature than the classic HEM. The HEM is generally regarded to produce reasonably accurate results at super-critical conditions. Furthermore, the simplicity of the HEM helps the stability of the model, which is important for use in optimization algorithms. However, this model fails to describe the ejector operation at lowpressure motive conditions. This effect is thought to be connected to relaxation effects. Two ways to deal with this non-equilibrium relaxation is presented in the literature; the relaxation approach (HRM) and the mixture approach. The HRM applies the equilibrium properties of R744, similarly to the HEM, but introduces a time-delay depending on the flow conditions. The studies show that a single relaxation parameter is not suitable for a wide range of operating conditions. By using three different zones for the relaxation time formulation fairly reasonable results can be achieved for most motive conditions. The mixture approach treats the phase change mechanisms explicitly and uses properties in meta-stable conditions for gas and liquid. While the mixture model can produce highly accurate motive MFR results for super-critical pressures, validation at low pressures has yet to be conducted. Taking into account the phase change mechanisms in the model, such as done in the mixture model and TFM, allows for a more realistic and physical description of the flow. This approach makes the model less dependent on experimental data, but it requires more extensive sub modeling. Still, this model suffers from numerical stability issues that can cause slow computations.

Turbulence
An accurate description of the mixing phenomenon and the flow structure inside the ejector is highly dependent on accurate modeling of the turbulent flow. Due to the complexity and lack of insight into multiphase turbulence, modeling turbulence in R744 ejectors is significantly more challenging than for single-phase ejectors. Multiphase turbulence has therefore often been left disregarded in current R744 ejector models.
To model compressible multiphase turbulence, fluctuations in velocity, density, and mass fraction, as well as interface effects must be accounted for. Thus, different averaging procedures of the flow fields are considered to reduce the need for excessive sub-modeling. A common approach is to use ensemble-and phase averaging. Model complexity can then be reduced to modeling to a few additional terms in the momentum and energy equations. Most prominant is the nonlinear momentum fluctuation u u i j , typically modelled by a Boussinesq approximation, Eq. (33), where turbulent fluctuations are treated as a diffusive turbulent viscosity. Fig. 7. Relative discrepancy of motive mass flow rate plotted at the corresponding motive condition. Downward triangle -Mixture model [62], Upward triangle -HEM [73], Circle -modified HRM [61].
Here, u is a turbulent fluctuation in the velocity field, and t is the turbulent viscosity. Commonly, a set of transport equation for turbulent kinetic energy, k Eq. (34), and the turbulent energy dissipation Eq. .
where the subscript, m, indicates mixture properties based on mass or volume weighted averaging. G k m , is the strain rate production of Interfacial turbulence interactions are included in the terms k m , and m , , which are further sub-modelled [59]. Additional multiphase effects, such as bubble induced turbulence, dispersion forces, and the influence of other neglected terms have so far not been discussed in any of the current models. For further details on multiphase turbulence see Ishii and Hibiki [91], Coutier-Delgosha et al. [107], and Morel [108].
The current state of the art R744 ejector models consider pseudofluid turbulence, i.e. assuming that the turbulent structures behave similarly to single-phase turbulence. Already, a large set of turbulence models have been studied for R744 ejector modeling, as shown in Table 7. The latest approaches recommend the use of k-SST models due to its better predictions of local and global flow parameters for single-phase flows [109][110][111].
Recently, a comparison of turbulence model performance for R744 ejectors was conducted by Haida et al. [112]. In this work, the authors compared four commonly used turbulence models, namely the k-Realizable model, the k-SST model, the Transition SST model, and the Reynolds stress model. It was found that the Transition SST model and the Reynold stress model were best at predicting the local wall temperatures. Furthermore, the global entrainment ratio was found to be best predicted by the k-Realizable and the Transition SST. The authors noted a significant dependence on the near-wall turbulence formulations for the Reynolds stress model and the k-Realizable model.
The disagreements in the literature indicate that further research is needed into appropriate turbulence models for two-phase ejector flow. Two-phase turbulence effects have largely been neglected in current literature and will require detailed investigations.
A more accurate, however computationally costly, turbulence model is the Large Eddy Simulation (LES). LES resolves more of the turbulent spectrum by using a very fine mesh and only modeling the effects of the smallest turbulent length scales (Kolmogorov length scale). For further literature on LES turbulence modelling for multiphase flows and atomization see the reviews by Sher et al. [33], Jiang et al. [44], Gorokhovski and Herrmann [115], Balachandar and Eaton [116], and Fox [117]. Current work is planned on LES simulation for R744 ejectors [68]. However, this model is considered immature for multiphase flows and not well established for ejectors. Furthermore, a good understanding of the bubble size distribution, as well as the Kolmogorov length scales, is critical for the appropriate use of the LES models [44]. Further work is needed for such models to be accurately applied for two-phase ejectors.
Looking forward, with the current developments in computational speed and massive parallelization, it will likely be possible to achieve full Direct Numerical Simulation (DNS) of the two-phase ejector within the next decade. DNS simulation a turbulence approach that fully resolves all turbulent length scales in a flow. This approach has already seen use for multiphase atomization [118].

Numerical solution methods
While model equations can provide highly accurate descriptions of the two-phase flow, the correct solution of these equations needs consideration. Among other things, achieving accurate numerical solutions involves choosing appropriate meshing, solver settings, and convergence strategies. Therefore, this section treats the applied numerical Table 6 Overview of currently available R744 two-phase ejector models, their limitations and advantages. Presented in previous work [105].

Model
Advantages Challenges and limitations

Mixture
Yazdani et al. [55], -Considers meta-stability, -Increased complexity, Giacomelli et al. [62] -Can more accurately evaluate the phase -Requires tuning of model parameters, fractions by mass transfer modeling -Less profound literature database -Highly accurate results for motive flows on R744 ejectors -Not yet tested for low motive pressure Table 7 Presentation of the different turbulence models used for R744 ejector modeling and the corresponding multiphase model.

Study Turbulence model Model
Colarossi et al. [63] k HRM Yazdani et al. [55] k SST Mixture Smolka et al. [45] k RNG HEM Banasiak et al. [113] k RNG HEM Lucas et al. [57] k SST HEM Palacz et al. [56] k Realizable HEM Palacz et al. [60] k Realizable HEM & HRM Giacomelli et al. [66] k RNG HEM Haida et al. [106] k SST HRM Giacomelli et al. [114] k Realizable HEM Giacomelli et al. [67] k SST Mixture & HEM Haida et al. [112] Comparative (4 models) HRM Fang et al. [68] Comparative (2 models) HEM K.E. Ringstad, et al. Thermal Science and Engineering Progress 20 (2020) 100647 approaches in modeling of R744 ejectors using CFD. Numerical solutions to fluid dynamic problems are actively researched. For an introduction to numerical fluid mechanics and heat transfer see e.g. the books by Patankar [119] and Pletcher et al. [120]. In the available literature on R744 ejectors, the most extensively used CFD frameworks are ANSYS Fluent and OpenFOAM. This section will therefore just discuss the numerical procedures used available in these software-frameworks. Solution algorithms in CFD have traditionally been divided into density-based and pressure-based solvers. The density-based solvers were initially developed for highly compressible flows and shock solutions, where the pressure-based were intended for low-speed incompressible flows. These solvers have since been extended to be able to handle most ranges of flow conditions. However, it is generally agreed that density-based solvers have an advantage for highly compressible flows, though good results for compressible multiphase flows with pressure-based solvers have been presented in the literature [121]. The pressure-based solvers typically use the coupled approach in ANSYS Fluent, combined with the PRESTO! pressure interpolation scheme and a second-order upwind [61] or a third-order QUICK [67] scheme.
So far, the R744 ejector models available in the literature have primarily been based on segregated or coupled pressure-based solvers [45,73,62,55]. This has been necessary as density-based solvers are incompatible with multiphase solvers in the ANSYS Fluent CFD software. Recently, an alternative numerical approach was considered by Fang et al. [68] for R744 converging-diverging nozzles. Their approach uses the rhoCentralFoam density solver which is based on the KT (Kurganov and Tadmor [122]) and KNP (Kurganov et al. [123]) secondorder spatial scheme. Boundary conditions based on in-and outgoing characteristics was applied for improved stability of the solution. Higher-order schemes for investigations of flashing R744 was conducted in the study by Gjennestad et al. [124]. Here, a higher-order CFD method based on the third-order Weighted Essentially Non-Oscillatory method and Span-Wagner EOS at equilibrium was presented.
The transient solution is often neglected in simulations for R744 ejectors. The assumption of steady-state R744 ejector flow has been considered in multiple studies [45,73,61,62]. This assumption has still not been fully investigated for R744 ejectors and further work is needed. Transient effects may be significant for other ejector geometries and CFD models and should be verified for each model case.

Model dimensionality
The simplicity of the ejector geometry makes a low dimensional model an attractive option. Ejector flow is commonly considered a 1D problem, as it primarily involves one flow direction. The one-dimensional models are typically divided into true 1D models and pseudo-1D models (also referred to as 0D models or thermodynamic models). The pseudo-1D models use empirical coefficients for the different sections of the ejector (nozzle, suction chamber, mixture chamber, diffuser) and connect these to predict performance. Such models have been extensively proposed in the available literature [24,125,126,18]. Among the first two-phase ejector models is that of Kornhauser [32], which has been extensively utilized and developed in later works [127,31]. Recently, Taslimi Taleghani et al. [54] presented a thermodynamic R744 ejector model able to reproduce both single choking and double choking conditions. Taslimitaleghani et al. [22] also used such a model to study the benefits of R744 ejectors in different cycle configurations. The authors concluded that the greatest benefits of ejectors were gained when used for expansion work recovery, noting an improvement of up to 23% in COP compared to throttling. The benefit of these simplified models is that design optimization comes at a much lower computational cost.
On the other hand, true 1D models treat the flow as two interacting 1D streams, one representing the motive flow and one representing the suction flow. Banasiak and Hafner [83] combined such an approach with multiphase models to account for meta-stable effects.
However, 1D models offer limited insight into the local flow field inside the ejector and are derived from assumptions that limit their range of validity. 3D models may have a significant effect on model fidelity as noted by Smolka et al. [128]. This is especially true for nonaxisymmetric suction chambers and high swirl flows. Mazzelli et al. [109] supported this by investigating 2D and 3D CFD models for ejectors, showing that for off-design conditions 2D models would fail to yield acceptable results. While full three-dimensional simulations of ejectors are fully feasible with CFD, the computational cost can be drastically reduced by assuming 2D flow. Ejector CFD models have extensively used a two-dimensional flow approach. The 2D flow assumption can give results with reasonable accuracy [67,62,57,68], which makes it suitable for model development as well as optimization algorithms where computational costs are limiting.

Experimental validation data
Experimental investigations on R744 ejector systems have been extensively reviewed in previous work [37,17]. In this section, some recent developments and future opportunities for model validation will be explored. Furthermore, a compiled data set used in some previous works for validation purposes will be presented. Extensive literature is available for ejectors with other working fluid. However, these are not applicable for validation purposes as R744 ejectors are significantly different from ejectors with other fluids. Specifically, in terms of ejector size, operating conditions (supercritical motive conditions), refrigerant properties, multiphase effects, and non-equilibrium flow, R744 ejectors deviate from other working fluids.
Ideally, to properly validate a model a large set of accurate experimental data for different flow variables, such as velocity profile, turbulence intensity, vapor quality, and pressure distribution, is required. Nonetheless, the typical R744 ejector dimensions are of magnitude millimeters, which heavily restrict the available experimental measurement techniques that can unobtrusively be applied. As an example, a common way to measure local velocity, even for supersonic flows, is using hotwire anemometers [129]. However, these devices typically still have a wire length of 1.2 mm [129], where the size of an R744 ejector mixing section can be 3 mm [113].
To date, the most common validation parameter is the mass flow rates (MFR) and the entrainment ratio. These have been extensively measured for different R744 ejectors at different operating conditions [73,65,18]. A compiled data set of some of the investigations conducted at SINTEF Energy Research (Trondheim, Norway) [60][61][62] that has been used for R744 ejector models in previous studies is presented in Table Appendix A for ease of use in future model validation and comparison to other models. The operating conditions of the different investigations are illustrated in a P-h diagram in Fig. 8 and include both near-critical and off-critical flows.
Even though the prediction of global parameters is an important validation parameter, proper validation of model results should be assessed based on local parameters. Local parameters give insight into the physical realism of the applied model, which is critical for reliable ejector models at large ranges of operating conditions and ejector designs. Some literature has been presented that considers local parameters are discussed below and summarized in Table 8.
An experimental study that has been extensively used for model validation is that of Nakagawa et al. [50]. A converging-diverging nozzle, similar to those used in ejectors, was equipped with pressure sensors. While this can be used to validate the expansion process in the nozzle, the full ejector was not considered. A visualization technique was used in the study by Berana et al. [51]. Here, the converging-diverging nozzle with transcritical R744 was investigated using direct photography and pressure measurements. This work followed up on the study by Nakagawa et al. [50] and further discussed the observed shock waves.
Haida et al. [61] presented experimentally obtained wall temperature data for an R744 ejector. Temperature probes were introduced into the ejector outer walls. However, validating a model using wall temperature data requires that the model includes wall heat transfer effects, adding additional model complexity.
Recently, visualization techniques have been used for imaging of the flow inside an R744 ejector. Such investigations were conducted by Zhu et al. [130] and Li et al. [52]. Using a direct photography visualization technique, the images identify the expansion angle of the primary flow. Zhu et al. [130] analyzed the different operating conditions and identified the expansion angle at different primary and secondary pressures. An example-image is included in Fig. 9. Furthermore, assuming an isentropic expansion with an irreversibility efficiency a liquid fraction in the suction chamber was estimated and compared with the observed grayscale in the images. The experiments showed that low primary flow pressure increased the expansion angle. The larger flow area caused by the increased expansion angle blocks the suction flow, limiting the entrainment.
Li et al. [52]131 used a high-speed camera as well as pressure sensors to investigate the phase change position for different operating conditions in an R744 ejector. Their findings indicate that the location of the phase change is largely dependent on the motive pressure. Furthermore, the obtained grayscale values were analyzed, and, under the assumption of a set bubble diameter, a qualitative bubble number density distribution was presented, shown in Fig. 10.
Looking forward, certain techniques not yet applied to two-phase ejectors may give more detailed insights into the ejector flow. As an example, Particle Image Velocimetry (PIV) is a velocity measurement technique that involves seeding particles into the flow and tracking these. The PIV method has already been used in single-phase- [133,134] and multiphase ejectors [135] with good results. Such an approach is, however, more challenging in the case of R744 ejectors due to their small size and high pressure. Currently, Haida et al. [132] are performing the first particle image velocimetry (PIV) measurements for a circular R744 ejector. Such an investigation could provide detailed insights into the local velocity field within a circular R744 ejector.

Generation of entropy
One primary interest in ejector modeling is the understanding of the origins of irreversibilities. The ejector is a theoretically isentropic expansion device, however, due to irreversible viscous forces, supersonic shocks, and turbulence, isentropic expansion is physically unobtainable. It is therefore of interest to examine how entropy is generated in ejector flows. Such an investigation was conducted by Banasiak et al. [113], where the CFD model of Smolka et al. [45] was used to identify zones of entropy generation in the R744 ejector. The results indicate that the largest production was located in the mixing section, primarily due to turbulent and shock losses. The authors argue that geometry    optimization should be conducted with respect to the interdependency of the geometric parameters, i.e. that optimization of one parameter in isolation does not in general yield optimal results. Further investigation of entropy can be found in the study by Sierra-Pallares et al. [136]. In this work entropy generation was evaluated for an R134a ejector by direct entropy method, solving a differential equation for entropy transport. This investigation supported that the primary contributor to irreversibilities in an ejector is turbulent viscous dissipation produced near the nozzle exit, as claimed by Banasiak et al. [113].

Shape design and optimization
Model-based design approaches have been extensively been used for design improvement for R744 ejectors. The studies based on modern CFD approaches are presented here. The HEM presented by Smolka et al. [45] has been extensively applied for R744 ejector shape design. One such study was conducted by Palacz et al. [60], where an R744 ejector was optimized using the EjectorPL software. This software was designed to automate parts of the CFD workflow, such as geometry generation and meshing, allowing for the rapid development of a large set of CFD results. This set can then be fed into an optimization algorithm to search for optimized geometry designs. In the work by Palacz et al. [60], six geometric dimensions were investigated to optimize for the maximal ejector efficiency. This was done using both an evolutionary and a genetic optimization algorithm. Later, Palacz et al. [75] furthered this study by investigating additional ejector geometries and geometric parameters, achieving up to 6% (percentage points) increase in ejector efficiency. The findings of these studies show that CFD models can be used to improve ejector efficiency robustly and reliably.
An interesting prospect for model design optimization is the implementation of adjoint methods [137,138]. Adjoint methods attempt to optimize an objective by looking at small changes in the design parameters. Adjoint methods then solve simultaneously the flow solution and the solution with different design parameters. However, adjoint methods are not sufficiently developed for complex multiphase flow, such as two-phase ejector flows.

Ejector flow control
An ejector can be used to control the high-side pressure in a refrigeration system and to regulate to achieve the optimal cooling capacity. The capacity control strategy of two-phase ejectors, i.e. the control of ejector pressure lift and entrainment, plays a pivotal role in guaranteeing high overall energy efficiency [20]. Therefore, different concepts have been proposed for flow control in ejectors. The methods which have gained to the most attention are the multi-ejector block [139], the adjustable needle design [31,18], and motive swirl control [140], see Fig. 11. CFD models have been a useful tool for investigating these solutions.
A comparative study between the fixed geometry ejector, used in a multi-ejector block solution, and the controllable needle ejector was conducted by Smolka et al. [128]. In this work, the two designs were modeled with the HEM CFD model and compared. The fixed geometry ejector showed good performance for all investigated operating conditions. The adjustable geometry design was able to outperform the fixed geometry design. The drawback of such an approach was that identifying the needle positions that yielded the best efficiency at a given flow rate were found difficult to asses, making the efficiency using such an ejector less predictable.
A numerical study of a full-scale multi-ejector module using the HEM approach was presented in Bodys et al. [141]. Their results indicated the high importance of properly designing outlet collectors. From this study, the authors considered the regulation capabilities of the multi-ejector block suitable for industrial refrigeration. Furthermore, Bodys et al. [142] looked into the swirl motion generation for multi-ejector applications. The study suggests that at high motive rotational speeds increased mass entrainment ratios could be achieved, independent of ejector size.
Bodys et al. [143] performed a numerical investigation on an R744 ejector with a bypass inlet [144], see Fig. 12. Different configurations of bypass placement and inlet angles were analyzed using the HEM approach. According to this study, up to 39% increase in mass entrainment ratio could be gained by using the suction flow bypass.

Model reduction
Besagni [42] emphasized the need for large reductions in model size to achieve coupling between the local scale modeling (such as CFD) to the large scale system simulations. The use of CFD models to approximate reasonable component scale Lumped Parameter Models (LPM) was indicated as an attractive prospect. However, large scale production of high accuracy CFD results is in many cases unfeasible due to the large computational and time costs. Nonetheless, advanced methods to achieve CFD-level fidelity at much lower computational costs have recently been presented in the literature. One such approach is the reduced-order models (ROM). ROM is a modeling approach where a set of results are processed to establish a reduced basis which optimally represents the previous results. Typically, high accuracy CFD and experimental results at a wide operating range are used as the input.
Such a ROM approach was applied for R744 ejector modeling by Haida et al. [106,145]. In these studies, the HEM model was used to create a set of results (referred to as snapshots) for the optimization of the reduced model. These snapshots were used to generate ejector performance maps [106] (such as those experimentally obtained by Banasiak et al. [20]) and develop accurate system models [145] based on the CFD model results. The model was implemented into a dynamic system scale simulation for different climate zones in Haida et al. [132]. This ROM was compared with simpler 0-and 1D models, commonly used in the literature. It was shown that the ROM predicted the ejector performance with superior accuracy to the thermodynamic models and could efficiently be implemented into system scale simulations. The accuracy of a ROM is however, limited by the accuracy of the analyzed input, and can therefore not be used to explore designs beyond the given input data. Therefore, more accurate CFD models will improve the ROM accuracy, which motivates further research into CFD modeling.

Conclusions and future developments
This paper has presented a complete overview of the current state of the art within two-phase R744 ejector modeling. The adoption of twophase ejector technology presents a substantial improvement in the energy efficiency of R744 vapor-compression systems, potentially leading the HVAC&R sector to have a future-proof solution for many applications. However, a rigorous design of the R744 two-phase ejector is necessary to achieve this target. As a consequence, the implementation of advanced models and tools for an accurate ejector design has become one of today's most important research topics. The implementation of computer-assisted tools, such as CFD, can improve R744 two-phase ejector design and operation. However, to the best of the authors' knowledge, no detailed review works have been carried out yet despite the importance of the research.
Therefore, the purpose of this investigation has been to bridge this knowledge gap by comprehensively evaluating all the aspects as well as challenges related to R744 two-phase ejector models, i.e. multiphase flow modeling, turbulence aspects, numerical solution methods, applications of models. This work has found that significant discrepancies can be observed when comparing ejector-model results with experimental measurements. The models are trending toward more advanced two-phase models that can treat motive flows at a wide range of operating conditions. The focus of recent research has been on the treatment of non-equilibrium expansion that occurs at low-pressure motive conditions. This review has shown that while the currently most advanced models can achieve reasonable accuracy at supercritical motive pressures, the accuracy quickly drops at lower pressures. Furthermore, current models still struggle to predict suction mass flow rates, likely due to its dependency on other flow features as well as models for turbulent momentum diffusion. Increasing model complexity does, however, require detailed knowledge about the ejector flow conditions, such as local pressures, velocity profile, turbulence intensity, and vapor quality. Currently, this knowledge is very limited due to the challenge of achieving accurate experimental measurements in the R744 ejector.
Future work should explore many aspects of ejector models. Special interest should be devoted to better the understanding of basic twophase effects, ie. meta-stability, velocity slip, atomization of jets, twophase turbulence, and multiphase thermodynamic properties. Additionally, exploration of the effects of various mixtures like oil-CO 2 in ejectors could uncover significant design improvements. More advanced modeling, such as two-fluid models, could enable better model accuracy while being less dependent on experimental tuning. Further development of thermodynamic and 1D models can be highly useful for optimization procedures and system modeling. Accounting for different non-equilibrium effects has been shown to improve accuracy in CFD models, and should be investigated further for thermodynamic and 1D models. 1D and thermodynamic models can also be supported by reduced-order models based on experiments and CFD modeling. Improving post-processing tools, such as exergy tubes [146], may allow for quantification of losses in ejectors and improve optimization tools. Artificial neural networks could be introduced to improve flow models by optimization of flow parameters [147], improving thermodynamic libraries [96], or predicting flow characteristics [148]. Improved design tools based on these models could provide an efficient way to explore and optimize ejector geometries at a limited computational cost. Modelling the near-wall flow is important for capturing important turbulent structures as well as flow separation. Such effects have been extensively discussed for single-phase ejectors [109,[149][150][151][152][153]152,151], however further work should be conducted for multiphase R744 ejectors. Achieving high accuracy local experimental measurements should also be highly prioritized as it can enable validation and tuning of future advanced models. Further experimental work using PIV could give insights into the local velocity-and turbulence fields. Additionally, further research using tools such as Phase Doppler Interferometry or high-resolution Schlieren imaging could help in understanding bubble size distributions and temperatures. Adapting such techniques may aid in achieving local experimental validation data, which will greatly improve future R744 ejector models.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.