Application of a fast transonic trajectory determination approach in 1-D modelling of steady-state two-phase carbon dioxide ﬂow

. An original generalised procedure of determination of the transonic trajectory has been proposed. The procedure is much faster than the commonly used Newton Critical Point approach. The approach was applied in modelling of a carbon dioxide transonic two-phase ﬂow through the converging-diverging nozzle by means of the Homogeneous Equilibrium Model and Delayed Equilibrium Model (DEM). The simulations concern ﬂows that were experimentally and theoretically investigated in the literature. DEM was previously used only in choked water ﬂow simulations. Its application in CO 2 ﬂow modelling and the supersonic trajectory part determination is a novel contribution. The adjusted for CO 2 version of the closure law was proposed. The investigation revealed that the applied Darcy friction factor determination approach has a signiﬁcant inﬂuence on the results. Moreover, the models are unable of producing physically acceptable solutions until the Lockhart-Martinelli approach is utilised. It was shown that the Friedel approach might be considered more proper for CO 2 ﬂows.


Introduction
Studies on implementation of the ejectors with converging-diverging nozzles in refrigeration, airconditioning and heat pumps demonstrated significant improvement of the energy efficiency of the discussed systems [1][2][3]. As a result, a growing interest in studying an operation and design of the ejector in mentioned systems can be observed [4][5][6]. The accurate prediction of a critical mass flow rate of the flow through the ejector motive nozzle is crucial for an appropriate design of the ejector [7]. The recent investigations of the CO 2 two-phase ejectors in refrigeration systems [4,8], heat pumps [6] and in other industry domains [9], provide experimental data in this field which gives an opportunity to study the prediction abilities of the available two-phase flow models. The schematic of a two-phase ejector with supercritical gas or subcooled liquid phase as a motive fluid and vapour phase as a secondary fluid is presented in Fig. 1a. The liquid phase due to depressurisation achieves saturation conditions and partly evaporates inside the nozzle so that two-phase flow emanates from the nozzle. Due to the momentum transfer between a motive jet and vapour, the vapour is entrained into the suction chamber and then flows into the mixing chamber. Further, as an e↵ect of the momentum transfer mixing process occurs along with formation of a two-phase shock wave with a rise of static pressure. Then the mixed homogeneous two-phase flow is additionally compressed in the diffuser achieving the discharge pressure pm. The partial evaporation due to depressurisation is called flashing and it is presented in Fig. 1b. The flashing process is ⇤ e-mail: w.angielczyk@pb.edu.pl characterised by high thermal none-equilibrium (e.g. temperature di↵erence across the interface). Despite that at a certain location, the equilibrium saturation pressure is achieved (the left vertical dashed line in Fig. 1), the nucleation process does not start and a metastable liquid flow occurs up to a location where the metastable liquid is sufficiently superheated to sustain the vapour nucleation (the right vertical dashed line) [10]. From that location, the two-phase flashing flow is formed. In this case, the growth rate of the vapour phase is limited by the interphase heat transfer rate rather instead of mechanical expansion. However, under real operation conditions, both mechanical and thermal non-equilibrium exist simultaneously making strong difficulties to accurately predict critical mass flow rate through the motive nozzle. This problem was widely investigated [7,[11][12][13] and it was proved that most physically consistent predictions may be obtained from models that take into account the thermal none-equilibrium e↵ects. Two of such kind models (HRM and DEM) together with Homogeneous Equilibrium Model and lumped parameter approach of Moody and Hanry-Fauske were investigated in [14]. It was shown that among tested models, DEM is the most accurate in terms of both the critical mass flow rate predictions and pressure distribution predictions. Therefore, the main aim of the present investigation was to check if the DEM with a closure law, originally developed for water, may be considered suitable for modelling of CO 2 transonic two-phase flows. HEM imposes the thermal and mechanical equilibrium between vapour phase and liquid phase (temperatures, pressures, and velocities of both phases are equal, the phases are perfectly mixed). The model consists of three conservation equations of mass, momentum, and energy, and two state equations, respectively: In the implementation of HEM used in this investigation, the independent variables are: the pressure p, the quality (the saturated gas mass fraction) x, and the velocity w.
Since the flow is assumed to be adiabatic then the heat flux q" is equal to 0. HEM is used as a reference case in this investigation.

Delayed Equilibrium Model (DEM)
DEM assumes the existence of three fractions: the metastable liquid phase (subscript ml), the saturated liquid phase, and the saturated vapour phase. The metastable fraction is assumed to have the same pressure as the saturated phases but higher temperature, as it is assumed to proceed an isentropic expansion. Therefore, DEM takes into account the thermal non-equilibrium e↵ects but it does not include the mechanical non-equilibrium e↵ects. DEM consists of the conservation equations (1-3) substituted with the following state equations: The vaporisation index y is the mass fraction of the saturated phases: In order to complete the model, a mass balance equation for saturated fractions has to be added to the system of equations (1-3, 6, 7). Recently, the closure equation for water has been presented in a new form [15]: The implementation of DEM used in this investigation treats the flow as adiabatic (q"=0). The independent variables are: p, x, y, w.

Solution procedure
Determination of the transonic solution of HEM or DEM requires a profound topological analysis. This analysis uses the theory of dynamical systems and was conducted by Bilicki et al. [12]. This section presents only crucial elements of the analysis which are necessary to describe the basis of the applied solution procedure. The proposed methods require some information that can be obtained by applying the Possible-Impossible Flow algorithm that was presented in [11,13].

Topological structure of the phase space
Practically all known one-dimensional models of a steadystate flow (including HEM and DEM) can be presented in a form of the following nonlinear ordinary first order differential equation system (Einstein summation convention has been applied): The size and elements of matrix A and vector b depend on the model type. The vector σ consists of n quantities describing the thermodynamic state of the fluid. The elements of matrix A depend only on vector σ components, and vector b elements additionally depend on the spatial coordinate z. The set of governing equations (10) supplied with vector σ B = [σ 1,B , σ 2,B , . . . , σ n,B ] (related to the flow inlet conditions) creates an initial-value problem. A solution to the problem is a trajectory σ(z) in n+1 dimensional phase space ⌦. Vector [dz, dσ 1 , dσ 2 , . . . , dσ n ] is tangent to σ(z) at each point. Therefore, vector V defined as: is also tangent to σ(z). Here D denotes a determinant of the A matrix, and N i are determinants of matrices that are created by replacing the i-th column of A with b. It is worth to notice that equation (11) is an autonomous form of the system (10). Fig. 2 presents a projection of trajectories corresponding to flows through a convergingdiverging nozzle on pressure p -spatial coordinate z plane (σ 1 = p). All trajectories related to the curves in Fig. 2 start form the same values of the inlet pressure p B , density ⇢ B , and specific enthalpy h B . The trajectory that corresponds to a fully subsonic flow (e.g. BL, BN, BP) or impossible flow (e.g. BF, BG, BH) can be determined through numerical forward-marching integration of (11) from the inlet conditions. The trajectory related to the curve passing through points B, S , E 1 is the sought transonic solution. Point S at which D = 0 and all N i = 0 is a singular point. The direction of V at this point is undefined because its magnitude is equal to 0. Thus, the numerical integration algorithms cannot neither "start from" nor "pass through" this kind of points. Therefore, another approach is required that is di↵erent than forward-marching integration. The conventional solution to this problem is based on the fact that when S is a singular saddle point then the eigenvectors of the Jacobian of V are tangent to the trajectories that pass through S (the Jacobian has to be calculated at the point S ) [12]. Let us denote the eigenvector tangent to the transonic trajectory by V 1 (Fig. 2). Making a step along V 1 determines a point that is eligible for starting numerical integration. The up-stream integration when it reaches z B , determines the inlet conditions corresponding to the considered singular saddle point. Similarly, down-stream integration determines the outlet conditions. Therefore, the described procedure is actually a solution to an initial-value problem of equations system (11) with initial values: z = z S and σ = σ S . In this conventional approach, the singular saddle point is chosen arbitrarily. Thus, in general case, the corresponding inlet conditions significantly di↵er from the given values. In order to fit into the given inlet conditions, a shooting method can be applied, e.g. method described in [16] that starts from a singular saddle point. A multidimensional, globally convergent Newton-Raphson algorithm can be applied for this purpose. In this case, the approach is called Newton Critical Point method (NCP) [17]. Each NCP iteration requires n integrations to determine how to change n − 1 singular saddle point parameters [16]. The enormous overall integration number makes this method relatively slow (especially if compared to Possible-Impossible Flow algorithm). The described issues and the laboriousness of solving the V Jacobian eigenvalue problem was a motivation for developing a faster approach.

Proposed generalised method of singular saddle point determination
The proposed method utilises the fact that the saddle singular point is an intersection point of the following curves: D(σ) = 0, N i (z, σ) = 0 ( Fig. 2 shows curve D = 0 and N p = 0). Thus, determination of the saddle singular point approximation requires the following steps (the PIF procedure has to be conducted previously): 1. Approximate the curve D(σ) = 0 by the straight line passing through point G and H.

Determine a point S i (z i , σ i ) that lies on the line. At this point z i = z H + ∆z. Assume thatṀ
3. Use the gradient descent method [16] (take S i as a starting point) to determine S i,0 (z i , σ i,0 ) where D = 0.
4. If all Ns at the point S i,0 have di↵erent signs than in point H, then repeat the previous steps with smaller ∆z until the required accuracy is reached. Otherwise, repeat the previous steps but replace H with S i,0 .
The proposed method is generalised in the sense that it may be applied for any model described by (10) regardless of the number of the equation. Fig. 2 shows the points L, N, P that are inflection points of the subsonic trajectories. Each of the vectors V L , V N , V P is tangent to a trajectory at a corresponding inflection point. The points M, O, R represent local pressure minima. It is worth to notice that the closer to the point S the subsonic trajectory lies, the closer the inflection point of this trajectory is to the pressure minimum point of this trajectory. Finally, at the transonic trajectory, those points merge together into point S . Consequently, the direction of V 1 can be approximated by the direction of V P . The direction of V P can be directly used to carry out the integration from the found saddle singular point (as it was described at the end of Subsect. 2.3.1) resulting with an approximation of the subsonic (curve BP) and supersonic (curve P 0 E1 0 ) parts of the sought transonic trajectory.

Results and discussion
DEM results compared with HEM predictions: The DEM correlation constants presented in eq. (9) were adjusted on the base of the Super Moby Dick experiment simulations described for instance in [13]. In the simulations, the model has been supplemented with Darcy friction factor calculated on the base of the Lockhart-Martinelli (LM) parameter and the flow was treated as adiabatic [13]. As a first step of the present investigation, calculations of all 10 runs of CO 2 flow cases presented in [8] were carried out with described DEM original setup. In addition, calculations with the use of HEM were carried out in order to obtain the reference data that are related to an instantaneous mass transfer between the considered fluid fractions. The calculated flow cases related to the described original DEM setup but conducted for CO 2 revealed the same tendency that here (for sake of the paper size limitation) is presented on the basis of only three flow cases. Pressure profiles for those representative cases are plotted in Fig. 3 -5 with red dashed lines and denoted by DEM. In the first of the representative cases, shown in Fig. 3 the nozzle outlet area (and consequently the angle of the nozzle divergent part ) is the highest one. Thus, a strong decompression was expected to occur in the nozzle divergent part. However, DEM generates a pressure profile located above the experimental pressure distribution despite the corresponding mass flow rate DEM is much higher than those calculated by means of any of HEM implementations (Fig.  3). The critical section, in this case, is located quite away from the throat and contains a singular saddle point. The shift of the critical section downstream from the throat is caused by a presence of thermal none-equilibrium and friction terms in the equations system (1-3, 6-9). This e↵ect was described, among others, by Bilicki et al. [12]. The second representative case is presented in Fig. 4. Again, the generated DEM pressure profile is located high above the experimental pressure distribution, despite the corresponding mass flow rate DEM is much higher than those calculated by means of any of HEM implementations. The singular saddle point is located almost in the middle of the nozzle divergent part. Consequently, the related trajectory is fully convex from this point until the nozzle outlet [12]. The last representative case is presented in Fig.  5. The generated DEM pressure distribution is still above the experimental data. The corresponding mass flow rate is still the greatest but the critical section, in this case, is located at the nozzle outlet and contains a turning point. The provided results reveal that the original DEM setup is not suitable to correctly predict CO 2 pressure profiles of Nakagawa et al. [8] experiment. The calculated profiles are passing above the experimental points. Thus, it seems that the dy/dz calculated by the original form of eq. (9) is too low for the case of CO 2 flows. An increase of dy/dz requires higher values of C 1 , C 2 in eq. (9). The results (presented in Fig. 3-5 by black solid lines and denoted by DEM-CO 2 -LM) confirmed the conclusion concerning dy/dz, since the curves calculated for the new set of the constants (which is: C 1 = 5.17, C 2 = 0.87 and C 3 = 0.25) are passing closer to the experimental points. However, the further increase of the constants causes that the calculated curves approach the HEM-LM (HEM using LM approach) pressure profiles (thus, they recede from the experimental points). This situation is caused by a fact that when DEM operates with significantly high dy/dz values, only few integration steps are required to fully convert the metastable fraction into saturated fractions. Some unexpected conclusions may be drawn from an analysis of HEM-LM curves (Fig. 3-5) black doted-dashed lines). Namely, the curves pass above the saturated pressure points (black dots) that are calculated for measured temperatures. However, this situation is not physically accurate. Namely, the measured temperature can be at most equal to T ml and at least equal to T l . Thus, the physically realistic HEM profiles always pass under the black dots. Only in such a case, an introduction of the non-equilibrium e↵ects could produce a temperature profile that realistically corresponds to the measured temperatures. Fig. 3-5 show pressure profiles calculated by means of two additional HEM implementations. Namely, HEM with ⌧ w = 0 (DEM-isentropic) and HEM using the Darcy friction fac-  tor calculated by means of the Friedel approach (DEM-Friedel). In all cases (also in those that are not depicted in the paper), the pressure distribution calculated with the LM approach passes above the saturated pressure points (black dots) that are calculated for measured temperature values. Consequently, the LM method commonly used for water flows, may be considered improper for CO 2 flows. It is worth to note that the corresponding mass flow rate is the lowest one in all of the cases. In most calculated cases, the HEM-Friedel curves pass under the experimental saturation pressure distributions. The corresponding mass flow rates are always between those related to HEM-isentropic and HEM-LM. The HEM-isentropic curves always repre-sent the steepest pressure drop and the highest mass flow rate. They always pass under the experimental distribution of saturation pressures (the black dots) but not always under points of the experimental pressure distributions. In all of the HEM-isentropic cases, the critical sections are located slightly beyond the throat.

Conclusions
The presented investigations revealed that the original DEM setup (the closure law developed for water, the friction factor calculated by means of the LM approach) is not suitable for CO 2 transonic flow calculations. More- over, the attempt of adjusting the closure law constants revealed that an increase in dy/dz causes an improvement of the pressure profiles predictions. However, this enhancement is restricted as a consequence of the LM approach application. The further investigation showed that the LM method is improper for CO 2 flows as it produces saturation pressure distributions that are passing above the saturated pressure points calculated for measured temperature values. In the case of Friedel approach, this happened only in three out of eight available experimental cases of two-phase flow. Thus, this approach can be considered more appropriate for CO 2 flows. However, among the investigated HEM implementations, only HEM-isentropic produces saturation pressure distributions that are always passing under the saturated pressure points. Therefore, DEM with ⌧ w = 0 is the most reasonable implementation for further adjusting attempts of eq. (9). This approach is planned to be used in the future works unless a more proper Darcy friction factor determination approach is found or developed. The conducted calculations confirmed that the proposed approach of transonic trajectory determination is reliable and practically useful. The combination of the PIF algorithm as a first step and the proposed approach (or the NCP method) as the second step is, to the authors' knowledge, the only available universal solution method. Moreover, the utilization of the proposed approach, instead of the NCP method, significantly shortens the computation time.