- Email: [email protected]
A computational study about the types of entropy generation in three different R134a ejector mixing chambers
Accepted Manuscript Title: A computational study about the types of entropy generation in three different r134a ejector mixing chambers Author: J. Sierra Pallares, J. García del Valle, P. García Carrascal, F. Castro Ruiz PII: DOI: Reference:
S0140-7007(15)00347-3 http://dx.doi.org/doi:10.1016/j.ijrefrig.2015.11.007 JIJR 3199
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
31-7-2015 5-11-2015 7-11-2015
Please cite this article as: J. Sierra Pallares, J. García del Valle, P. García Carrascal, F. Castro Ruiz, A computational study about the types of entropy generation in three different r134a ejector mixing chambers, International Journal of Refrigeration (2015), http://dx.doi.org/doi:10.1016/j.ijrefrig.2015.11.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A computational study about the types of entropy generation in three different R134a ejector mixing chambers Sierra Pallares, J.a,*, García del Valle, J.b, García Carrascal, P.a, Castro Ruiz, F.a a Departamento de Ingeniería Energética y Fluidomecánica - Grupo de Ingeniería de Fluidos Escuela de Ingenierías Industriales, Universidad de Valladolid Paseo del Cauce 59, 47010 Valladolid (Spain) b Área de Energía, Facultad de Ingeniería Civil y Mecánica Universidad Técnica de Ambato Avd. Los Chasquis y Rio Payamino s/n, Ambato (Ecuador) *Corresponding author Email addresses: [email protected] (Sierra Pallares, J.), [email protected] (García del Valle, J.), [email protected] (García Carrascal, P.), [email protected] (Castro Ruiz, F.) Highlights: • • • •
CFD based second law analysis for three ejector geometries. Division of the entropy generation into four basics machanims. Discussion of flow features based on CFD local entropy generation maps. Comparison of the compression mechanisms inside the ejector.
Abstract The purpose of this investigation is to study the efficiency of a R134a ejector operating with three different mixing chambers by means of a CFD based entropy generation analysis. With the aid of the differential equation for entropy, the local entropy generation is pursued. Using this method, the areas where the irreversibilities occur are identified and geometry improvements suggested. The novelty of the numerical procedure presented herein is that the bulk entropy generation is analysed by means four sources related with viscous dissipation and heat transfer, divided in mean and fluctuating terms. Entropy generation within the boundary layer has also been considered. The latter has proved to be small, as well as the heat transfer contribution. The fluctuating viscous dissipation accounts for more than 75% of the added entropy, being its main sources the shear layer after the nozzle exit and the shock wave trains, independently of their position. Keywords: Ejector, entropy, CFD, R-134a, shock wave Nomenclature A area [m2] M mach number p pressure [bar] s specific entropy [J·kg−1·K−1] sg volume rate of entropy generation [W·K−1·m−3] Sgen rate of entropy generation [W·K−1] T temperature [K] u velocity component in the x direction [m·s−1] v velocity component in the y direction [m·s−1] w velocity component in the z direction [m·s−1] q heat flux [W·m−2] v velocity vector [m·s−1] c control volume surface
Page 1 of 17
control volume Greek symbols: α thermal diffusivity [m2·s−1] ϵ dimensionless performance parameter ε turbulent dissipation rate [W·kg−1] κ conductivity [W·m−1·K−1] µ viscosity [Pa·s] ρ density [kg·m−3] σ entropy flux [W·m−2·K−1] τ stress tensor [J·m−3] ω entrainment ratio Subscripts: heat transfer term h m motive flow e entrained flow f flow cond condenser µ viscous dissipation term τ indicates friction velocity mean mean terms contribution to entropy fluc fluctuating terms contribution to entropy pro production rate Superscripts: ~ mean term (Favre averaged) – mean term ʺ fluctuating term T turbulent + indicates wall coordinate ʹ per unit area Acronyms: SKE Standard k − ϵ turbulence model RKE Realizable k − ϵ turbulence model RNG RNG k − ϵ turbulence model SST SST k − ω turbulence model AAD Average Absolute Deviation CFL Courant - Friedrichs - Lewy CFD Computational Fluid Dynamics 1. Introduction Supersonic ejectors have been used widely for more than a century. In refrigeration systems they may contribute to reducing energy consumption because they can be powered by solar energy or by waste heat (García del Valle et al., 2014). An ejector is a simple device in which a high pressure stream (primary fluid) is used to compress a low pressure stream (secondary fluid) to a higher pressure (discharge pressure). In jet-ejector cooling systems, a refrigerant stream of low velocity and low pressure is evaporated to generate a cooling capacity, and compressed to an adequate condensation pressure level by the interaction with a high velocity stream generated by heat input and adiabatic near-reversible expansion.
Page 2 of 17
The current trends in ejector research are focused in minimizing the ejector internal irreversibility, considering the flow deceleration the most important factor affecting entropy generation. In order to design and develop a high performance ejector, a clear understanding of the flow and mixing inside the ejector is first needed. Thus, modelling seems essential to afford this task. Although one dimensional models are useful for global state-space models, (Xue et al., 2015), it is commonly accepted that these models are not sufficient to understand the physics of a supersonic ejector (Bartosiewicz et al., 2005). For instance, (Sriveerakul et al., 2007) remarks that these models only can be used when the ejector is operated at its design conditions. Thus, the use of CFD codes is highly recommendable for the understanding of ejector flow physics, and in turn, once a figure of merit is defined, they may improve the ejector design, (Zhang and Mohamed, 2015). The figure of merit used for ejector geometry optimization, independently of its formulation, is ultimately related to the entropy generation minimization in the ejector. The pitfalls encountered in CFD computations are related with the selection of the turbulence model (Gagan et al., 2014), which makes a CFD validation mandatory in order to asses the uncertainties attained in the computed results. Entropy generation minimization is the method of modelling and optimization of real devices that owe their thermodynamic imperfection to heat transfer, mass transfer, and fluid flow irreversibilities (Bejan, 1996). In order to minimize entropy generation, its sources should be studied. Entropy production in turbulent flows can be determined using direct or indirect methods (Herwig and Kock, 2006). In the direct method, the differential entropy transport equation is employed to calculate the entropy production terms in the post-processing phase of a CFD calculation. That means they are determined by using the known field quantities (time-mean) velocity, temperature and composition. Using this method, different local parts of entropy production (due to viscous dissipation and heat and mass transfer) can be determined. Integration of these field quantities over the whole flow domain results in the overall entropy production rate. In the indirect method, entropy generation is computed using a control volume, in a black - box approach, with the integral entropy transport equation, (Kock and Herwig, 2004). Though this equation is exact, its application for practical problems involves several approximations, specially regarding the entropy flux terms across the inlet and outlet boundaries, normally neglected since unknown. Further approximations are involved when inlet and outlet quantities are assumed to have uniform values in the cross sections, neglecting their actual profiles. Some studies have addressed the problem of entropy generation in ejectors for refrigeration applications by means of the latter approach. (Arbel et al., 2003), using the indirect method, performed an analysis of the ejector as a whole and of each of its internal processes. They showed that entropy production is equivalent to performance losses, and that minimizing entropy production could be used as a criterion for optimization. They also indicated that the ejector irreversibility is due to three main factors: “pure mixing”, “kinetic energy losses”, and the normal shock wave produced in the mixing chamber. (Liu and Groll, 2013) presented a method to determine the internal-ejector component efficiencies by combining zero-dimensional modelling with measurements performed for a R744 ejector refrigeration cycle. Measured and simulation data were used to determine the isentropic efficiencies of the motive nozzle, the suction nozzles, and the mixing section efficiency. On the CFD side, (Banasiak et al., 2014) investigated about the performance of R744 ejectors. In that work, a new approach was introduced to assess the contribution of the local irreversibilities to the overall entropy increase
Page 3 of 17
inside the ejector. This analysis was based on close examination of the flow variable profiles, as well as on the energy and the entropy balances formulated for individual ejector passages using the indirect method. This approac enabled the precise location of intensified irreversibility sites and delivered detailed information on the irreversibilities origin. To the best of the author knowledge, all of the studies described thus far have been limited to the application of the indirect method. Therefore, our aim in this work is to apply the direct method to the analysis of different ejector chambers for refrigeration applications. It is clear that the direct method is superior in the analysis of complex flow applications, such as ejectors: from the direct method we obtain information about the total entropy production and how this production is distributed between the single mechanisms conforming the flow pattern; information the indirect method cannot provide. Therefore, the present study extends the scope of the existing literature by providing relations to compute the entropy generation using CFD computations combining the ideas of other authors (Bejan, 1982; Kock and Herwig, 2004) and allowing us to discriminate between the different entropy production types in an easy and accurate way. In this way, we shall provide information that cannot be obtained using the indirect method, thus adding value to the existing literature on the topic. The rest of the paper is developed in the following way: The computational model is first described, along with detailed information about the computation of the entropy production rate in ejectors both direct and indirect methods. Secondly, the three mixing chambers under study will be detailed geometrically, along with the boundary conditions employed in the analysis. The last Section is related to the comparison between experimental results and CFD simulations and a discussion regarding the entropy generation rate in the three chambers under study. It is important to note that the model validation was already addressed in a previous work, (García del Valle et al., 2015). A short account of the results found will be included in the description of the computational model. 2. Computational model 2.1. Fluid dynamic and thermodynamic models. Numerical implementation The commercial software ANSYS Fluent 14.5 has been used to compute the compressible steady-state, 2D axis-symmetric turbulent flow governing equations which apply inside the ejector. This software is based on the finite volume approach, a technique that discretizes all governing equations by first dividing the physical space into a number of control volume grid elements. When dealing with turbulent flow, a suitable averaging method must be used. The density averaging technique or Favre averaging (Versteeg and Malalasekera, 1996) is appropriate for compressible flows and will be used in this work. The total energy equation is included and coupled to the mass and momentum equations and the thermodynamic model. Turbulence modelling is based on the Standard k - ϵ model (Launder and Spalding, 1972b). Real gas thermodynamic models must be considered due to the departure from ideal gas found in the region of superheated vapour, close to the saturation line, a region in which the ejector commonly operates. The REFPROP library, (Lemmon et al., 2002), integrated into ANSYS Fluent, has been used as the source of thermodynamic data. For the case of R-134a, this formulation is based in the multiparameter equation of state of Tillner-Roth and Baer (Tillner-Roth and Baehr, 1994). This equation constitutes a 21 parameter model that represents nearly all selected experimental data within their estimated accuracy with the exception of heat capacities and speed of sound in the region close to the critical point. Typical accuracies are 0.05% for density, 0.02% for the vapor pressure, or 0.5 and 1% for the heat capacity. A second order
Page 4 of 17
upwind scheme was selected for the advective terms of each equation except for the pressure, where a Least Squares Cell Based approach is chosen. High order term relaxation was also applied. The density - based algorithm was selected for the pressure - velocity coupling problem. A strict Courant-Friedrichs-Lewi (CFL) condition is needed here in order to prevent the solution from diverging during the first integration steps. This is due mainly to the thermodynamic state falling outside of the vapor (both stable and metastable) regions when using the REFPROP library. A CFL of 0.8 is used for the first 500 steps, and then gradually increased to reach a value of 2.5. 2.2. Computation of the entropy generation rate 2.2.1. Direct method The entropy in a flow field is a state variable that in its specific form has a balance equation that for a single phase, non-reacting and single component flow reads (Sciacovelli et al., 2015):
1 1 q = 2 q T τ : v Dt T T T
Ds
(1)
which can recast into
Ds
= σ s g
Dt
(2) Where σ is the entropy flux vector and sg is the entropy generation rate per unit volume. A comparison between (1) and (2) leads to σ =
q T
(3) 1
sg =
q T
2
T
1
τ : v
T
(4) It can be observed that the total entropy generation rate can be splitted in two contributions, each one strictly related to a specific transport phenomenon: a contribution due to heat transfer (6) and a viscous dissipation term (7): sg = sg
sg
h
(5)
where sg sg
h
=
1 T
=
2
q T
1
(6) τ : v
T
(7) In a turbulent flow, Equation (1) has to be time averaged in order to find an equation for the mean entropy. For a compressible flow, Favre averaged is employed to obtain
u s v s w s = sg Dt y z x
Ds
h
sg
(8) Unclosed terms are modelled using the method proposed by Herwig and Kock (Herwig and Kock, 2006): the entropy generation rate is divided between the contribution of the mean gradients and the contribution of the fluctuating components:
Page 5 of 17
sg
= sg
h
sg
h
= sg
sg
mean
h
sg
(9)
fluc
(10) where the contibutions to the entropy generation rate due to mean gradients are: sg
= h
m ean
2 T
T
m ean
fluc
2 2 T T y z x
2
u
sg
m ean
2
2
v w 2 = z T x y
2
(11) u v 2 y x
v w u w z x z y
(12)
2
2
and the fluctuating ones are: sg
=
h
sg
fluc
T
sg
h
(13)
mean
= ~ T
(14) Although equations (11) - (14) can be used in every point of the computational domain, special care is needed when using wall functions in the computation of the turbulent flow. In this work, standard wall functions are used and in order to obtain the correct estimation of entropy generation in the cells adjacent to the walls, the method proposed by Walsh and McEligot (Walsh and McEligot, 2008) is used. Under boundary layer approximations for an un-heated two-dimensional flow, the pointwise entropy generation rate may be approximated (using wall coordinates) as:
s
fluc
2
g
u uv y
u y
(15)
where
uv
u 1 y
(16) The value of pointwise entropy is desired to identify regions where most losses occur and to deduce the integrated entropy generation per unit surface area
s y
g
=
T s g
u
3
s y
=
0
g
dy
=
y
0
u y
s y
g
. It can be written as
dy
(17) In this way, the total entropy generation is easily obtained given the wall function employed and the entropy generation per unit area due to the wall function computed as:
Page 6 of 17
u 3 y u s g = 0 y T
dy
(18) Thus, given the area of the computational cell and the value of y+, the total entropy generated in the boundary layer can be computed. In a standard computation, the continuity, momentum and energy equations are first solved. Then, in a postprocessing step, Equations (11), (12), (13), (14) are computed at cell centers except in the first row of cells adjacent to the wall, where Eq. (18) is applied as commented. All of these apportations to the total entropy generation rate constitutes the entropy production rate per unit volume, that can be plotted as map over the computational mesh. 2.2.2. Indirect method From the entropy transport equation (2) is possible to derive its integral form, which constitute the basis of the indirect method (please see the appendix for a detailed derivation): d
dt
s d c
c
q n d c c T
s v n d c =
sg d f
(19) Recasting the last equation, it is possible to obtain the entropy production rate S p ro for a given control volume, S pro =
sg d = f
d dt
s d c
s v n d c c
q c T
n d c
(20) Therefore, the local rate of entropy production can be computed by determining the right hand side of (20). As explained by Kock and Herwig (Kock and Herwig, 2004), this method has the drawbak of the necessary computation of the integral of the entropy flux term:
q n d c c T
(21) which is always neglected in the literature, or even unmentioned (Arbel et al., 2003; Banasiak et al., 2014). However, the relative value of this term is of importance in some situations. If the control volume selected cuts down the geometry under study in a region of big temperature gradients, the error produced in the computation of the entropy production may be significant due to the high value of the heat flux accross the control surface. Validation of the CFD model 2.3. In a previous study, (García del Valle et al., 2015), the CFD tool was subjected to a four step validation procedure, related with the real gas implementation, turbulence modelling, the shock wave-boundary layer interaction and the ejector simulation itself using different validation cases extracted from the literature. The purpose of the validation cases was to independently address the different flow features which appear in the ejector. The general results were as follows: 1. The real gas implementation was compared with the shock tube (Riemman problem) solution for a R-134a vapour flow. Albeit some diffusion was found in the shock and expansion wave, their location and the jump in thermophysical properties were correctly predicted. 2. Four turbulence models, according to Table 1, were considered. Experimental data of a supersonic air jet expanding into quiescent air, (Eggers, 1966), was used as a test case for this purpose. It was found that, although the supersonic region of the jet was accurately
Page 7 of 17
predicted, the numerical results in the subsonic region exhibit higher diffusive rates than the experimental case. The SKE and SST model had the lower overall error. 3. The shock wave-boundary layer interaction was compared against the experimental work of (Simeonides, 1994), who studied the Mach 5 air flow over a hollow cylinder flare. Both the pressure bump before the shock wave and the reattachment point were predict with accuracy with the SKE and RKE, while the RNG and SST were poorer in performance. 4. For the ejector validation, both the entrainment ratio and the pressure profiles along the wall were considered for four cases, one operating in the mixed regime and the rest in the double choked regime. The ejector geometries are the same as those studied in this investigation. The results are best understood if divided into the mixed regime and the double choked regime. For the former, a discharge pressure dependent value for the entrainment ratio was found both experimentally and numerically. The average value over a time span was used to compute the entrainment ratio. It was found that, due to the nature of the flow, where the turbulence modelling is of paramount relevance, the entrainment ratio was poorly predicted by all turbulence models. The pressure profile, however was accurately predicted. Further numerical studies shown that the pressure profile is decoupled from the entrainment ratio for low values of the latter, which explain the incongruence describe previously. For the double choked regime, both the SKE and SST model yielded the best results, with an overall error of 15% for the entrainment ratio and 5% for the pressure profile. Based on the validation phase, the SKE turbulence model will be used in the present investigation. 3. Ejector geometries under study The R-134a refrigerant ejector was exhaustively described in (García del Valle et al., 2014). The reader is referred to the original source, as only a short summary will be included here focused on the experimental uncertainties. The experimental installation, rated at 10.5 kW of boiler and evaporator power, was designed to accurately control and measure the mass flow, pressure and degree of superheating of the primary and secondary mass flow and the discharge pressure and temperature of a number of ejector geometries. For this work, the results from three geometries as shown in Fig. 1 have been employed. The main characteristics of each geometry are as follows: 1. Chamber “A” is the classical shape with a 30° half-angle conical entry and a constant area zone with a length to diameter ratio of 8.5. 2. Chamber “B”, promotes a supersonic compression in a convergent region rather than a shock system in a constant area region. That is why the present design has a 1.14° total conical entrance and a short constant area region. 3. Chamber “C”, pursues a supersonic compression through a sharp contraction. This geometry has two constant area regions with diameters of 5.5 mm and 4.8 mm respectively. Each geometry is tested with the same primary and secondary conditions, only the discharge pressure is modified in each run. The discharge pressure values range from 6.8 bar up to the point in which the entrainment ratio tends to zero. It is worth mentioning that operation close to zero mass entrainment is unstable, although the secondary pressure remains approximately stationary due to the mass and thermal inertia of the systems, oscillations were found for the secondary mass flow rate. Furthermore, due to the poor regulation both of the
Page 8 of 17
resistor and of the electronic expansion valve at low mass flow rates, slight deviations from the set point are found. All the relevant experimental data, in which the numerical investigation is based, are included in Table 3. Further columns have been added and will be used later on. The nozzle exit is positioned at -5.38 mm in all measurements. The uncertainties of the different measured variables are shown in Table 2. The following comments are needed in order to fully understand this table. The column “Range” represents the calibration range of the sensor. The term “Dev. FS” represents the deviation over full scale. It should be noted that two terms are provided, the first is attributed to the sensor and the second to the data acquisition system. The column “Dev. abs. max” is the maximum absolute deviation, computed as (“Max Range”-“Min Range”)*“Dev. FS”. An uncertainty analysis, performed according to (Moffat, 1988), reveals that the entrainment ratio error is equal to ± 0.0061. Each point corresponds to the average value over a time span of 5 minutes (which corresponds to 60 logged data points, since the data acquisition rate is equal to 0.2 Hz), in which the deviation of selected measured variables is within a predefined limit to ensure the steady state condition. 3.1. Computational details The SKE turbulence model is employed to predict the performance of the ejector using the different mixing chambers described. The flow was considered in steady state, 2D axis-symmetric and compressible. Second order upwind schemes are used for the flow variables except for pressure, where a Least Squares Cell Based approach is employed. As an example, the computational domain along with boundary conditions for Chamber B is sketched in Figure 2. Pressure inlets are modelled using the values of static pressure and temperature as total pressure and temperature values respectively, due to the low value of the velocity in these inlets in relation to the velocity in the chamber. Turbulence intensity is 5% in each inlet, and the other boundary conditions (pressure and temperature at inlets and discharge pressure) are specified in Table 3. Different grid sizes have been employed for the computational domain, looking for mesh independence, ranging from 40000 to 150000 cells. The solution showed independence form the mesh using 50000, 90000 and 60000 cells for chambers A, B and C respectively. Mesh adaptation has been carried out based on Mach number gradients and pressure gradients. CFL condition is fixed to 1 for the first hundreds of iterations, and then progressively increased to a value of 2.5 to prevent the solution from diverging during the first integration steps. This is due mainly to the thermodynamic state falling outside of the vapor (both stable and metastable) regions when using the REFPROP library. The mean y+ value for all the chambers under study is around 125. Results 4. The numerical results are structured in two subsections. First, the entrainment ratio and the overall entropy generation rate are addressed. Experimental and numerical results will be compared at this point. Second, contour plots charts for the local entropy generation rate will be reported. Based on these results, a comparison of the three geometries will be addressed. 4.1. Prediction of the entrainment ratio. Total entropy production This section is useful for validation purposes, since both the entrainment ratio and entropy generation is compared against experimental values. The former is shown in Figure 3 whilst the latter is included in Figure 4. Fig. 3 shows the comparison between the experimental and numerical entrainment ratio for the three mixing chambers under study for a range of discharge pressures. A non dimensional performance parameter ϵ is used in this section to discuss the performance of the different
Page 9 of 17
geometries. This parameter was introduced in (García del Valle et al., 2014) and is defined as the quotient between the entrainment ratio and the entrainment ratio of a isentropic ejector with the same primary and secondary flows stagnation conditions and discharge pressure. Therefore, the parameter ϵ is inversely related with the entropy generation. Its advantage over the latter is that its value is bounded between 0 and 1, and its meaning its similar to that of the isentropic efficiency. Numerical results compare better to experimental data for mixing chamber “A”, where the double choked regime is accurately predicted by the model. For chamber “B”, a mean deviation of 14% is found, close to the previously reported value in (García del Valle et al., 2015). Mixing chamber “C” operates in the mixed regime in the selected discharge pressure range. While the critical pressure is successfully computed for the mixing chamber “B”, for mixing chamber “A” it is overestimated. Fig. 4 shows the entropy production in each chamber computed using the indirect method, Eq. (20). It is worth to mention here that to avoid the computation of the entropy flux term, Eq. (21) the control surface is located far away from regions were temperature gradients are found. In this way, it is safe to neglect it in the computation of the entropy production. From Figure 3 and 4 the aforementioned inverse relation between the entropy generated and the performance parameter ϵ is clearly seen. The peak in the performance parameter match the minimum in entropy generation. The experimental values of entropy generation (see last column of Table 3) are lower than those obtained numerically. For instance, the minimum numerical values are 2.3 W·K−1 and 2.65 W·K−1, against the minimum experimental values of 1.67 W·K−1 and 2.1 W·K−1 for mixing chamber “A” and “B” respectively. Although it is a deviation close to 40%, is worth noting the strong correlation which exists between the variation of the entropy generation with the discharge temperature. Simple computations with the control volume approach yield a relation of 0.2 W·K−1 per K of discharge temperature increment. So the variation of numerical to experimental discharge temperature is in the range of 3 to 4 K. Fig. 5 reports the comparison of the direct and the indirect methods in computing the entropy production of the ejector. When using the direct method, the entropy generation flowfield is first computed and then integrated for the whole computational domain in order to obtain the entropy production for the ejector. Theoretically, both methods should give the same results of total entropy generation, however, according to Fig. 5, a deviation within ± 10% of the latter over the former is found. This difference can be attributed the poor convergence of some cases under study, specially in the mixed regime, where a non stady flow is found (the procedure used to handle these cases is reported in (García del Valle et al., 2015)). 4.2. Entropy generation analysis In this section, regions of the ejector where entropy generation occurs are identified, along with the mechanism which produce every type of irreversibility. Figs. 6 shows the contour plots of Mach number, mixture fraction and entropy generation for the mixing chamber “A”. The mixture fraction is computed as a passive scalar, with boundary conditions equal to one for the primary flow and 0 for the secondary flow, and provides a measure of the extent of mixing at the macro-scale. From the Mach number distribution, the geometry “A” operates with supersonic flow up to the end of the constant area region, where a shock wave train adjust the pressure of the mixture to the discharge pressure. Mixture fraction maps reveals that the primary and secondary flows do not completely mix over the constant area region, but rather the mixing
Page 10 of 17
is finished through the train of shock waves, were high turbulence is found. Provided the moderate rate of mixing in the constant area region shown in Fig. 6, the Mach number remains approximately constant through this region. Due to the increment of the shear layer between the primary and secondary flows through the constant area region, the rate of entropy generation slightly falls along this region. As expected, the irreversibility is concentrated within the shock wave train. Fig. 7 depicts the entropy production per unit length, decomposed in its different sources. This figure is produced integrating the entropy generation map at selected surfaces where the axial distance is constant along the computational domain. The sharp peak depicted in the figure shows how the main sources of irreversibility are concentrated in the shock train region. Besides, the rate of entropy generation within the boundary layer is small compared with the bulk entropy generation. This fact is repeated for all geometries. The long tapered mixing chamber introduced in geometry “B”, changes the flow pattern described for geometry “A” (See Fig. 8). With this geometry, a compression occurs through the convergent region, which translates into a decrement of the Mach number. Furthermore, the mixture between the primary and secondary flow is almost completed at the entrance of the constant area region. In this short region, the flow is transonic. Through the diffuser the mixture decelerates up to the discharge pressure. Fig. 9 shows the entropy production per unit length for mixing chamber “B”. Through the long tapered region, both the compression and the shear layer contribute to the entropy generation. But, since the velocity gradient is smaller than for geometry “A”, the entropy generation is slightly lower. As for the former case, a peak of irreversibility is found at the end of the constant area region, though this time is smaller, mainly because the bulk velocity is lowered. The geometry “C” (see Fig. 10) resembles the operation of an unadapted two throat system. A shock wave occurs at the nozzle exits in order to obtain a subsonic flow at the entrance of the “secondary” throat, i.e. the constant area region, in which the mixture accelerates again up to sonic flow. All the entropy generation is thus concentrated immediately after the nozzle exit, as is clearly visible from Fig. 11. In order to analyze the different entropy generation mechanisms, as described in Section 2.2, the Fig. 12 has been included. It is a multiple pie chart in which the different sources of irreversibility are identified for the different geometries and ejector regions. One of the superior aspects of the direct method over the indirect one is the capability of determining the source of irreversibilities. From Fig. 12, it can be seen that the main entropy generation mechanism is viscous dissipation for all the chambers under study. Irreversibilities induced by heat transfer (due to mean gradients and fluctuating ones) merely account for less than 2% of the total amount. Furthermore, fluctuating viscous dissipation is the main entropy source, accounting for a share of more than 75%. This figure is obtained using the SKE model, but can be reproduced using other RANS approaches. For example, using the SST model, small differences are found for the cases under study. If we pay attention to the entropy generation per zone, we observe that in the region previous to the mixing chamber (what we call “premixing” region) entropy generation due to the mean dissipation Eq. (12) is predominant. This is because the main source of this type of dissipation is located close to the walls, completely in the boundary layer of the supersonic flow inside the nozzle. There, the velocity gradients are very important and contribute the most to the
Page 11 of 17
entropy production. In this region, fluctuating dissipation is unimportant because the mixing process has not even started. This effect is common to all geometries. However, in the mixing chamber, we observe different behaviours for each geometry. In this region, most of the dissipation is due to the fluctuating one. In mixing chamber “A”, most of the entropy generated due to the fluctuating dissipation is located in the diffuser, but in the other geometries the behaviour is the opposite. This is connected with the location of the normal shocks and the mixing processes that occur inside the ejector. Depending on the location of the shock and the mixing intensity, this type of entropy generation is placed mostly in a region or another. For mixing chamber “A”, the wake of the normal shock is responsible of most of the irreversibility generated in the diffuser (see the entropy generation map of Fig. 6) but since mixing is relatively slow, the mixing chamber do not present a significant entropy generation of this type due to mixing. For mixing chambers “B” and “C”, due to the increased mixing rate, most of this irreversiblity is located in the mixing chamber itself, although the shock is located downwards the constant area region in the “B” geometry. For mixing chamber “C”, both mixing and normal shock happen prior to the reduction of diameter, and that is why up to 88% of the fluctuating dissipation is located backwards. Regarding the entropy generation due to mean and fluctuating temperature gradients (i.e. heat transport inside the ejector), the comments of the above paragraph hold. The temperature difference of the primary and secondary flow is dissipated by means of mixing, and there is a strong temperature gradient associated to the normal shock. That is why this source of irreversibility more or less matches with the entropy generation due to fluctuating dissipation. 4.3. Design optimization by means of entropy generation minimization The end point of this investigation should provide some guidelines to the design of high performance mixing chambers. The entropy generation should be taken with caution, because, although for an operation condition one configuration can exhibit the lower entropy generation rate, over different conditions this is not true. Furthermore, although the entropy is closely related with the performance parameter ϵ, it does not mean that this relation is inversely linear. From Fig. 4, mixing chamber “B” perform better in the lower discharge pressure range against mixing chamber “A”, which does better in the upper discharge range. Not only that but, the differences in the minimum entropy generation 2.3 and 2.65 W·K−1 for mixing chambers “A” and “B” respectively are not directly translated to the performance parameter ϵ, which is almost similar from one to another. The key point in the discussion is to analyze why the geometry “B” does not improve the performance over the mixing chamber “A”. Our vision is somehow skewed by the description given earlier for the double choked regime, i.e. that the entropy generation after the constant are region is smaller for geometry “B” (softer shock wave train) and that the irreversibility lost into the shear layer is smaller or of the same order as that of geometry “A”. However, as the critical point is reached, the entropy generation in the sock wave train in geometry “A” is decreased while the entropy generated in the shear layer will remain approximately constant. At the critical point, the pressure increment in mixing chamber “A” occurs solely through the sock wave train, for mixing chamber “B” the compression occurs through the convergent region and through a softer sock wave train. From both the experimental and numerical data, the performance of both procedures is similar, although the mechanism described for geometry “B” it is “a priori” more favourable. However, there is an additional element missed in the former description.
Page 12 of 17
As for the efficiency of a subsonic diffuser (divergent area) is negatively influenced by low boundary layer velocities and boundary layer detachment, the efficiency of a supersonic diffuser (convergent area) is also negatively influenced by the low boundary layer velocities created by the secondary fluid. Coupled with the negative pressure gradient, the thickness of the boundary layer should grow way more rapidly than for geometry “A”, as is clearly visible from the mixture fraction of Figs. 6 and 8. The authors suggest that this behaviour is the cause that negligible performance increments are found for long tapered ejectors in relation with long constant are one. Thus, the possible gain induced by geometry variations in the peak ejector performance parameter ϵ is marginal. For an specific application and operating conditions, the critical feature is to get a standard geometry (say “A” or “B”) under which the operating point is closed to the critical point, this is entropy generation minimization. 5. Conclusions The mechanisms of entropy generation have been studied in this investigation based on experimental and computational results for a R134a ejector, covering three different geometries and a number of operating conditions. Using the differential balance for entropy, the four mechanisms relevant for entropy production are identified inside the ejector chambers under analysis. Between all the mechanisms analyzed, the fluctuating viscous dissipation, which bears a share of more than 75% of the total entropy generation for all geometries, is the single most important dissipative effect which drags ejector performance. Furthermore, a special treatment has been applied for dealing with the irreversibilities inside the boundary layer. Not only the total entropy generation has been identified, but rather the volumetric rate of entropy generation through the domain has been pursued. Based on these results, and for the sizes of cross-sectional areas studied herein, it is has been concluded that the efficiency of the compression through a train of shock waves is similar to that occurred in a long convergent conduit. The authors suggest that limited improvement may be gain upon geometry variation of cases “A” and “B” because the physical behaviour will keep being the same as in those geometries. An operation close to the critical point may provide better results in this regard. Acknowledgements The authors gratefully acknowledge the financial support given to the investigation by the “Junta de Castilla y León” under grant 04/09/VA/0025 and “Secaderos y Frigoríficos Industriales S.L.”. References Arbel, a., Shklyar, a., Hershgal, D., Barak, M., Sokolov, M., 2003. Ejector Irreversibility Characteristics. Journal of Fluids Engineering 125 (1), 121. Banasiak, K., Palacz, M., Hafner, A., Buliski, Z., Smoka, J., Nowak, A. J., Fic, A., Apr. 2014. A CFD-based investigation of the energy performance of two-phase R744 ejectors to recover the expansion work in refrigeration systems: Anirreversibility analysis. International Journal of Refrigeration 40, 328–337. Bartosiewicz, Y., Aidoun, Z., Desevaux, P., Mercadier, Y., Feb. 2005. Numerical and experimental investigations on supersonic ejectors. International Journal of Heat and Fluid Flow 26 (1), 56–70. Bejan, A., 1982. Entropy generation through heat and fluid flow. John Wiley & Sons Inc. Bejan, A., 1996. Entropy generation minimization: The new thermodynamics of finite-size devices and finite-time processes. Journal of Applied Physics 79 (3), 1191.
Page 13 of 17
Eggers, J., 1966. Velocity profiles and eddy viscosity distributions downstream of a mach 2.22 nozzle exhausting to quiescent air. NASA TN D-3601. Gagan, J., Smierciew, K., Butrymowicz, D., Karwacki, J., Apr. 2014. Comparative study of turbulence models in application to gas ejectors. International Journal of Thermal Sciences 78, 9–15. García del Valle, J., Saíz Jabardo, J., Castro Ruiz, F., San José Alonso, J., Oct. 2014. An experimental investigation of a R-134a ejector refrigeration system. International Journal of Refrigeration 46, 105–113. García del Valle, J., Sierra Pallares, J., García Carrascal, P., Castro Ruiz, F., 2015. An experimental and computational study of the flow pattern in a refrigerant ejector. Validation of turbulence models and real-gas effects. Applied Thermal Engineering 00, 000–000. Herwig, H., Kock, F., Mar. 2006. Direct and indirect methods of calculating entropy generation rates in turbulent convective heat transfer problems. Heat and Mass Transfer 43 (3), 207–215. Kock, F., Herwig, H., May 2004. Local entropy production in turbulent shear flows: a high-Reynolds number model with wall functions. International Journal of Heat and Mass Transfer 47 (10-11), 2205–2215. Launder, B., Spalding, D., 1972a. Lectures in Mathematical Models of Turbulence. Academic Press. Launder, B. E., Spalding, D. B., 1972b. Lectures in Mathematical Models of Turbulence. Academic Press. Lemmon, E. W., McLinden, M. O., Huber, M. L., 2002. REFPROP: Reference Fluid Thermodynamic and Transport Properties. NIST Standard Database 23. Liu, F., Groll, E. A., Apr. 2013. Study of ejector efficiencies in refrigeration cycles. Applied Thermal Engineering 52 (2), 360–370. Menter, F., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal 32, 1598–1605. Moffat, R., 1988. Describing the uncertainties in experimental results. Experimental Thermal and Fluid Science 1, 3–17. Sciacovelli, a., Verda, V., Sciubba, E., 2015. Entropy generation analysis as a design toolA review. Renewable and Sustainable Energy Reviews 43, 1167–1181. Shih, T., Liou, W., Shabbir, A., Yang, Z., 1995. A new k − ϵ eddy viscosity model for high reynolds. Computers and Fluids 24, 227. Simeonides, G., 1994. Test cases for CFD code validation for inlets and nozzles of airbreathing propulsion systems. Doc. EWP/1784, ESA/ESTEC. Sriveerakul, T., Aphornratana, S., Chunnanond, K., Aug. 2007. Performance prediction of steam ejector using computational fluid dynamics: Part 2. Flow structure of a steam ejector influenced by operating pressures and geometries. International Journal of Thermal Sciences 46 (8), 823–833. Tillner-Roth, R., Baehr, H. D., 1994. An international standard formulation for the thermodynamic properties of 1, 1, 1, 2-tetrafluoroethane (hfc-134a) for temperatures from 170 k to 455 k and pressures up to 70 mpa. Journal of Physical and Chemical Reference Data 23 (5), 657–729. Versteeg, H. K., Malalasekera, W., 1996. An introduction to computational fluid dynamics. Longman. Walsh, E. J., McEligot, D., 2008. Relation of entropy generation to wall laws for turbulent flows. International Journal of Computational Fluid Dynamics 22 (March 2015), 649–657.
Page 14 of 17
Xue, B., Cai, W., Wang, X., Jan. 2015. State-space modelling for the ejector-based refrigeration system driven by low grade energy. Applied Thermal Engineering 75, 430–444. Yakhot, V., Orszag, S., 1986. Renormalization group analysis of turbulence. i. basic theory. Journal of Scientific Computing 1, 3–51. Zhang, T., Mohamed, S., 2015. Conceptual Design and Analysis of Hydrocarbon-Based Solar Thermal Power and Ejector Cooling Systems in Hot Climates. Journal of Solar energy engineering 137 (April), 1–9. Appendix A. Derivation of the integral transport equation for entropy From the entropy transport equation (2) is possible to derive its integral form, which constitute the basis of the indirect method:
Ds
d
=
Dt
f
σd
f
sg d
(A.1)
f
since
Ds
d
d
=
Dt
f
dt
sd
(A.2)
f
and taken into account the transport theorem, we obtain a transport equation for a control volume d
dt
sd
d
=
dt
f
sd
c
s v nd
(A.3)
c
c
therefore, d dt
sd
c
s v nd
c
=
c
σd c
sg d
(A.4)
c
Applying the divergence theorem, the volume integral of the entropy flux σ can be written
σd
=
c
c
q nd T
c
(A.5)
Introducing the above equation in (25) one obtains d
dt
sd c
s v nd c
c
q = nd c T
c
sg d
(A.6)
f
Figure 1: Ejector geometries under study. Note the origin of the coordinate system for each of the three cases. Figure 2: Sketch of the computational domain and details of the meshing in selected areas for the ejector geometry under study (Chamber B) Figure 3: Comparison between the experimental and numerical entrainment ratio for the three mixing chambers. The performance parameter ϵ, based on CFD data, has also been included. Figure 4: Numerical entropy generation for the three mixing chambers. The performance parameter ϵ, based on CFD data, has also been included. Figure 5: Comparison between the “direct” and “indirect” entropy generation numerical procedures.
Page 15 of 17
Figure 6: Numerical Mach number, mixture fraction and rate of entropy generation for the mixing chamber “A”, turbulence model SKE and conditions. Case A-2 Figure 7: Rate of entropy generation as a function of the axial coordinate for mixing chamber “A”. Case A-2 Figure 8: Numerical Mach number, mixture ratio and rate of entropy generation for the mixing chamber “B”, turbulence model SKE and conditions. Case B-2 Figure 9: Rate of entropy generation as a function of the axial coordinate for mixing chamber “B”. Case B-2 Figure 10: Numerical Mach number, mixture ratio and rate of entropy generation for the mixing chamber “C”, turbulence model SKE and conditions. Case C-2 Figure 11: Rate of entropy generation as a function of the axial coordinate for mixing chamber “C”. Case C-2 Figure 12: Multiple pie chart in which the entropy generating mechanisms as well as its share in the different regions is reported for all mixing chambers. Table 1: Turbulence models considered during tha validation stage. Model Standard k − ϵ
Realizable k−ϵ RNG k − ϵ SST k − ω
Reference (Launder and Spalding, 1972a) (Shih et al., 1995) (Yakhot and Orszag, 1986) (Menter, 1994)
Main application Fully turbulent flows
Acronym used in this work SKE
Planar and round jets, rotation and separated flows High streamline curvature flows Adverse pressure gradient flows
RKE RNG SST
Table 2: Measurement uncertainties for the relevant physical properties. Magnitude Position Range Dev. FS Dev. abs. max. Pressure Primary 0-40 bar ±(0.5%+0.25%) ±0.30 bar fluid inlet Secondary 0-30 bar ±(1%+0.25%) ±0,375 bar fluid inlet Discharge Temperature Primary 0-150 °C ±(0.1%+0.25%) ±0.52 °C fluid inlet Secondary -50-50 °C ±0.35 °C fluid inlet
Page 16 of 17
Mass flow
Discharge Primary Secondary
0-100 °C 0-0.125 kg·s−1
±(0.18%+0.25%)
±0.35 °C ±5.37·10−4 kg·s−1
Table 3: Relevant measured physical properties for each mixing chamber, along with the entropy produced. Entropy generation is computed here from experimental measurements “A ”
“B”
“C”
Cas e 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 1 2 3 4 5
p0mf[bar] T0mf[°C] p0ef[bar] 29.15 29.17 29.08 29.16 29.08 29.27 29.19 29.15 29.05 29.09 29.18 29.17 29.09 29.12 29.04 29.16 29.10 29.12 28.91 29.01 29.01 28.98 28.95
95.10 94.95 94.93 95.09 94.97 95.19 95.00 95.01 94.95 95.01 95.04 95.03 94.98 95.08 94.91 95.11 95.00 95.09 94.94 95.01 94.99 95.05 94.95
4.14 4.12 4.12 4.14 4.14 4.14 4.14 4.14 4.13 4.17 4.20 4.24 4.11 4.13 4.11 4.12 4.15 4.34 4.11 4.14 4.15 4.17 4.14
T0ef[°C ] 19.96 20.00 20.00 20.02 20.00 20.01 20.01 20.03 20.02 20.26 21.97 24.04 20.00 20.00 20.00 20.01 20.86 23.47 20.10 20.01 20.07 20.22 23.34
pcond[bar ] 6.82 7.00 7.17 7.38 7.58 7.81 7.99 8.19 8.39 8.58 8.81 9.00 6.79 7.03 7.17 7.36 7.80 8.04 6.83 7.00 7.19 7.29 7.59
Tcond[°C] me[kg·s−1] mm[kg·s−1] Sgen[W·K−1] 46.96 47.39 47.11 47.66 48.45 48.85 48.83 49.14 50.68 52.86 54.95 56.32 45.83 46.33 46.58 47.01 52.73 55.82 47.70 47.88 49.04 49.94 53.88
1.46E-2 1.44E-2 1.46E-2 1.47E-2 1.47E-2 1.47E-2 1.46E-2 1.46E-2 1.19E-2 7.93E-3 3.99E-3 1.33E-3 1.77E-2 1.79E-2 1.78E-2 1.78E-2 6.50E-3 8.16E-4 1.39E-2 1.35E-2 1.11E-2 9.66E-3 2.92E-3
3.65E-2 3.66E-2 3.65E-2 3.67E-2 3.66E-2 3.67E-2 3.68E-2 3.66E-2 3.65E-2 3.64E-2 3.63E-2 3.65E-2 3.65E-2 3.66E-2 3.65E-2 3.66E-2 3.64E-2 3.63E-2 3.63E-2 3.64E-2 3.63E-2 3.62E-2 3.63E-2
2.38 2.32 2.11 2.05 2.01 1.92 1.79 1.67 1.69 1.81 1.86 1.87 2.29 2.17 2.09 2.01 2.22 2.26 2.38 2.31 2.25 2.25 2.31
Page 17 of 17
S0140-7007(15)00347-3 http://dx.doi.org/doi:10.1016/j.ijrefrig.2015.11.007 JIJR 3199
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
31-7-2015 5-11-2015 7-11-2015
Please cite this article as: J. Sierra Pallares, J. García del Valle, P. García Carrascal, F. Castro Ruiz, A computational study about the types of entropy generation in three different r134a ejector mixing chambers, International Journal of Refrigeration (2015), http://dx.doi.org/doi:10.1016/j.ijrefrig.2015.11.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A computational study about the types of entropy generation in three different R134a ejector mixing chambers Sierra Pallares, J.a,*, García del Valle, J.b, García Carrascal, P.a, Castro Ruiz, F.a a Departamento de Ingeniería Energética y Fluidomecánica - Grupo de Ingeniería de Fluidos Escuela de Ingenierías Industriales, Universidad de Valladolid Paseo del Cauce 59, 47010 Valladolid (Spain) b Área de Energía, Facultad de Ingeniería Civil y Mecánica Universidad Técnica de Ambato Avd. Los Chasquis y Rio Payamino s/n, Ambato (Ecuador) *Corresponding author Email addresses: [email protected] (Sierra Pallares, J.), [email protected] (García del Valle, J.), [email protected] (García Carrascal, P.), [email protected] (Castro Ruiz, F.) Highlights: • • • •
CFD based second law analysis for three ejector geometries. Division of the entropy generation into four basics machanims. Discussion of flow features based on CFD local entropy generation maps. Comparison of the compression mechanisms inside the ejector.
Abstract The purpose of this investigation is to study the efficiency of a R134a ejector operating with three different mixing chambers by means of a CFD based entropy generation analysis. With the aid of the differential equation for entropy, the local entropy generation is pursued. Using this method, the areas where the irreversibilities occur are identified and geometry improvements suggested. The novelty of the numerical procedure presented herein is that the bulk entropy generation is analysed by means four sources related with viscous dissipation and heat transfer, divided in mean and fluctuating terms. Entropy generation within the boundary layer has also been considered. The latter has proved to be small, as well as the heat transfer contribution. The fluctuating viscous dissipation accounts for more than 75% of the added entropy, being its main sources the shear layer after the nozzle exit and the shock wave trains, independently of their position. Keywords: Ejector, entropy, CFD, R-134a, shock wave Nomenclature A area [m2] M mach number p pressure [bar] s specific entropy [J·kg−1·K−1] sg volume rate of entropy generation [W·K−1·m−3] Sgen rate of entropy generation [W·K−1] T temperature [K] u velocity component in the x direction [m·s−1] v velocity component in the y direction [m·s−1] w velocity component in the z direction [m·s−1] q heat flux [W·m−2] v velocity vector [m·s−1] c control volume surface
Page 1 of 17
control volume Greek symbols: α thermal diffusivity [m2·s−1] ϵ dimensionless performance parameter ε turbulent dissipation rate [W·kg−1] κ conductivity [W·m−1·K−1] µ viscosity [Pa·s] ρ density [kg·m−3] σ entropy flux [W·m−2·K−1] τ stress tensor [J·m−3] ω entrainment ratio Subscripts: heat transfer term h m motive flow e entrained flow f flow cond condenser µ viscous dissipation term τ indicates friction velocity mean mean terms contribution to entropy fluc fluctuating terms contribution to entropy pro production rate Superscripts: ~ mean term (Favre averaged) – mean term ʺ fluctuating term T turbulent + indicates wall coordinate ʹ per unit area Acronyms: SKE Standard k − ϵ turbulence model RKE Realizable k − ϵ turbulence model RNG RNG k − ϵ turbulence model SST SST k − ω turbulence model AAD Average Absolute Deviation CFL Courant - Friedrichs - Lewy CFD Computational Fluid Dynamics 1. Introduction Supersonic ejectors have been used widely for more than a century. In refrigeration systems they may contribute to reducing energy consumption because they can be powered by solar energy or by waste heat (García del Valle et al., 2014). An ejector is a simple device in which a high pressure stream (primary fluid) is used to compress a low pressure stream (secondary fluid) to a higher pressure (discharge pressure). In jet-ejector cooling systems, a refrigerant stream of low velocity and low pressure is evaporated to generate a cooling capacity, and compressed to an adequate condensation pressure level by the interaction with a high velocity stream generated by heat input and adiabatic near-reversible expansion.
Page 2 of 17
The current trends in ejector research are focused in minimizing the ejector internal irreversibility, considering the flow deceleration the most important factor affecting entropy generation. In order to design and develop a high performance ejector, a clear understanding of the flow and mixing inside the ejector is first needed. Thus, modelling seems essential to afford this task. Although one dimensional models are useful for global state-space models, (Xue et al., 2015), it is commonly accepted that these models are not sufficient to understand the physics of a supersonic ejector (Bartosiewicz et al., 2005). For instance, (Sriveerakul et al., 2007) remarks that these models only can be used when the ejector is operated at its design conditions. Thus, the use of CFD codes is highly recommendable for the understanding of ejector flow physics, and in turn, once a figure of merit is defined, they may improve the ejector design, (Zhang and Mohamed, 2015). The figure of merit used for ejector geometry optimization, independently of its formulation, is ultimately related to the entropy generation minimization in the ejector. The pitfalls encountered in CFD computations are related with the selection of the turbulence model (Gagan et al., 2014), which makes a CFD validation mandatory in order to asses the uncertainties attained in the computed results. Entropy generation minimization is the method of modelling and optimization of real devices that owe their thermodynamic imperfection to heat transfer, mass transfer, and fluid flow irreversibilities (Bejan, 1996). In order to minimize entropy generation, its sources should be studied. Entropy production in turbulent flows can be determined using direct or indirect methods (Herwig and Kock, 2006). In the direct method, the differential entropy transport equation is employed to calculate the entropy production terms in the post-processing phase of a CFD calculation. That means they are determined by using the known field quantities (time-mean) velocity, temperature and composition. Using this method, different local parts of entropy production (due to viscous dissipation and heat and mass transfer) can be determined. Integration of these field quantities over the whole flow domain results in the overall entropy production rate. In the indirect method, entropy generation is computed using a control volume, in a black - box approach, with the integral entropy transport equation, (Kock and Herwig, 2004). Though this equation is exact, its application for practical problems involves several approximations, specially regarding the entropy flux terms across the inlet and outlet boundaries, normally neglected since unknown. Further approximations are involved when inlet and outlet quantities are assumed to have uniform values in the cross sections, neglecting their actual profiles. Some studies have addressed the problem of entropy generation in ejectors for refrigeration applications by means of the latter approach. (Arbel et al., 2003), using the indirect method, performed an analysis of the ejector as a whole and of each of its internal processes. They showed that entropy production is equivalent to performance losses, and that minimizing entropy production could be used as a criterion for optimization. They also indicated that the ejector irreversibility is due to three main factors: “pure mixing”, “kinetic energy losses”, and the normal shock wave produced in the mixing chamber. (Liu and Groll, 2013) presented a method to determine the internal-ejector component efficiencies by combining zero-dimensional modelling with measurements performed for a R744 ejector refrigeration cycle. Measured and simulation data were used to determine the isentropic efficiencies of the motive nozzle, the suction nozzles, and the mixing section efficiency. On the CFD side, (Banasiak et al., 2014) investigated about the performance of R744 ejectors. In that work, a new approach was introduced to assess the contribution of the local irreversibilities to the overall entropy increase
Page 3 of 17
inside the ejector. This analysis was based on close examination of the flow variable profiles, as well as on the energy and the entropy balances formulated for individual ejector passages using the indirect method. This approac enabled the precise location of intensified irreversibility sites and delivered detailed information on the irreversibilities origin. To the best of the author knowledge, all of the studies described thus far have been limited to the application of the indirect method. Therefore, our aim in this work is to apply the direct method to the analysis of different ejector chambers for refrigeration applications. It is clear that the direct method is superior in the analysis of complex flow applications, such as ejectors: from the direct method we obtain information about the total entropy production and how this production is distributed between the single mechanisms conforming the flow pattern; information the indirect method cannot provide. Therefore, the present study extends the scope of the existing literature by providing relations to compute the entropy generation using CFD computations combining the ideas of other authors (Bejan, 1982; Kock and Herwig, 2004) and allowing us to discriminate between the different entropy production types in an easy and accurate way. In this way, we shall provide information that cannot be obtained using the indirect method, thus adding value to the existing literature on the topic. The rest of the paper is developed in the following way: The computational model is first described, along with detailed information about the computation of the entropy production rate in ejectors both direct and indirect methods. Secondly, the three mixing chambers under study will be detailed geometrically, along with the boundary conditions employed in the analysis. The last Section is related to the comparison between experimental results and CFD simulations and a discussion regarding the entropy generation rate in the three chambers under study. It is important to note that the model validation was already addressed in a previous work, (García del Valle et al., 2015). A short account of the results found will be included in the description of the computational model. 2. Computational model 2.1. Fluid dynamic and thermodynamic models. Numerical implementation The commercial software ANSYS Fluent 14.5 has been used to compute the compressible steady-state, 2D axis-symmetric turbulent flow governing equations which apply inside the ejector. This software is based on the finite volume approach, a technique that discretizes all governing equations by first dividing the physical space into a number of control volume grid elements. When dealing with turbulent flow, a suitable averaging method must be used. The density averaging technique or Favre averaging (Versteeg and Malalasekera, 1996) is appropriate for compressible flows and will be used in this work. The total energy equation is included and coupled to the mass and momentum equations and the thermodynamic model. Turbulence modelling is based on the Standard k - ϵ model (Launder and Spalding, 1972b). Real gas thermodynamic models must be considered due to the departure from ideal gas found in the region of superheated vapour, close to the saturation line, a region in which the ejector commonly operates. The REFPROP library, (Lemmon et al., 2002), integrated into ANSYS Fluent, has been used as the source of thermodynamic data. For the case of R-134a, this formulation is based in the multiparameter equation of state of Tillner-Roth and Baer (Tillner-Roth and Baehr, 1994). This equation constitutes a 21 parameter model that represents nearly all selected experimental data within their estimated accuracy with the exception of heat capacities and speed of sound in the region close to the critical point. Typical accuracies are 0.05% for density, 0.02% for the vapor pressure, or 0.5 and 1% for the heat capacity. A second order
Page 4 of 17
upwind scheme was selected for the advective terms of each equation except for the pressure, where a Least Squares Cell Based approach is chosen. High order term relaxation was also applied. The density - based algorithm was selected for the pressure - velocity coupling problem. A strict Courant-Friedrichs-Lewi (CFL) condition is needed here in order to prevent the solution from diverging during the first integration steps. This is due mainly to the thermodynamic state falling outside of the vapor (both stable and metastable) regions when using the REFPROP library. A CFL of 0.8 is used for the first 500 steps, and then gradually increased to reach a value of 2.5. 2.2. Computation of the entropy generation rate 2.2.1. Direct method The entropy in a flow field is a state variable that in its specific form has a balance equation that for a single phase, non-reacting and single component flow reads (Sciacovelli et al., 2015):
1 1 q = 2 q T τ : v Dt T T T
Ds
(1)
which can recast into
Ds
= σ s g
Dt
(2) Where σ is the entropy flux vector and sg is the entropy generation rate per unit volume. A comparison between (1) and (2) leads to σ =
q T
(3) 1
sg =
q T
2
T
1
τ : v
T
(4) It can be observed that the total entropy generation rate can be splitted in two contributions, each one strictly related to a specific transport phenomenon: a contribution due to heat transfer (6) and a viscous dissipation term (7): sg = sg
sg
h
(5)
where sg sg
h
=
1 T
=
2
q T
1
(6) τ : v
T
(7) In a turbulent flow, Equation (1) has to be time averaged in order to find an equation for the mean entropy. For a compressible flow, Favre averaged is employed to obtain
u s v s w s = sg Dt y z x
Ds
h
sg
(8) Unclosed terms are modelled using the method proposed by Herwig and Kock (Herwig and Kock, 2006): the entropy generation rate is divided between the contribution of the mean gradients and the contribution of the fluctuating components:
Page 5 of 17
sg
= sg
h
sg
h
= sg
sg
mean
h
sg
(9)
fluc
(10) where the contibutions to the entropy generation rate due to mean gradients are: sg
= h
m ean
2 T
T
m ean
fluc
2 2 T T y z x
2
u
sg
m ean
2
2
v w 2 = z T x y
2
(11) u v 2 y x
v w u w z x z y
(12)
2
2
and the fluctuating ones are: sg
=
h
sg
fluc
T
sg
h
(13)
mean
= ~ T
(14) Although equations (11) - (14) can be used in every point of the computational domain, special care is needed when using wall functions in the computation of the turbulent flow. In this work, standard wall functions are used and in order to obtain the correct estimation of entropy generation in the cells adjacent to the walls, the method proposed by Walsh and McEligot (Walsh and McEligot, 2008) is used. Under boundary layer approximations for an un-heated two-dimensional flow, the pointwise entropy generation rate may be approximated (using wall coordinates) as:
s
fluc
2
g
u uv y
u y
(15)
where
uv
u 1 y
(16) The value of pointwise entropy is desired to identify regions where most losses occur and to deduce the integrated entropy generation per unit surface area
s y
g
=
T s g
u
3
s y
=
0
g
dy
=
y
0
u y
s y
g
. It can be written as
dy
(17) In this way, the total entropy generation is easily obtained given the wall function employed and the entropy generation per unit area due to the wall function computed as:
Page 6 of 17
u 3 y u s g = 0 y T
dy
(18) Thus, given the area of the computational cell and the value of y+, the total entropy generated in the boundary layer can be computed. In a standard computation, the continuity, momentum and energy equations are first solved. Then, in a postprocessing step, Equations (11), (12), (13), (14) are computed at cell centers except in the first row of cells adjacent to the wall, where Eq. (18) is applied as commented. All of these apportations to the total entropy generation rate constitutes the entropy production rate per unit volume, that can be plotted as map over the computational mesh. 2.2.2. Indirect method From the entropy transport equation (2) is possible to derive its integral form, which constitute the basis of the indirect method (please see the appendix for a detailed derivation): d
dt
s d c
c
q n d c c T
s v n d c =
sg d f
(19) Recasting the last equation, it is possible to obtain the entropy production rate S p ro for a given control volume, S pro =
sg d = f
d dt
s d c
s v n d c c
q c T
n d c
(20) Therefore, the local rate of entropy production can be computed by determining the right hand side of (20). As explained by Kock and Herwig (Kock and Herwig, 2004), this method has the drawbak of the necessary computation of the integral of the entropy flux term:
q n d c c T
(21) which is always neglected in the literature, or even unmentioned (Arbel et al., 2003; Banasiak et al., 2014). However, the relative value of this term is of importance in some situations. If the control volume selected cuts down the geometry under study in a region of big temperature gradients, the error produced in the computation of the entropy production may be significant due to the high value of the heat flux accross the control surface. Validation of the CFD model 2.3. In a previous study, (García del Valle et al., 2015), the CFD tool was subjected to a four step validation procedure, related with the real gas implementation, turbulence modelling, the shock wave-boundary layer interaction and the ejector simulation itself using different validation cases extracted from the literature. The purpose of the validation cases was to independently address the different flow features which appear in the ejector. The general results were as follows: 1. The real gas implementation was compared with the shock tube (Riemman problem) solution for a R-134a vapour flow. Albeit some diffusion was found in the shock and expansion wave, their location and the jump in thermophysical properties were correctly predicted. 2. Four turbulence models, according to Table 1, were considered. Experimental data of a supersonic air jet expanding into quiescent air, (Eggers, 1966), was used as a test case for this purpose. It was found that, although the supersonic region of the jet was accurately
Page 7 of 17
predicted, the numerical results in the subsonic region exhibit higher diffusive rates than the experimental case. The SKE and SST model had the lower overall error. 3. The shock wave-boundary layer interaction was compared against the experimental work of (Simeonides, 1994), who studied the Mach 5 air flow over a hollow cylinder flare. Both the pressure bump before the shock wave and the reattachment point were predict with accuracy with the SKE and RKE, while the RNG and SST were poorer in performance. 4. For the ejector validation, both the entrainment ratio and the pressure profiles along the wall were considered for four cases, one operating in the mixed regime and the rest in the double choked regime. The ejector geometries are the same as those studied in this investigation. The results are best understood if divided into the mixed regime and the double choked regime. For the former, a discharge pressure dependent value for the entrainment ratio was found both experimentally and numerically. The average value over a time span was used to compute the entrainment ratio. It was found that, due to the nature of the flow, where the turbulence modelling is of paramount relevance, the entrainment ratio was poorly predicted by all turbulence models. The pressure profile, however was accurately predicted. Further numerical studies shown that the pressure profile is decoupled from the entrainment ratio for low values of the latter, which explain the incongruence describe previously. For the double choked regime, both the SKE and SST model yielded the best results, with an overall error of 15% for the entrainment ratio and 5% for the pressure profile. Based on the validation phase, the SKE turbulence model will be used in the present investigation. 3. Ejector geometries under study The R-134a refrigerant ejector was exhaustively described in (García del Valle et al., 2014). The reader is referred to the original source, as only a short summary will be included here focused on the experimental uncertainties. The experimental installation, rated at 10.5 kW of boiler and evaporator power, was designed to accurately control and measure the mass flow, pressure and degree of superheating of the primary and secondary mass flow and the discharge pressure and temperature of a number of ejector geometries. For this work, the results from three geometries as shown in Fig. 1 have been employed. The main characteristics of each geometry are as follows: 1. Chamber “A” is the classical shape with a 30° half-angle conical entry and a constant area zone with a length to diameter ratio of 8.5. 2. Chamber “B”, promotes a supersonic compression in a convergent region rather than a shock system in a constant area region. That is why the present design has a 1.14° total conical entrance and a short constant area region. 3. Chamber “C”, pursues a supersonic compression through a sharp contraction. This geometry has two constant area regions with diameters of 5.5 mm and 4.8 mm respectively. Each geometry is tested with the same primary and secondary conditions, only the discharge pressure is modified in each run. The discharge pressure values range from 6.8 bar up to the point in which the entrainment ratio tends to zero. It is worth mentioning that operation close to zero mass entrainment is unstable, although the secondary pressure remains approximately stationary due to the mass and thermal inertia of the systems, oscillations were found for the secondary mass flow rate. Furthermore, due to the poor regulation both of the
Page 8 of 17
resistor and of the electronic expansion valve at low mass flow rates, slight deviations from the set point are found. All the relevant experimental data, in which the numerical investigation is based, are included in Table 3. Further columns have been added and will be used later on. The nozzle exit is positioned at -5.38 mm in all measurements. The uncertainties of the different measured variables are shown in Table 2. The following comments are needed in order to fully understand this table. The column “Range” represents the calibration range of the sensor. The term “Dev. FS” represents the deviation over full scale. It should be noted that two terms are provided, the first is attributed to the sensor and the second to the data acquisition system. The column “Dev. abs. max” is the maximum absolute deviation, computed as (“Max Range”-“Min Range”)*“Dev. FS”. An uncertainty analysis, performed according to (Moffat, 1988), reveals that the entrainment ratio error is equal to ± 0.0061. Each point corresponds to the average value over a time span of 5 minutes (which corresponds to 60 logged data points, since the data acquisition rate is equal to 0.2 Hz), in which the deviation of selected measured variables is within a predefined limit to ensure the steady state condition. 3.1. Computational details The SKE turbulence model is employed to predict the performance of the ejector using the different mixing chambers described. The flow was considered in steady state, 2D axis-symmetric and compressible. Second order upwind schemes are used for the flow variables except for pressure, where a Least Squares Cell Based approach is employed. As an example, the computational domain along with boundary conditions for Chamber B is sketched in Figure 2. Pressure inlets are modelled using the values of static pressure and temperature as total pressure and temperature values respectively, due to the low value of the velocity in these inlets in relation to the velocity in the chamber. Turbulence intensity is 5% in each inlet, and the other boundary conditions (pressure and temperature at inlets and discharge pressure) are specified in Table 3. Different grid sizes have been employed for the computational domain, looking for mesh independence, ranging from 40000 to 150000 cells. The solution showed independence form the mesh using 50000, 90000 and 60000 cells for chambers A, B and C respectively. Mesh adaptation has been carried out based on Mach number gradients and pressure gradients. CFL condition is fixed to 1 for the first hundreds of iterations, and then progressively increased to a value of 2.5 to prevent the solution from diverging during the first integration steps. This is due mainly to the thermodynamic state falling outside of the vapor (both stable and metastable) regions when using the REFPROP library. The mean y+ value for all the chambers under study is around 125. Results 4. The numerical results are structured in two subsections. First, the entrainment ratio and the overall entropy generation rate are addressed. Experimental and numerical results will be compared at this point. Second, contour plots charts for the local entropy generation rate will be reported. Based on these results, a comparison of the three geometries will be addressed. 4.1. Prediction of the entrainment ratio. Total entropy production This section is useful for validation purposes, since both the entrainment ratio and entropy generation is compared against experimental values. The former is shown in Figure 3 whilst the latter is included in Figure 4. Fig. 3 shows the comparison between the experimental and numerical entrainment ratio for the three mixing chambers under study for a range of discharge pressures. A non dimensional performance parameter ϵ is used in this section to discuss the performance of the different
Page 9 of 17
geometries. This parameter was introduced in (García del Valle et al., 2014) and is defined as the quotient between the entrainment ratio and the entrainment ratio of a isentropic ejector with the same primary and secondary flows stagnation conditions and discharge pressure. Therefore, the parameter ϵ is inversely related with the entropy generation. Its advantage over the latter is that its value is bounded between 0 and 1, and its meaning its similar to that of the isentropic efficiency. Numerical results compare better to experimental data for mixing chamber “A”, where the double choked regime is accurately predicted by the model. For chamber “B”, a mean deviation of 14% is found, close to the previously reported value in (García del Valle et al., 2015). Mixing chamber “C” operates in the mixed regime in the selected discharge pressure range. While the critical pressure is successfully computed for the mixing chamber “B”, for mixing chamber “A” it is overestimated. Fig. 4 shows the entropy production in each chamber computed using the indirect method, Eq. (20). It is worth to mention here that to avoid the computation of the entropy flux term, Eq. (21) the control surface is located far away from regions were temperature gradients are found. In this way, it is safe to neglect it in the computation of the entropy production. From Figure 3 and 4 the aforementioned inverse relation between the entropy generated and the performance parameter ϵ is clearly seen. The peak in the performance parameter match the minimum in entropy generation. The experimental values of entropy generation (see last column of Table 3) are lower than those obtained numerically. For instance, the minimum numerical values are 2.3 W·K−1 and 2.65 W·K−1, against the minimum experimental values of 1.67 W·K−1 and 2.1 W·K−1 for mixing chamber “A” and “B” respectively. Although it is a deviation close to 40%, is worth noting the strong correlation which exists between the variation of the entropy generation with the discharge temperature. Simple computations with the control volume approach yield a relation of 0.2 W·K−1 per K of discharge temperature increment. So the variation of numerical to experimental discharge temperature is in the range of 3 to 4 K. Fig. 5 reports the comparison of the direct and the indirect methods in computing the entropy production of the ejector. When using the direct method, the entropy generation flowfield is first computed and then integrated for the whole computational domain in order to obtain the entropy production for the ejector. Theoretically, both methods should give the same results of total entropy generation, however, according to Fig. 5, a deviation within ± 10% of the latter over the former is found. This difference can be attributed the poor convergence of some cases under study, specially in the mixed regime, where a non stady flow is found (the procedure used to handle these cases is reported in (García del Valle et al., 2015)). 4.2. Entropy generation analysis In this section, regions of the ejector where entropy generation occurs are identified, along with the mechanism which produce every type of irreversibility. Figs. 6 shows the contour plots of Mach number, mixture fraction and entropy generation for the mixing chamber “A”. The mixture fraction is computed as a passive scalar, with boundary conditions equal to one for the primary flow and 0 for the secondary flow, and provides a measure of the extent of mixing at the macro-scale. From the Mach number distribution, the geometry “A” operates with supersonic flow up to the end of the constant area region, where a shock wave train adjust the pressure of the mixture to the discharge pressure. Mixture fraction maps reveals that the primary and secondary flows do not completely mix over the constant area region, but rather the mixing
Page 10 of 17
is finished through the train of shock waves, were high turbulence is found. Provided the moderate rate of mixing in the constant area region shown in Fig. 6, the Mach number remains approximately constant through this region. Due to the increment of the shear layer between the primary and secondary flows through the constant area region, the rate of entropy generation slightly falls along this region. As expected, the irreversibility is concentrated within the shock wave train. Fig. 7 depicts the entropy production per unit length, decomposed in its different sources. This figure is produced integrating the entropy generation map at selected surfaces where the axial distance is constant along the computational domain. The sharp peak depicted in the figure shows how the main sources of irreversibility are concentrated in the shock train region. Besides, the rate of entropy generation within the boundary layer is small compared with the bulk entropy generation. This fact is repeated for all geometries. The long tapered mixing chamber introduced in geometry “B”, changes the flow pattern described for geometry “A” (See Fig. 8). With this geometry, a compression occurs through the convergent region, which translates into a decrement of the Mach number. Furthermore, the mixture between the primary and secondary flow is almost completed at the entrance of the constant area region. In this short region, the flow is transonic. Through the diffuser the mixture decelerates up to the discharge pressure. Fig. 9 shows the entropy production per unit length for mixing chamber “B”. Through the long tapered region, both the compression and the shear layer contribute to the entropy generation. But, since the velocity gradient is smaller than for geometry “A”, the entropy generation is slightly lower. As for the former case, a peak of irreversibility is found at the end of the constant area region, though this time is smaller, mainly because the bulk velocity is lowered. The geometry “C” (see Fig. 10) resembles the operation of an unadapted two throat system. A shock wave occurs at the nozzle exits in order to obtain a subsonic flow at the entrance of the “secondary” throat, i.e. the constant area region, in which the mixture accelerates again up to sonic flow. All the entropy generation is thus concentrated immediately after the nozzle exit, as is clearly visible from Fig. 11. In order to analyze the different entropy generation mechanisms, as described in Section 2.2, the Fig. 12 has been included. It is a multiple pie chart in which the different sources of irreversibility are identified for the different geometries and ejector regions. One of the superior aspects of the direct method over the indirect one is the capability of determining the source of irreversibilities. From Fig. 12, it can be seen that the main entropy generation mechanism is viscous dissipation for all the chambers under study. Irreversibilities induced by heat transfer (due to mean gradients and fluctuating ones) merely account for less than 2% of the total amount. Furthermore, fluctuating viscous dissipation is the main entropy source, accounting for a share of more than 75%. This figure is obtained using the SKE model, but can be reproduced using other RANS approaches. For example, using the SST model, small differences are found for the cases under study. If we pay attention to the entropy generation per zone, we observe that in the region previous to the mixing chamber (what we call “premixing” region) entropy generation due to the mean dissipation Eq. (12) is predominant. This is because the main source of this type of dissipation is located close to the walls, completely in the boundary layer of the supersonic flow inside the nozzle. There, the velocity gradients are very important and contribute the most to the
Page 11 of 17
entropy production. In this region, fluctuating dissipation is unimportant because the mixing process has not even started. This effect is common to all geometries. However, in the mixing chamber, we observe different behaviours for each geometry. In this region, most of the dissipation is due to the fluctuating one. In mixing chamber “A”, most of the entropy generated due to the fluctuating dissipation is located in the diffuser, but in the other geometries the behaviour is the opposite. This is connected with the location of the normal shocks and the mixing processes that occur inside the ejector. Depending on the location of the shock and the mixing intensity, this type of entropy generation is placed mostly in a region or another. For mixing chamber “A”, the wake of the normal shock is responsible of most of the irreversibility generated in the diffuser (see the entropy generation map of Fig. 6) but since mixing is relatively slow, the mixing chamber do not present a significant entropy generation of this type due to mixing. For mixing chambers “B” and “C”, due to the increased mixing rate, most of this irreversiblity is located in the mixing chamber itself, although the shock is located downwards the constant area region in the “B” geometry. For mixing chamber “C”, both mixing and normal shock happen prior to the reduction of diameter, and that is why up to 88% of the fluctuating dissipation is located backwards. Regarding the entropy generation due to mean and fluctuating temperature gradients (i.e. heat transport inside the ejector), the comments of the above paragraph hold. The temperature difference of the primary and secondary flow is dissipated by means of mixing, and there is a strong temperature gradient associated to the normal shock. That is why this source of irreversibility more or less matches with the entropy generation due to fluctuating dissipation. 4.3. Design optimization by means of entropy generation minimization The end point of this investigation should provide some guidelines to the design of high performance mixing chambers. The entropy generation should be taken with caution, because, although for an operation condition one configuration can exhibit the lower entropy generation rate, over different conditions this is not true. Furthermore, although the entropy is closely related with the performance parameter ϵ, it does not mean that this relation is inversely linear. From Fig. 4, mixing chamber “B” perform better in the lower discharge pressure range against mixing chamber “A”, which does better in the upper discharge range. Not only that but, the differences in the minimum entropy generation 2.3 and 2.65 W·K−1 for mixing chambers “A” and “B” respectively are not directly translated to the performance parameter ϵ, which is almost similar from one to another. The key point in the discussion is to analyze why the geometry “B” does not improve the performance over the mixing chamber “A”. Our vision is somehow skewed by the description given earlier for the double choked regime, i.e. that the entropy generation after the constant are region is smaller for geometry “B” (softer shock wave train) and that the irreversibility lost into the shear layer is smaller or of the same order as that of geometry “A”. However, as the critical point is reached, the entropy generation in the sock wave train in geometry “A” is decreased while the entropy generated in the shear layer will remain approximately constant. At the critical point, the pressure increment in mixing chamber “A” occurs solely through the sock wave train, for mixing chamber “B” the compression occurs through the convergent region and through a softer sock wave train. From both the experimental and numerical data, the performance of both procedures is similar, although the mechanism described for geometry “B” it is “a priori” more favourable. However, there is an additional element missed in the former description.
Page 12 of 17
As for the efficiency of a subsonic diffuser (divergent area) is negatively influenced by low boundary layer velocities and boundary layer detachment, the efficiency of a supersonic diffuser (convergent area) is also negatively influenced by the low boundary layer velocities created by the secondary fluid. Coupled with the negative pressure gradient, the thickness of the boundary layer should grow way more rapidly than for geometry “A”, as is clearly visible from the mixture fraction of Figs. 6 and 8. The authors suggest that this behaviour is the cause that negligible performance increments are found for long tapered ejectors in relation with long constant are one. Thus, the possible gain induced by geometry variations in the peak ejector performance parameter ϵ is marginal. For an specific application and operating conditions, the critical feature is to get a standard geometry (say “A” or “B”) under which the operating point is closed to the critical point, this is entropy generation minimization. 5. Conclusions The mechanisms of entropy generation have been studied in this investigation based on experimental and computational results for a R134a ejector, covering three different geometries and a number of operating conditions. Using the differential balance for entropy, the four mechanisms relevant for entropy production are identified inside the ejector chambers under analysis. Between all the mechanisms analyzed, the fluctuating viscous dissipation, which bears a share of more than 75% of the total entropy generation for all geometries, is the single most important dissipative effect which drags ejector performance. Furthermore, a special treatment has been applied for dealing with the irreversibilities inside the boundary layer. Not only the total entropy generation has been identified, but rather the volumetric rate of entropy generation through the domain has been pursued. Based on these results, and for the sizes of cross-sectional areas studied herein, it is has been concluded that the efficiency of the compression through a train of shock waves is similar to that occurred in a long convergent conduit. The authors suggest that limited improvement may be gain upon geometry variation of cases “A” and “B” because the physical behaviour will keep being the same as in those geometries. An operation close to the critical point may provide better results in this regard. Acknowledgements The authors gratefully acknowledge the financial support given to the investigation by the “Junta de Castilla y León” under grant 04/09/VA/0025 and “Secaderos y Frigoríficos Industriales S.L.”. References Arbel, a., Shklyar, a., Hershgal, D., Barak, M., Sokolov, M., 2003. Ejector Irreversibility Characteristics. Journal of Fluids Engineering 125 (1), 121. Banasiak, K., Palacz, M., Hafner, A., Buliski, Z., Smoka, J., Nowak, A. J., Fic, A., Apr. 2014. A CFD-based investigation of the energy performance of two-phase R744 ejectors to recover the expansion work in refrigeration systems: Anirreversibility analysis. International Journal of Refrigeration 40, 328–337. Bartosiewicz, Y., Aidoun, Z., Desevaux, P., Mercadier, Y., Feb. 2005. Numerical and experimental investigations on supersonic ejectors. International Journal of Heat and Fluid Flow 26 (1), 56–70. Bejan, A., 1982. Entropy generation through heat and fluid flow. John Wiley & Sons Inc. Bejan, A., 1996. Entropy generation minimization: The new thermodynamics of finite-size devices and finite-time processes. Journal of Applied Physics 79 (3), 1191.
Page 13 of 17
Eggers, J., 1966. Velocity profiles and eddy viscosity distributions downstream of a mach 2.22 nozzle exhausting to quiescent air. NASA TN D-3601. Gagan, J., Smierciew, K., Butrymowicz, D., Karwacki, J., Apr. 2014. Comparative study of turbulence models in application to gas ejectors. International Journal of Thermal Sciences 78, 9–15. García del Valle, J., Saíz Jabardo, J., Castro Ruiz, F., San José Alonso, J., Oct. 2014. An experimental investigation of a R-134a ejector refrigeration system. International Journal of Refrigeration 46, 105–113. García del Valle, J., Sierra Pallares, J., García Carrascal, P., Castro Ruiz, F., 2015. An experimental and computational study of the flow pattern in a refrigerant ejector. Validation of turbulence models and real-gas effects. Applied Thermal Engineering 00, 000–000. Herwig, H., Kock, F., Mar. 2006. Direct and indirect methods of calculating entropy generation rates in turbulent convective heat transfer problems. Heat and Mass Transfer 43 (3), 207–215. Kock, F., Herwig, H., May 2004. Local entropy production in turbulent shear flows: a high-Reynolds number model with wall functions. International Journal of Heat and Mass Transfer 47 (10-11), 2205–2215. Launder, B., Spalding, D., 1972a. Lectures in Mathematical Models of Turbulence. Academic Press. Launder, B. E., Spalding, D. B., 1972b. Lectures in Mathematical Models of Turbulence. Academic Press. Lemmon, E. W., McLinden, M. O., Huber, M. L., 2002. REFPROP: Reference Fluid Thermodynamic and Transport Properties. NIST Standard Database 23. Liu, F., Groll, E. A., Apr. 2013. Study of ejector efficiencies in refrigeration cycles. Applied Thermal Engineering 52 (2), 360–370. Menter, F., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal 32, 1598–1605. Moffat, R., 1988. Describing the uncertainties in experimental results. Experimental Thermal and Fluid Science 1, 3–17. Sciacovelli, a., Verda, V., Sciubba, E., 2015. Entropy generation analysis as a design toolA review. Renewable and Sustainable Energy Reviews 43, 1167–1181. Shih, T., Liou, W., Shabbir, A., Yang, Z., 1995. A new k − ϵ eddy viscosity model for high reynolds. Computers and Fluids 24, 227. Simeonides, G., 1994. Test cases for CFD code validation for inlets and nozzles of airbreathing propulsion systems. Doc. EWP/1784, ESA/ESTEC. Sriveerakul, T., Aphornratana, S., Chunnanond, K., Aug. 2007. Performance prediction of steam ejector using computational fluid dynamics: Part 2. Flow structure of a steam ejector influenced by operating pressures and geometries. International Journal of Thermal Sciences 46 (8), 823–833. Tillner-Roth, R., Baehr, H. D., 1994. An international standard formulation for the thermodynamic properties of 1, 1, 1, 2-tetrafluoroethane (hfc-134a) for temperatures from 170 k to 455 k and pressures up to 70 mpa. Journal of Physical and Chemical Reference Data 23 (5), 657–729. Versteeg, H. K., Malalasekera, W., 1996. An introduction to computational fluid dynamics. Longman. Walsh, E. J., McEligot, D., 2008. Relation of entropy generation to wall laws for turbulent flows. International Journal of Computational Fluid Dynamics 22 (March 2015), 649–657.
Page 14 of 17
Xue, B., Cai, W., Wang, X., Jan. 2015. State-space modelling for the ejector-based refrigeration system driven by low grade energy. Applied Thermal Engineering 75, 430–444. Yakhot, V., Orszag, S., 1986. Renormalization group analysis of turbulence. i. basic theory. Journal of Scientific Computing 1, 3–51. Zhang, T., Mohamed, S., 2015. Conceptual Design and Analysis of Hydrocarbon-Based Solar Thermal Power and Ejector Cooling Systems in Hot Climates. Journal of Solar energy engineering 137 (April), 1–9. Appendix A. Derivation of the integral transport equation for entropy From the entropy transport equation (2) is possible to derive its integral form, which constitute the basis of the indirect method:
Ds
d
=
Dt
f
σd
f
sg d
(A.1)
f
since
Ds
d
d
=
Dt
f
dt
sd
(A.2)
f
and taken into account the transport theorem, we obtain a transport equation for a control volume d
dt
sd
d
=
dt
f
sd
c
s v nd
(A.3)
c
c
therefore, d dt
sd
c
s v nd
c
=
c
σd c
sg d
(A.4)
c
Applying the divergence theorem, the volume integral of the entropy flux σ can be written
σd
=
c
c
q nd T
c
(A.5)
Introducing the above equation in (25) one obtains d
dt
sd c
s v nd c
c
q = nd c T
c
sg d
(A.6)
f
Figure 1: Ejector geometries under study. Note the origin of the coordinate system for each of the three cases. Figure 2: Sketch of the computational domain and details of the meshing in selected areas for the ejector geometry under study (Chamber B) Figure 3: Comparison between the experimental and numerical entrainment ratio for the three mixing chambers. The performance parameter ϵ, based on CFD data, has also been included. Figure 4: Numerical entropy generation for the three mixing chambers. The performance parameter ϵ, based on CFD data, has also been included. Figure 5: Comparison between the “direct” and “indirect” entropy generation numerical procedures.
Page 15 of 17
Figure 6: Numerical Mach number, mixture fraction and rate of entropy generation for the mixing chamber “A”, turbulence model SKE and conditions. Case A-2 Figure 7: Rate of entropy generation as a function of the axial coordinate for mixing chamber “A”. Case A-2 Figure 8: Numerical Mach number, mixture ratio and rate of entropy generation for the mixing chamber “B”, turbulence model SKE and conditions. Case B-2 Figure 9: Rate of entropy generation as a function of the axial coordinate for mixing chamber “B”. Case B-2 Figure 10: Numerical Mach number, mixture ratio and rate of entropy generation for the mixing chamber “C”, turbulence model SKE and conditions. Case C-2 Figure 11: Rate of entropy generation as a function of the axial coordinate for mixing chamber “C”. Case C-2 Figure 12: Multiple pie chart in which the entropy generating mechanisms as well as its share in the different regions is reported for all mixing chambers. Table 1: Turbulence models considered during tha validation stage. Model Standard k − ϵ
Realizable k−ϵ RNG k − ϵ SST k − ω
Reference (Launder and Spalding, 1972a) (Shih et al., 1995) (Yakhot and Orszag, 1986) (Menter, 1994)
Main application Fully turbulent flows
Acronym used in this work SKE
Planar and round jets, rotation and separated flows High streamline curvature flows Adverse pressure gradient flows
RKE RNG SST
Table 2: Measurement uncertainties for the relevant physical properties. Magnitude Position Range Dev. FS Dev. abs. max. Pressure Primary 0-40 bar ±(0.5%+0.25%) ±0.30 bar fluid inlet Secondary 0-30 bar ±(1%+0.25%) ±0,375 bar fluid inlet Discharge Temperature Primary 0-150 °C ±(0.1%+0.25%) ±0.52 °C fluid inlet Secondary -50-50 °C ±0.35 °C fluid inlet
Page 16 of 17
Mass flow
Discharge Primary Secondary
0-100 °C 0-0.125 kg·s−1
±(0.18%+0.25%)
±0.35 °C ±5.37·10−4 kg·s−1
Table 3: Relevant measured physical properties for each mixing chamber, along with the entropy produced. Entropy generation is computed here from experimental measurements “A ”
“B”
“C”
Cas e 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 1 2 3 4 5
p0mf[bar] T0mf[°C] p0ef[bar] 29.15 29.17 29.08 29.16 29.08 29.27 29.19 29.15 29.05 29.09 29.18 29.17 29.09 29.12 29.04 29.16 29.10 29.12 28.91 29.01 29.01 28.98 28.95
95.10 94.95 94.93 95.09 94.97 95.19 95.00 95.01 94.95 95.01 95.04 95.03 94.98 95.08 94.91 95.11 95.00 95.09 94.94 95.01 94.99 95.05 94.95
4.14 4.12 4.12 4.14 4.14 4.14 4.14 4.14 4.13 4.17 4.20 4.24 4.11 4.13 4.11 4.12 4.15 4.34 4.11 4.14 4.15 4.17 4.14
T0ef[°C ] 19.96 20.00 20.00 20.02 20.00 20.01 20.01 20.03 20.02 20.26 21.97 24.04 20.00 20.00 20.00 20.01 20.86 23.47 20.10 20.01 20.07 20.22 23.34
pcond[bar ] 6.82 7.00 7.17 7.38 7.58 7.81 7.99 8.19 8.39 8.58 8.81 9.00 6.79 7.03 7.17 7.36 7.80 8.04 6.83 7.00 7.19 7.29 7.59
Tcond[°C] me[kg·s−1] mm[kg·s−1] Sgen[W·K−1] 46.96 47.39 47.11 47.66 48.45 48.85 48.83 49.14 50.68 52.86 54.95 56.32 45.83 46.33 46.58 47.01 52.73 55.82 47.70 47.88 49.04 49.94 53.88
1.46E-2 1.44E-2 1.46E-2 1.47E-2 1.47E-2 1.47E-2 1.46E-2 1.46E-2 1.19E-2 7.93E-3 3.99E-3 1.33E-3 1.77E-2 1.79E-2 1.78E-2 1.78E-2 6.50E-3 8.16E-4 1.39E-2 1.35E-2 1.11E-2 9.66E-3 2.92E-3
3.65E-2 3.66E-2 3.65E-2 3.67E-2 3.66E-2 3.67E-2 3.68E-2 3.66E-2 3.65E-2 3.64E-2 3.63E-2 3.65E-2 3.65E-2 3.66E-2 3.65E-2 3.66E-2 3.64E-2 3.63E-2 3.63E-2 3.64E-2 3.63E-2 3.62E-2 3.63E-2
2.38 2.32 2.11 2.05 2.01 1.92 1.79 1.67 1.69 1.81 1.86 1.87 2.29 2.17 2.09 2.01 2.22 2.26 2.38 2.31 2.25 2.25 2.31
Page 17 of 17