Internati onal seminar on ejector/jetpump technol ogy and applicati ons Louvain-la-Neuve, Belgium September 7-9, 2009 Paper No. 34 MULTIDIMENSIONAL MODELING OF CONDENSING TWO-PHASE EJECTOR FLOW Davi d Schmi dt and Michael Colarossi University of Massachusetts Amherst Mark J. Bergander Magnetic Development Inc. ABSTRACT Condensing ejectors utilize the beneficial thermodynamics of condensation to produce an exit ing static pressure that can be in excess of either entering static pressure. The phas e change process is driven by both mixing and interphase heat transfer. The co mp lete details of the flu id flow are not well understood and cannot be calculated directly. However, semi emp irical models could be used in conjunction with computational flu id dynamics (CFD) to gain some understanding of how condensing ejectors should be designed and operated. The current work describes the construction of a preliminary mu ltid imensional simu lation capability that is built around an Eulerian pseudo -fluid approach. The transport equations for mass and mo mentum t reat the two phases as a continuus mixture. The rate of mass transfer is assumed to be limited by heat transfer, and a modified form of the Ho mogenous Relaxation Model (HRM) is employed. This model was originally intended for representing flash-boiling, but the authors hypothesize that with suitable modification, the same ideas could be used for condensing flow. The co mputational flu id dynamics code is constructed using the open -source OpenFOAM lib rary, wh ich facilitates object-oriented fluid flow calculations. The resulting code is a transient, finite-volu me code that can operate in either two or three dimensions, though two -dimensional calculat ions appear to be much more nu merically stable. The simulat ions are parallelized using a standard message passing interface (MPI) co mbined with domain decomposition. The simulat ion can use general polyhedral cells, for flexibility in geo metry. Fluid properties are evaluated using the REFPROP database from NIST , which includes many common fluids and refrigerants. Development is currently underway for handling multi-co mponent fluids. The results of the simulat ions will allow the study of the flow within the two-phase condensing ejectors. Of particular interest are the stability of the condensation process and the permissible boundary conditions. The model can also be used to study flow separation and consequent diffuser stall that may occur in the diverging portion of the ejector. Due to the strong adverse pressure gradients that coincide with condensation, extreme care is warranted in avoiding such recirculations. NOMENCLATURE Θ Timescale (s) ρ Density (kg/ m3 ) h Enthalpy (J/ kg) p Pressure (Pa) τ Stress tensor (Pa) p crit Critical pressure (Pa) φ Mass flu x (kg/ m2 *s) p sat Saturated pressure (Pa) ψ Dimensionless pressure difference t Time (s) U Velocity (m/s) x Vapor mass fraction Equilibriu m vapor mass fraction α Vapor volu me fraction INTRODUCTION The main objective of this project is to model the flow inside a condensing ejector, including the condensation shock that occurs, when the vapor phase is quickly condensed onto the liquid stream, producing rapid transformation fro m two-phase into single-phase flow with a resulting rise in pressure. Figure 1. Princip le of a condensing ejector The condensing ejector is a two-phase jet device in which a sub-cooled fluid in a liquid state is mixed with its vapor phase, producing a liquid stream with a pressure that is higher than the pressure of either of the two inlet streams. Figure 1 above presents the predicted behavior of a condensing ejector. The high relative velocity between vapor and liquid streams produces a high value of heat transfer and the rate of condensation of the vapor phase is high. This condensed vapor adds to the mo mentum flu x of liquid stream. The remaining vapor and liquid jet (original plus condensed vapor) enter the constant area mixing section where a condensation shock may occur such that a completely liquid state exists downstream of the shock. The pressure rises due to lowering the velocity of the stream and the mass flow balance is maintained by a sudden change of density (condensation). A computational fluid dynamics (CFD) analysis was emp loyed to overcome the limitat ions of a conventional control volume analysis and to provide an input to design methodology of the condensing ejector. The simu lation method is needed to predict the behavior of the condensing ejector in various operating conditions and for different refrigerants. This provides the basis for establishing fundamentally new engineering and design methods for a broad spectrum of equip ment operating on two-phase flow, not only condensing ejectors but also pumps, heat exchangers, jet compressors, supersonic nozzles and diffusers. MODELING APPROACH For a condensing flow simulat ion, the basic laws of conservation are followed as for a single fluid. This is known as the “pseudo-fluid” approach. The equations for the conservation of mass, mo mentum, and energy are given below. In these equations,  is the mass flu x and  is the stress tensor.     0 t  U  (U )  p   t   h  U  p  (h)  t t (1) (2) (3) These three equations are not a closed system. For a flu id in single -phase, an equation of state would be necessary. In the case of two-phase flow, the condensation process involves mixing and heat transfer between the phases, and there is no direct method to calculate this. Also, the two-phase mixture is not in thermodynamic equilibriu m. To bring the mixture towards equilibriu m, a modified form o f the Ho mogenous Relaxation Model (HRM) is used. The model was originally used for flash-boiling cases, but with some modificat ions the same idea can be applied for condensing flow. The Homogeneous Relaxation Model is based on a linearized expansion proposed by Bilicki and Kestin [1]. The equation for the total derivative of quality, the mass fraction of vapor, is shown below. The derivative shows the relaxation of the quality, x, co mpared to the equilibriu m quality, , over the timescale Θ. Dx x  x  Dt  (4) The formula suggested by Downar-Zapolski et al. [2] for the timescale, Θ, is shown below in Eqn. 5.    o a b (5) The values for Θ0 , a, and b were determined by best-fit values to flashing experiments [2]. The values are Θ0 = 6.51 10-4 [s], a = -0.257, and b = -2.24. The variable α is the vapor volume fraction and the variable  is a dimensionless pressure difference represented in Eqn. 6. Though the current application is far beyond the intended use of this empiricism, it serves as a starting point for further investigation.   p sat  p p crit  p sat x ,h Dp   Dt x p ,h Dx   Dt h p,x Dh Dt (7) As shown in the above equation, density is a function of pressure, quality, and enthalpy [1]. The first and last terms on the right-hand side of Eqn. 7 are currently not included. The HRM model is an empiricis m that has been tuned to represent flash boiling. The current application to rapid condensation represents a stretch of the model beyond its original intent, and quite possibly, beyond its validity. One should certainly expect that some adjustment of the model is required. In the current application, the following changes were made to the HRM model:  A turbulent mixing model based on a single length scale was added to the transport equations for mo mentum, energy, and density (Eqn. 7) [3]  An extra term was added to Eqn. 5 to provide symmetry of phases. According to Eqn. 5, the timescale of phase change becomes infinite at zero void fraction. This behavior represents the effect of nucleation. A similar behavior should be present for the case of unity void fraction.  The timescale was bounded, both above and below, in order to avoid numerical instability and floating point exceptions. By subtracting the conservation of mass, Eqn. 1, fro m Eqn. 7, the left side gives an expression for the velocity divergence at the new time step.    U   p x ,h Dp   Dt x a pU p  H U   p (6) The HRM model solves for the pressure that satisfies the chain rule and utilizes the continuity equation. Because a non-equilibriu m fluid is being considered, the chain rule takes the place of an equation of state. In the chain rule, the pressure reacts to the change in compressibility and density caused by the phase change. Eqn. 7 describes the chain rule as a total derivative of density. D   Dt p To make an equation for pressure, a discretized form of the momentu m equation is coupled with the chain rule. A generalized form of the mo mentum equation is shown below in Eqn. 9, where a p is the coefficient term of the matrix o f velocity equations and H represents convection and diffusion as discretized equation coefficients multip lied by neighboring velocities. The pressure equation, Eqn. 10, does not include the neglected terms fro m Eqn. 7. p ,h Dx Dt (8) H    1 p    ap x  ap     p ,h (9) Dx 0 Dt (10) A benefit of this pressure equation is that most of the terms are linear in relation to pres sure. In the case of constant density, an incompressible form is established. Schmidt et. al [4] used this idea in a two-step projection method on a staggered mesh approach. This method was used for their two-d imensional structured grid solver. In order to apply this model to a threedimensional solution and for an unstructured, general polyhedral mesh, a collocated variable approach is used. The numerical details can be found in Gopalakrishnan and Schmidt [5]. Flu id properties are obtained using REFPROP (Reference Flu id Thermodynamic and Transport Properties Database) fro m the NIST (National Institute of Standards and Technology) database. Given two thermodynamic inputs, REFPROP calls its necessary subroutines to calculate the thermodynamic and transport properties requested by the user. The user selects the two thermodynamic properties to use as inputs and selects an appropriate range and step size. The properties calculated at these specified points are arranged in a table format. Properties can be calculated for many flu ids, including mixtures; in this case the working fluid is water. Typically, five to ten PISO/secant iterations are emp loyed, and the solution of the pressure equation is necessary for every iteration. For the full pressure equation, a diagonal incomp lete LU preconditioned biconjugate gradient is used. Once the pressure equation is solved, the pressure field is used to correct the flu xes and the time step is fin ished. The structure for solving these equations is provided by OpenFOAM, wh ich allows rapid construction of CFD codes in an object-oriented framework [6]. OpenFOAM also provides a message passing interface (MPI) that can run decomposed cases in parallel. LEVY-B ROWN CONDENS ING EJ ECTOR The nozzle being modeled in the CFD simu lations is one that has been used in experiments by E.K. Levy and G. A. Brown [7], and is shown below in Figure 2. Wall Gas inlet Li qui d inlet Outflow Centerline axis Figure 2: Two-dimensional representation of nozzle used by Levy and Brown [7] This nozzle acts as a condensing ejector for twodimensional flo w of a flu id. In the figure above, the image can be mirrored around the centerline axis. For these simulations, saturated vapor enters through the annular channel and saturated liquid enters near the centerline. As the saturated liquid and saturated vapor flow past the splitter, the two streams should start to mix. Around the diverging part of the nozzle there is a shock wave that leads to more mixing and circulation of the two streams. The location and size of the shock wave is difficu lt to determine at this time. Due to the unstableness of the flow that occurs during condensation, the validity and accuracy of the mixing prediction both before and after the throat are in question. Choosing appropriate boundary conditions are essential to running a meaningful simu lation. Because of the strong pressure gradients and general instability of the flow, find ing acceptable boundary conditions can be difficult. The velocity at the gas inlet is set to a value much lower than what was used in Levy and Brown [7]; 75 meters per second compared to approximately 472 meters per second in the experiments. Zero pressure gradient boundary conditions and specified densities , 0.5 kg/m3 at the gas inlet and 996.7 kg/ m3 at the liquid inlet, are used at the inlets. At the outlet, a specified pressure is imposed at an arbitrary distance, creating a partially non-reflection boundary condition. The outlet boundary condition on velocity is zero gradient. A ll walls are adiabatic. RES ULTS The simulat ions were run for over 0.3 s of simu lated time. During this period, the flow never achieved a steady state due to diffuser stall. Vortices were created due to flow separation, as shown in Fig. 3. These vortices create oscillations in pressure and convolute the condensation shock. The upper half of Fig. 3 shows that, by the beginning of the diffuser section, the high enthalpy of the gas mixes with the lo w enthalpy of the liquid. This mixing predisposes the vapor to condense by lowering its temperature. In Fig. 4, the mass fraction of vapor declines just after the end of the high enthalpy zone. Pressure shows a modest increase, as shown in the upper half of Fig. 4. However, the vortices created by the flow separation cause locally low pressures that reduce the effectiveness of the condensing ejector. This separation is exacerbated by the strong adverse pressure gradient in this part of the flow. Figure 3. The enthalpy near the diffuser is shown in the top half of the figure. In a mirror image, the streak lines, indicating the instantaneous velocity field, are shown in the bottom half of the figure. Figure 4. The pressure field in the co mputational domain is show in the top half of the figure wh ile the mass fraction of vapor is mirrored in the bottom half o f the figure. Figure 5. The velocity magnitude in the computational do main is show in the top half of the figure while the time derivative of mass fraction of vapor is mirrored in the bottom half of the figure. Figure 5 shows the deceleration that accompanies the pressure increase and condensation shock. As can be seen in the upper half of the figure, the shock is oblique. Locally low pressures, corresponding to the locations of vortices, are visible in the diffuser section. The bottom half of Fig. 5 shows the high rate of condensation in the shock. Though not clearly visible in the figure, small amounts of liquid can vaporize in the high-shear zone near the inlet. Though, in these images, the condensing ejector is producing the desired static pressure increase, the strongly oscillatory nature of the flow caused very large transient variations in pressure. Depending on the time considered, the static pressure at the exit could be below the upstream pressures. Further investigations will include time-averaging of performance and investigations of the acoustic reflect ivity of various exit boundary conditions. CONCLUS IONS The idea embodied in the Ho mogenous Relaxation Model was used to construct a CFD model of condensing ejectors. Though further model adjustment and experimental validation is likely required, it is an encouraging early result. To date, the authors are unaware of any existent CFD simu lations of condensing ejectors. The current imp lementation benefits from the flexib ility of the OpenFOAM libraries and is fully parallelized. The preliminary results show the potential for a condensation shock to occur at the beginning of the diffuser section of the condensing ejector, triggering diffuser stall. The separated flo w causes a highly unstable flow field and interferes with ejector efficiency. The shed vortices produce locally low pressures and contribute to the oscillations. Continuing work will use experimental validation to adjust the model to improve its accuracy and check the overall applicability to predicting condensation shocks. It is also likely that a more sophisticated model of interphase turbulent mixing will be required, since the mixing rate strongly affects the performance of the ejector. These areas will be exp lored in future work. ACKNOWLEDGEMENT This material is based upon work supported by the National Science Foundation, STTR Phase II Project No. 0822525 and by the Depart ment of Energy, under Award Nu mber DE-FG36-06GO16049 DISCLAIMER This report was prepared as an account of work sponsored by an agency of the United States Govern ment. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liab ility or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific co mmercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Govern ment or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Govern ment or any agency thereof. REFERENCES [1] Bilicki Z., Kestin J.: Physical Aspects of the Relaxation Model in Two-Phase Flo w, Proc. Roy. Soc. London A., 428(1990), 379-397. [2] Downar-Zapolski P., Bilicki Z., Bo lle L., Franco J.: The Non-Equilibriu m Relaxation Model for One-Dimensional Flashing Liquid Flow, IJMF, 22(1996), No. 3, 473-483. [3] Prandtl, L.: About Fully-Developed Turbulence (Über d ie ausgebildete Turbulenz), Z. Angew. Math. Mech., 5(1925), 136-139. [4] Schmidt D.P., Corradin i M.L., Rut land C.J.: A Two-Dimensional, Non-Equilibriu m Model of Flashing Nozzle Flo w, 3rd ASME/JSM E Joint Flu ids Engineering Conference, 208(1999), No. 616. [5] Gopalakrishnan S., Sch midt D.P.: A Co mputational Study of Flashing Flow in Fuel Injector Nozzles, SA E Int. J. Engines 1(2008), No. 1, 160-170. [6] Weller H.G., Tabor G., Jasak H., Fureby C.: A Tensorial Approach to Computational Continuum Mechanics Using Object-Oriented Techniques, Computers in Physics, 12(1998), No. 6, 620. [7] Levy E.K., Bro wn G.A.: Liquid-Vapor Interactions in a Constant-Area Condensing Ejector, Journal of Basic Engineering, 94(1972), No. 1, 169-180.