Design of Radial-Inflow Turbines for Low-Temperature Organic Rankine Cycle

Abstract

This study presents the development of a design method that has been extended to the design of radial-inflow turbines operating in organic Rankine cycles (ORC). Both the conventional design method and the circulation method available in the literature have been reviewed. The two main limitations of the current circulation method that make it not suitable for the ORC turbine design are the lack of real gas capability and 3D blades with high stresses. Using the circulation method, the flow field is decomposed into a potential part and a rotational part. The mean velocity field and the periodic velocity field are solved separately. To model the thermodynamic properties of the real gas, NIST REFPROP or CoolProp are used. The blade geometry is then solved iteratively by assuming that the velocity vector is parallel to the blade surface. The blade boundary condition is modified to force the blade camber to be radial-fibred, which is helpful to reduce the centrifugal bending stress on the blade. All the formulations are derived step by step, and the numerical treatments, including grid generation, numerical differentiation, computational scheme, and convergence, are discussed in detail. This method is validated by designing a R245fa ORC turbine rotor. The performance of the rotor design is predicted by CFD and FEA simulations, and it is compared to the results using other methodologies in the literature.

1. Introduction

The organic Rankine cycle is a special variation of the Rankine cycle, which is named for its use of organic, high-molecular-mass fluids whose vaporization temperature is lower than that of water. The main advantage of ORC is that it is suitable for extracting energy from low-grade heat sources compared to the conventional Rankine cycle. One of the earliest research projects about ORC was carried out by Tabor and Bronicki in the 1960s [1]. They designed and tested several small solar ORC systems with monochlorobenzene as the working fluid, and the system power outputs were 2 kW and 10 kW. Since then, ORC has been widely used for heat recovery from low-temperature sources, including biomass combustion [2,3,4], geothermal heat [5,6], industrial waste heat [7,8], and solar thermal heat [9,10,11,12]. The working fluid used in the ORC system plays a very important role in the overall cycle efficiency. Hung et al. [13] compared the thermophysical properties of wet and dry working fluids used in the ORC, while the fluids were characterized based on their slope of saturation vapor diagram. Dai et al. [14] compared different working fluids for high-temperature ORC systems, and Usman et al. [15] performed a comparative study on the low-medium-temperature geothermal ORC system using R245fa and R1233zde. Pang et al. [16] carried out an experiment to compare using R245fa, R123, and their mixtures with varying mass fractions in the ORC systems. Eyerer et al. [17] experimentally investigated modern ORC working fluids R1224yd(Z) and R1233zd(E) as replacements for R245fa.

Turbines or expanders are another important part of the ORC system. Radial-inflow turbines, scroll expanders, and screw expanders are typically used. The most favorable or widely used turbine nowadays is the radial-inflow turbine due to its light weight and high efficiency. However, the cost of radial-inflow turbines will be higher than others. Different researchers have tried to design and use radial-inflow turbines in the ORC [18,19,20,21,22]. The method they used for the design of radial-inflow turbines is called the conventional method. The initial leading edge (LE) and trailing edge (TE) blade angles are calculated based on 1D meanline analysis and velocity triangles. A smooth cubic spline or Bezier curve would be used to define the blade angle or wrap angle distribution along the streamwise direction. The wrap angle at other spanwise locations (from hub to shroud) will be obtained through radial filament, which minimizes the blade’s centrifugal bending stress. The 3D blade geometry is constructed by lofting all the sections together, and CFD simulations or experiments will be used to evaluate the turbine’s performance.

A different method called the circulation method (or inverse design method) was first proposed by Hawthorne et al. [23], and they used it to design infinitely thin 2D cascades for inviscid and incompressible flow. The cascade blade is modeled as a single vorticity sheet using this method, and the blade shape is calculated iteratively by using the boundary condition that the flow is assumed to be parallel to the blade surface. Tan et al. [24] extended this method to the design of 3D cascades with constant hub and shroud diameters. One limitation of their method is that the blade is still assumed to be infinitely thin (zero blade thickness). Borges [25] extended Tan’s method further to the design of low-speed turbomachinery with arbitrary meridional geometry, which makes it possible to design centrifugal and mixed-flow pumps using this method. Ghaly [26] added flow compressibility to this method and used it to design radial-inflow turbines with compressible and subsonic flow. The blade thickness effect was then introduced by Zangeneh [27] using a blockage factor in the mean flow continuity equation. Since then, the circulation method has been widely used for the design of a very wide range of turbomachinery blades, including compressors [28], fans [29], pumps [30], and hydraulic turbines [31].

To apply the existing circulation method to the design of radial-inflow turbines for ORC, there are two main limitations. The first is that the flow has been assumed to be incompressible or compressible ideal gas, while the ORC working fluid is real gas. Therefore, the ideal gas assumption does not hold any more. Colonna et al. [32] investigated the real gas effect in ORC turbines by comparing the polytropic ideal gas law, the Peng–Robinson–Stryjek cubic equation of state (EoS) and the Span–Wagner EoS. The second limitation is that the blade designed using the circulation method has quite high stress due to its complex 3D blade shape. To control the stress within allowable material strength, post-modifications, including modifying the blade thickness or wrap angle distribution (radial filament modification), are required to convert the blade to be radial-filament or radial-fibred blading.

In this paper, a design method that can be used to design radial-inflow turbines with real gas (organic) fluids and generate radial-filament geometries directly is proposed. The performance of the radial-inflow turbine will be evaluated and validated through 3D CFD and FEA simulations.

2. Methodology and Governing Equations

2.1. Basic Assumptions

First, the following assumptions are made:

• The flow is assumed to be adiabatic, steady, and inviscid;
• The fluid is compressible and a real gas;
• The blade is assumed to be infinitely thin with zero blade thickness;
• The flow near TE is free vortex, which means there is no vortex shedding after TE.

Using the first assumption with Kelvin’s theorem, it can be concluded that the only vorticity in the flow field is the blade-bound vorticity, which is bounded on all the blade surfaces. Because of the assumption of infinitely thin blades, one vorticity sheet can be used to represent the blade geometry. A blockage factor is used in the mean flow continuity equation to account for the blockage effect caused by the blade thickness.

The blade camber (or wrap angle) can be defined by:

$$\alpha \left(r,\theta ,z\right)=\theta -f\left(r,z\right)=\pm n\frac{2\pi }{B}$$

(1)

where α is all the blade surfaces, r, θ, and z are cylindrical coordinates, f is the blade camber, and B is the blade number. It should be noted that f represents the wrap angle of one single blade, while $\alpha$ represents the wrap angles of all the blades.

The overbar symbol ‘−’ is used to define a way of averaging along the circumferential direction, for example:

$$\overline{A}\left(r,z\right)=\frac{2\pi }{B}{\int }_0^{\frac{2\pi }{B}}A\left(r,\theta ,z\right)d\theta$$

(2)

To model the periodic flow properties, a Dirac delta function is used:

$${\sum }_{n=-\infty }^{\infty }\delta \left(\alpha -n\frac{2\pi }{B}\right)=\left\{\begin{array}{c}+\infty \alpha =n\frac{2\pi }{B}\\ 0 \alpha \ne n\frac{2\pi }{B}\end{array}\right.$$

(3)

It can be easily seen that the value of the delta function is infinite on the blade and remains zero anywhere else. Equation (3) can be rewritten in the form of a Fourier series (Lighthill [33]):

$${\delta }_P\left(\alpha \right)=\frac{2\pi }{B}{\sum }_{n=-\infty }^{\infty }\delta \left(\alpha -n\frac{2\pi }{B}\right)={\sum }_{n=-\infty }^{\infty }{e}^{inB\alpha }$$

(4)

The vorticity field can then be defined by Equation (5). It should be noted that the vorticity is periodic along the circumferential direction and only holds non-zero values on all the blade surfaces.

$$\overrightarrow{\Omega }\left(r,\theta ,z\right)=\overline{\overrightarrow{\Omega }}\left(r,z\right){\delta }_P\left(\alpha \right)$$

(5)

2.2. Monge–Clebsch Decomposition

The Monge–Clebsch decomposition is used to decompose the velocity field into the summation of a potential term and a rotational term (Wu et al. [34]).

$$\overrightarrow{V}\left(r,\theta ,z\right)=\nabla \Phi \left(r,\theta ,z\right)+\lambda \left(r,\theta ,z\right)\nabla \mu \left(r,\theta ,z\right)$$

(6)

where ∇Φ is the potential term and λ∇μ is the rotational term. Taking the curl of Equation (6) to obtain the vorticity:

$$\overrightarrow{\mathsf{\Omega }}\left(r,\theta ,z\right)=\nabla \lambda \left(r,\theta ,z\right)\times \nabla \mu \left(r,\theta ,z\right)$$

(7)

One of the Clebsch variables can be assumed to be α. Therefore, Equation (7) can be rewritten as:

$$\overrightarrow{\Omega }\left(r,\theta ,z\right)=\left[\frac{\partial \lambda \left(r,\theta ,z\right)}{\partial r},\frac{\partial \lambda \left(r,\theta ,z\right)}{r\partial \theta },\frac{\partial \lambda \left(r,\theta ,z\right)}{\partial z}\right]\times \left[\frac{-\partial f\left(r,z\right)}{\partial r},\frac{1}{r},\frac{-\partial f\left(r,z\right)}{\partial z}\right]$$

(8)

The mean vorticity field $\overrightarrow{\Omega }$can be obtained by taking the circumferential average of both sides of Equation (8).

$$\overline{\overrightarrow{\Omega }}\left(r,z\right)=\left[-\frac{\partial \overline{\lambda }\left(r,z\right)}{r\partial z},\frac{\partial \overline{\lambda }\left(r,z\right)}{\partial r}\frac{\partial f\left(r,z\right)}{\partial z}-\frac{\partial \overline{\lambda }\left(r,z\right)}{\partial z}\frac{\partial f\left(r,z\right)}{\partial r},\frac{\partial \overline{\lambda }\left(r,z\right)}{r\partial r}\right]$$

(9)

Taking the circumferential average of the velocity curl to obtain the mean vorticity $\overline{\overrightarrow{\Omega }}$:

$$\overline{\overrightarrow{\Omega }}\left(r,z\right)=\left[-\frac{\partial {\overline{rV}}_{\theta } \left(r,z\right)}{r\partial z},\frac{\partial {\overline{V}}_r \left(r,z\right)}{\partial z}-\frac{\partial {\overline{V}}_z \left(r,z\right)}{\partial r},\frac{\partial {\overline{rV}}_{\theta } \left(r,z\right)}{r\partial r}\right]$$

(10)

By comparing Equations (9) and (10), it can be easily seen that:

$$\overline{\lambda }\left(r,z\right)={\overline{rV}}_{\theta } \left(r,z\right)$$

(11)

Equations (5) and (9) can then be rewritten as:

$$\overrightarrow{\Omega }\left(r,\theta ,z\right)=\left[\nabla {\overline{rV}}_{\theta }\left(r,z\right)\times \nabla \mathsf{\alpha }\left(r,\theta ,z\right)\right]{\delta }_P\left(\alpha \right)$$

(12)

$$\overline{\overrightarrow{\mathsf{\Omega }}}\left(r,z\right)=\nabla {\overline{\mathrm{r}\mathrm{V}}}_{\mathsf{\theta }}\left(r,z\right)\times \nabla \mathsf{\alpha }\left(r,\theta ,z\right)$$

(13)

A periodic sawtooth function $S\left(\alpha \right)$ can be defined by Equations (14) and (15).

$$S^{\prime }\left(\alpha \right)={\delta }_P\left(\alpha \right)-1={\sum }_{n=-\infty }^{\infty }{e}^{inB\alpha }-1=2{\sum }_{n=1}^{\infty }{e}^{inB\alpha }$$

(14)

$$\mathrm{S}\left(\alpha \right)=2\sum _{n=1}^{\mathrm{\infty }}\frac{e^{inB\alpha }}{inB}+C$$

(15)

Using Equations (12), (14) and (15), the absolute vorticity field $\overrightarrow{\Omega }$ can be written in the following form:

$$\begin{gathered} \nabla \times \overrightarrow{V}\left(r,\theta ,z\right)=\overrightarrow{\Omega }\left(r,\theta ,z\right) \\ =\nabla \times \left[\nabla \Phi \left(r,\theta ,z\right)+{\overline{rV}}_{\theta }\left(r,z\right)\nabla \mathsf{\alpha }\left(r,\theta ,z\right)-S\left(\alpha \right)\nabla {\overline{\mathrm{r}\mathrm{V}}}_{\mathsf{\theta }}\left(r,z\right)\right] \end{gathered}$$

(16)

The velocity field $\overrightarrow{V}$ can then be expressed by the following equation.

$$\overrightarrow{V}\left(r,\theta ,z\right)=\nabla \Phi \left(r,\theta ,z\right)+{\overline{rV}}_{\theta }\left(r,z\right)\nabla \mathsf{\alpha }\left(r,\theta ,z\right)-S\left(\alpha \right)\nabla {\overline{\mathrm{r}\mathrm{V}}}_{\mathsf{\theta }}\left(r,z\right)$$

(17)

The velocity field can also be expressed by the summation of a mean term and a periodic term.

$$\overrightarrow{\mathrm{V}}\left(r,\theta ,z\right)=\overline{\overrightarrow{V}}\left(r,z\right)+\tilde{\overrightarrow{V}}\left(r,\theta ,z\right)$$

(18)

Comparing Equations (17) and (18), the equations for the mean flow velocity $\overline{\overrightarrow{V}}$ and the periodic flow velocity $\tilde{\overrightarrow{V}}$ can be obtained.

$$\overline{\overrightarrow{V}}\left(r,z\right)=\nabla \overline{\Phi }\left(r,z\right)+{\overline{rV}}_{\theta }\left(r,z\right)\nabla \mathsf{\alpha }\left(r,\theta ,z\right)$$

(19)

$$\tilde{\overrightarrow{V}}\left(r,\theta ,z\right)=\nabla \tilde{\Phi }\left(r,\theta ,z\right)-S\left(\alpha \right)\nabla {\overline{\mathrm{r}\mathrm{V}}}_{\mathsf{\theta }}\left(r,z\right)$$

(20)

where $\overline{\Phi }$ is the mean flow potential function and $\tilde{\Phi }$ is the periodic flow potential function. The last terms on the RHS of Equations (19) and (20) are equal to zero in the inlet and outlet regions since the flow is irrotational.

2.3. Continuity Equations

The typical steady and compressible flow continuity equation can be written in the following form:

$$\begin{gathered} \nabla ·\left[\rho \left(r,\theta ,z\right)\overrightarrow{W}\left(r,\theta ,z\right)\right]=0 \\ \nabla ·\overrightarrow{W}\left(r,\theta ,z\right)=-\overrightarrow{W}\left(r,\theta ,z\right)·\nabla \mathrm{l}\mathrm{n}\rho \left(r,\theta ,z\right) \end{gathered}$$

(21)

To obtain the mean flow continuity equation, we can take the circumferential average on both sides of Equation (21).

$$𝛻·\overline{\overrightarrow{W}}\left(r,z\right)=-\overline{\overrightarrow{W}\left(r,\theta ,z\right)·𝛻ln\rho \left(r,\theta ,z\right)}$$

(22)

The relative velocity $\overrightarrow{W}$ can be written in the following form:

$$\overrightarrow{W}\left(r,\theta ,z\right)=\overline{\overrightarrow{W}}\left(r,z\right)+\tilde{\overrightarrow{V}}\left(r,\theta ,z\right)$$

(23)

To obtain the periodic flow continuity equation, we can take the divergence of Equation (23) and compare it with Equations (21) and (22).

$$𝛻·\tilde{\overrightarrow{V}}\left(r,\theta ,z\right)=-\overrightarrow{W}\left(r,\theta ,z\right)·𝛻ln\rho \left(r,\theta ,z\right)+\overline{\overrightarrow{W}\left(r,\theta ,z\right)·𝛻ln\rho \left(r,\theta ,z\right)}$$

(24)

2.4. Mean Flow Equations

2.4.1. Blockage Effect

The blade thickness at the hub section of a typical radial-inflow turbine is quite large, and its blockage effect should not be neglected. A blockage factor $B_f$ is introduced into the continuity equation for the mean flow to account for the blockage effect caused by the finite blade thickness.

$$𝛻·\left[{{\rho }_m\left(r,z\right)B}_f\left(r,z\right)\overline{\overrightarrow{W}}\left(r,z\right)\right]=0$$

(25)

where ${\rho }_m$ is a special artificial density, which will be discussed in the next subsection. $B_f$ is defined as:

$$B_f\left(r,z\right)=1-\frac{t_{\theta }\left(r,z\right)B}{2\pi r}$$

(26)

where $t_{\theta }$ is the tangential thickness distribution for the blade and B is the blade number. The relation between the blade’s normal thickness and the tangential thickness is defined by Equation (27). It is more common to directly define the blade’s normal thickness for stress considerations.

$$t_{\theta }\left(r,z\right)=t_N\left(r,z\right)\sqrt{1+r^2{\left[\frac{\partial f\left(r,z\right)}{\partial m}\right]}^2}$$

(27)

2.4.2. Artificial Density

A special artificial density ${\rho }_m$ is defined (Equation (25)) so that the mean flow velocities can be solved through a Stokes stream function for 3D axisymmetric flow. ${\rho }_m$ can be calculated through Equations (22) and (25):

$$\overline{\overrightarrow{W}}\left(r,z\right)·𝛻ln{\rho }_m\left(r,z\right)=\overline{\overrightarrow{W}\left(r,\theta ,z\right)·𝛻ln\rho \left(r,\theta ,z\right)}-\overline{\overrightarrow{W}}\left(r,z\right)·𝛻ln{B}_f\left(r,z\right)$$

(28)

To simplify equations later, ${\rho }_m$ is normalized by a constant reference density ${\rho }_i$ (e.g., inlet stagnation density) to obtain ${\rho }_a$:

$${\rho }_a\left(r,z\right)=\frac{{\rho }_i}{{\rho }_m\left(r,z\right)}$$

(29)

2.4.3. Stokes Stream Function

To solve the mean flow continuity equation (Equation (25)), a stream function Ψ for axisymmetric flow can be defined:

$$\begin{gathered} \overline{W_r}\left(r,z\right)=\overline{V_r}\left(r,z\right)=-\frac{{\rho }_a\left(r,z\right)}{r{B}_f\left(r,z\right)}\frac{\partial \Psi \left(r,z\right)}{\partial z} \\ \overline{W_z}\left(r,z\right)=\overline{V_z}\left(r,z\right)=\frac{{\rho }_a\left(r,z\right)}{r{B}_f\left(r,z\right)}\frac{\partial \Psi \left(r,z\right)}{\partial r} \end{gathered}$$

(30)

The stream function Ψ satisfies the mean flow continuity equation Equation (25) automatically. Another equation is required to solve the unknown stream function Ψ. Comparing Equations (9) and (10), it can be concluded that:

$$\frac{\partial {\overline{V}}_r \left(r,z\right)}{\partial z}-\frac{\partial {\overline{V}}_z \left(r,z\right)}{\partial r}=\frac{\partial {\overline{rV}}_{\theta }\left(r,z\right)}{\partial r}\frac{\partial f\left(r,z\right)}{\partial z}-\frac{\partial {\overline{rV}}_{\theta }\left(r,z\right)}{\partial z}\frac{\partial f\left(r,z\right)}{\partial r}$$

(31)

Substituting Equation (30) into Equation (31) with some rearrangements, the stream function Ψ can be expressed in the following form:

$$\begin{gathered} \frac{{\partial }^2\Psi }{\partial r^2}-\frac{\partial \Psi }{r\partial r}+\frac{{\partial }^2\Psi }{\partial z^2}+\frac{\partial \Psi }{\partial r}\left(\frac{\partial ln{\rho }_a}{\partial r}-\frac{\partial ln{B}_f}{\partial r}\right)+\frac{\partial \Psi }{\partial z}\left(\frac{\partial ln{\rho }_a}{\partial z}-\frac{\partial ln{B}_f}{\partial z}\right) \\ =\frac{r{B}_f}{{\rho }_a}\left(\frac{\partial {\overline{rV}}_{\theta }}{\partial z}\frac{\partial f}{\partial r}-\frac{\partial {\overline{rV}}_{\theta }}{\partial r}\frac{\partial f}{\partial z}\right) \end{gathered}$$

(32)

where the RHS of Equation (32) remains zero since the flow is irrotational in the upstream and downstream regions of the blade region.

2.4.4. Boundary Conditions

The stream function equation (Equation (32)) is an elliptic partial differential equation that can be solved subject to boundary conditions applied at all the boundaries.

There is no flow passing through the hub and the shroud surface, which means that the flow streamlines align with the hub and the shroud walls. Therefore, the following boundary conditions can be applied at the hub and the shroud:

$$\Psi \left(r,z\right)=const$$

(33)

The upstream (inlet) boundary condition is defined by the known velocity at the upstream boundary:

$$-\frac{{\rho }_a\left(r,z\right)}{r}\frac{\partial \Psi \left(r,z\right)}{\partial s\left(r,z\right)}={\overline{\overrightarrow{V}}}_{-\infty }·\overrightarrow{n}$$

(34)

The velocity at the downstream boundary will be uniform and parallel in the absence of the shed vorticity. This condition is defined as:

$$\frac{\partial \Psi \left(r,z\right)}{\partial \overrightarrow{n}}=0$$

(35)

2.5. Periodic Flow Equations

To solve the periodic flow equation (Equation (24)), we would first like to obtain the equation for the potential function for the periodic flow $\tilde{\Phi }$. Take the divergence of Equation (20) and substitute Equation (24), and we can get:

$$\begin{gathered} {𝛻}^2\tilde{\Phi }\left(r,\theta ,z\right)=S\left(\alpha \right){𝛻}^2{\overline{rV}}_{\theta } \left(r,z\right)+S^{\prime }\left(\alpha \right)𝛻\alpha ·𝛻{\overline{rV}}_{\theta } \left(r,z\right)-\overrightarrow{W}\left(r,\theta ,z\right) \\ ·𝛻ln\rho \left(r,\theta ,z\right)+\overline{\overrightarrow{W}\left(r,\theta ,z\right)·𝛻ln\rho \left(r,\theta ,z\right)} \end{gathered}$$

(36)

where the first two items on the RHS remain zero in the inlet and outlet regions, the last two items on the RHS require the calculation of the full 3D velocity field and density field. Due to the fact that the flow is periodic along the circumferential direction, it is possible to express the potential function for the periodic flow $\tilde{\Phi }$ in the form of a complex Fourier series.

$$\tilde{\Phi }\left(r,\theta ,z\right)=\sum _{n=-\infty ,n\ne 0}^{\infty }{\tilde{\Phi }}_n\left(r,z\right)e^{inB\theta }$$

(37)

where ${\tilde{\Phi }}_n$ is the Fourier coefficients of the potential function for the periodic flow $\tilde{\Phi }$. The periodic potential function $\tilde{\Phi }$ can be further approximated in the form of a truncated IDFT:

$$\begin{gathered} \tilde{\Phi }\left(r,\theta ,z\right)\approx \sum _{n=-\frac{N}{2}}^{\frac{N}{2}-1}{\tilde{\Phi }}_n\left(r,z\right)e^{inB\theta } \\ \tilde{\Phi }\left(r,{\theta }_j,z\right)\approx \sum _{n=-N/2}^{N/2-1}{\tilde{\Phi }}_n\left(r,z\right)e^{i{k}_n{\theta }_j} \end{gathered}$$

(38)

where ${\theta }_j$ is defined by Equation (39) and its value is the tangential coordinate of grid points in the tangential direction. The total number of grids is N and the spacing between grids is ∆θ.

$${\theta }_j=j∆\theta \left(j=\mathrm{0,1},2,\cdots ,N\right)$$

(39)

$k_n$ has been chosen to make the function $e^{i{k}_n{\theta }_j}$ space equally and periodically.

$$k_n=nB=\frac{2\pi n}{N∆\theta } \left(n=-\frac{N}{2},\cdots ,-\mathrm{1,0},1,\cdots ,\frac{N}{2}-1\right)$$

(40)

Using Equations (36), (38)–(40), the periodic flow potential function can be expressed in the following form:

$$\begin{gathered} \frac{{\partial }^2{\tilde{\Phi }}_n}{\partial r^2}+\frac{\partial {\tilde{\Phi }}_n}{r\partial r}+\frac{{\partial }^2{\tilde{\Phi }}_n}{\partial z^2}-\frac{n^2{B}^2}{r^2}{\tilde{\Phi }}_n=\frac{e^{-inBf}}{inB}{𝛻}^2{\overline{rV}}_{\theta }-e^{-inBf}𝛻f·𝛻{\overline{rV}}_{\theta }+F_n \\ \left(n=-\frac{N}{2},\cdots ,-\mathrm{1,0},1,\cdots ,\frac{N}{2}-1\right) \end{gathered}$$

(41)

where $F_n$ is the Fourier coefficient of $-\overrightarrow{W}·𝛻ln\rho +\overline{\overrightarrow{W}·𝛻ln\rho }$. The terms when n = 0 need to be neglected for the solution of the periodic flow since they represent the Monge–Clebsch formulation of the mean flow equation. Therefore, Equation (41) becomes:

$$\begin{gathered} \frac{{\partial }^2{\tilde{\Phi }}_n}{\partial r^2}+\frac{\partial {\tilde{\Phi }}_n}{r\partial r}+\frac{{\partial }^2{\tilde{\Phi }}_n}{\partial z^2}-\frac{n^2{B}^2}{r^2}{\tilde{\Phi }}_n=\frac{e^{-inBf}}{inB}{𝛻}^2{\overline{rV}}_{\theta }-e^{-inBf}𝛻f·𝛻{\overline{rV}}_{\theta }-R_n \\ \left(n=-\frac{N}{2},\cdots ,-\mathrm{1,1},\cdots ,\frac{N}{2}-1\right) \end{gathered}$$

(42)

where $R_n$ is defined as:

$$R_n\left(r,z\right)=\frac{1}{N} \sum _{j=1}^N\left[\overrightarrow{W}\left(r,{\theta }_j,z\right)·𝛻ln\rho \left(r,{\theta }_j,z\right)\right] e^{-\frac{i2\pi nj}{N}}$$

(43)

Since the periodic flow potential function is a real function $\tilde{\Phi }$, its Fourier transform satisfies the following equation:

$${\tilde{\Phi }}_{-n}\left(r,z\right)={\tilde{\Phi }}_n^{*}\left(r,z\right) \left(n=\mathrm{1,2},\cdots ,\frac{N}{2}\right)$$

(44)

where * denotes the complex conjugation and only half of the frequency spectra in Equation (42) need to be solved.

In order to solve Equation (42) for ${\tilde{\Phi }}_n$, boundary conditions need to be applied to the four boundaries of the physical domain, including the hub wall, the shroud wall, and the upstream (inlet) and downstream (outlet) boundaries.

The periodic velocity $\tilde{\overrightarrow{V}}$ should align with the hub and the shroud walls, which can be expressed by:

$$\tilde{\overrightarrow{V}}\left(r,\theta ,z\right)·\overrightarrow{n}=0$$

(45)

where $\overrightarrow{n}$ is the unit vector that is perpendicular to the hub or the shroud.

Substituting Equation (20) into Equation (45) and using Equations (15) and (37), the hub and the shroud boundary conditions for the periodic potential function can be obtained:

$$\frac{\partial {\tilde{\Phi }}_n\left(r,z\right)}{\partial \overrightarrow{n}}=\frac{e^{-inBf\left(r,z\right)}}{inB}\frac{\partial {\overline{rV}}_{\theta }\left(r,z\right)}{\partial \overrightarrow{n}}$$

(46)

The flow is uniform at the far upstream and far downstream boundaries, which implies that the periodic velocity vanishes as the flow approaches the upstream/downstream boundaries. This condition can then be expressed in the following form:

$$\frac{\partial {\tilde{\Phi }}_n\left(r,z\right)}{\partial \overrightarrow{n}}=0$$

(47)

The periodic velocity $\tilde{\overrightarrow{V}}$ is computed by using the IDFT (inverse discrete Fourier transform) of Equation (20).

$$\begin{gathered} \widetilde{\overrightarrow{V_r}}\left(r,\theta ,z\right)=\sum _{n=-\frac{N}{2},N\ne 0}^{\frac{N}{2}-1}\left[\frac{\partial {\tilde{\Phi }}_n\left(r,z\right)}{\partial r}-\frac{\partial {\overline{rV}}_{\theta }\left(r,z\right)}{\partial r}\frac{e^{-i{k}_{n}f\left(r,z\right)}}{i{k}_n}\right] \\ \widetilde{\overrightarrow{V_{\theta }}}\left(r,\theta ,z\right)=\sum _{n=-\frac{N}{2},N\ne 0}^{\frac{N}{2}-1}\frac{{i{k}_n\tilde{\Phi }}_n\left(r,z\right)}{r} \\ \widetilde{\overrightarrow{V_z}}\left(r,\theta ,z\right)=\sum _{n=-N/2,N\ne 0}^{N/2-1}\left[\frac{\partial {\tilde{\Phi }}_n\left(r,z\right)}{\partial z}-\frac{\partial {\overline{rV}}_{\theta }\left(r,z\right)}{\partial z}\frac{e^{-i{k}_{n}f\left(r,z\right)}}{i{k}_n}\right] \end{gathered}$$

(48)

2.6. Calculation of Thermophysical Properties

The thermophysical properties of any real gas can be obtained by calling the NIST REFPROP Fortran subroutine [35] or CoolProp C++ library [36] with any two known thermodynamic properties. The inlet total enthalpy $h_{01}$ and inlet entropy $s_1$ can be obtained from the given inlet total pressure $P_{01}$ and inlet total temperature $T_{01}$. The entropy remains constant throughout the whole flow domain since the flow is inviscid and isentropic.

The rothalpy, which is defined by Equation (49), is conserved along any streamline in a steady adiabatic flow.

$$I=h+\frac{1}{2}\left(W^2-U^2\right)=h_0-\varpi r{V}_{\theta }=const$$

(49)

where $h$ is the static enthalpy, $W$ is the relative velocity, $U$ is the blade speed, and $\varpi$ is the rotational speed. The inlet rothalpy can be calculated based on the inlet total enthalpy $h_{01}$ and the inlet ${\overline{rV}}_{\theta }$ value. The total enthalpy value $h_0$ and static enthalpy $h$ at other locations rather than the inlet can be calculated from the constant rothalpy value I and known ${\overline{rV}}_{\theta }$ or relative velocity $W$ at that location. Using the known entropy $s$ and enthalpy $h_0$($h$) with REFPROP or CoolProp, other fluid properties, including density $\rho$, pressure $P$ and temperature $T$, at any location can be obtained. It should be noted that all the properties discussed in this section are circumferentially averaged values, which means they are functions of $r$ and $z$. The static enthalpy jump across the blade can be calculated from Equation (50).

$$h^{+}\left(r,z\right)-h^{-}\left(r,z\right)=\frac{2\pi }{B}{\overrightarrow{W}}_{bl}\left(r,z\right)·𝛻{\overline{rV}}_{\theta }\left(r,z\right)$$

(50)

where $h^{+}$ is the static enthalpy on the blade pressure surface, $h^{-}$ is the static enthalpy on the blade suction surface and ${\overrightarrow{W}}_{bl}$ is the mean blade surface velocity. The mean static enthalpy is known, and it is equal to the mean of $h^{+}$and $h^{-}$. Using Equations (50) and (51), $h^{+}$and $h^{-}$ can be obtained. Similarly, using the constant entropy, the blade surface pressures $P^{+}$ (static pressure on the pressure surface) and $P^{-}$ (static pressure on the suction surface) will be obtained.

$$h\left(r,z\right)=\frac{1}{2}\left[h^{+}\left(r,z\right)+h^{-}\left(r,z\right)\right]$$

(51)

The 3D enthalpy field $h\left(r,\theta ,z\right)$can be calculated using Equation (49) using rothalpy I and the 3D relative velocity field $\overrightarrow{W}\left(r,\theta ,z\right).$ The 3D density field $\rho \left(r,\theta ,z\right)$ can then be obtained from the constant entropy and $h\left(r,\theta ,z\right)$using REFPROP or CoolProp. The Fourier transform of the natural logarithm of the density field is computed by the following equation once the density field $\rho \left(r,\theta ,z\right)$is available.

$${\rho }_n\left(r,z\right)=\frac{1}{N}\sum _{j=1}^N\left[\mathit{ln}{\rho }_n\left(r,{\theta }_j,z\right)\right]e^{\frac{i2\pi jn}{N}} { \theta }_j=j\Delta \theta$$

(52)

Taking the IDFT of $\frac{i{k}_n}{r}{\rho }_n\left(r,z\right)$ we can get the derivative of $\mathit{ln}{\rho }_n\left(r,{\theta }_j,z\right)$ in the tangential direction. Similarly, to get the axial and radial derivatives of $\mathit{ln}{\rho }_n\left(r,{\theta }_j,z\right)$ we need to take the IDFT of the corresponding derivatives of ${\rho }_n\left(r,z\right).$

$\overrightarrow{W}\left(r,\theta ,z\right)·𝛻ln\rho \left(r,\theta ,z\right)$ can be computed based on the derivatives of $\mathit{ln}{\rho }_n\left(r,{\theta }_j,z\right)$ and the relative velocity field $\overrightarrow{W}\left(r,\theta ,z\right)$. The Fourier transform of $\overrightarrow{W}\left(r,\theta ,z\right)·𝛻ln\rho \left(r,\theta ,z\right)$ is calculated using Equation (43) to obtain $R_n$. The 0^th component of Equation (43) which is equivalent to $\overline{\overrightarrow{W}\left(r,\theta ,z\right)·𝛻ln\rho \left(r,\theta ,z\right)}$ will be used in Equation (28) to calculate the artificial density ${\rho }_m$.

2.7. Calculation of Blade Shape

The blade shape is calculated based on the boundary condition that the relative velocity is parallel to the blade surface.

$${\overrightarrow{W}}_{bl}\left(r,\theta ,z\right)·\nabla \mathsf{\alpha }\left(r,\theta ,z\right)=0$$

(53)

Substituting Equations (1) and (23) into Equation (53) gives:

$$\left({\overline{V}}_r+{\tilde{V}}_{r,bl}\right)\frac{\partial f}{\partial r}+\left({\overline{V}}_z+{\tilde{V}}_{z,bl}\right)\frac{\partial f}{\partial z}=\frac{{\overline{rV}}_{\theta } }{r^2}+\frac{{\tilde{V}}_{\theta ,bl}}{r}-\varpi$$

(54)

To solve Equation (54), the initial value of $f$ along a quasi-orthogonal needs to be specified, and this is called the stacking condition. To make the blade radial-fibred, the gradient of the camber $f$ along the radial direction is forced to be zero using Equation (55).

$$\frac{\partial f}{\partial r}=0$$

(55)

3. Numerical Techniques

3.1. Grid Generation

To numerically solve all the PDEs above, a computational mesh is needed. In the case of irregular-shaped domains (e.g., curved hub and shroud profiles for a turbine rotor), transformation from the physical domain (r–z plane) to a regular (mostly rectangle) computational domain (ξ–η plane) is required. It should be noted that the numerical grid in the θ direction is not needed since the flow equations in the θ direction are solved using the discrete Fourier Transform. Figure 1 below shows an example of this coordinate transformation. It should be noted that the physical domain boundaries/edges (ABCD) do not have to be straight and can be curved or have irregular shapes. The newly generated ξ–η coordinate is called a body-fitted curvilinear coordinate, and its main advantage is that it makes the discretization of PDEs and imposing the boundary conditions more straightforward and simpler compared to solving PDEs on the r–z plane.

Using Equation (56), all the mesh point coordinates can be obtained through transfinite interpolation [37]. The curves ${\overrightarrow{c}}_1\left(\xi \right)$ and ${\overrightarrow{c}}_3\left(\xi \right)$are one pair of opposite sides of the physical domain, and ${\overrightarrow{c}}_2\left(\eta \right)$ and ${\overrightarrow{c}}_4\left(\eta \right)$are another pair shown in Figure 1. ${\overrightarrow{P}}_{\mathrm{1,2}}$ is the point A where curves ${\overrightarrow{c}}_1$ and ${\overrightarrow{c}}_2$ intersect.

$$\begin{array}{c}\overrightarrow{S}\left(\xi ,\eta \right)=\left(1-\eta \right){\overrightarrow{c}}_1\left(\xi \right)+\eta {\overrightarrow{c}}_3\left(\xi \right)+\left(1-\xi \right){\overrightarrow{c}}_2\left(\eta \right)+\xi {\overrightarrow{c}}_4\left(\eta \right)\\ -\left[\left(1-\xi \right)\left(1-\eta \right){\overrightarrow{P}}_{\mathrm{1,2}}+\xi \eta {\overrightarrow{P}}_{\mathrm{3,4}}+\xi \left(1-\eta \right){\overrightarrow{P}}_{\mathrm{1,4}}+\eta \left(1-\xi \right){\overrightarrow{P}}_{\mathrm{3,2}}\right]\end{array}$$

(56)

A typical meridional geometry of a radial-inflow turbine rotor is shown in Figure 2 and it is divided into three regions, including the inlet region, the blade region, and the outlet region by blade LE and TE. The body-fitted curvilinear mesh for each region can then be generated using the method explained above. The number of streamwise grid points (along ξ direction) in each region is proportional to its averaged meridional length, which will result in a relatively uniform grid size along the streamwise direction.

3.2. Transformed Equations in Curvilinear Coordinates

Once the numerical grid is ready, the equations described in Section 2 need to be transformed from (r,z) coordinates to (ξ,η) coordinates using the method described by Thompson et al. [38], which is shown in Appendix A.

Equation (32) is transformed into:

$$\begin{gathered} \alpha {\mathsf{\Psi }}_{\xi \xi }-2\beta {\mathsf{\Psi }}_{\xi \eta }+\gamma {\mathsf{\Psi }}_{\eta \eta }+\left\{\tau +\frac{z_{\eta }J}{r}+\alpha \left[{\left(\ln {\rho }_a\right)}_{\xi }-{\left(\ln B_f\right)}_{\xi }\right]-\beta \left[{\left(\ln {\rho }_a\right)}_{\eta }-{\left(\ln B_f\right)}_{\eta }\right]\right\}{\mathsf{\Psi }}_{\xi } \\ +\left\{\sigma -\frac{z_{\xi }J}{r}+\gamma \left[{\left(\ln {\rho }_a\right)}_{\eta }-{\left(\ln B_f\right)}_{\eta }\right]-\beta \left[{\left(\ln {\rho }_a\right)}_{\xi }-{\left(\ln B_f\right)}_{\xi }\right]\right\}{\mathsf{\Psi }}_{\eta } \\ =Jr{B}_f\left[{\left({\overline{rV}}_{\theta }\right)}_{\xi }{f}_{\eta }-{\left({\overline{rV}}_{\theta }\right)}_{\eta }{f}_{\xi }\right]/{\rho }_a \end{gathered}$$

(57)

Mean velocities $\overline{V_{\xi }}$ and $\overline{V_{\eta }}$ can be calculated using the following equations:

$$\begin{gathered} \overline{V_{\xi }}\left(\xi ,\eta \right)=\frac{J{\rho }_a\left(\xi ,\eta \right)}{r{B}_f\left(\xi ,\eta \right)}\frac{\partial \Psi \left(r,z\right)}{\partial \eta } \\ \overline{V_{\eta }}\left(\xi ,\eta \right)=-\frac{J{\rho }_a\left(\xi ,\eta \right)}{r{B}_f\left(\xi ,\eta \right)}\frac{\partial \Psi \left(r,z\right)}{\partial \xi } \end{gathered}$$

(58)

Equation (42) is transformed into:

$$\begin{gathered} \alpha {\left({\Phi }_n\right)}_{\xi \xi }-2\beta {\left({\Phi }_n\right)}_{\xi \eta }+\gamma {\left({\Phi }_n\right)}_{\eta \eta }+\left(\tau -\frac{z_{\eta }J}{r}\right){\left({\Phi }_n\right)}_{\xi }+\left(\sigma +\frac{z_{\xi }J}{r}\right){\left({\Phi }_n\right)}_{\eta }-\frac{J^2{n}^2{B}^2}{r^2}{\Phi }_n \\ =J^2\frac{e^{-inBf\left(\xi ,\eta \right)}}{inB}\left[\alpha {\left({\overline{rV}}_{\theta }\right)}_{\xi \xi }+{\gamma \left({\overline{rV}}_{\theta }\right)}_{\eta \eta }-2\beta {\left({\overline{rV}}_{\theta }\right)}_{\xi \eta }+\tau {\left({\overline{rV}}_{\theta }\right)}_{\xi }+\sigma {\left({\overline{rV}}_{\theta }\right)}_{\eta }\right] \\ -{J^{2}e}^{-inBf\left(\xi ,\eta \right)}\left\{f_{\xi }\left[\alpha {\left({\overline{rV}}_{\theta }\right)}_{\xi }-\beta {\left({\overline{rV}}_{\theta }\right)}_{\eta }\right]+f_{\eta }\left[\gamma {\left({\overline{rV}}_{\theta }\right)}_{\eta }-\beta {\left({\overline{rV}}_{\theta }\right)}_{\xi }\right]\right\} \\ -J^2{R}_n\left(\xi ,\eta \right) \end{gathered}$$

(59)

where $R_n\left(\xi ,\eta \right)$is calculated by the following equation.

$$R_n\left(\xi ,\eta \right)=\frac{1}{N} \sum _{j=1}^N\left[\overrightarrow{W}\left(\xi ,{\theta }_j,\eta \right)·𝛻ln\rho \left(\xi ,{\theta }_j,\eta \right)\right] e^{-\frac{i2\pi nj}{N}}$$

(60)

Blade boundary conditions, Equations (54) and (55), are transformed into:

$$\begin{gathered} \left({\overline{V}}_{\xi }+{\tilde{V}}_{\xi ,bl}\right)\frac{\partial f}{\partial \xi }+\left({\overline{V}}_{\eta }+{\tilde{V}}_{\eta ,bl}\right)\frac{\partial f}{\partial \eta }=J^2\left(\frac{{\overline{rV}}_{\theta } }{r^2}+\frac{{\tilde{V}}_{\theta ,bl}}{r}-\varpi \right) \\ \left(z_{\xi }\frac{\partial f}{\partial \eta }-z_{\eta }\frac{\partial f}{\partial \xi }\right)=0 \end{gathered}$$

(61)

where the periodic blade surface velocities ${\tilde{V}}_{\xi ,bl}$, ${\tilde{V}}_{\eta ,bl}$ and ${\tilde{V}}_{\theta ,bl}$can be calculated by Equation (62) and are shown below.

$$\begin{gathered} {\tilde{V}}_{\xi ,bl}\left(\xi ,\eta \right)={\sum }_{n=-\frac{N}{2},N\ne 0}^{\frac{N}{2}-1}\left\{\left[\alpha {\left({\Phi }_n^c\right)}_{\xi }-\beta {\left({\Phi }_n^c\right)}_{\eta }\right]\cos nBf-\left[\alpha {\left({\Phi }_n^s\right)}_{\xi }-\beta {\left({\Phi }_n^s\right)}_{\eta }\right]\sin nBf\right\} \\ {\tilde{V}}_{\eta ,bl}\left(\xi ,\eta \right)={\sum }_{n=-\frac{N}{2},N\ne 0}^{\frac{N}{2}-1}\left\{\left[\gamma {\left({\Phi }_n^c\right)}_{\eta }-\beta {\left({\Phi }_n^c\right)}_{\xi }\right]\cos nBf-\left[\gamma {\left({\Phi }_n^s\right)}_{\eta }-\beta {\left({\Phi }_n^s\right)}_{\xi }\right]\sin nBf\right\} \\ {\tilde{V}}_{\theta ,bl}\left(\xi ,\eta \right)={\sum }_{n=-\frac{N}{2},N\ne 0}^{\frac{N}{2}-1}-\frac{nB}{r}\left({\Phi }_n^c\sin nBf+{\Phi }_n^s\cos nBf\right) \end{gathered}$$

(62)

It should be noted that only the real part of the function is kept since the velocity is a real function. Therefore, the superscript ‘c’ indicates the real part of the function, and the superscript ‘s’ indicates the imaginary part of the function.

3.3. Numerical Differentiation

The elliptic Equations (57) and (59) are discretized by using central difference at the internal points and second-order upwind/downwind at the boundary points (inlet, outlet, hub, and shroud). Assuming a uniform grid in the transformed computational domain ($\Delta \xi =\Delta \eta$), for any point (j,k), where j is the point index along the ξ direction and k is the point index along the η direction, the derivatives of a function F(ξ,η) can be expressed in the following forms:

$$\begin{gathered} {\left(F_{\xi }\right)}_{j,k}=\frac{\partial F\left(\xi ,\eta \right)}{\partial \xi }=\frac{1}{2}\left(F_{j+1,k}-F_{j-1,k}\right) \\ {\left(F_{\eta }\right)}_{j,k}=\frac{\partial F\left(\xi ,\eta \right)}{\partial \eta }=\frac{1}{2}\left(F_{j,k+1}-F_{j,k-1}\right) \\ {\left(F_{\xi \xi }\right)}_{j,k}=\frac{{\partial }^{2}F\left(\xi ,\eta \right)}{\partial {\xi }^2}=F_{j+1,k}-2F_{j,k}+F_{j-1,k} \\ {\left(F_{\eta \eta }\right)}_{j,k}=\frac{{\partial }^{2}F\left(\xi ,\eta \right)}{\partial {\eta }^2}=F_{j,k+1}-2F_{j,k}+F_{j,k-1} \\ {\left(F_{\xi \eta }\right)}_{j,k}=\frac{{\partial }^{2}F\left(\xi ,\eta \right)}{\partial \xi \eta }=\frac{1}{4}\left(F_{j+1,k+1}-F_{j+1,k-1}+F_{j-1,k-1}-F_{j-1,k+1}\right) \end{gathered}$$

(63)

The above discretization equations are second-order accurate. For the boundary points (inlet, outlet, hub, and shroud), second-order upwind/downwind difference schemes are used similarly.

The first-order hyperbolic Equation (61) is discretized by using the central difference at the fictitious point (j+1/2,k) with second-order accuracy and the derivatives of a function F(ξ,η) can be expressed in the following forms:

$$\begin{gathered} {\left(F_{\xi }\right)}_{j+\frac{1}{2},k}=\frac{\partial F\left(\xi ,\eta \right)}{\partial \xi }=F_{j+1,k}-F_{j-1,k} \\ {\left(F_{\eta }\right)}_{j+1/2,k}=\frac{\partial F\left(\xi ,\eta \right)}{\partial \eta }=\frac{1}{4}\left(F_{j,k+1}-F_{j,k-1}+F_{j+1,k+1}-F_{j+1,k-1}\right) \end{gathered}$$

(64)

Equation (28) that solves the special artificial density ${\rho }_m$ is very similar to Equation (61). It can be discretized and solved using the same method.

All the discretized equations will be solved iteratively using the conventional Gauss–Seidel method with some initial conditions and the specified boundary conditions. The root mean square values of the difference between the current step solution and the previous step solution will be calculated and monitored.

3.4. Computation Scheme

The iterative computational scheme consists of the following steps:

• Step 1: Specify the design inputs, including the meridional geometry, the inlet volume flow rate (or the inlet velocity), the inlet total pressure, the inlet total temperature, the fluid name, the stacking condition, the RPM, the blade normal thickness distribution, and the blade number;
• Step 2: Generate the curvilinear-body-fitted mesh on the (ξ,η) plane (see Section 3.1);
• Step 3: The inlet entropy, the inlet enthalpy, and the inlet total density can be obtained from REFPROP or CoolProp based on the given inlet total pressure and inlet total temperature;
• Step 4: The density and velocity fields are initialized with the inlet values. The initial values for the blade camber f are set to zero;
• Step 5: Calculate the blade tangential thickness based on the blade camber and the normal thickness using Equation (27);
• Step 6: Calculate the blockage B_f using Equation (26);
• Step 7: Calculate the 3D velocity field using Equation (18);
• Step 8: Calculate the 3D enthalpy field using Equation (49);
• Step 9: Calculate the 3D density field based on the 3D enthalpy field and the inlet entropy using REFPROP or CoolProp;
• Step 10: Calculate the special artificial density using Equation (28);
• Step 11: Calculate the mean flow stream function $\Psi$ using Equation (57);
• Step 12: Calculate the mean flow velocity using Equation (58);
• Step 13: Calculate the potential function of the periodic flow $\tilde{\Phi }$using Equation (59);
• Step 14: Calculate the periodic flow velocity field using Equation (48) and the blade surface periodic velocity using Equation (62);
• Step 15: Calculate the blade camber using Equation (61);
• Step 16: If the maximum difference between the blade camber $f^n$ in the current iteration and the blade camber $f^{n-1}$from the previous iteration is less than the specified tolerance (1.0 × 10⁻⁵ radian), the solution is considered converged. Otherwise, go to step 5 and repeat the process.

3.5. Convergence

To improve the convergence of the blade updating process, it is found useful to use some underrelaxation [39]. This is achieved according to the expression:

$$f^n=f^{n-1}+RF\times \left(f^n-f^{n-1}\right)$$

(65)

where RF is the relaxation factor with a typical value between 0.1 and 0.8, $f^n$ is the new blade camber, the solution of Equation (61), and $f^{n-1}$ is the blade camber obtained in the previous iteration.

A similar approach is used to update the blade surface periodic velocities ${\tilde{V}}_{\xi ,bl}$, ${\tilde{V}}_{\eta ,bl}$ and ${\tilde{V}}_{\theta ,bl}$, where RF is gradually increased from 0.5 to 1.0 during the initial 10 iterations and remains constant after that.

$${\tilde{V}}_{bl}^{new}=RF\times {\tilde{V}}_{bl}^n$$

(66)

It is also found that the calculation of the blade camber is very easy to diverge quickly during the initial several iterations. This is caused by the simple zero camber values that are used as initialization. To obtain more reasonable initial blade camber values with a robust method, only the mean flow equation is solved (steps 13 and 14 are omitted) and the periodic flow equation is neglected, which means a zero periodic velocity field for the first 10–20 iterations. The blade camber is calculated only based on the mean flow velocities and will be used as a starting point for the next step of calculation involving all the equations. The resulting solution is the so-called actuator duct or asymmetrical through flow, which assumes there are an infinite number of blades. This assumption forces the flow to be uniform in the tangential direction (zero periodic velocities) [40].

4. Validation and Results

4.1. Design Specifications and Main Inputs

To validate the design method described above, it is applied to the design of an organic Rankine cycle radial-inflow turbine rotor. The turbine design specifications and the meridional geometry are taken from the work of da Silva et al. [21], where they performed a preliminary design, optimization, and CFD analysis of a R245fa radial turbine rotor. The rotor design specifications are summarized in

Table 1. Turbine rotor design specifications.

Parameters	Value	Unit
Working fluid	R245fa	-
Inlet total pressure P01	987.53	kPa
Inlet total temperature T01	95.89	°C
Inlet absolute flow angle a1	69.15°	-
Outlet static pressure P2	166.99	kPa
Mass flow rate $\dot{m}$	8.62	kg/s
RPM	9000	rev/min
Blade number	14	-

, and its main meridional geometry parameters are summarized in

Table 2. Turbine rotor main meridional parameters.

Parameters	Value	Unit
Tip radius r1	192.7	mm
Tip width b1	5.3	mm
Outlet hub radius r2h	57.8	mm
Outlet shroud radius r2s	97.3	mm
Outlet width b2	39.5	mm
Axial length L	65.7	mm

. The meridional geometry parameters are illustrated in Figure 3. The inlet and outlet domains are added to the rotor meridional geometry, and the $100\times 30$ body-fitted curvilinear mesh is generated using the method described in Section 3.1 is shown in Figure 4.

The blade thickness at the hub and the shroud sections are defined by cubic splines. The thickness between the hub and the shroud sections is obtained through linear interpolation, which is explained in Figure 5. The rotor blade LE ${\overline{rV}}_{\theta }$ (circulation) can be calculated based on the inlet velocity triangles. The TE ${\overline{rV}}_{\theta }$ is assumed to be zero to minimize the exit kinetic energy and maximize the turbine power. The ${\overline{rV}}_{\theta }$ contour on the meridional plane is obtained through linear extrapolation by assuming constant circulation along the spanwise direction, which is shown in Figure 6. It should be noted that ${\overline{rV}}_{\theta }$ remains constant inside the inlet and outlet regions due to the free-vortex assumption. The zero-stacking condition is applied at the LE. The blade wrap angle distribution can then be iteratively solved using the circulation method described above, and the resulting wrap angle and blade angle distribution are shown in Figure 7. The blade angle is calculated based on the wrap angle using Equation (67).

$$\beta ={\mathit{tan}}^{-1}\left(r\frac{\partial \theta }{\partial m}\right)$$

(67)

4.2. CFD Simulation

The structured mesh is generated using ANSYS TurboGrid for the rotor. The tip clearance of the rotor is assumed to be 0.5 mm, and at least 20 layers of grids are used to resolve the flow inside the tip gap. The near-wall first element size is set to 0.001 mm so that the resulting y+ is small enough (around 1.0) to resolve the boundary layers well when using the SST k–ω turbulence model [41]. Four different levels of grids are generated and their CFD results (mass flow rate and torque) are summarized in

Table 3. Mesh independence study.

Grid Number	No. of Elements	Mass Flow Rate (kg/s)	$\mathbf{Torque} (\mathit{N} \times \mathit{m})$
1	377,260	8.5605	260.44
2	715,010	8.5724	260.61
3	1,057,730	8.6186	262.53
4	2,079,220	8.6395	263.50

. By comparing the mesh independence study results in Table 3, it can be seen that the mass flow and the torque difference between Grid 3 and Grid 4 (the finest grid) are less than 0.4%. Therefore, Grid 3 is chosen for the CFD simulation and analysis. The mesh from Grid 3 is shown in Figure 8.

The computational domain consists of two parts, including one rotational part and one stationary part. The rotational part (blade region) is rotating with an RPM of 9000 rev/min, while the blade shroud wall is stationary with a contra-rotating wall velocity. A Frozen-Rotor interface is set between the blade passage and the outlet passage. The total pressure and the total temperature (937.53 kPa and 95.89 °C) are used as the inlet boundary conditions, and the static pressure (166.99 kPa) is applied at the outlet boundary. The inlet flow angle is specified as 69.15°. The periodic boundary condition is applied to all the periodic surfaces, and the SST k–ω turbulence model is used [41]. A R245fa RGP (real gas property) file is used to accurately simulate the thermodynamic properties of R245fa. ANSYS CFX is used to perform a steady-state RANS CFD analysis. The target residue for all the equations is set as 1.0 × 10⁻⁵. The inlet/outlet mass flow rate and the torque value measured on the blade surface and the blade hub surface are monitored during the CFD iterations. The solution is considered as converged when both mass flow rate and torque values converge to a constant value.

Once the CFD simulation is converged, the turbine rotor’s performance parameters are calculated based on CFD solutions and compared with the design from reference [21] in

Table 4. Comparison of CFD results.

Design	Mass Flow Rate (kg/s)	ηt-t	ηt-s	$\mathbf{Power} \left(\mathbf{k}\mathbf{W}\right)$
Use circulation method	8.6186	94.10%	87.28%	247.428
From reference [21]	8.62	94.89%	-	221.674

. As it can be seen, with the same inlet and outlet boundary conditions (pressure and temperature), almost the same mass flow rate (0.02% difference), and the same RPM (9000 rev/min), the rotor designed using the circulation method has good efficiency and 11.6% higher power.

Figure 9 shows the static pressure contour on the meridional plane. The static pressure is gradually decreased from the inlet (412 kPa) to the outlet (160 kPa). Figure 10 shows the Mach number contour on the rotor meridional plane. The Mach number inside the rotor is between 0.32 and 0.62, and the flow remains subsonic. The minimum Mach number occurs near the rotor outlet shroud region. Figure 11 shows the relative velocity contours in the blade-to-blade view at three different spanwise locations (10%, 50%, and 90%). In general, the flow is well-behaved and attached to the blade surface. Some positive LE incidence can be observed, and the low-momentum flow region near the shroud (90% span) is caused by the tip leakage vortex flow.

Figure 12 shows the blade loading distribution at three different (10%, 50%, and 90%) spanwise locations. The loading distribution of all three locations is similar and fore loaded with positive incidence. No obvious separation can be observed since the blade surface pressure is distributed smoothly. Figure 13 shows the comparison of the blade loading distributions at 50% spanwise location between the design using the circulation method and the design from reference [21]. It can be clearly seen that the blade loading of the design using the circulation method is much higher than that of the reference design, which explains the 11.6% higher power shown in Table 4. Figure 14 shows the inlet and outlet velocity triangles. It should be noted that the outlet velocity triangle is measured at the mid-span of TE, where R = 74.55 mm. As can be seen, there is still some tangential velocity component left in the exit of the turbine rotor, which means the blade’s turning is not enough to remove all the swirl/tangential velocities. This can be improved by reducing the TE ${\overline{rV}}_{\theta }$ in Figure 6 to a negative value and generating a new design.

4.3. FEA Simulation

The FEA simulation is performed using ANSYS Mechanical. A one-sector model of the turbine rotor is created by adding a 3 mm radius fillet at the blade root, and an unstructured mesh with around 510,000 elements is generated. Both the CAD model and the mesh are shown in Figure 15. The material used for the FEA analysis is AMS4133, which is a type of aluminum alloy, and its material properties are shown in

Table 5. Material properties for AMS4133.

Parameters	Value	Unit
Density	2.8	g/cm3
Modulus of elasticity	72.3	GPa
Poisson’s ratio	0.33	-
Tensile yield strength	410	MPa

. The rotor is rotating at a rotational speed of 9000 rev/min. A cylindrical support is applied to the backplate where the rotor is connected to the shaft. The resulting von-Mises stress contours on the blade surface and the back plate are shown in Figure 16. As can be seen, the maximum stress is around 20–30 MPa, which is well below the material yield strength. Therefore, it can be concluded that this turbine rotor is a structurally sound design.

5. Conclusions

This paper presents a turbomachinery blade design method that can be applied to the design of radial-inflow turbines for ORC applications. Using Monge–Clebsch decomposition, the velocity field is decomposed into a potential part and a rotational part. The vorticity field can be expressed by the cross product of the gradients of two Clebsch scalars, which are the circumferentially averaged circulation on the blade (${\overline{rV}}_{\theta }$) and the blade camber surface ($\alpha$). A blockage factor (B_f) is used in the mean flow equation to represent the blockage effect caused by the blade thickness, and the mean flow velocity is solved through the Stokes stream function. The periodic velocity is solved through a discrete Fourier transform in the frequency domain. The flow model is expanded to include compressible real gas from incompressible fluid and compressible ideal gas by using NIST REFPROP or CoolProp. This makes it suitable for modeling the organic fluid used in the ORC. The turbine blade geometry can be calculated iteratively by using the blade boundary condition. The governing equations are derived step by step, and the numerical schemes are discussed in detail. This design method is validated and used for the design of a R245fa turbine rotor, which is available in the literature. The resulting design shows similar t–t efficiency and 11.6% higher power values compared to the design obtained using the conventional method in the steady-state CFD simulation, and its structural integrity is confirmed by the FEA simulation.

In future work, this method will be expanded to properly model the blade thickness effect by using two vorticity sheets and be able to generate splitter turbine blades, which may be helpful to reduce the mass of the turbine wheel and increase its flow capacity.