Turbulent dense gas flow characteristics in swirling conical diffuser

Turbulent dense gas flow characteristics in swirling conical diffuser

Computers and Fluids 149 (2017) 100–118 Contents lists available at ScienceDirect Computers and Fluids journal homepage: www.elsevier.com/locate/com...

4MB Sizes 0 Downloads 11 Views

Computers and Fluids 149 (2017) 100–118

Contents lists available at ScienceDirect

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

Turbulent dense gas flow characteristics in swirling conical diffuser C.S. From a, E. Sauret a,∗, S.W. Armfield 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 flows 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 efficiency. This increase of efficiency is critical for the overall cycle efficiency of renewable power cycles based on low temperature renewable resources. Optimising the performance of a conical diffuser in renewable power cycles using high-density fluids 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 firstly 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 flow behaviour for IG while being ineffective on the RG. The final application analysed the diffuser performance using the inlet conditions extracted directly from a potential radial-inflow turbine working with R143a. The change of conditions highlighted that radial components can be reduced, and thus the swirling number too. By implementing the first numerical study on real gas swirling conical diffuser, it was established that real gas flow regimes differ from the ones previously established for ideal gas, and thus preliminary flow regimes for R143a, specifically, 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 fluids to power, in most cases, a radial-inflow turbine. The increase in turbine efficiency can dramatically improve the overall cycle efficiency, 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 flow through the conical diffuser that must now be investigated. The turbine component is critical in achieving high cycle efficiency 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 flow through the conical diffuser. The diffuser is responsible for the pressure recovery at the exit of the turbine and thus can significantly contribute to a gain of



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

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

overall efficiency. However, this pressure recovery can be affected by the separation or recirculation generated by the flow 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 flow inside the diffuser and the fluid 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. Armfield and Fletcher [2,3] established that the effect of swirling fluid flows compared to non-swirling flows in conical diffusers is to produce a positive radial pressure gradient in the core flow at the inlet, preventing flow detaching from the walls as

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

the flow diffuses within the conical diffuser. As the flow 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 flow constraint, the axial velocity near the wall is accelerating, preventing separation. Hence swirling flows aim to ensure sufficient 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] defined the boundary separation to be driven by the mostly viscous layer and recirculation in the core flow as an effect of the essentially inviscid flow. 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 identified 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 flow separation. Turbulent velocity profiles 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 flow 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 flow regime that affects the performance of the diffuser. Senoo et al. [47] experimentally studied five diffusers with different free-vortex inlet swirling flow 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 flow 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 flow reversal. Mathematical simplifications have been made to axisymmetric Navier Stokes equations, known as the Bragg–Hawthorne equations [9] to describe the inviscid process in swirling flows at high Reynolds numbers. These analytical methods, commonly known as vorticity-streamfunction methods were originally limited to the application of inviscid and incompressible flows [7,9]. That is until the recent extension of the BH equations to compressible axisymmetric swirling flows, made by Maicke and Majdalani [29,30]. Despite these advances, these analytical approaches are limited to inviscid flows and thus can’t reproduce boundary layer separation without the addition of corrective terms. This is an important limitation for dense gas flows in swirling conical diffusers. Modelling swirling and rotating turbulent flows also remains challenging as stated by Jakirlic´ et al. [20] due to the complex flow structures which cannot be accurately reproduced by turbulence models commonly established and validated for simple flows. Lee et al. [26] also recently added that conical diffusers under adverse pressure gradient are among the most difficult 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 flow features. For example, Armfield 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 configurations 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 flows and recommended the use of an Explicit Algebraic Reynolds Stress Model (EARSM) instead for three-dimensional rotating flows. 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 fluid 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 inflow 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 flow structures encountered in swirling turbulent conical diffusers remains a challenge. In addition, the inflow conditions are critical for the accuracy of the numerical model. Furthermore, all the existing studies consider “common” fluids in liquid or gas phases such as air or water while ORC turbines work using high-density fluids, commonly referred to as dense or real gases, such as refrigerants. These dense fluids 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” fluids [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 flows to date [34]. This is most likely due to the large variety of available real fluid mixtures and their limitations as well as feasible applicability to various practical applications that potentially benefit 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 fluid flow [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 configurations compared to ideal gas? Will the inflow conditions, e.g. turbulence and swirl, affect the real gas flow 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 configurations?

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 inflow conditions on the performance of the ERCOFTAC swirling conical diffuser using air and refrigerant gas R143a (real gas). Air used within this study satisfies 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-Trifluoroethane). NIST REFPROP model is considered as the international standard for the fluid properties of R143a [27]. This model utilises high-fidelity 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-inflow 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 perhydrofluorene (PP10) flow 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 flow. 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 flows are implemented in the finite 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 β coefficients are functions of the following five 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 coefficients 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 flows, 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 coefficients 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 β coefficients 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 specific 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 defined 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 fluid, 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 flows 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 flow 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 fluid domain based on length scale, maximum velocity scale in the domain, density, dynamic viscosity, mass flow, 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 coefficients 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 refined. Furthermore, in order to predict transport properties due to compressible effects of the real gas, the energy within the fluid

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 clarification 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 profiles 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 defined as the axis along and the axis normal to the diffuser wall respectively in which the experimental profiles 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 coefficient, CP defined 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 defined 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 fluid density and Uo the mean reference inflow velocity.

The construction of the computational mesh followed the work of Bounous [8], where the mesh was divided into five 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 flow-separation at the walls. The mesh was thus refined 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 five longitudinal sections with nx# grid points in the longitudinal direction for each section. In total three grids with increasing refinement 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 coefficient 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 refined enough to provide accurate results for this configuration, both in the viscous and inviscid regions of the flow. Grid 2 is presented in Fig. 4 and has a total of 920,0 0 0 nodes with the y+ < 2 requirement satisfied 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 refinements in the near-wall region. Firstly grid 2 was modified by varying the progression ratio in order to decrease the first 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 refinements were performed by firstly 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 first near-wall distance to yw = 0.5e−6 m. Results indicated no improvements to the simulation by either of the refined meshes. The standard grid which has a total of 920,0 0 0 nodes with modified progression ratio and yw = 1.1e−6 m, was thus considered as a suitable grid for the study of dense gas turbulence swirling flows. 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 flow 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 inflow 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 fixed 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] defined 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, defined 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 define 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 profiles, 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 inflow turbulence conditions were not provided in the experimental data and thus the inlet turbulence is defined and controlled by specifying the turbulence intensity, T¯u , where T¯u is taken as 5% based on Armfield 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 specific dissipation √ 3 rate, ω, defined 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 justifies the method of defining 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 profiles (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). Specifically, this indicates

that the BSLkω EARSM accurately predicts interactions between the viscous near-wall flow and the essentially inviscid core flow, throughout the diffuser as the vorticity is (inviscidly) stretched. This accurate prediction of near-recirculating flows is an essential component for turbulent swirling flows in conical diffusers since this will proportionally predict flow 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 profiles are considerably more accurate with the BSLkω EARSM model for S4 and S7. The turbulent kinetic energy profiles 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 significant impact on the flow. 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 flow characteristics through the diffuser. This also confirms 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 first 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 flow characteristics throughout the ERCOFTAC conical diffuser, providing good agreement with the experimental data for both the near-wall (viscous) and centreline profiles (inviscid) for all parameters at S0, S4 and S7. Despite far-wall kinetic energy profiles 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 flow, despite the sensitivity. This further supports the importance of near-wall parameters to accurately predict turbulent swirling flow 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 flows 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 fluid has been shown to be a suitable candidate for the development of efficient 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 fluid at the rotor exit of the turbine, the inlet temperature, T, and outlet pressure, Pout , are then modified. 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 five 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 profiles 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 profiles 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 first 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 flow 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 confirmed 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 magnified 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 profiles RG1 and RG2. This suggests that the effect of density does not directly correspond to different flow characteristics in conical diffuser where a solid-body rotation is applied. It is, however, obvious that the flow conditions differ with the different fluid 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 flows 4.2. Evolution of the boundary layer parameters in the diffuser In a first attempt to establish a direct comparison between the two fluids (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 influences the flow 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 fluids 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 fluids. 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 confirming 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 fluids, 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 confirm 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 profiles differ between the two fluids from S0, leading us to investigate the effect of the parameters at the diffuser inlet(i) in the next Section 4.3 and the inflow turbulence in the later Section 4.4 in order to establish the influence of inflow components.

4.3. Characteristics of the boundary layer and velocity profiles 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 fluid flow 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 influence of the flow characteristics developed within the diffuser for both fluids, flow 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 flow parameters, including the swirl number, Sn , as defined in Eq. (24) in order to quantify the strength of swirling flows 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 flow generates a positive radial profile 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 profiles of RG at S0 establishes that there is dependence on the inlet, i, flow conditions. To compare how IG, RG1 and RG2 differ at the inlet(i), Mach number Ma , mass flow rate M˙ and boundary layer components are detailed in Table 4. Comparing IG and RG2 indicates that the density ρ and similarly the mass flow 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 fluid 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 flow 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 profile, 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



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 profiles, 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 flow 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 flow 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 flow rate for both IG and RG2 will be presented in Section 5. However, varying the flow rate will affect the swirl through a change of inlet axial velocity and thus is not considered here.

Fig. 15. Velocity profiles 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 inflow turbulence Investigations at the entrance of the diffuser (inlet(i) ) are conducted in an attempt to explain the large variance in flow characteristics seen so far between IG and RG. It is known that the turbulence quantities have an important impact on the swirling flows conditions within the diffuser. As detailed by Armfield and Fletcher [3], turbulent velocity profiles 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 flows with high turbulent intensity against flows 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 inflow conditions. However, the actual production throughout this rotational domain is not known. Therefore, to investigate the effect of the turbulence production for both fluids 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 profiles towards the centreline, as seen in Fig. 14. Within the diffuser, where the flow 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 flow profiles 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 sufficient quantity to benefit from any variation of turbulence intensity in the fluid inflow. Looking at the axial velocity profile 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 flow 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 benefits 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 profiles 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 inflow conditions. This is potentially due to the excessive increase in density - induced rotational momentum which hence leads to no changes when the inflow 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 flow 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 flow 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 fluid properties, dynamic viscosity and density, modifies the characteristics of swirling fluid flows within the conical diffuser. The difference in the radial velocity component Ur for both fluids appears to be critical, since evidence indicates that the significant increase in the radial components (e.g. turbulence components, radial profiles, momentum) is due to the density-induced rotational momentum and hence the radial component is representative of that effect. More robust wall flows are produced by dense gas R143a, although consequent to this is a large reduction in the centreline axial velocities to the extent that recirculating flows 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 flow 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 sufficient

113

impact on the fluid, with a density beyond that of RG1, to have any real significant effect on the fluid flow. This can also explain the absence of a change in turbulence intensity (T¯u ) of the flow of RG2 compared to IG. This indicates that conical diffusers utilized in ORC applications may have less sensitivity to inflow 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. Inflow swirl effects 5.1. Flow characteristics within the diffuser at various Swirl numbers As discussed, the use of high density fluids 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 fluid could potentially help reduce the tendency for recirculation at the centreline. It is known from both literature [19,36,45] and findings 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 profiles due to the density-induced radial components discussed in Section 4. As such, the possibility to establish an inlet(i) radial profile nearly identical to that of the IG case is desired at the same flow rate. In order to match the inlet(i) radial profiles 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) profiles for IG and RG2, at standard angular velocities,  = 56.1 rads−1 and modified with the lower angular velocity setting of  = 36.1 rads−1 . Hence the modified real gas case RG∗ matches the inlet(i) radial velocity profile Ur of IG and the modified IG∗ now has a considerably lower Ur profile than the standard IG. The respective swirl number Sn of standard cases are; SnIG = 0.3, SnRG2 = 0.278 and modified 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 profile, where the peak of IG∗ is scaled slightly more towards the wall. Near-identical kinetic energy k profiles are displayed between the standard and modified cases for the respective fluids at S0. However, due to the lower strength induced by the solid body-rotation, the peaks of k have lower magnitudes for the modified cases. The radial velocity profile 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 profiles 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 confirms 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 , profiles of Standard and Modified, 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 flow 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 profiles that the increase in Ur from S0 to S4 is still sufficiently 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 confirms 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 coefficient in Fig. 18, defined by Eq. (23). It is clear that the modified IG∗ has the lowest performance, which is most likely due to the separated flow 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 Modified∗ 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 modified 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. Modified IG∗ has a shift in profile at approximately xw = 0.13 m, which corresponds to the re-circulating flow at the wall (e.g. separation occurs between S3 and S4). The recirculating flow 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 flow detaching from the wall entirely. The modified 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 modified cases do not benefit from this recovering wall shear past S4, which is the result of lower solid-body rotation decreasing wall-components as the flows diffuses. Formulation for the kinetic energy production term, Pk is estimated by Eq. (32), which follows the simplified 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 modified: −−, 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 fluids 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 modified cases are displayed in Fig. 21(a) and (b), respectively. The momentum thickness indicates a continuous increase for both modified 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 modified cases are overall lower, however they do not display any distinctive difference in characteristics from their corresponding standard cases, essentially mirroring the profiles 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 fluid, 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 flows. 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 profile at the diffuser even at lower rotational speed resulting from the higher density. This indicates that real gas and ideal gas flow regimes are different in swirling conical diffusers as detailed in the following section. 6. Flow regimes Armfield and Fletcher [3] characterized the swirling flows in the conical diffusers into three categories: weak swirl, moderate swirl and strong swirl depending on the swirl number (defined by Eq. (24)), conical half-angle,  and flow type which consists

116

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

Fig. 21. Standard and modified∗ 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 flow quality between standard and modified∗ with half-cone angle of  = 10◦ . ERCOFTAC (IG)

IG∗

RG2

RG∗

Ideal flow with some reduction at centreline

Strong separation at wall

Ideal flow, with negligible recirculation at the exit

Ideal flow, with near separation at wall

Table 7 Preliminary proposal of R143a swirling flow 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 modified conditions are considered elliptic flows, parabolic flows 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 flow regimes of the table detailed in Armfield and Fletcher [3]. i.e. no separation at the walls but near-reversing flow 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 flow 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 flow, firstly 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 fluid presented greater wall shear strength, despite the lower friction velocity uf . Considering these factors further confirmed the effect of higher-density fluid 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 flow 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 flow characteristics for RG. This could potentially also provide reasoning as to how RG1 and RG2 display such near-identical profiles despite the dramatic difference in density. These speculations are yet to be assessed further in greater detail, by using various real-gas fluids 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 flows. If future findings are consistent with the results presented here that will confirm that the inlet turbulence intensity is not one of the main factors affecting the swirling flow characteristics for turbulent high-density fluids, unlike ideal gas swirling flows. More detailed investigations are needed to fully understand the turbulent flow 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 influence of inflow swirl. It was shown that a large decrease of the angular rotation was required to match the radial velocity profile 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 flow 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 confirmed that higher-density fluids induce greater positive radial velocity components even when subject to solid-body rotation, and thus benefits from lower swirl intensities when compared at the ideal gas assumption. This confirms previous speculations that overall the high-density fluids provide greater and more robust wall components at lower swirl intensities. As such it is evident that real gas flows 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 defining the optimum working conditions, such as the geometry, inlet velocities, turbulence intensities and other factors for high-density swirling flows in conical diffuser for specific density and dynamic viscosity ranges is critical to develop robust guidelines for conical diffuser applied to ORC turbines. Thus, the improved understanding of the flow 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 sufficient swirl level through an increase of the axial velocity in order to establish favourable flow characteristics for the conical diffuser. In return, this will potentially improve the overall ORC turbine efficiency. In conclusion this paper presents a preliminary numerical investigation of real gas swirling turbulent flows 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 flows. The core issue, and one of the most important aspects to address is whether the turbulent real gas flow phenomena can be accurately, or even sufficiently, predicted with numerical models coupled with EOS models that are based on high-fidelity 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 financial support (DE130101183). References [1] Armfield SW, Cho C-H, Fletcher CAJ. Prediction of turbulence quantities for swirling flow in conical diffusers. AIAA J 1990;28:453–60. [2] Armfield SW, Fletcher CAJ. Numerical simulation of swirling flow in diffusers. Int J Numer Methods Fluids 1986;6(8):541–56. doi:10.10 02/fld.1650 060805. [3] Armfield SW, Fletcher CAJ. Comparison of k-epsilon and Algebraic Reynolds Stress models for swirling diffuser flow. Int J Numer Methods Fluids 1989;9:987–1009. [4] Azad RS. Turbulent flow in a conical diffuser: a review. Exp Therm Fluid Sci 1996;13:318–37. [5] Azad RS, Kassab SZ. Turbulent flow 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 fluid 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 flow 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 flows 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 flows 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 flows. 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 flow 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.compfluid.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 flow. 23rd IAHR Symposium, Yokohama, Japan; 2006. [20] Jakirlic´ S, Hanjalic´ K, Tropea C. Modeling rotating and swirling turbulent flows: a perpetual challenge. AIAA J 2002;40:1984–96. [21] Jeyachandran K, Ganesan V. Numerical modelling of turbulent flow through conical diffusers with uniform and wake velocity profiles 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 flow in a conical diffuser. J Turbul 13, N30. [27] Lemmon EW, Huber ML, McLinden MO. NIST reference fluid 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-trifluoroethane (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 flowfield 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 flow 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 flows. 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 quantification of thermodynamic uncertainties in dense gas flows. Reliab Eng Syst Saf 2015;134:305–23. [35] Okwuobi PAC, Azad RS. Turbulence in a conical diffuser with fully developed flow 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-inflow 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 flux-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 flow 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-inflow 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 quantification 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-inflow turbines and high-density working fluids for geothermal power systems. Energy 2011;36:4460–7. [47] Senoo Y, Kawaguchi N, Nagata T. Swirl flow 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 flows. 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.