Numerical Solution of Flow through Turbomachinery with the SIMPLEN Algorithm

 

 

G. D. Thiart^I; T. W. von Backström^II

^ISenior Lecturer, Department of Mechanical Engineering University of Stellenbosch 7600 Stellenbosch
^IIProfessor, Department of Mechanical Engineering University of Stellenbosch 7600 Stellenbosch

 

 

---

SUMMARY

The SIMPLEN algorithm for the numerical solution of convection-diffusion problems on nonstaggered grids is applied to two turbomachinery through flow problems. Incompressible flow only is considered, and the k-ε model is used to simulate turbulence effects. The first problem is that of flow through a propeller. Both a uniform and a shear approach flow are considered. The thrust and torque exerted by the propeller on the fluid are modelled as prescribed body forces. Good agreement between the numerical results and the corresponding results available in the literature is obtained. The second problem is that of the flow through an axialflow fan. Both an axisymmetric and a distorted inflow is considered. Blade element theory is used to model the thrust and torque exerted by the fan blades on the air. Reasonable correspondence with experimental results is obtained.

---

 

 

Introduction

In recent years, methods based on the SIMPLE (Semi-Implicit Method for Pressure-Linked equations) algorithm of Patankar and Spalding [1] have been used successfully to solve a wide range of flow problems. In the original SIMPLE algorithm separate grids, staggered in space, are used for the calculation of scalar properties such as pressure, and each of the velocity components. Thus four grids are used for a three-dimensional flow problem, which complicates the implementation of boundary conditions. If all variables are located at the same grid points, on the other hand, a checkerboard pressure field may result. This has been described by, amongst others, Patankar [2], Rhie and Chow [3] and Peric [4].

Several methods that have been proposed to overcome the checkerboard problem on nonstaggered grids, are described by Patankar [5] in a recent review paper. Thiart [6] subsequently presented a new finite difference scheme that also eliminates the need for staggered grids. The scheme was incorporated in an algorithm called SIMPLEN: the "N" attached to SIMPLE being an indication that a nonstaggered grid is used.

The finite difference scheme on which SIMPLEN was originally based, has subsequently been refined by Thiart [7], It was shown that the new scheme is capable of producing numerical solutions of comparable or superior accuracy to that of other numerical methods. The differencing scheme was originally developed specifically for a Cartesian coordinate system; Thiart and Von Backström [8] recently extended the scheme to cylindrical polar coordinates, which is particularly useful for turbomachinery applications.

It was also shown by Thiart [6] that the SIMPLEN methodology is able to simulate accurately the one-dimensional flow through a thrust-producing actuator disc. In this paper, the actuator disc simulation is extended to two typical turbomachinery throughflow problems: a propeller operating with uniform inflow and with a sheared inflow, and an axial flow fan operating with axisymmetric inflow and with circumferentially distorted inflow.

 

Governing Equations

The governing equations for the flow problems that are considered in this paper, are the well-known Navier-Stokes equations for incompressible flow (see for example Schlichting [9]), augmented by the k-ε turbulence model. These equations are those for the conservation of mass, i.e.

and for the conservation of axial momentum, radial momentum, moment of momentum, turbulent kinetic energy (k) and the decay rate of turbulent kinetic energy (ε), which can be written in the general form

where φ denotes a general variable representing either the axial velocity µ, the radial velocity v, the angular velocity ω, k or ε, and where m is equal to one for φ = ωand zero otherwise. Eg. (2) can be regarded as a transport equation for the general variable φ, with transient and convective terms on the left hand side being balanced by a source term and diffusion terms on the right hand side. The source terms are as follows:

The dynamic viscosity μin Eqs. (2)-(7) is replaced by the eddy viscosity for turbulent flow problems:

The k and ε-equations are standard except for the last term on the right hand side of Eq. (7), which is the modification of Launder et al. [10] to account for the destabilizing effect of swirl on turbulence. The empirical constants that appear in Eqs. (6)-(8) are assigned the following values:

T and Q represent respectively the external thrust and torque on the fluid, while G_k is the turbulence kinetic energy production term, given by (Gupta and Lilley [11]):

 

Solution Procedure

The computational domain is subdivided into a number of non-overlapping control volumes, one control volume per grid point. The governing equations are discretized by integration over these control volumes; these discretised equations are solved using the SIMPLEN algorithm. The numerical procedure was incorporated in a computer code called FLO WAX, written in single precision FORTRAN 77. Boundary conditions for k and εat solid walls are determined by means of wall functions in the manner described by Burns and Wilkes [12]. Details about the discretization, the code itself and the verification of the code are given by Thiart [13]. All computations reported in this paper were performed on a CONVEX 120 computer.

 

Flow through a Propeller

The three-dimensional turbulent flow through a propeller was simulated by Pelletier [14] (see also Pelletier and Schetz [15]), using a finite element method with a relatively simple turbulence model (compared to the k-ε model): the eddy viscosity is assumed to vary only in the axial direction, the variation being described by a first order ordinary differential equation. The propeller and flow characteristics used for the simulation correspond to those of the wind-tunnel experiments of Kotb [16]. The propeller was modeled as an actuator disk of constant thickness Δ and diameter D. The radial distributions of thrust and swirl force per unit volume, i.e. (dQ/dV)/r, per unit volume was assumed to be trapezoidal, increasing linearly from zero at r/D = 0,125 to a maximum at r/D = 0,35 and decreasing linearly again from r/D = 0,425 to zero at r/D = 0,5.

Two cases were considered by Pelletier and Kotb: one for uniform flow into the propeller, and one for a shear flow. These two cases were solved with FLOWAX using the same computational domain and, as far as possible, the same boundary conditions as for the calculations of Pelletier.

(a) Uniform Flow

The propeller and flow characteristics are as follows: Propeller diameter: D = 0,492 m

Propeller thickness: Δ = 0,0407 D = 0,020 m

Air density: p = 1,177 kg/m³

Relative air velocity: U = 8,52 m/s

Total thrust: T = 3,80 Ν

Total torque: Q = 0,348 Nm

The computational grid is depicted in Fig. 1. In the axial direction, the computational domain extends two diameters upstream of the propeller face and three diameters downstream of the rear face of the propeller. The rear face is located at ζ = 0, thus locating the inlet plane at z = -2,0407D and the outlet plane at z = 3D. Radially, the computational domain extends from the propeller axis tor = 1,2D. The computational grid has 34 stations in the axial direction and 22 stations in the radial direction, giving a total of 748 grid points. The grid has been designed so that expansion ratios (that is, the ratio of successive δz's or δr's) never exceed 1,3. Other constraints were that one radial grid line had to coincide with each of the following axial stations within the computational domain: z = - 0,0407D (propeller face), z = 0 (rear face for propeller), and z = 0,025D (location where velocity components and static pressure were measured by Kotb). Similarly, in the radial direction, one grid line had to coincide with each of the radial stations r = 0,425D (location for which axial distributions of velocity components and static pressure are given by Pelletier) and r = 0,5D (outer radius of propeller).

 

 

Boundary conditions for the numerical solution (other than those of the zero-gradient type) were as follows:

z = -2,0407D: u = U, ν = 0, ω = 0, k = 0,0035U²,

ε = 2000pC_µk²;

z = 3D: p = p₀ = 0; r = 1,2D: u = U

The free-stream value of k is one commonly used with the k-ε model, while the free-stream value of εwas chosen so that the free-stream value of eddy viscosity corresponds to that used by Pelletier. Some results are presented in Fig. 2.

The axial distribution of the axial velocity component at r = 0,425D (Fig. 2a) appears to have a lower peak value compared to that of Pelletier, and the decay of axial velocity component also seems to be much slower. This is not the case, however; at r = 0,425D the thrust and torque distributions and therefore the velocity components change rapidly with radius so that small differences in radial profiles appear as large differences in axial distributions.

The axial velocity profile at z = 0,025/) is compared with those of Pelletier and Kotb in Fig. 2b (values shown at a negative radius in Figs. 2b, 2d and 2f correspond to the measurements of Kotb on a vertical plane, and values at a positive radius to measurements on a horizontal plane). The numerical procedure does of course yield identical results for the vertical and horizontal planes. The predicted axial velocity profile is almost identical for the two numerical solutions, and shows good correspondence with the measured velocity profile.

The corresponding axial distributions of the azimuthal velocity component are shown in Fig. 2b. Here also a lower peak value and a slower decay with axial distance is observed. In this case, however, there is a real difference, as borne out by the azimuthal velocity profiles, Fig. 2d. The current calculations seems to yield better correspondence with the experimental results. Note also in Fig. 2c the absence of "wiggles" upstream of the propeller compared to the substantial oscillations produced by the finite-element calculations.

Finally, axial distributions and profiles of the pressure coefficient C_p = (p-p₀)/0,5pU² are shown in Fig. 2e-f. Evident in Fig. 2e is the pressure jump associated with thrust provided by the propeller. The axial distributions are almost identical, with the current calculations showing a slightly higher peak value on the upstream face and a slightly lower peak value at the downstream face of the propeller. Correspondence with the experimental results is again reasonable at z = 0,025D.

(b) Shear Flow

The propeller dimensions and air density are the same as for the uniform flow case, the other characteristics being as follows:

Relative air velocity: U = U(0) = 8,52 - O,67cos0 m/s

Total thrust: T = 4,00 Ν

Total torque: Q = 0,332 Nm

The computational grid is the same in the axial and radial coordinate directions, but now has twelve equi-spaced azimuthal planes instead of only one. The boundary conditions used here are similar to those for the uniform flow case:

z = -2,0407D: u = U(θ), ν = 0, ω= 0,

k = 0,0035U², ε = 2000pC_µk²;

z = 3D: p = p₀ = 0; r = 1,2D: u = U(θ)

Note that the inlet boundary condition on k (and therefore also on ε) is based on the average inlet velocity. Some results are presented in Figs. 3-4.

Axial distributions of the axial velocity component at r/D = 0,425 are shown in Fig. 3a at four azimuthal sections through the computational domain. Pelletier's results differ from these in the same way as for the uniform flow case and are therefore not reproduced here. The nonsymmetrical behaviour of the 0° and 180° axial velocity curves reported by Pelletier is clearly evident in Fig. 3a.

Corresponding axial distributions of the aximuthal velocity component are presented in Fig. 3b. Once again no "wiggles" can be observed upstream of the propeller, and the near-perfect axial symmetry reported by Pelletier is also evident.

Radial profiles of axial and azimuthal velocity components are compared with those of Pelletier and Kotb in Fig. 4. Values shown at negative radii in these figures correspond to results at θ= 180° and θ = 270°, and values at positive radii to results at θ = 0° and θ= 90°. The same comments as for the uniform flow case apply generally.

 

Flow through an Axial Flow Fan

Thiart [13] recently carried out wind tunnel experiments with an axial flow fan operating under distorted inflow conditions; the corresponding numerical simulations will now be described. The idealized computational domain corresponding to the experimental configuration is shown in Fig. 5. In the experiment the test fan was mounted in a round duct, the intake of which was located in a hole in the floor of the subsonic windtunnel of the Departement of Mechanical Engineering at the University of Stellenbosch. The windtunnel was used to simulate cross-wind conditions. Measurements of cross-flow velocity, volumetric flow rate, fan static pressure, fan speed and shaft power were taken for two cases: one with no cross-flow and one with crossflow of approximately 10 m/s, which is of the same magnitude as the average axial velocity through the duct. Detailed measurements of velocity and static pressure were also taken with a five-hole probe at a location approximately one and a half duct diameters downstream from the inlet.

 

 

Several idealizations are incorporated in the computational domain. The wind tunnel test section is modelled as two parallel plates; the spacing between the plates is equal to the height of the test section. No account is taken of the sidewalls of the test section. The duct section extends from the test section to the five-hole probe measuring position. The fan rotor hub and fan motor is modelled as a round cylinder, the length and diameter of which correspond to the combined length of the fan hub and motor and the hub diameter respectively (the motor is actually slightly smaller in diameter than the hub). The small gap between the hub and the motor is not modelled; neither are the cooling fins on the motor surface. The fan blades are not modelled directly as rotating solid surfaces, but indirectly: the influence of the fan blades are modelled as body forces by means of the blade element theory commonly employed for aircraft and ship propeller calculations. Full details are given by Thiart [13].

The computational grid for axisymmetric inflow is shown in Fig. 6. There are 67 axial stations and 45 radial stations, for a total of 2301 grid points. Expansions ratios of 1,08 in the axial direction and 1,15 in the radial direction are used throughout. The computational grid for distorted inflow is the same as for axisymmetric inflow, but now has sixteen equispaced azimuthal planes instead of only one. The total number of grid points in this case is 36,036.

The boundary conditions for the fan performance predictions are as follows: at the solid surfaces, the no-slip conditions (ι/ = ν = ω = 0 except at the rotor hub, were ω= Ω) and wall functions for k and ε are used; at the duct outlet plane zero-gradient boundary conditions are used for all variables except static pressure, for which the radial equilibrium boundary condition dp/dr = prω² is used; at the inlet and outlet planes between the parallel surfaces zero-gradient boundary conditions are used for k and ε, while the velocity components are fixed to provide the correct flow rate through the fan.

The velocity profiles to be used as boundary conditions between the parallel surfaces (representing the wind tunnel top and bottom walls) posed a problem: these boundaries are not sufficiently far removed from the fan for uniform flow velocity profiles to be used. To overcome this problem, the inlet velocity profiles (µ and ν only - ωis equal to zero) were determined from a potential flow solution on part of the computational domain, excluding the duct section and assuming uniform inflow into the duct section. For the distorted inflow case, the required cross-flow velocity is simply added to the potential flow solution.

The potential flow solution was also used to provide initial conditions for the fan performance predictions; static pressure conditions are readily obtained from Bernoulli's equation. In the duct section the initial values for ν and ωare taken as zero, while the axial velocity component µ is taken as uniform everywhere between the inlet and outlet planes of the duct.

(a) Axisymmetric Inflow

The main experimental data for the case without cross-flow are as follows:

Air density: p = 1,172 kg/m³

Volume flow rate through fan: q = 2,627 m³/s

Increase in static pressure through fan: Δp = 79 N/m²

Fan speed: Ν= 1409 r.p.m.

Shaft power: Ρ = 981 W

Several solutions were obtained for the axisymmetric inflow case in order to establish reasonable values for some of the constants in the blade element model. For the chosen set of constants, the predicted values of fan static pressure and fan power differed by -5% and ~0% respectively compared to the experimental results. A similar calculation was done for uniform inflow boundary conditions - in this case the predictions differed from the experimental results by -8% and +4% respectively. The introduction of the inflow boundary conditions calculated from the potential flow solution is therefore considered worthwhile.

The general features of the flow field are depicted in Fig. 7, which shows the streamlines in an azimuthal plane. The reverse flow at the outer part of the fan blades, which was also observed during the experiment, is clearly evident in the figure. Another zone of recirculatory flow is predicted downstream of the fan motor, extending beyond the outflow boundary. The reverse flow at the outflow boundary was not observed during the experiment; possible reasons for the discrepancy are inadequacies of the turbulence model and the influence of a right angle bend and turning vanes just downstream of the outflow plane which are not accounted for in the numerical predictions.

The effect of the reverse flow at the outer part of the fan blades is depicted in Fig. 8a-c. Fig. 8a shows the local angles of flow onto the fan blades: separated flow (beyond the blade stall angle, which was taken as 15°) prevails over the outer quarter of the blade span and attached flow over the inner three-quarters. The corresponding thrust and torque distributions, non-dimensionalized by means of the mean thrust and torque respectively, are shown in Fig. 8b-c. There is a sharp drop in thrust, and a corresponding but smaller drop in torque (and hence power requirement) at the outer quarter of the blade span.

Predicted profiles of velocity and static pressure at the outflow boundary are compared with experimental results in Fig. 8d-f (measurements were taken at eight equi-spaced azimuthal sections). As can be expected from the foregoing discussion, the predicted axial velocity profile (Fig. 8d) does not compare well with the measured one: the reverse flow region at the centreline has to be compensated for by a higher maximum value of axial velocity in the outer region. Also, the predicted radial velocities at the outflow boundary are close to zero, in accordance with the zero-gradient boundary condition on the axial velocity component, whereas the experimental results indicate a significant inward radial flow. Another calculation with a computational domain two duct diameters longer than the original one was also performed: the results were almost identical, indicating that the discrepancies cannot be ascribed to the zero-gradient outflow boundary conditions.

The results for the azimuthal velocity and the static pressure (Figs. 8e-f) are better; both these variables are fairly well predicted over the region excluding the reverse flow zone. The good prediction of w is to be expected to some extent, because the amount of swirl is essentially proportional to the power input by the fan, and the constants in the blade element model have been chosen so as to give the correct power. Also, the static pressure at the outflow boundary is related to w through the radial equilibrium boundary condition, so that if w is predicted correctly, the static pressure must also be predicted correctly. It is nevertheless encouraging that the shapes of the predicted profiles are also essentially correct.

(b) Circumferentially Distorted Inflow

The main experimental data for the case with 10 m/s crossflow are as follows:

Air density: p= 1,176 kg/m³

Volume flow rate through fan: q = 2,590 m³/s

Increase in static pressure through fan: Δp = 102 N/m²

Fan speed: Ν= 1369 r.p.m.

Shaft power: Ρ = 1165 W

In this case, both the fan static pressure and power were underpredicted by -11,5% with respect to the experimental values. Two recirculation zones are predicted at the blade tips and roots, as can be seen in Fig. 9, which shows the local angles of attack (the crossflow is from left to right; the fan rotation is clockwise). These angles of attack are defined locally as the angle between the relative velocity vector and the blade chord, so that an angle of attack larger than the blade chord angle (which is the angle between the plane of rotation and the blade chord) indicates that the axial velocity is negative and hence that reverse flow exists at that position. The blade tip chord angle is approximately 30°, therefore there must be reverse flow between θ = 135° and θ= 180° near the duct surface (0 is measured counterclockwise relative to the crossflow direction). Similarly, because the blade root chord angle is approximately 55°, there must be reverse flow between θ = 225° and θ = 315° at the hub surface.

 

 

It is evident from Fig. 9 that the flow onto the fan blades is highly distorted, both radially and azimuthally, and that separated flow prevails over a large fraction of the fan blade area. This means that the blade loading is very much non-uniform, and that the stresses in a blade can vary considerably during a revolution. The non-di-mensionalized distributions of thrust and torque over the fan blade area are depicted in Figs. 10-11: both distributions show maxima of more than twice their mean values.

 

 

 

 

The numerically predicted and experimentally determined axial velocity distributions at the duct outlet are compared in Figs. 12-13. The predicted value and position of the maximum axial velocity at the duct outlet compare very well with the experimental values: approximately 2.2 times the mean axial velocity at θ = 225° and r/D = 0,4. Comparisons for the azimuthal velocity component and static pressure at the duct outlet (not shown here due to space restrictions) are similar to those for the axisymmetric inflow case.

 

 

 

 

Concluding Remarks

The application of the SIMPLEN algorithm to problems involving distorted flow through turbomachinery has been presented in this paper. It is evident from the two case studies that it is possible to simulate practical turbo-machinery flow problems with the SIMPLEN methodology incorporating the modified k-ε turbulence model.

The general performance of the numerical procedure is good, considering the complexity of some of the problems, but there is one major deficiency: that of excessive computing time requirements for especially the three-dimensional simulations. Investigations are currently under way to speed up the computer code by means of multigrid schemes, in order to make the methodology more useful as an engineering tool. Further work that is envisaged include the generalization of the numerical procedure to a general curvilinear coordinate system, and the calculation of the flow through multistage turbo-machinery.

Nomenclature

C₁,C₂,C_µ - constants in the k-ε model

C_p - pressure coefficient

D - diameter

G_k - turbulence kinetic energy production term

k - turbulence kinetic energy

m - exponent

Ν - rotational speeds

p,pο,Δp - static pressure, static pressure at outlet, static pressure difference Ρ - power

q - volume flow rate Q - torque

r,δr - radial direction, radial distance between grid points

S - source term

T - thrust

u - axial velocity component

U - free stream velocity for propeller flow

ν - radial velocity component

V - volume

w - azimuthal velocity component

z,δz - axial direction, axial distance between grid points

Δ - actuator disc thickness

ε - turbulence kinetic energy decay rate

μ,μ₁ - dynamic viscosity, turbulent dynamic viscosity

φ - general dependent variable

p - density

σ_k,σ_E - Schmidt numbers in the k-ε model

θ - azimuthal direction

w - angular velocity of fluid particles

 

References

1. Patankar, S. V., and Spalding, D. B. - "A Calculation Procedure for Heat, Mass and Momentum Transfer in Three-Dimensional Parabolic Flows", Int. J. Heat Mass Transfer, Vol. 15, pp. 1787-1806, 1972.         [ Links ]

2. Patankar, S. V. - "Numerical Heat Transfer and Fluid Flow", Hemisphere, Washington D.C., 1980.

3. Rhie, C. M., and Chow, W. L. - "A Numerical Study of the Turbulent Flow Past an Isolated Airfoil with Trailing Edge Separation", AI A A J., vol. 21, pp. 1525-32, 1983.         [ Links ]

4. Peric, M. - "A Finite Volume Method for the Prediction of Three-Dimensional Fluid Flow in Complex Ducts", Ph.D. thesis, University of London, 1985.         [ Links ]

5. Patankar, S. V. - "Recent Developments in Computational Heat Transfer", ASMEJ. Heat Transfer, Vol. 110, pp. 1037-45, 1988.         [ Links ]

6. Thiart, G. D. - "Finite Difference Scheme for the Numerical Solution of Fluid Flow and Heat Transfer Problems on Nonstaggered Grids", Num. Heat Transfer, Vol. 17(B), pp. 43-63, 1990.         [ Links ]

7. Thiart, G. D. - "Improved Finite Difference Scheme for the Solution of Convection-Diffusion Problems with the SIMPLEN Algorithm", Num. Heat Transfer, Vol. 18(B), pp. 81-95, 1990.         [ Links ]

8. Thiart, G. D., and Von Backström, T. W. - "Extension of the SIMPLEN Algorithm Differencing Scheme to Cylindrical Polar Coordinates", to be published in Num. Heat Transfer, 1992.

9. Schlichting, H. - "Boundary Layer Theory", 6th ed., McGraw-Hill, New York, 1968.

10. Launder, B. E., Pridden, C. H., and Sharma, B. I. - "The Calculation of Turbulent Boundary Layers on Spinning and Curved Surfaces", ASMEJ. Fluids Engrng, vol. 99, pp. 231-9, 1977.         [ Links ]

11. Gupta, A. K., and Lilley, G. D. - "Flowfield Modeling and Diagnostics", Abacus, Tunbridge Wells, Kent, 1985.

12. Burns, A. D., and Wilkes, N. S. - "A Finite Difference Method for the Computation of Fluid Flows in Complex Three-Dimensional Geometries", AERE Report No. 12342, Harwell Laboratory, 1987.

13. Thiart, G. D. - "A Numerical Procedure for Predicting the Effects of Distorted Inflow Conditions on the Performance of Axial Flow Fans", Ph.D. thesis, University of Stellenbosch, 1990.         [ Links ]

14. Pelletier, D. H. - "Finite Element Solutions of the Navier-Stokes Equations for 3-D Turbulent Shear Flows", Ph.D. thesis, Virginia Polytechnic Institute and State University, 1984.         [ Links ]

15. Pelletier D. H., and Schetz, J. A. - "Finite Element Navier-Stokes Calculation of Three-Dimensional Turbulent Flow Near a Propeller", Aí A A J., Vol. 24, pp. 1409-16, 1986.

16. Kotb, M. A. - "Experimental Investigation of 3-D Turbulent Free Shear Flows Past Propellers and Windmills", Ph.D. thesis, Virginia Polytechnic Institute and State University, 1984.         [ Links ]

 

 

First received December 1991
Final Form December 1991