- Email: [email protected]

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Turbulent dense gas ﬂow characteristics in swirling conical diffuser C.S. From a, E. Sauret a,∗, S.W. Armﬁeld b, S.C. Saha a, Y.T. Gu a a

School of Chemistry, Physics and Mechanical Engineering, Science and Engineering Faculty Queensland University of Technology, Brisbane, QLD 4000, Australia b School of Aerospace, Mechanical and Mechatronic Engineering, University of Sydney, NSW 2006, Australia

a r t i c l e

i n f o

Article history: Received 8 August 2016 Revised 14 March 2017 Accepted 19 March 2017 Available online 20 March 2017 Keywords: Turbulence Swirling ﬂows Conical diffuser Dense gas Explicit algebraic Reynolds stress model

a b s t r a c t Diffusers placed at the exit of turbines are essential to recover pressure and increase turbine eﬃciency. This increase of eﬃciency is critical for the overall cycle eﬃciency of renewable power cycles based on low temperature renewable resources. Optimising the performance of a conical diffuser in renewable power cycles using high-density ﬂuids can be established by examining the turbulence characteristics of both air considered as an ideal gas (IG) and R143a, a refrigerant with high-density in a non-ideal state, considered as a real gas (RG). Turbulence was ﬁrstly modelled and validated against experimental results from the ERCOFTAC swirling conical diffuser database and previous numerical results. The real gas thermodynamic and transport properties of refrigerant R143a were then obtained from the NIST REFPROP database. Investigating both RG and IG revealed that general trends remain, where the stronger wall components in RG help improve the diffuser performance. Furthermore, investigations regarding turbulence intensities indicated a clear effect on the ﬂow behaviour for IG while being ineffective on the RG. The ﬁnal application analysed the diffuser performance using the inlet conditions extracted directly from a potential radial-inﬂow turbine working with R143a. The change of conditions highlighted that radial components can be reduced, and thus the swirling number too. By implementing the ﬁrst numerical study on real gas swirling conical diffuser, it was established that real gas ﬂow regimes differ from the ones previously established for ideal gas, and thus preliminary ﬂow regimes for R143a, speciﬁcally, are proposed. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Low-to-medium temperature power cycles are being examined to provide better renewable solutions that reduce fossil fuel consumption and CO2 emissions. Organic Rankine Cycles (ORCs) as one solution use high-density ﬂuids to power, in most cases, a radial-inﬂow turbine. The increase in turbine eﬃciency can dramatically improve the overall cycle eﬃciency, critical for energy conversion cycles based on low temperature renewable resources, e.g. waste, geothermal. While the turbine optimisation has been studied extensively, it is the exhaust ﬂow through the conical diffuser that must now be investigated. The turbine component is critical in achieving high cycle eﬃciency and a lot of work has been done in this area [23,37,39,44,46,50]. However, one aspect that is easily neglected when considering the system is the exhaust ﬂow through the conical diffuser. The diffuser is responsible for the pressure recovery at the exit of the turbine and thus can signiﬁcantly contribute to a gain of

∗

Corresponding author. E-mail address: [email protected] (E. Sauret).

http://dx.doi.org/10.1016/j.compﬂuid.2017.03.021 0045-7930/© 2017 Elsevier Ltd. All rights reserved.

overall eﬃciency. However, this pressure recovery can be affected by the separation or recirculation generated by the ﬂow conditions and the geometry of the diffusers. One solution to mitigate this issue is to introduce some swirl at the inlet of the diffuser. McDonald et al. [31] experimental studies discovered that the reduction of separation due to the introduction of swirl resulted in increased performance. They found that the more extreme the separation, the more swirl was needed. However, they also concluded that the introduction of swirl could only provide increase performance when the diffuser was stalled, e.g. when separation is present, and if that was not the case, then the introduction of swirl would have almost no effect. Thus to optimise the diffuser, the two attributes decreasing the conical diffuser’s ability to recover pressure have to be considered, i.e. the recirculation of ﬂow inside the diffuser and the ﬂuid separation from the walls of the diffuser. Senoo et al. [47] also experimentally demonstrated that appropriate inlet swirl for a given diffuser can eliminate partially or even totally separation and recirculation and thus lead to optimum pressure recovery. Armﬁeld and Fletcher [2,3] established that the effect of swirling ﬂuid ﬂows compared to non-swirling ﬂows in conical diffusers is to produce a positive radial pressure gradient in the core ﬂow at the inlet, preventing ﬂow detaching from the walls as

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

the ﬂow diffuses within the conical diffuser. As the ﬂow advances through the conical diffuser the swirl-induced angular momentum rapidly decays and which as a result rapidly decreases the positive radial pressure gradient. This dramatic decrease in the rotational component in turn results in an increase in the positive axial pressure gradient. The larger positive axial pressure gradient is thus reducing the axial velocity at the centreline. To compensate for the mass ﬂow constraint, the axial velocity near the wall is accelerating, preventing separation. Hence swirling ﬂows aim to ensure suﬃcient strength of positive axial pressure gradient and axial near-wall velocity to be maintained throughout the diffuser so that both boundary separation and recirculation at the centreline can be prevented. Clausen et al. [14] deﬁned the boundary separation to be driven by the mostly viscous layer and recirculation in the core ﬂow as an effect of the essentially inviscid ﬂow. This explains why swirl number is critical to balance both of these competing effects, wall separation and centreline recirculation. Furthermore, from the 1930s, inlet conditions and notably turbulence have been identiﬁed as playing a key role in the performance of diffusers [38]. Since then, a large amount of work has been published on the effects of inlet conditions on the performance of conical diffuser as reviewed by Klein [22], and turbulence, as well as swirl, have been demonstrated to play a major role in preventing ﬂow separation. Turbulent velocity proﬁles exhibit a thicker boundary layer when fully developed and thus a larger momentum thickness better suited to resist the adverse pressure gradient [3]. Many experimental and theoretical works have been carried out on diffusers as summarized by Okwuobi and Azad [35] and Azad [4] for both swirling [15,31,47], and non-swirling conical diffusers [5,41]. As discussed above, behaviour and characteristics of turbulent swirling ﬂow in conical diffuser have been extensively studied and reported. McDonald et al. [31] investigated twenty-four different diffusers and established a relationship between inlet swirl number and the ﬂow regime that affects the performance of the diffuser. Senoo et al. [47] experimentally studied ﬁve diffusers with different free-vortex inlet swirling ﬂow while Clausen et al. [14] conducted experiments in an attempt to eliminate wall separation and recirculation within the diffuser region. To achieve a ’perfect’ operating scenario in a conical diffuser, they effectively ’tuned’ the inlet ﬂow conditions and notably the swirl using a solid body rotation so that no recirculation or wall separation was present in the diffuser region. Their case has been added to the ERCOFTAC database and is commonly referred to as the ’ERCOFTAC conical diffuser’. This case has attracted lots of attention for modelling as it can be considered as extreme [1] since any small variation of inlet swirl number for this diffuser can lead to ﬂow reversal. Mathematical simpliﬁcations have been made to axisymmetric Navier Stokes equations, known as the Bragg–Hawthorne equations [9] to describe the inviscid process in swirling ﬂows at high Reynolds numbers. These analytical methods, commonly known as vorticity-streamfunction methods were originally limited to the application of inviscid and incompressible ﬂows [7,9]. That is until the recent extension of the BH equations to compressible axisymmetric swirling ﬂows, made by Maicke and Majdalani [29,30]. Despite these advances, these analytical approaches are limited to inviscid ﬂows and thus can’t reproduce boundary layer separation without the addition of corrective terms. This is an important limitation for dense gas ﬂows in swirling conical diffusers. Modelling swirling and rotating turbulent ﬂows also remains challenging as stated by Jakirlic´ et al. [20] due to the complex ﬂow structures which cannot be accurately reproduced by turbulence models commonly established and validated for simple ﬂows. Lee et al. [26] also recently added that conical diffusers under adverse pressure gradient are among the most diﬃcult cases to numerically reproduce using turbulence models. Adding rotation effects through

101

Coriolis and centrifugal forces requires sophisticated turbulence models to predict turbulence anisotropy, and streamline curvature effects. It is well established now that the standard k − model or other similar two-equation models are not capable to accurately reproduce such complex ﬂow features. For example, Armﬁeld and Fletcher [3] demonstrated that two k − models underperformed compared to their proposed Algebraic Stress Model (ASM), especially for the prediction of the near-wall velocity peak and turbulence quantities for swirling and non-swirling conical diffusers from the experimental studies by Clausen et al. [14] and Okwuobi and Azad [35], respectively. Similar conclusions were drawn by Cho and Fletcher [10] on the same conﬁgurations using an improved ASM model. However, their model over predicted the Reynolds shear stress towards the exit of the diffuser. Gatski and Speziale [18] demonstrated that the ASM models are not suitable for this type of ﬂows and recommended the use of an Explicit Algebraic Reynolds Stress Model (EARSM) instead for three-dimensional rotating ﬂows. More recently, conducted a numerical study on the “internal layers” which are the result of a sudden change in streamwise pressure gradients due to the diverging geometry in diffusers. Wu et al. [52] applied Large-Eddy Simulations (LES) to an asymmetric planar diffuser, known as the Obi diffuser [5], to elucidate the ﬂuid mechanics of the internal layer of the boundary layer within the diffuser region. Duprat [17] also performed LES but on the ERCOTAC swirling conical diffuser. The overall agreement with the experimental data was not satisfactory even though it predicted the experimental trends. This was mainly explained by the strong requirement for the simulations to use highly accurate inﬂow boundary conditions. This requirement is also highlighted by the works of Gyllenram and Nilsson [19], Payette [36] and Bounous [8] who all worked on the ERCOFTAC swirling conical diffuser. To sum up, lots of existing literature has been published on swirling and non-swirling turbulent conical diffusers from experimental, analytical and numerical point of view. Limitations of numerical models have been highlighted and the development of turbulence models capable of accurately reproducing the complex ﬂow structures encountered in swirling turbulent conical diffusers remains a challenge. In addition, the inﬂow conditions are critical for the accuracy of the numerical model. Furthermore, all the existing studies consider “common” ﬂuids in liquid or gas phases such as air or water while ORC turbines work using high-density ﬂuids, commonly referred to as dense or real gases, such as refrigerants. These dense ﬂuids operate in conditions for which compressible effects are present. In addition to this is the added molecular complexity of dense gases, when compared to “common” ﬂuids [12]. To account for these complexities, transport properties of dense gases are calculated through the use of Equation of State (EOS), which are based on high accuracy empirical correlations. Experimental studies of real gas have to this date only been conducted to obtain their thermodynamic properties at various phases. However, there has not been any experimental study of real gas ﬂows to date [34]. This is most likely due to the large variety of available real ﬂuid mixtures and their limitations as well as feasible applicability to various practical applications that potentially beneﬁt from their transport properties. Many numerical studies of dense gases have been conducted using various types of EOS to calculate the transport properties of dense gases in ﬂuid ﬂow [12,13,16,24,34]. So far, to the best of the authors’ knowledge no studies either numerical or experimental have been carried out on swirling and non-swirling turbulent conical diffusers using real gas. Thus, many questions remain open. For example, what will be the real gas behaviour in such conﬁgurations compared to ideal gas? Will the inﬂow conditions, e.g. turbulence and swirl, affect the real gas ﬂow in swirling conical diffuser as in the ideal case? How the turbulence will be changed by using real gas instead of ideal gas in these conﬁgurations?

102

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

The purpose of this paper is therefore to numerically compare the effects of inﬂow conditions on the performance of the ERCOFTAC swirling conical diffuser using air and refrigerant gas R143a (real gas). Air used within this study satisﬁes the ideal gas law assumption. For readers convenience we will hereinafter refer to this component as the Ideal Gas case, denoted IG. For R143A (Real Gas) case will be denoted hereinafter RG. Firstly in Section 2, the k − ω based Baseline - Explicit Algebraic Reynolds Stress Model, denoted as BSLkω EARSM will be presented and detailed. This turbulence model is then validated against the ERCOFTAC swirling conical diffuser database including both experimental and numerical results in Section 3. Section 4 will detail how the National Institute of Standards and Technology (NIST) Reference Fluid Thermodynamic and Transport Properties (REFPROP) database is applied to account for the real gas behaviour of refrigerant R143a (1,1,1-Triﬂuoroethane). NIST REFPROP model is considered as the international standard for the ﬂuid properties of R143a [27]. This model utilises high-ﬁdelity empirical correlations in the Helmholtz Energy Equation of State (HEEOS) to calculate the transport properties of dense gas [28]. The effects of inlet conditions and in particular turbulence and swirl intensities are studied in details for optimal ORC conditions in Section 5. Finally, an investigation is proposed based on inlet conditions extracted directly from a potential radial-inﬂow turbine working with R143a.

where Dk is the total turbulent diffusion. Following the EARSM assumptions:

ui uj Dai j = 0 , Di j = Dk , Dt k Eq. (3) is then reduced to:

ui uj = k ai j +

(1)

where δ ij is the Kronecker operator. The Reynolds stress from Eq. (1) is then introduced in the Reynolds stress transport equation below:

∂ ui uj ∂ ui uj + uk = Pi j + i j + Di j − i j , ∂t ∂ xk

(2)

where Pij is the production term, ij is the pressure-strain, Dij is the diffusion and ij is the dissipation. Eq. (2) thus becomes:

Dai j k − Dt

Di j −

ui uj k

D

k

= Pi j + i j − i j −

ui uj k

( P − ),

(3)

k

( P − ).

(5)

This is an implicit form for the anisotropy tensor aij in which the pressure-strain ij and turbulent dissipation ij terms need to be modelled. The isotropic dissipation assumption is then used for the dissipation leading to:

2 3

i j = δi j ,

(6)

while the pressure-strain term is decomposed into the slow term si j from Rotta [43] and the rapid term ri j from Launder et al. [25] as:

i j = si j + ri j = −c1 ai j − −

2 δi j , 3

ui uj

Pi j + i j − i j =

2. Governing equations and turbulence model The numerical investigation in Cinnella and Congedo [12] of incompressible perhydroﬂuorene (PP10) ﬂow over an airfoil at Re = 500, compared and validated the model against the Blasius solution for incompressible boundary layers. They claimed that compressible effects were negligible and as such validated their numerical approach. Here, the ERCOFTAC ideal gas test case has a Mach number of M = 0.04 and is therefore incompressible. The cases with refrigerant R143a do not exceed Mach numbers of M = 0.07, which correspond to an incompressible ﬂow. However, due to the rapid change in viscosity and density in these dense gases, which vary spatially throughout the domain, a fully viscous compressible solver adapted for low Mach numbers was chosen to account for any of these characteristic changes. The Reynolds-Averaged Navier Stokes (RANS) equations for viscous compressible ﬂows are implemented in the ﬁnite volume solver adapted for low Mach numbers to accommodate for dense gas simulations. Following Gatski and Speziale [18] results, an extension of the Explicit Algebraic Reynolds Stress Model (EARSM) model based on the initial work from Wallin and Johansson [49] is implemented as the turbulence model in order to accurately capture the streamline curvature and rotational effects through the explicit addition of the rotation ij and anisotropy aij tensors. The following relationship (Eq. (1)) is thus used as the starting point of the EARSM linking the Reynolds stresses tensor ui uj to the anisotropy tensor aij and the turbulent kinetic energy k:

(4)

c2 + 8 2 30c2 − 2 ∂ ui ∂ u j Pi j − P δi j − ρk + 11 3 55 ∂ x j ∂ xi

8c2 − 2 2 i j − P δi j , 11 3

(7)

ri j is then re-formulated as a function of the anisotropy aij , strain rate Sij and vorticity ij tensors. Eq. (5) is then expressed as a function of a, S and only from which the anisotropy tensor can be expressed as a projection onto the tensorial integrity basis initially proposed by Pope [40] as: a= β n Tn , (8) n

for which the β coeﬃcients are functions of the following ﬁve invariants:

IIS = tr S2 ; II = tr

V = tr S 2

2

2

; IIIS = tr S3 ; IV = tr S

2

;

.

(9)

The rotation and strain rate tensor are then normalized by appropriate turbulent time scale τ depending on the two-equation model chosen to be linked to the EARSM, either τ = k/ or τ = 1/(Cμ ω ) with Cμ =0.09. In the rotational frame of reference, the non-dimensional form of the vorticity and strain rate tensors are then written:

=

1 τ 2

∂ ui ∂ u j − ∂ x j ∂ xi

− τ i j rot ,

S=

1 τ 2

∂ ui ∂ u j + . ∂ x j ∂ xi

(10)

Where rot is the rotational vector component. Following Wallin and Johansson [49] work using the nondimensional form of the vorticity and strain rate tensors (Eq. (10)), Eq. (8) thus becomes:

Na = −A1 S + (a − a ) − A2 aS − Sa −

2 tr{aS}I . 3

(11)

with N = A3 + A4 (Pk / ) and the Ai coeﬃcients are taken from Wallin and Johansson [49]; A1 = 1.245, A2 = 0, A3 = 1.8 and A4 = 2.25. In order to get the function N in three-dimensional ﬂows, a sixth-order equation needs to be solved for which no explicit solution can be derived. Thus, the analytic solution recommended by Wallin and Johansson [49] is used. Once the N function is established, the Reynolds stress anisotropy tensor aij can be

C.S. From et al. / Computers and Fluids 149 (2017) 100–118 Table 1 All relevant coeﬃcients for Eqs. (16) and (17).

is calculated by the Total Energy Eq. (22). Based on the total enthalpy, htotal :

β

α1

β1

σ k1

σ ω1

α2

β2

σ k2

σ ω2

0.09

5/9

0.075

2

2

0.44

0.0828

1

1 0.856

∂ (ρ htotal ) ∂ P ∂ − + (ρ Ux,y,z htotal ) ∂t ∂ t ∂ xi ∂ ∂ ∂T λ = + (U τ ), ∂ xi ∂ xi ∂ xi x,y,z i j

determined through Eq. (8) as:

ai j =

1 β1 Si j + β3 ik k j − kl lk δi j + β4 Sik k j − ik Sk j 3 2 +β6 Sik kl l j + ik kl Sl j − Skl lm mk δi j − kl lk Si j . 3

(12) the β coeﬃcients being given in Eq. (9). Finally the Reynolds stress tensor is obtained by:

2 kδi j − 2νe Si j + kaextra , ij 3

ui uj =

(13)

with ν e the effective viscosity is given in two dimensions by:

1 2

νe = − β1 kτ ,

(14)

and the additional anisotropy is provided by:

aextra = β2 (S − S ). ij

(15)

In this paper, the BSL − kω model proposed by Menter [33] is added to complete the two-equation EARSM. The complete BSLkω EARSM is shown in Eqs. (16) and (17). This model combines the turbulent kinetic energy dissipation rate ( ) equation from Menter et al. [32] with the speciﬁc rate of dissipation ω equation from Wilcox [51] so that the near-wall turbulence is modelled using ω for higher accuracy and the far wall regions are modelled using .

∂ (ρ k ) ∂ ∂ μt ∂ k + (ρ U j k ) = + Pk − β ρ kω, (16) μ+ ∂t ∂xj ∂xj σk 3 ∂ x j

∂ (ρω ) ∂ ∂ μt ∂ω + (ρ U j ω ) = μ+ ∂t ∂xj ∂xj σω 3 ∂ x j 1 ∂ k ∂ω ω +(1 − F1 )2ρ + α3 Pk − β3 ρω2 . σω 2 ω ∂ x j ∂ x j k

2 ∂ Uk ∂ Ui ∂ Ui ∂ U j ∂ Uk + − 3μt + ρk . ∂ x j ∂ x j ∂ xi 3 ∂ xk ∂ xk

(18)

A blended formulation for the near-wall treatment as described by Eqs. (19)–(21) is implemented in the solver. The blending function in the transport equation (Eq. (17)) is deﬁned by:

F1 = tanh arg41 , where

arg1 = min max with:

CDkω = max 2ρ

(19)

√

k 500ν , βωy y2 ω

,

4ρ k

CDkω σω2 y2

∂k ∂ω , 1 × 10−10 . σω 2 ω ∂ x j ∂ x j 1

(22)

where ρ , is the density, P is the pressure, λ is the thermal conductivity of the ﬂuid, T is the static temperature, τ ij is the stress tensor and Ux, y, z is the vector function of all three velocity components. A non-staggered, cell-centred, Pressure-Velocity coupled compressible solver was used, solving the three momentum equations in addition to pressure, at each integration point. The solver computes compressible ﬂows at any Mach number through the implicit discretization of the product of the density and the masscarrying advecting velocity. The implicit discretization is achieved through Newton–Raphson linearization between the current and the new iterations, where the current iteration of density is linearized in terms of pressure. The pressure is interpolated using the Rhie and Chow [42] approach, where the mass ﬂow terms are discretized to avoid the decoupling of pressure and velocity at adjacent cells. The solution method consists of linearisation of the non-linear equations, implemented into a matrix solution. The discrete system of linearised equations is then solved using an Algebraic Multi-Grid approach. A second-order advection scheme is used based on the boundedness principles detailed by Barth and Jespersen [6] to compute the variable ’β ’ in a non-linear fashion so that it is as close as possible to one (β ≈ 1) for each individual mesh node. The spatial discretization is second-order accurate. The Courant Friedrichs Lewy (CFL) number is based on the time step, which is related to the estimate of time scale of the ﬂuid domain based on length scale, maximum velocity scale in the domain, density, dynamic viscosity, mass ﬂow, total mass and impact of compressive time scale. Here the CFL number is limited to not exceed 5, CF LMax = 5. The time scale is factored with a set constant value of one for simplicity. The CFL number for the ideal case was measured at 2.75, and real gas case at 4.25. 3. Ideal gas swirling conical diffuser: ERCOFTAC case

(17)

The coeﬃcients in Eqs. (16) and (17) are provided in Table 1. The viscous production of turbulence kinetic energy, Pk , is modelled as:

Pk = μt

103

,

(20)

(21)

This formulation automatically switches from wall functions to a low-Reynolds near wall formulation as the mesh is reﬁned. Furthermore, in order to predict transport properties due to compressible effects of the real gas, the energy within the ﬂuid

3.1. Test case The ERCOFTAC swirling conical diffuser has been extensively studied in the literature [1,3,8,10,17,19,36]. The geometric details of the diffuser are provided by Clausen et al. [14]. The diffuser length is 510 mm with a half-cone angle, , of 10° and an area ratio, AR, of 2.84. A swirl generator which consisted of a pipe with an internal honeycomb, had a rotating velocity of 550 rpm, placed at the inlet of the diffuser to generate the turbulence and solid body rotation type swirl, refer to Fig. 1. All the dimensions are provided in Fig. 1, and for clariﬁcation purposes the computational rotational domain inlet and diffuser inlet are denoted as ii and i respectively. As stated earlier, this case is particularly interesting as the swirl number has been set so that neither separation at the wall nor recirculation occurs in the diffuser section. An extension was added at the end of the diffuser domain to account for the experimental outlet condition, i.e. the air exhausted into open space at atmospheric pressure. This extension, commonly referred to a tailpipe, has a length of 1.02 m which corresponds to twice the length of the diffuser as illustrated in Fig. 2(c).Bounous [8] showed that this extended geometry was the best suited to accurately predict the proﬁles towards the exit of the diffuser compared to no extension or a dump section followed by a restriction. The coordinate system in this paper is set in terms of Cylindrical (x, y, z) coordinates, where the x-axis is aligned with the

104

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

Fig. 1. Sketch of the ERCOFTAC conical diffuser adapted from Clausen et al. [14] experiments.

(b)

(a)

(c)

S0

Fig. 2. Graphical representation of; (a) the O-H grid layout on circumferential faces, (b) the corresponding O-H grid nodal allocation and (c) allocation of longitudinal nodes, nx# .

centreline and y-axis directed radially. Corresponding axial, radial and circumferential velocity components are denoted as U, V and W respectively. The radial velocity component along any probe-line (which are normal to the diffuser wall), is denoted as Ur , to avoid confusion with velocity component V. Cartesian coordinate (xw , yn ) is deﬁned as the axis along and the axis normal to the diffuser wall respectively in which the experimental proﬁles have been measured. Both mean velocities and Reynolds stresses were measured at seven different locations (S0-S7 in Fig. 1) within the diffuser as well as the pressure recovery coeﬃcient, CP deﬁned by Clausen et al. [14] as:

Cp =

|P| − Po , 1 ρ Uo2 2

Clausen et al. [14] estimated measurements’ errors for the various parameters as being approximately ± 2% for the mean velocities, ± 5% for the axial velocity at the inlet which is assumed to be uniform, and as large as ± 10% for the Reynolds stresses while the wall shear stress at the inlet was estimated to be within ± 5%. Flow induced by rotation has an intensity which is referred to a dimensionless parameter known as the Swirl Number , commonly denoted as Sn , is deﬁned in this work by the following Cho [11]:

Rin

Sn =

UW r 2 dr . Rin Rin 0 U 2 r dr 0

(24)

3.2. Grid

(23)

where P is the absolute pressure at each given point along the diffuser wall, Po the diffuser outlet pressure, ρ the ﬂuid density and Uo the mean reference inﬂow velocity.

The construction of the computational mesh followed the work of Bounous [8], where the mesh was divided into ﬁve sections as seen in Fig. 2(c). The mesh was made up entirely of threedimensional hexahedron elements, with a standard O-H grid on

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

105

Fig. 3. Grid convergence test comparison; (a) Cp along the diffuser wall, (b) Velocity component U at S0 and S7, (c) Close up view of Cp , and (d) Near-wall U at S0 and S7: Grid {· · · , 1 ; −−, 2; − − −, 3}; Experiment[14] {+, Wal l ; ×, S0; , S7}.

Table 2 Grid parameters. Grid

r#

l#1,2

nx1

nx2

nx3

nx4

nx5

Total nodes

1. 2. 3.

20 50 55

25 50 55

8 11 13

4 5 10

3 8 12

25 35 45

10 20 20

625,0 0 0 987,500 1,250,0 0 0

the circumferential faces, illustrated in Fig. 2(a) and (b) and in Fig. 4. The near-wall components are of high importance due to the associated high possibility, as well as, sensitivity of ﬂow-separation at the walls. The mesh was thus reﬁned in order to capture nearwall components with high accuracy following the previous works of Sauret et al. [45] in agreement with Payette [36] and Bounous [8]. The O-H grid was kept consistent throughout all of the circumferential faces along the computational domain. The offset distance from the outer radius as shown in Fig. 2 was set to 50 mm. Plane view of the grid in the longitudinal direction is illustrated in Fig. 2 (c), showing the ﬁve longitudinal sections with nx# grid points in the longitudinal direction for each section. In total three grids with increasing reﬁnement were constructed using this methodology to carry out grid convergence tests and determine mesh dependency. Details of the boundary conditions are given in Section 3.3. The characteristics of the grids are presented in Table 2. The pressure recovery coeﬃcient Cp and the axial velocity component, U, at S0 and S7 obtained from the three grids are presented in Fig. 3(a), (c) and (b), (d) respectively. From Fig. 3(c), it is clear that grids 2 and 3 give very close results for Cp while the difference with mesh 1 is noticeable. At the inlet of the diffuser (S0), there appears to be negligible difference between the three grids for U, as seen in Fig. 3(b). Towards the end of the diffuser, at S7, it can be seen that grid 1 slightly overestimates the axial velocity U towards the centreline. Grids 2 and 3 are almost identical and accurately predict the centreline velocity. Near the

diffuser wall, the three grids provide similar results for the axial velocity at (S0), while differences are noticed towards the end of the diffuser at (S7) as shown in Fig. 3(d). Even for a highly sensitive parameter such as the axial velocity U at S7, grid 2 shows very close agreement with grid 3 and is thus considered reﬁned enough to provide accurate results for this conﬁguration, both in the viscous and inviscid regions of the ﬂow. Grid 2 is presented in Fig. 4 and has a total of 920,0 0 0 nodes with the y+ < 2 requirement satisﬁed for the wall nearest nodes [36]. Results presented were considered iteratively converged when the Root Mean Squared (RMS) residuals for the mass, momentum and turbulence variables were decreased below 1 × 10−6 . By using the same grid dependency test the same standard grid from IG (grid 2 in Table 2) proved to be also valid for the study of real-gas with the same conditions as ideal-gas, denoted as RG1. However, for the study of real-gas in near critical conditions, denoted as RG2 and presented in Section 4, the grid dependency test showed the grid required reﬁnements in the near-wall region. Firstly grid 2 was modiﬁed by varying the progression ratio in order to decrease the ﬁrst near-wall distance to yw = 1.1e−6 m, compared to 7e−6 m for ideal gas. This was done to account for the increase in density and thus satisfy the y+ w < 2 requirement [36]. Further reﬁnements were performed by ﬁrstly adding 10 extra nodes along the radial direction r# , which allowed for a grid to have similar progression ratio as the IG standard grid and secondly by decreasing the ﬁrst near-wall distance to yw = 0.5e−6 m. Results indicated no improvements to the simulation by either of the reﬁned meshes. The standard grid which has a total of 920,0 0 0 nodes with modiﬁed progression ratio and yw = 1.1e−6 m, was thus considered as a suitable grid for the study of dense gas turbulence swirling ﬂows. 3.3. Boundary conditions The boundary conditions for the BSLkω EARSM, set at the inlet of the rotational frame (ii), are based on the experimental conditions of Clausen et al. [14]. At the end of the tailpipe extension outlet pressure Po , is set to atmospheric with zero-gradient turbulence

106

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

Fig. 4. Computational mesh, top circumferential and bottom longitudinal (not scaled).

treatment. No-slip condition is applied to all walls, including the extension wall and the rotating wall. Opening boundary condition on all the walls of the extension has been tested. However, the results were similar to the case without extension, leading to an increased ﬂow recirculation at the centreline compared to the experiments. Consequently, the no-slip condition was maintained. The inlet, ii, axial velocity component Uo is set to the experimental reference mean inﬂow velocity, 11.6 m.s−1 , with the Ur velocity component set to zero. The temperature is set to room temperature of T = 298.15 Kelvin (K). In order to generate the swirl component a ﬁxed rotational velocity, , is applied on the rotating domain (inlet, ii, and rotating wall), as such the inlet circumferential velocity, W, component is equivalent of W = × Radius. The following section details how the rotation was established. In the Clausen et al. [14] experiment the honeycomb rotates at 550 rpm, an angular velocity of = 57.6 rad.s−1 . However, it was found that using this value in the numerical model results in a circumferential velocity, W, larger than that given by Clausen et al. [14] at the diffuser entrance, S0. As a result, the simulated near-wall components were stronger, and subsequently the centreline velocity was lower than that found experimentally. Clausen et al. [14] deﬁned the swirl in terms of the maximum circumferential velocity, Wmax , divided by the axial reference velocity Uo , where it was also noted that the swirl at the inlet of the diffuser was Wmax /Uo = 0.59. This equates to an angular velocity, deﬁned as W = × r (where r is the radius of rotating domain), of = 52.64 rad.s−1 , which is lower than the set rotation by the swirl generator. This is mainly explained by the losses in the honey comb section. Previous works, and in particular that provided by Bounous [8] and Gyllenram and Nilsson [19], have utilized this swirl relationship to deﬁne the angular velocity, , of the computational rotating domain through the relationship for angular velocity, W. However, this does not account for the decay of rotational velocity in the 100 mm stationary between the domain inlet and the diffuser entrance, i, which resulted in a lower swirl at the entrance of the diffuser. Therefore the angular velocity of the computational domain rotational frame is assumed to be within the bounds of 52.64 < < 57.6 rad.s−1 .

Following Payette [36], who investigated and discussed the importance of accuracy of the inlet radial proﬁles, the angular velocity of the computational rotating domain is estimated to ≈ 56.1 rad.s−1 to match the Clausen et al. [14] experimental data for W along S0. This value was found to give an inlet swirl of Wmax /Uo = 0.61 at the diffuser inlet, i, which is close to the 0.59 value obtained by Clausen et al. [14]. This therefore validates the method of translating the system rotation of the honeycomb to the computational rotating domain. The inﬂow turbulence conditions were not provided in the experimental data and thus the inlet turbulence is deﬁned and controlled by specifying the turbulence intensity, T¯u , where T¯u is taken as 5% based on Armﬁeld et al. [1]. As proposed by Gyllenram and Nilsson [19], the turbulent length scale lt is equal to 0.0032 m at the inlet of the rotational frame which is equal to the cell diameter of the rotating honeycomb as detailed by Clausen et al. [14]. It is possible to approximate kinetic energy, k:

k=

3 (Uo × T¯u )2 , 2

(25)

and from the k − ω based model using the speciﬁc dissipation √ 3 rate, ω, deﬁned as ω = 2 /(k β ) = k/(β lt ) [19], the turbulent dissipation rate, , can be calculated from: 3

=

k2 . lt

(26)

Using Eqs. (25) and (26) with T¯u level of 5% and assuming lt = 0.0032 m gives theoretical values of k = 0.5 m2 .s−2 , and = 112 m3 .s−2 , at the inlet of the diffuser (S0), while the proposed BSLkω EARSM with T¯u = 5% gives numerical values of k and of approximately 0.46 m2 .s−2 and 137 m3 .s−2 , respectively. This excellent agreement with the theoretical values justiﬁes the method of deﬁning the inlet boundary turbulence condition purely by its intensity (T¯u ). Table 3 summaries the boundary conditions used in the BSLkω EARSM.

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

107

Fig. 5. Validation graphs comparing the proposed model against experimental data and previous works by Bounous [8] in terms of dimensionless parameters U, W and k at probe-lines; (S0) Left, (S4) Middle and (S7) Right: − − • − −, BSLkω EARSM; ×, Experiment [14]; −−, k − [8]. Table 3 Boundary conditions. Unit

Value

m/s %

11.6 5

Pa

101,325

rad/s

56.1

Walls (No Slip) Ux, y, z Temperature (T)

m/s K

0 298.15

InletParameters Uo Tu OutletParameters Pout RotationalDomain

3.4. Validation of computational model The results obtained with the BSLkω EARSM are compared in Fig. 5 against the experimental results [14] and the numerical results from Bounous [8] who used a k − model, at locations S0, S4 and S7 (refer to Fig. 1). Three normalised parameters are compared including; axial velocity, U, tangential velocity, W, and kinetic energy, k, at all three locations. These parameters present increasing sensitivity along the three locations, with S7 being the most sensitive. It should be noted that Bounous [8] mentioned that the prediction of the velocity components (U, W) were identical between the k − and the SST models. However, the SST model had considerably worse prediction for the turbulent kinetic energy. Overall, the BSLkω EARSM performs better than the k − model. It is clear that the peaks of mean velocities both axial and tangential and the turbulent kinetic energy near the wall are correctly positioned and estimated with the right intensity (Fig. 5). Both models accurately predict the axial velocity proﬁles (Fig. 5) at different Slocations along the diffuser. However, the BSLkω EARSM better positions the near-wall velocity peak and more importantly accurately predicts the reduction of the axial velocity, U, at the centreline towards the exit of the diffuser (S7, Fig. 5). Speciﬁcally, this indicates

that the BSLkω EARSM accurately predicts interactions between the viscous near-wall ﬂow and the essentially inviscid core ﬂow, throughout the diffuser as the vorticity is (inviscidly) stretched. This accurate prediction of near-recirculating ﬂows is an essential component for turbulent swirling ﬂows in conical diffusers since this will proportionally predict ﬂow parameters near the wall. The tangential velocity W is also better estimated by the BSLkω EARSM all along the diffuser, both in terms of the near-wall intensity and the position of the velocity peak. Far-wall proﬁles are considerably more accurate with the BSLkω EARSM model for S4 and S7. The turbulent kinetic energy proﬁles in Fig. 5 are also slightly improved with the BSLkω EARSM compared to the k − model with a net increase of the peak magnitudes which closely approach the experimental values. It is clear that the k − model better predicts the far-wall kinetic energy. However, given the results presented in Fig. 5 it is evident that the greater far-wall scaling of k does not produce a signiﬁcant impact on the ﬂow. It should also be noted that the k − model initially predicts a greater peak at S0, which then dissipates too much, and consequently under predicts the peaks at S4 and S7. The BSLkω EARSM is also clearly capable of capturing the decrease of the near-wall turbulent kinetic energy, k, due to the inlet solid-body rotation swirl between Sections S4 and S7. This demonstrates the ability of the BSLkω EARSM to correctly capture the main ﬂow characteristics through the diffuser. This also conﬁrms the validity of the boundary conditions applied to the tailpipe downstream the diffuser section, which were implemented to emulate the open outlet. The pressure recovery is also accurately reproduced by the BSLkω EARSM as shown in Fig. 6. Furthermore, the wall shear τ wx and τ wz components are obtained based on the wall shear as:

τw =

2 + τ2 + τ2 = μ τwx wy wz

du dy

|y=0 .

(27)

108

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

Fig. 6. Validated model prediction of pressure recovery, CP along the diffuser, xw ; − • −, BSLkω EARSM; ×, Experiment[14].

τ wx and τ wz are then compared against the experimental data in Fig. 7, where each experimental data point corresponds to a probe-line along the diffuser wall, e.g. the ﬁrst data point is measured at S0 and the last one at S7. Hence collectively there are only eight data-points (S0–S7). Despite this limitation in the available data, the BSLkω EARSM results for the wall shear components are in relatively good agreement with the experimental data, and within the error bounds given by Clausen et al. [14], ± 5%. 3.5. Closing remarks The BSLkω EARSM allowed for accurate simulations of the swirling turbulent ﬂow characteristics throughout the ERCOFTAC conical diffuser, providing good agreement with the experimental data for both the near-wall (viscous) and centreline proﬁles (inviscid) for all parameters at S0, S4 and S7. Despite far-wall kinetic energy proﬁles slightly more accurately predicted by the k − model, the overall prediction are not improved. This is potentially due to the relatively low value of k resulting in negligible effects in ﬂow, despite the sensitivity. This further supports the importance of near-wall parameters to accurately predict turbulent swirling ﬂow characteristics in conical diffusers. The BSLkω EARSM therefore displayed crucial advantages to accurately predict wall separation, and/or, centreline recirculation in swirling conical diffusers, hence establishing solid grounds for the investigation of real-gas swirl ﬂows in conical diffusers. 4. Real gas conical diffuser In order to evaluate the effects of real gas on the performance of a swirling turbulent conical diffuser, the refrigerant R143a is chosen. This ﬂuid has been shown to be a suitable candidate for the development of eﬃcient ORC turbine in low-to-medium renewable power cycles [44,46]. The NIST REFPROP database is directly coupled to the CFD solver to provide the real gas properties. This database models the thermodynamic and transport properties of R143a using the formulation developed by Lemmon and Jacobsen [28], which is based entirely on experimental data. It is then explicitly applied to the ancillary equations and Equation of State (EOS) for Helmholtz energy [27]. The R143a properties produced by this database are considered the international standard and are valid for the triple-point temperature ranging from 161.34 K to 450 K [28].

maintained unchanged compared to the ideal case described in Section 3.3. This case will hereinafter be referred in the text as ’RG1’. At these operating conditions, the thermodynamic state of RG1 is vapour and would result the ORC cycle to be within sub-critical conditions. In order to approach a realistic state of the ﬂuid at the rotor exit of the turbine, the inlet temperature, T, and outlet pressure, Pout , are then modiﬁed. The temperature is thus set to 367.1 K and the outlet pressure to 1.845 MPa based on the R143a turbine investigated by Sauret and Gu [44]. This case is referred to as ’RG2’ in the paper. Here it should be noted that the turbulence intensity, T¯u , is set to ﬁve percent (5%) to emulate the production of turbulence for both RG1 and RG2 cases. A direct comparison between air (IG) and R143a (RG1 and RG2) is demonstrated in Figs. 8 and 9. Similar velocity proﬁles are obtained, although near-wall peaks of the axial U and circumferential W velocities occur much closer to the wall for both RG1 and RG2 than for IG downstream of the diffuser. Intensities of the peaks for U and W are relatively similar. The peaks’ intensity of the kinetic energy k is much smaller for both RG cases from S0 . Whilst the kinetic energy proﬁles are lower, they also appear to reach their peaks closer to the wall compared to IG. However, the evolution of k in the diffuser for both RG cases is consistent with the observations for IG when a solid-body rotation is added to the conical diffuser, i.e., the turbulent kinetic energy increases ﬁrst with peaks pushed near the walls and further downstream of the diffuser the peaks of k decrease [1]. With lower near-wall turbulent kinetic energy and peaks of velocity and k closer to the walls, the RG cases also present lower velocities at the centreline compared to the IG case. This slow down of U velocity component results in a higher tendency for ﬂow reversal at the centreline for RG cases. Towards the exit of the diffuser, at S7, this reduction of the centreline velocity also corresponds to a slight increase of the velocity peaks near the walls. It is noted that very small recirculation right at the exist of the diffuser (located longitudinally approximately 40 mm past the end point of probe-line S7) is present for the RG2 case. Despite being located so close to the end of the diffuser, there seems to be no apparent effect on the pressure recovery performance Cp , as demonstrated in later Section 5 in Fig. 18. The slight increase of near-wall velocity to compensate for the reduction of centreline velocity will, contrary to the IG case, reduce the tendency for wall separation for the RG cases. This is conﬁrmed by the higher shear wall value τ w for RG than for IG as presented in Fig. 10. The trends seen for RG1 in Fig. 8 are only slightly magniﬁed in RG2. For RG1 and RG2, their respective density, ρ , is approximately 3.5 kg/m3 and 60.012 kg/m3 and corresponding Reynolds number, Re, approximately 9.55 × 105 and 1.3 × 107 respectively while the density of the IG (air) is 1.185kg/m3 with Re = 2 × 105 . There is a far greater difference in density between RG1 and RG2 than there is between RG1 and IG. However, this magnitude of difference is not mirrored when comparing the various component proﬁles RG1 and RG2. This suggests that the effect of density does not directly correspond to different ﬂow characteristics in conical diffuser where a solid-body rotation is applied. It is, however, obvious that the ﬂow conditions differ with the different ﬂuid properties and as such different turbulence is produced by the swirl generator and different boundary layer characteristics are expected. This will be discussed in the following sections. However, due to the thermodynamic state of vapour-gas, in a sub critical conditions for RG1, this case will be mainly disregarded in the later sections as this is an undesirable state for the ORC turbomachinery application and RG2 will be the main focus of the following sections.

4.1. Initial investigation of swirling real gas ﬂows 4.2. Evolution of the boundary layer parameters in the diffuser In a ﬁrst attempt to establish a direct comparison between the two ﬂuids (ideal IG and real gas RG), the geometry and the basic boundary conditions (Uo , Tu , Pout and , T) in Table 3 are

The evolution and the characteristics of the boundary layer along the diffuser wall for both IG and RG are investigated. One

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

109

Fig. 7. Validation of Wall Shear components (a) τ wx and (b) τ wz , normalised by the density along the wall ρ w and the velocity squared : − − • − −, BSLkω EARSM; ×, Experiment [14].

Fig. 8. Dimensionless parameters U, W and k comparison between RG1, RG2 and BSLkω EARSM for at probe-lines; (S0) Left, (S4) Middle and (S7) Right: − − • − −, IG; − · − · −, RG1; − − −, RG2.

of the characteristic parameters which inﬂuences the ﬂow at the walls is the wall shear stress, τ w , as calculated using Eq. (27). The wall shear, τ w , presented in Fig. 10, is clearly larger for RG2 than IG. This is the result of RG2 having a greater velocity gradient, dU /dy [s−1 ], at the walls, as the dynamic viscosity measures μ = 1.416 × 10−5 Pa.s, lower than IG (μ = 1.813 × 10−5 Pa.s). Therefore, with a greater wall shear compared to IG, RG2 has a lower tendency to separate at the walls. The wall shear for IG drops faster just after the inlet of the diffuser and reaches its lower value before the middle of the diffuser. The decrease in wall shear for the RG2 case is less pronounced and remains relatively stable after the initial drop while the IG increases after reaching its lowest value. This can be explained by the velocity peak near the wall that does not drop much from S4 to S7 (Fig. 8) for RG2 while the velocity magnitude keeps decreasing for the IG case due to a relatively larger decrease of the turbulent kinetic energy near the wall (Fig. 8). The boundary layer thickness δ 99% is estimated at each probeline from S0 to S7, refer to Fig. 1, by measuring the distance (yn )

between the wall and the point for which 99%Umax is reached for a particular probe-line. In order to understand why RG better resists wall separation than air (IG), another interesting parameter to consider is the friction velocity, which is calculated as per Eq. (28) below:

uf =

τw . ρ

(28)

The normalized friction velocity, uf /Uo , is shown in Fig. 11 (a). The evolution for both ﬂuids display similar characteristics. However, RG2, despite a higher wall shear seen in Fig. 10 presents an overall lower friction velocity than IG. This is due to the much higher density and viscosity of RG2 compared to IG, again resulting in a better capacity to resist wall separation. Additionally, the boundary layer thickness, δ 99% for RG2 as seen in Fig. 11 (b), is approximately half that of IG along the diffuser wall. Thus, it is clear that the boundary layer develops slower for RG2 than for IG. The thinner boundary layer δ 99% seen for RG2 is potentially due to the lower

110

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

Fig. 9. Mid-plane Velocity, U, contour plot: Comparison between (a) IG and (b) RG2.

before suddenly decreasing at the last section, S7. This behavior is similar for both ﬂuids. However, the drop in θ between S6 and S7 is accentuated for RG2. This drop also corresponds to the drop in the peak magnitude of the near-wall turbulent kinetic energy at S7 as seen in Fig. 8. It can also be seen that the momentum thickness is overall larger for RG2 than for IG conﬁrming again its better tendency to resist wall separation than IG and explaining the slower development of the boundary layer for RG2 as shown in Fig. 11(b). Similarly with the shape factor, H in Fig. 12(b), there is a steady increase throughout the diffuser wall even past S6 for both ﬂuids, and thus there is a greater opposition to adverse pressure gradients in RG2, and hence a more stable near-wall boundary layer. The results presented in this section conﬁrm that RG2 has less tendency to wall separation compared to IG due to higher wall shear and a thinner boundary layer which also develops at a slower pace due to more momentum θ in the boundary layer. This, in turn, results in a greater shape factor, H, and thus increased resistance to adverse pressure gradient. It was also noticed that the boundary layer and velocity and kinetic energy proﬁles differ between the two ﬂuids from S0, leading us to investigate the effect of the parameters at the diffuser inlet(i) in the next Section 4.3 and the inﬂow turbulence in the later Section 4.4 in order to establish the inﬂuence of inﬂow components.

4.3. Characteristics of the boundary layer and velocity proﬁles at the entrance of the diffuser

Fig. 10. Wall Shear τ w comparison between IG and RG2: −−, IG; − − −, RG2;.

friction velocity uf produced at the wall thus allowing for greater ﬂuid ﬂow at the wall as shown in Fig. 8. The displacement, δ ∗ , and momentum, θ , thicknesses are calculated using the relationships in Eqs. (29) and (30) respectively, following both [48] and [21].

δ = ∗

θ=

yn

U (r ) 1− Umax

0

0

yn

r dr, yn

U (r ) U (r ) 1− Umax Umax

r dr, yn

(29)

(30)

where, r, is the distance from the wall along yn , and the velocity component U(r) is the axial velocity at each corresponding distance r, extracted at each individual velocity probe-line S0 through to S7, refer to Fig. 1. Both equations are then solved by utilizing the trapezoidal integration technique. Additionally the shape factor H is determined using Eq. (31):

H = δ ∗ /θ .

(31)

The momentum thickness, θ , for both IG and RG2, is presented in Fig. 12(a). θ increases throughout the diffuser until location S6

In order to identify the primary inﬂuence of the ﬂow characteristics developed within the diffuser for both ﬂuids, ﬂow characteristics at the inlet(i) of the diffuser are further investigated. The review of literature indicates that the application of swirl for ideal gas induces an increased positive radial pressure gradient, and thus increases the diffuser wall components. This is also what is observed so far for the dense gas. As such this section investigates the boundary layer characteristics at the entrance of the diffuser (i) and the inlet ﬂow parameters, including the swirl number, Sn , as deﬁned in Eq. (24) in order to quantify the strength of swirling ﬂows produced by both IG and RG2 and the resulting velocity components. This swirl is produced by the computational rotational domain angular velocity set to = 56.1 rads−1 . So far, for the conditions presented, IG has an inlet(i) swirl strength of Sn = 0.3, with both RG1 and RG2 cases performing similarly with swirl number of approximately 0.28, calculated using Eq. (24). The swirling ﬂow generates a positive radial proﬁle towards the wall at S0, where comparatively both RG cases were found to have a greater strength in the radial velocity Ur , as shown in Fig. 13. The increase in the positive radial velocity proﬁles of RG at S0 establishes that there is dependence on the inlet, i, ﬂow conditions. To compare how IG, RG1 and RG2 differ at the inlet(i), Mach number Ma , mass ﬂow rate M˙ and boundary layer components are detailed in Table 4. Comparing IG and RG2 indicates that the density ρ and similarly the mass ﬂow M˙ are 50 times greater for RG2 than IG and thus could explain the 14% increase in the peak of Ur as well as the changes in the boundary layer parameters. However, when compared to RG1, for which the density is only increased 3 times compared to IG, but is still 17.2 times lower than RG2, it is clear that the swirl level is not directly linked to the ﬂuid density as Sn decreases by almost the same quantity for both RG cases compared to IG. Furthermore, the maximum of Ur increases much more between RG1 and IG than between RG2 and RG1 despite a much larger increase in density, probably due to similar swirl levels for both RG cases. It appears that this density change and thus ﬂow rate variation mainly affect the boundary layer thickness and momentum but have almost no effect on the shape factor.

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

111

Fig. 11. IG and RG2 comparison: (a) friction velocity, uf , and (b) boundary layer thickness estimated at S0 − S7: (a){ − − • − −, IG; − − −, RG2 } : (b){ −−, IG; − − −, RG2 }.

Fig. 12. IG and RG comparison in; (a) Momentum thickness θ and (b) Shape-factor H = δ ∗ /θ estimated at each probe-line (S0 − S7): −−, IG; − − −, RG2.

Fig. 13. Dimensionless inlet radial proﬁle, r/Rin comparison against real gas cases, at probe-line S0 (a) Overview, (b) Close up: −−, IG; − · −·, RG1; − − −, RG2.

Table 4 Inlet parameters post-processed at S0.

ρ [units]: IG: RG1: RG2:

Ur(max) −3

[kgm ] 1.185 3.489 60.012

−1

[ms ] 0.248 0.274 0.283

Wrot −1

[J.s ] 2.514e−3 7.286e−3 0.1236

M˙

Ma

δ

HS0

θ

Sn

[kg.s−1 ] 9.453e−5 3e−4 4.65e−3

[−] 0.033 0.065 0.066

[m] 0.0139 0.0104 0.0085

[−] 1.053 1.059 1.064

[m] 2.8e−3 3.1e − 3 3.4e−3

[−] 0.3 0.28 0.278

112

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

Fig. 14. Near-wall dimensionless turbulence kinetic energy proﬁles, k/Uo2 , at various T¯u %, (a) S0 located near the inlet(i) and (b) within the diffuser domain at S7: −− ,IG {, T¯u 10%; •, T¯u 5%; , T¯u 1%}: − − − ,RG2 {, T¯u 10%; •, T¯u 5%; , T¯u 1%}.

The slightly lower swirl for RG cases also indicates that there is some reduction in velocity components U and W. It is obvious that the higher density, ρ , and the near-wall circumferential velocity, W, of RG2 results in an increased mass ﬂow rate, M˙ , and therefore produces greater rotational work Wrot , thus a greater force subjected radially towards the wall and higher adverse pressure gradient. This additionally provides a reason for the increased wall shearing component, resulting in RG2 producing more stronger wall components capable to better resist separation. However, should this be true, then this effect should be linearly proportional to the increase in density. This is not the case given the results of RG1 (Table 4). This suggests that the boundary conditions have potential impact on the ﬂow characteristics in an exponential fashion rather than linear. To clearly differentiate the effects of density from the swirl strength, a variation of swirl at constant ﬂow rate for both IG and RG2 will be presented in Section 5. However, varying the ﬂow rate will affect the swirl through a change of inlet axial velocity and thus is not considered here.

Fig. 15. Velocity proﬁles U/Uo at various T¯u % for RG2 and IG, at S0 and S7: −− ,IG {, T¯u 10%; •, T¯u 5%; , T¯u 1%}: − − − ,RG2 {, T¯u 10%; •, T¯u 5%; , T¯u 1%}.

4.4. Effect of inﬂow turbulence Investigations at the entrance of the diffuser (inlet(i) ) are conducted in an attempt to explain the large variance in ﬂow characteristics seen so far between IG and RG. It is known that the turbulence quantities have an important impact on the swirling ﬂows conditions within the diffuser. As detailed by Armﬁeld and Fletcher [3], turbulent velocity proﬁles are fuller than laminar in the boundary layer, and thus a larger momentum θ at the wall is created. The same analogy may be applied when considering ﬂows with high turbulent intensity against ﬂows with less intensity. In the ERCOFTAC case some turbulence is produced by the swirl generator, referred to as the computational rotating domain, where additional turbulence originates from the inﬂow conditions. However, the actual production throughout this rotational domain is not known. Therefore, to investigate the effect of the turbulence production for both ﬂuids and its corresponding effects in the diffusing domain, three inlet turbulence intensities, T¯u = 1%, T¯u = 5% and T¯u = 10% were applied at the inlet of the rotational computational domain (ii). It was found that for IG, the variation of T¯u do not appear to have much effect on the boundary layer at the inlet of the diffuser, although they do have an important impact on the k proﬁles towards the centreline, as seen in Fig. 14. Within the diffuser, where the ﬂow is subject to diffusion, these same observations are also valid for the medium (T¯u = 5%) and high (T¯u = 10%) turbulence in-

tensity settings, where the low (T¯u = 1%) setting looses strength at the peak as a result of the diffusion. Also surprisingly, it appeared that the results for RG2 are almost identical for the different turbulence intensities. In fact, it was possible to set a turbulence intensity as high as 10% and still maintain the same ﬂow proﬁles throughout the diffuser, as seen in Fig 14(b). This result indicates that the turbulence in the RG2 is primarily produced by the set rotating domain, meaning that the inlet boundary conditions for the computational rotating domain in combination with a reference velocity, Uo , do not have a suﬃcient quantity to beneﬁt from any variation of turbulence intensity in the ﬂuid inﬂow. Looking at the axial velocity proﬁle U for RG2 and IG at the various turbulence intensities T¯u at probe-line S0 and S7 in Fig. 15, echoes the discussed observations for kinetic turbulence energy in Fig. 14. Additional effects for IG result in increased near-wall velocity subsequent to the increase in T¯u causing an increase in near-wall momentum and a slight increase at the centreline. It is noted that the 10% intensity for IG helps increase the robustness of the ﬂow characteristics, as generally recirculation occurs within the extension component (past-diffuser exit), however at 10% the recirculation is completely absent. If this result was reproducible for RG2 it would potentially allow easing the tendency of recirculation in the RG cases. However, it is evidenced that RG2 is completely unaffected by an increase of T¯u , and thus the question of the beneﬁts of varying the intensity T¯u for R143a remains open and requires further research.

C.S. From et al. / Computers and Fluids 149 (2017) 100–118 Table 5 Inlet parameters at various turbulence intensities T¯u post-processed at S0.

δ θ H

[unit]

IG1%

IG5%

IG10%

RG21%

RG25%

RG210%

[mm] [mm] [-]

11.47 3 1.0561

13.9 2.8 1.053

22.17 2.1 1.042

8.5 3.4 1.0633

8.5 3.4 1.064

8.5 3.3 1.064

4.5. Discussion As stated earlier T¯u provides some advantages in diffuser performance when using ideal gas. However, the same cannot be concluded for RG. In Fig. 14, the proﬁles for all T¯u intensity levels for RG2 present almost no change at either the inlet or near the diffuser exit probe-line S7. Considering the boundary layer parameters in Table 5, it can be noted that the boundary layer is a lot thinner for RG than IG. The boundary layer thickness remains constant with the variation of T¯u for RG while it keeps increasing for IG. The momentum and the shape factor H are also larger and almost constant for RG while they decrease with an increase of T¯u for IG. This suggests that there is a high possibility that the swirl generator develops greater turbulence that outweighs any inﬂow conditions. This is potentially due to the excessive increase in density - induced rotational momentum which hence leads to no changes when the inﬂow condition has a high or low turbulence intensity. From the results presented in Section 4.3, it was initially speculated that the boundary conditions do potentially propose some limitation in how high-density swirling ﬂow behave. It is possible that given the excessive density, and in turn high Reynolds number, an extremely high quantity of turbulence is produced, where most of this turbulence is then directed towards the wall. These preliminary observations can be further utilized to establish guidelines as to how real-gas ﬂow effects should be further investigated. It is stressed that these results for RG2 are limited to the set boundary conditions and geometry considered here. 4.6. Closing remarks In this section, it was demonstrated that the change in ﬂuid properties, dynamic viscosity and density, modiﬁes the characteristics of swirling ﬂuid ﬂows within the conical diffuser. The difference in the radial velocity component Ur for both ﬂuids appears to be critical, since evidence indicates that the signiﬁcant increase in the radial components (e.g. turbulence components, radial proﬁles, momentum) is due to the density-induced rotational momentum and hence the radial component is representative of that effect. More robust wall ﬂows are produced by dense gas R143a, although consequent to this is a large reduction in the centreline axial velocities to the extent that recirculating ﬂows occur at the very exit of the diffuser domain for RG2. This recirculation is located close to the exit of the diffuser and has no apparent impact on the diffuser pressure recovery. Furthermore, the higher density causes a stronger wall-shear, despite a lower friction velocity. A stronger wall shear associated with a larger near-wall turbulence momentum and greater shape factor leads the RG cases to have more robust wall behaviour and thus less tendency for wall separation. Moreover, the presented results suggest that the increase in density alone does not have a proportional linear effect on the development of ﬂow characteristics within the diffuser. Rather, there is almost an inverse exponential effect for which the change in density only provides dramatic changes until a certain point, thus, the larger variance between cases IG and RG1, rather than RG1 and RG2 as demonstrated in Table 4. This is particularly the case for the given conditions, as the magnitude of angular velocity combined with the mean inlet velocity does not have a suﬃcient

113

impact on the ﬂuid, with a density beyond that of RG1, to have any real signiﬁcant effect on the ﬂuid ﬂow. This can also explain the absence of a change in turbulence intensity (T¯u ) of the ﬂow of RG2 compared to IG. This indicates that conical diffusers utilized in ORC applications may have less sensitivity to inﬂow conditions when subject to solid-body rotation. More reliability in the consistency of performance, even when subjected to disturbances, may also be a great advantage for this application. In order to investigate this aspect further, the following section will address the effect of a lower swirl setting for both IG and RG2.

5. Inﬂow swirl effects 5.1. Flow characteristics within the diffuser at various Swirl numbers As discussed, the use of high density ﬂuids results in more robust wall components due to the density-induced radial components leading to a lesser tendency to separation. This comes at the expense of lower centreline velocities making the RG more susceptible to recirculation. It is therefore suggested that a reduction in swirl intensity for the high-density ﬂuid could potentially help reduce the tendency for recirculation at the centreline. It is known from both literature [19,36,45] and ﬁndings within this paper, that the diffuser inlet conditions are critical. Given that U, W and Sn are all near identical at the inlet(i), we will consider the inlet radial velocity Ur proﬁles due to the density-induced radial components discussed in Section 4. As such, the possibility to establish an inlet(i) radial proﬁle nearly identical to that of the IG case is desired at the same ﬂow rate. In order to match the inlet(i) radial proﬁles of RG to the IG case, the computational rotating domain is decreased to an angular velocity of = 36.1 rads−1 for the RG2 case. This can be seen in Fig. 16, which directly compares the radial inlet(i) proﬁles for IG and RG2, at standard angular velocities, = 56.1 rads−1 and modiﬁed with the lower angular velocity setting of = 36.1 rads−1 . Hence the modiﬁed real gas case RG∗ matches the inlet(i) radial velocity proﬁle Ur of IG and the modiﬁed IG∗ now has a considerably lower Ur proﬁle than the standard IG. The respective swirl number Sn of standard cases are; SnIG = 0.3, SnRG2 = 0.278 and modiﬁed cases both at approximately Sn = 0.18. All parameters, U, W, k, and Ur are compared in Fig. 17, at the three main probe-lines S0, S4 and S7. Focusing on inlet probe-line S0 in Fig. 17 it can be seen that all cases, except for IG∗ , have the same inlet U proﬁle, where the peak of IG∗ is scaled slightly more towards the wall. Near-identical kinetic energy k proﬁles are displayed between the standard and modiﬁed cases for the respective ﬂuids at S0. However, due to the lower strength induced by the solid body-rotation, the peaks of k have lower magnitudes for the modiﬁed cases. The radial velocity proﬁle Ur in Fig. 17 is as per discussed previously in Fig. 16. While decreasing Ur through a decrease of Sn , there is a continuous increase of the turbulent kinetic energy from S0 to S7 for IG∗ . This is counter to the behaviour of IG for which the solid-body rotation is strong enough to reduce the turbulent kinetic energy towards the exit of the diffuser at S7. This can also be seen on the radial velocity proﬁles where Ur keeps increasing for IG while for IG∗ , Ur increases from S0 to S4 but to a lesser extent and then noticeably drops from S4 to S7. This conﬁrms that the solid-body rotation is not strong enough to maintain its effect all the way through the diffuser. Due to this continuous increase of k for IG∗ , and its maximum peak displaced away from the wall, the peak of axial velocity near the wall is shifted towards the centreline with a dramatic decrease of the near-wall velocity leading to separation. This also corresponds to Ur reversing near the wall at S7. To compensate for this decrease of velocity near the wall, the centreline velocity

114

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

Fig. 16. Dimensionless inlet radial, r/Rin , proﬁles of Standard and Modiﬁed, IG and RG cases measured at S0; (a) Overview, (b) Close up: −−, IG; − − −−, IG∗ ; − − −, RG2; , RG∗ .

Fig. 17. Comparison between RG2=36.1 , RG2 and BSLkω EARSM dimensionless parameters U, W, k and Ur , at probe-lines; (S0) Left, (S4) Middle and (S7) Right: −−, IG; − − −−, IG∗ ; − − −, RG2; , RG∗ .

dramatically increases for mass ﬂow conservation, which suggests a stronger inviscid core and weaker viscous boundary layer. For the RG cases, the effect of the decrease of Sn is not as critical as for IG. Indeed, the turbulent kinetic energy increases for RG∗ from S0 to S4 but decreases from S4 to S7 as observed for RG2 and IG. It is clear from the radial velocity proﬁles that the increase in Ur from S0 to S4 is still suﬃciently large to almost reach the level achieved by the IG case. Because of that, the decrease of Ur from S4 to S7 is not as dramatic as for IG∗ , ending at a maximum close to the one reached by IG∗ at S4. As for IG∗ , the peak of axial velocity for RG∗ is thus moved away from the wall with a

decreased intensity and the centreline velocity increases but to a lesser extent than IG∗ . Therefore, for RG∗ , no separation nor recirculation occur in the diffuser, as opposed to IG∗ which shows some wall separation. This conﬁrms that the real gas is more resistant to adverse pressure gradient than the ideal gas and can better handle a reduction of inlet solid-body rotation. Furthermore the diffuser performance of each case is displayed in terms of pressure recovery coeﬃcient in Fig. 18, deﬁned by Eq. (23). It is clear that the modiﬁed IG∗ has the lowest performance, which is most likely due to the separated ﬂow at the wall. All other cases perform similarly, although it should be noted that the

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

Fig. 18. Standard and Modiﬁed∗ cases Diffuser pressure recovery performance CP comparison: −−, IG; − − −−, IG∗ ; − − −, RG2; , RG∗ .

Fig. 19. Wall Shear τ w comparison between standard and modi f ied ∗ : −−, IG; − − −−, IG∗ ; − − −, RG2; , RG∗ .

standard RG2 presents the best performance overall, despite the slight recirculation present at the diffuser exit. 5.2. Boundary layer development within the diffuser The diffuser wall shear, in Fig. 19, indicates that the modiﬁed and standard cases initially have same shearing at the start of the diffuser, where this shear then dramatically drops. This drop is even more pronounced for IG compared to RG. Modiﬁed IG∗ has a shift in proﬁle at approximately xw = 0.13 m, which corresponds to the re-circulating ﬂow at the wall (e.g. separation occurs between S3 and S4). The recirculating ﬂow also appears to induce wall shear as it increases shear until xw = 0.3 m before slightly loosing strength. This means that any additional length to this diffuser case could potentially result in the recirculating ﬂow detaching from the wall entirely. The modiﬁed high-density case RG∗ displays a continuous decrease right after the diffuser entry. At approximately half way through the diffuser the shearing strength becomes lower than IG, and at the end shares the same wall shear as IG∗ . Unlike the standard cases, the modiﬁed cases do not beneﬁt from this recovering wall shear past S4, which is the result of lower solid-body rotation decreasing wall-components as the ﬂows diffuses. Formulation for the kinetic energy production term, Pk is estimated by Eq. (32), which follows the simpliﬁed BSLkω EARSM assumption detailed by Menter et al. [32]. Pk = −uν

dU dU = Cμ k , where along the wall xw , y = yn = 0. (32) dy dy

115

Fig. 20. Kinetic energy production Pk , along the diffuser wall comparison between standard IG and RG2, and modiﬁed: −−, IG; − − −−, IG∗ ; − − −, RG2; , RG∗ .

where the BSLkω constant Cμ = β ∗ = 0.09, and k is the kinetic energy at the given point of measurement, where dU/dy is the velocity gradient along the diffuser wall. The near-wall production Pk , is then plotted along the diffuser wall xw in Fig. 20. This indicates that at the wall the maximum peak production is located at the entrance of the diffuser, with the high-density ﬂuids having greater production due to greater velocity gradients, linked to higher density effects. The standard cases have a steeper decrease than the RG cases, suggesting that the gradients are maintained more strongly for RG than for IG. The production term, Pk , of the standard cases decreases until around xw = 0.25m, before slightly increasing, which corresponds to the behaviour observed in the wall shear in Fig. 19. The momentum θ thickness and shape factor H, for standard and modiﬁed cases are displayed in Fig. 21(a) and (b), respectively. The momentum thickness indicates a continuous increase for both modiﬁed cases, RG∗ and IG∗ , while the standard IG and RG2 decrease at the last section. This continuous increase in momentum thickness θ explains the continuous increase in k all the way through to S7, seen in Fig. 17. The shape factor H for the modiﬁed cases are overall lower, however they do not display any distinctive difference in characteristics from their corresponding standard cases, essentially mirroring the proﬁles of the standard cases. It is however rather unexpected that the shape factor H is greater for IG∗ than RG∗ . Although the point where H for IG∗ starts having a greater value corresponds exactly to where the wall recirculation initiates as shown in Fig. 19. At this point, the displacement thickness remains stable with a continuously low momentum thickness leading to a greater shape factor. It is clear that for the same ﬂuid, the effect of decreasing the swirl level and thus reducing the radial inlet velocity leads towards weaker near-wall components, velocity and turbulence, leading towards separation, as expected by the viscous-inviscid interaction in diffusing swirling ﬂows. However, due to the effect of viscosity and density, the RG presents a greater resistance to separation than ideal gas when subject to diffusion. As presented in Table 6, a large separation at the wall is observed for IG∗ but no separation occurs for RG∗ . This is due to a larger radial velocity proﬁle at the diffuser even at lower rotational speed resulting from the higher density. This indicates that real gas and ideal gas ﬂow regimes are different in swirling conical diffusers as detailed in the following section. 6. Flow regimes Armﬁeld and Fletcher [3] characterized the swirling ﬂows in the conical diffusers into three categories: weak swirl, moderate swirl and strong swirl depending on the swirl number (deﬁned by Eq. (24)), conical half-angle, and ﬂow type which consists

116

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

Fig. 21. Standard and modiﬁed∗ IG and RG comparison in; (a) Momentum thickness θ and (b) Shape-factor H = δ ∗ /θ estimated at each probe-line (S0 − S7): − − , IG; − − ∗, IG∗ ; − − −, RG2; − − − + , RG∗ .

Table 6 Summary of diffuser ﬂow quality between standard and modiﬁed∗ with half-cone angle of = 10◦ . ERCOFTAC (IG)

IG∗

RG2

RG∗

Ideal ﬂow with some reduction at centreline

Strong separation at wall

Ideal ﬂow, with negligible recirculation at the exit

Ideal ﬂow, with near separation at wall

Table 7 Preliminary proposal of R143a swirling ﬂow regimes for diffusers with half-cone angle of = 10◦ . Description

IG

RG2

Weak Swirl - Separation at the wall and small reduction in centreline axial velocity Moderate Swirl - No-separation, near-reversing centreline axial velocity

(Sn < 0.3) (0.3 < Sn < 0.6)

(Sn 0.18) (0.18 < Sn < 0.3)

of either elliptic or parabolic. In this paper, the IG and RG cases presented in both the standard and the modiﬁed conditions are considered elliptic ﬂows, parabolic ﬂows being generally associated with small half angle diffusers. The results presented so far found that the near-wall velocity components are far greater for the high-density refrigerant R143a (RG), where a subsequent velocity reduction is seen at the centreline, which thus potentially places RG within the moderate swirl ﬂow regimes of the table detailed in Armﬁeld and Fletcher [3]. i.e. no separation at the walls but near-reversing ﬂow at the centreline. However, considering the slightly weaker swirl places RG within the sub-weak swirl region (

BSLkω model proposed by Menter et al. [32]. Validation was established by comparison against experimental data, which concluded that the BSLkω EARSM, in conjunction with the boundary conditions provided highly accurate prediction of the ﬂow characteristics throughout the diffuser. The validated CFD model was then coupled with REFPROP to model the thermodynamic and transport properties of dense gas R143a, to simulate the real gas swirling ﬂow, ﬁrstly at the same temperature and pressure settings as the ERCOFTAC case (RG1) then at its near-critical state (RG2). It was found that the increase in density of the real gas caused an increase in the radial velocity components, as a result of density induced rotational momentum. Consequently all near-wall velocity components have dramatically increased, which in return reduced the centreline axial velocity. Further investigations found that near the wall the turbulence production Pk was greater for RG, as well as the boundary layer momentum thickness θ and shape factor, H, and a thinner boundary layer, δ , subsequent to a lower friction velocity uf . It was also found that at the wall, the high-density ﬂuid presented greater wall shear strength, despite the lower friction velocity uf . Considering these factors further conﬁrmed the effect of higher-density ﬂuid in increasing the density-induced rotational momentum. Investigations also found that surprisingly RG is not affected by changes in the turbulence intensities T¯u at the inlet of the computational domain, in contrast to IG that displayed sensitivity in kinetic energy produced at the walls and centreline when subject to various turbulence intensities T¯u . This unexpected result suggests that the turbulent ﬂow components are entirely driven by the rotating domain, such that there is almost no transfer in turbulence from before to after the rotating domain. It is hypothesized that this is due to the boundary conditions which are not within a range to have any potential effect on the

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

turbulence characteristics for RG. Even for turbulence intensities up to 10%, the effect is negligible on the ﬂow characteristics for RG. This could potentially also provide reasoning as to how RG1 and RG2 display such near-identical proﬁles despite the dramatic difference in density. These speculations are yet to be assessed further in greater detail, by using various real-gas ﬂuids with different densities. Solving this question will allow for a greater understanding of the effect of the inlet turbulence components such as k and in swirling high-density ﬂows. If future ﬁndings are consistent with the results presented here that will conﬁrm that the inlet turbulence intensity is not one of the main factors affecting the swirling ﬂow characteristics for turbulent high-density ﬂuids, unlike ideal gas swirling ﬂows. More detailed investigations are needed to fully understand the turbulent ﬂow behaviour in the rotating domain and a straight pipe subject to solid-body rotation may be a judicious case to study this aspect. The present study also highlighted the importance of the radial velocity component at the entrance of the diffuser through the investigation of the inﬂuence of inﬂow swirl. It was shown that a large decrease of the angular rotation was required to match the radial velocity proﬁle of the RG at the inlet of the diffuser with IG values. In doing so, the swirl number was reduced but the RG case was still capable of preventing the ﬂow from separation and recirculation. The reduction of swirl level in the IG case on the other hand leads to lower radial velocity and to a large separation at the wall due to a dramatic decrease of the radial velocity component through the diffuser. These preliminary investigations of real gas when subject to different swirl strengths conﬁrmed that higher-density ﬂuids induce greater positive radial velocity components even when subject to solid-body rotation, and thus beneﬁts from lower swirl intensities when compared at the ideal gas assumption. This conﬁrms previous speculations that overall the high-density ﬂuids provide greater and more robust wall components at lower swirl intensities. As such it is evident that real gas ﬂows with higher-density could also potentially lead to longer diffusers and/or diffusers with a greater diffusing angle. Future work should now investigate the sensitivity of these geometric factors on the performance of real gas swirling conical diffusers to improve the understanding of the density-induced rotational momentum. From this study, it is apparent that deﬁning the optimum working conditions, such as the geometry, inlet velocities, turbulence intensities and other factors for high-density swirling ﬂows in conical diffuser for speciﬁc density and dynamic viscosity ranges is critical to develop robust guidelines for conical diffuser applied to ORC turbines. Thus, the improved understanding of the ﬂow characteristics in high-density swirling conical diffusers may inform the design of the rotor blades of the ORC turbine. For example, the radial velocity component at the exit of the rotor could be reduced while maintaining a suﬃcient swirl level through an increase of the axial velocity in order to establish favourable ﬂow characteristics for the conical diffuser. In return, this will potentially improve the overall ORC turbine eﬃciency. In conclusion this paper presents a preliminary numerical investigation of real gas swirling turbulent ﬂows in a conical diffuser coupled with the most accurate transport properties for R143a derived from empirical correlations, with the most realistic RANS model for turbulent swirling ﬂows. The core issue, and one of the most important aspects to address is whether the turbulent real gas ﬂow phenomena can be accurately, or even suﬃciently, predicted with numerical models coupled with EOS models that are based on high-ﬁdelity empirical correlations. Future advancements to take place in real gas applications such as renewable ORC turbines, would therefore require experimental validation and DNS simulations’.

117

Acknowledgments The authors would like to acknowledge the High Performance Computing (HPC) facilities at Queensland University of Technology (QUT) for its support to carry out the numerical simulations. Dr. E. Sauret would like to thank the Australian Research Council for its ﬁnancial support (DE130101183). References [1] Armﬁeld SW, Cho C-H, Fletcher CAJ. Prediction of turbulence quantities for swirling ﬂow in conical diffusers. AIAA J 1990;28:453–60. [2] Armﬁeld SW, Fletcher CAJ. Numerical simulation of swirling ﬂow in diffusers. Int J Numer Methods Fluids 1986;6(8):541–56. doi:10.10 02/ﬂd.1650 060805. [3] Armﬁeld SW, Fletcher CAJ. Comparison of k-epsilon and Algebraic Reynolds Stress models for swirling diffuser ﬂow. Int J Numer Methods Fluids 1989;9:987–1009. [4] Azad RS. Turbulent ﬂow in a conical diffuser: a review. Exp Therm Fluid Sci 1996;13:318–37. [5] Azad RS, Kassab SZ. Turbulent ﬂow in a conical diffuser: overview and implications. Phys Fluids A 1989;1:564–73. doi:10.1063/1.857425. [6] Barth T, Jespersen D. The design and application of upwind schemes on unstructured meshes. Aerospace Sciences Meetings; American Institute of Aeronautics and Astronautics; 1989. [7] Batchelor GK. An introduction to ﬂuid dynamics:. Cambridge Mathematical Library; Cambridge University Press; 20 0 0. ISBN 978051180 0955. doi:10.1017/ CBO9780511800955. [8] Bounous O. Studies of the ercoftac conical diffuser with openfoam. Tech. Rep.. Division of Fluid Dynamics, Department of Applied Mechanics, Chalmers University of Technology, Göteborg, Sweden; 2008. [9] Bragg SL, Hawthorne WR. Some exact solutions of the ﬂow through annular cascade actuator discs. J Aeronaut Sci 1950;17(4):243–9. doi:10.2514/8.1597. [10] Cho C-H, Fletcher CAJ. Computation of turbulent conical diffuser ﬂows using a non-orthogonal grid system. Comput Fluids 1991;19:347–61. doi:10.1016/ 0 045-7930(91)90 060-U. [11] Cho N-H. Computation of turbulent swirling ﬂows in conical diffusers with tailpipes. Ph.D. thesis. Department of Mechanical Engineering, The University of Sydney; 1990. [12] Cinnella P, Congedo PM. Numerical solver for dense gas ﬂows. AIAA J 2005;43(11):2458–61. doi:10.2514/1.16335. [13] Cinnella P, Congedo PM. Inviscid and viscous aerodynamics of dense gases. J Fluid Mech 2007;580:179–217. doi:10.1017/S0 0221120 070 05290. [14] Clausen P, Koh S, Wood D. Measurements of a swirling turbulent boundary layer developing in a conical diffuser. Exp Therm Fluid Sci 1993;6:39–48. [15] Clausen P, Wood D. Some measurements of swirling ﬂow through an axisymmetric diffuser. In: Proceedings of 6th Symposium on Turbulent Shear Flows, Toulouse, France; 1987. [16] Congedo P, Corre C, Cinnella P. Numerical investigation of dense-gas effects in turbomachinery. Comput Fluids 2011;49(1):290–301. http://dx.doi.org/10.1016/ j.compﬂuid.2011.06.012. [17] Duprat C. Simulation numérique instationnaire des écoulements turbulents dans les diffuseurs des turbines hydrauliques en vue de l’amélioration des performances. Ph.D. thesis. Laboratoire des Ecoulements Géophysiques et Industriels (LEGI, UMR 5519, CNRS-UJF-INPG), Université de Grenoble, France; 2010. [18] Gatski TB, Speziale CG. On explicit algebraic stress models for complex turbulent model. J Fluid Mech 1993;254:59–78. [19] Gyllenram W, Nilsson H. Very large eddy simulation of draft tube ﬂow. 23rd IAHR Symposium, Yokohama, Japan; 2006. [20] Jakirlic´ S, Hanjalic´ K, Tropea C. Modeling rotating and swirling turbulent ﬂows: a perpetual challenge. AIAA J 2002;40:1984–96. [21] Jeyachandran K, Ganesan V. Numerical modelling of turbulent ﬂow through conical diffusers with uniform and wake velocity proﬁles at the inlet. Math Comput Model 1988;10(2):87–97. http://dx.doi.org/10.1016/0895-7177(88) 90070-2. [22] Klein A. Review: effects of inlet conditions on conical-diffuser performance. J Fluids Eng 1981;103:250–7. [23] Klonowicz P, Borsukiewicz-Gozdur A, Hanausek P, Kryllowicz W, Bruggemann D. Design and performance measurements of an organic vapour turbine. Appl Therm Eng 2014;63:297–303. [24] Kluwick A, Meyer G. Shock regularization in dense gases by viscousâ;;inviscid interactions. J Fluid Mech 2010;644:473–507. doi:10.1017/S0 0221120 09992618. [25] Launder BE, Reece GJ, Rodi W. Progress in the development of a reynolds-stress turbulence closure. J Fluid Mech 1975;68:537–66. doi:10.1017/ S0 0221120750 01814. [26] Lee J., Jang S.J., Sung H.J.. Direct numerical simulations of turbulent ﬂow in a conical diffuser. J Turbul 13, N30. [27] Lemmon EW, Huber ML, McLinden MO. NIST reference ﬂuid thermodynamic and transport properties-REFPROP, version 8.0. Tech. Rep.. Physical and Chemical Properties Division, National Institute of Standards and Technology; 2007. [28] Lemmon WW, Jacobsen RT. An international standard formulation for the thermodynamic properties of 1,1,1-triﬂuoroethane (hfc-143a) for temperatures from 161 to 450 k and pressures to 50 mpa.. J Phys Chem Ref Data 20 0 0;29 (4):521–52.

118

C.S. From et al. / Computers and Fluids 149 (2017) 100–118

[29] Maicke B, Majdalani J. On the compressible bidirectional vortex. part 1: a bragg-hawthorne stream function formulation. In: 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition; 2012a. p. 1103. [30] Maicke B, Majdalani J. On the compressible bidirectional vortex. part 2: a beltramian ﬂowﬁeld approximation. In: 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition; 2012b. p. 1104. [31] McDonald AT, Fox RW, Van Dewoestine RV. Effects of swirling inlet ﬂow on pressure recovery in conical diffusers. AIAA J 1971;9:2014–18. [32] Menter F, Garbaruk A, Egorov Y. Explicit algebraic reynolds stress models for anisotropic wall-bounded ﬂows. EUCASS Proc Ser 2012;3:89–104. [33] Menter FR. Two equation eddy-viscosity turbulence models for engineering applications. AIAA J 1994;32, No. 8:1598–605. [34] Merle X, Cinnella P. Bayesian quantiﬁcation of thermodynamic uncertainties in dense gas ﬂows. Reliab Eng Syst Saf 2015;134:305–23. [35] Okwuobi PAC, Azad RS. Turbulence in a conical diffuser with fully developed ﬂow at entry. J Fluid Mech 1973;57:603–22. [36] Payette F-A. Simulation de l’écoulement turbulent dans les aspirateurs de turbines hydrauliques : Impact des paramètres de modélisation. Master’s thesis. Faculté des Sciences et de Génie, Université de Laval, Québec; 2008. [37] Persky R, Sauret E, Ma L. Optimisation methods for coupled thermodynamic and 1D design of radial-inﬂow turbines. In: Proceedings of the ASME 2014 Joint US-European Fluids Engineering, Division Summer Meeting and 11th International Conference on Nanochannels, Microchannels, and Minichannels, FEDSM 2014, August 3–7, 2014, Chicago, Illinois, USA; 2014. [38] Peters H. Energieumsetzung in querschnittserweiterungen bei verschiedenen zulaufbedingungen. Ingenieur-Archiv 1931;2(1):92–107. doi:10.1007/BF02079816. [39] Pini M, Persico G, Casati E, Dossena V. Preliminary design of a centrifugal turbine for organic rankine cycle applications. J Eng Gas Turbine Power 2013;135:042312–1-042312-9.

[40] Pope SB. A more general effective-viscosity hypothesis. J Fluid Mech 1975;72:331–40. doi:10.1017/S0 0221120750 03382. [41] Pordal H, Khosla P, Rubin S. Pressure ﬂux-split viscous solutions for swirl diffusers. Comput Fluids 1993;22(4–5):663–83. Special Issue Egon Krause Honour Issue. http://dx.doi.org/10.1016/0 045-7930(93)90 031-4 [42] Rhie C, Chow W. Numerical study of the turbulent ﬂow past an airfoil with trailing edge separation. AIAA J 1983;21(11):1525–32. [43] Rotta JC. Statistische theorie nichthomogener turbulenz. J Fluid Mech 1951;129:547–72. [44] Sauret E, Gu Y. Three-dimensional off-design numerical analysis of an organic rankine cycle radial-inﬂow turbine. Appl Energy 2014;135:202–11. doi:10.1016/ j.apenergy.2014.08.076. [45] Sauret E, Persky R, Chassaing J-C, Lucor D. Uncertainty quantiﬁcation applied to the performance analysis of a conical diffuser. In: 19th Australasian Fluid Mechanics Conference, Melbourne, Australia, 8–11 December; 2014. [46] Sauret E, Rowlands AS. Candidate radial-inﬂow turbines and high-density working ﬂuids for geothermal power systems. Energy 2011;36:4460–7. [47] Senoo Y, Kawaguchi N, Nagata T. Swirl ﬂow in conical diffusers. Bull Japan Soc Mech Eng JSME 1978;21:112–19. [48] Senoo Y, Nishi M. Improvement of the performance of conical diffusers by vortex generators. J Fluids Eng 1974;96(1):4. [49] Wallin S, Johansson AV. An explicit algebraic reynolds stress model for incompressible and compressible turbulent ﬂows. J Fluid Mech 20 0 0;403:89–132. doi:10.1017/S0 0221120990 070 04. [50] Wheeler APS, Ong J. The role of dense gas dynamics on organic rankine cycle turbine performance. J Eng Gas Turbines Powerngineering 2013;135:102603–1-102603-9. [51] Wilcox D. Turbulence modelimg for cfd. third ed. DCW Industries; 2006. [52] Wu X, Schlüter J, Moin P, Pitsch H, Iaccarino G, Ham F. Computational study on the internal layer in a diffuser. J Fluid Mech 2006;550:391–412.