Nội dung

8 Numerical flow calculations Real flows are described by partial differential equations which cannot be solved analytically in the general case. By dividing a complex flow domain into a multitude of small cells, these equations can be solved in an approximate manner by numerical methods. Because of their wide range of application, numerical flow calculations (“computational fluid dynamics” or “CFD” for short) have become a special discipline of fluid dynamics. The information given in this chapter is intended to help the understanding and interpretation of numerical flow calculations of centrifugal pumps. The focus is on viscous methods since these provide the best chance of describing realistically the boundary layers and secondary flows in decelerated flow on curved paths as encountered in radial and semi-axial pumps. This statement is not meant to preclude that simpler methods may be applied meaningfully to make a first design. The following discussion of CFD methods and possibilities focuses also on limits and uncertainties of CFD modeling as well as on quality criteria and issues in CFD applications, Chaps. 8.3.2, 8.3.3, 8.8 and 8.10. 8.1 Overview Because of the complex flow phenomena in centrifugal pumps, the design of impellers, diffusers, volutes and inlet casings is frequently based on empirical data for determining the flow deflection in the impeller and estimating performance and losses. The design of flow channels and blades then relies on experience and coefficients which originate from test data, Chaps. 3 and 7. The availability of relatively inexpensive computers with high computing powers has fostered the development of numerical methods which are able to solve the 3-dimensional Navier-Stokes-equations in complex components with reasonable effort. Therefore, numerical methods are used also in the pump industry with the object of optimizing the hydraulic components, to increase the reliability of performance prediction and thus to reduce testing costs. Computational fluid dynamics are still in a phase of development. Depending on the task at hand and the available resources, different types of CFD programs can be applied (see also Chap. 8.10 in this context): 1. Design programs for impellers or diffusers where the blades and the meridional section are modified based on simplified flow calculations until the designer 440 2. 3. 4. 5. 6. 7. 8. 9. 8 Numerical flow calculations considers the component as acceptable or optimum. Quasi-3D-calculations with computation times of a few minutes are often used for this purpose. Separate calculations of impellers, diffusers, volutes and inlet casings using 3D-Navier-Stokes-methods (“3D-NS”). If the results of such an analysis are unsatisfactory, the components are modified until an optimum design is reached. The success depends on the experience of the designer in using the programs and in assessing the results. However, the effort necessary to modify the geometry and the considerable amount of computation times limit the number of optimization cycles which are feasible in practice. Stage calculations which model inlet casing, impeller and collector as an entity, e.g. [8.4]. Stage computations can be conducted as unsteady calculations by the “sliding mesh” method (long computer times, high memory requirements), or as steady calculation by the “frozen rotor” method or by introducing “mixing planes” between stationary and rotating component, [8.28, 8.66]. Calculation of hydraulic forces. Inverse methods which determine that geometry which is able to generate a prescribed pressure distribution. Such methods may be used in special cases, for example for designing blades with favorable cavitation characteristics in two-dimensional flows, [8.2]. Other procedures are developed for optimizing blade shapes which produce less secondary flows, [8.67]. Automatic optimization of impellers where the program varies the geometry until the losses are minimized and/or a prescribed pressure distribution or any other target values or criteria are reached. Expert systems for design and analysis of the components and for interpretation of the results. They may be based on the methods according to items 1) and 2) (possibly item 3) and comprise programs for design and geometry generation in 3D. Furthermore, programs can be integrated into the system for generating the data for numerically controlled machining (e.g. for producing patterns and components) and for stress analysis. Calculation of two-phase flows. Modeling cavitating flows which capture the impact of the cavities on the flow field in the impeller. The methods mentioned are available in many variations. In addition methods for non-viscous flows are in use, such as 3D-Euler calculations or singularity methods for 2-dimensional flows – for example for developing blades of axial pumps with favorable cavitation characteristics. 3D-Euler programs can be supplemented by boundary layer models in order to include friction losses and flow separation. However, the losses caused by turbulent exchange of momentum (dominating in turbulent flows) cannot be captured that way. These methods will not be discussed in detail because they are being replaced more and more by 3D-Navier-Stokes programs which receive most of the development efforts. Three mathematical theories for the numerical calculation of incompressible flow shall be mentioned: 1. Irrotational non-viscous (potential) flows are handled by the Laplace equation. 8.2 Quasi-3D-procedures and 3D-Euler-calculations 441 2. Rotational non-viscous flows are handled by the Euler equations. 3. Rotational viscous flows are handled by the Navier-Stokes equations: laminar flows by means of Prandtl’s equation; turbulent flows by means of the Reynolds equations and turbulence models; modeling of large vortices (LES: “large-eddy-simulation”); direct numerical simulation of the turbulence (DNS: “direct-numerical-simulation”). There is a wealth of literature on numerical flow calculations, only a fraction of which can be quoted. An overview of newer developments and basic concepts, theories and procedures can be found in [8.9 to 8.14, 8.40 to 8.43, 8.58 to 8.60, 8.76 to 8.77]. Details on turbulence modeling are given in [8.8], [8.26, 8.33, 8.61, 8.79 and 8.80]. The performance of some commercial 3D-Navier-Stokes programs was compared using the test case of a turbulent flow through a 180°-bend, [8.27]. It was found that the patterns of secondary flow in the bend strongly depend on the turbulence model – even qualitatively. The velocity profiles in the bend were not at all well predicted. The evaluation of three 3D-Navier-Stokes programs and the comparison between CFD calculation and velocity measurements in a pump by laser anemometry was reported in [8.22 and 8.23]. The attainable accuracy of 3DNavier-Stokes calculations and their sensitivity to the influences of the calculation grid and numerical parameters was investigated in [8.35]. Technically speaking, 3D-Navier-Stokes calculations are basically appropriate for all applications because they offer the chance – if used correctly – of getting the most accurate solution. This does not imply that their use would be in every case economically justifiable or technically sensible. Many a case attacked by CFD could presumably be handled with a fraction of the effort by an engineering assessment and engineering common sense. 8.2 Quasi-3D-procedures and 3D-Euler-calculations 8.2.1 Quasi-3D- procedures The quasi-3D-method (“Q-3D”) was developed by Wu [8.5] at a time when highperformance computers were not yet available. In the Q-3D method the 2-D Euler-equations are solved in cylindrical coordinates in that the flow between the front and rear shrouds (“S2-surfaces”) and between the blades (“S1-surfaces”) of an impeller is iteratively superimposed. In doing so, rotational symmetry is assumed for the flow on the S2-surfaces in the meridional section. At the impeller blade trailing edge equal static pressures are imposed on the suction and pressure surfaces. The continuity equation is solved via a stream function. Frequently only one average stream surface is used in the meridional section in order to accelerate convergence. In this way efficient design programs may be developed, for example the “real-time-procedure” described in [8.6]. In this program system the first draft of an impeller can be optimized by interactive modification of the blades, 442 8 Numerical flow calculations while the resulting pressure and velocity distributions are instantly (“in real time”) available. The flow calculations can be done either by the Q-3D-method or by 3DEuler calculations. Quasi-3D-procedures are unable to predict losses and secondary flows. Within limits they may be used for: • Impeller calculations at the best efficiency point (BEP) if the influence of secondary flows is deemed small. Q-3D is therefore less appropriate for radial impellers of high specific speeds and/or large outlet widths. • Impeller design systems within the limits stipulated above. • Determination of cavitation inception based on the calculated pressure distributions on the blades – provided the approach flow conditions can be specified correctly and the grid resolution is made sufficiently high. • Q-3D-methods (in combination with empirical 1D-design procedures) were used extensively for the first draft of an impeller and its pre-optimization, [8.7]. The procedure needs extensive calibration by comparisons between calculation and test data in order to adapt the calculation parameters and to ease interpretation of the calculations. Essential geometric parameters must be defined, based on tested impellers with good performance, prior to being able to make efficient use of such programs. • Q-3D-methods can be supplemented by a boundary layer model in order to estimate friction losses. Quasi-3D-methods are not suitable for the calculation of diffusers, volutes, inlet casings and curved channels. They are of little use for impeller calculations at partload. As discussed extensively in Chap. 5.2, secondary flows in curved and/or rotating channels are influenced by non-uniform velocity distributions generated by boundary layers. Therefore velocity distributions calculated for non-viscous flow always deviate (more or less) markedly from test data. Consequently, nonviscous calculations of an inlet casing or a return channel yield incorrect incidences which falsify the pressure distributions. With axial approach flow or when losses are deemed insignificant, non-viscous calculations can predict the theoretical head (i.e. the flow deflection in the impeller) reasonably well at the best efficiency point. However, the velocity distribution over the impeller outlet width usually exhibits a trend which is contradicted by measurements; examples may be found in [8.45]. It may therefore appear difficult to draw relevant conclusions for improving the design from such velocity distributions (which does not preclude achieving good designs, as is indeed also possible without resorting to CFD). 8.2.2 Three-dimensional Euler-procedures The Euler equations represent the non-viscous terms of the Navier-Stokes equations discussed in Chap. 8.3. Consequently, 3D-Euler methods do not predict losses, but the balance of centrifugal, Coriolis and pressure forces is captured correctly with the exception of the impact of shear stresses and boundary layers. 3D- 8.3 Basics of Navier-Stokes calculations 443 Euler methods fail in applications where boundary layers play an essential role, such as local flow separation and secondary flows which are dominated by boundary layer flow. The theoretical head is predicted correctly only if the velocity distribution at the impeller outlet is not affected by boundary layer effects and secondary flow to the extent that the integral of u×cu over the exit is falsified. 3D-Euler procedures are better suited for accelerated flows (turbines) and for the same applications as listed above for Quasi-3D. Accuracy and relevance are significantly superior to with Quasi-3D. In fully separated flow, the recirculation in the impeller is determined to a large extent by centrifugal and Coriolis forces. Fully developed recirculation (but not the onset of recirculation) can therefore be modeled rather well by 3D-Euler calculations. The application of Euler methods is not recommended for diffusers and volutes. With respect to velocity distributions the same limits apply as for Quasi-3D. 8.3 Basics of Navier-Stokes calculations 8.3.1 The Navier-Stokes equations Because of the widespread use of 3D-Navier-Stokes programs, this topic is discussed here in detail in order to show the possibilities and limits of numerical calculations and to give recommendations for assessing CFD results. Consider a 3-dimensional, incompressible flow with the relative velocities wx, wy, wz in a Cartesian x,y,z-coordinate system; the rotation is around the z-axis. Since the present discussion focuses on the basic aspects, conservation of momentum is written only for the x-direction (complete equations and derivations can be found in many text books, e.g. [1.11]): ∂w x ∂w x ∂w x ∂w x 1 ∂p + − ω2 x + 2ωw y = + wx + wy + wz ∂t ∂x ∂y ∂z ρ ∂x 2 2 § ∂2w x + ∂ wx + ∂ wx ν¨ 2 2 ¨ ∂x ∂y ∂z 2 © · § ∂σ' ∂τ' ¸ + ¨ x + xy + ∂τ'xz ¸ ¨© ∂x ∂y ∂z ¹ ·1 ¸ ¸ρ ¹ (8.1) On the left side of the equation are the substantial acceleration terms (Chap. 1.4.1), the effect of the pressure and the body forces in the rotating system (i.e. centrifugal and Coriolis acceleration, gravity being negligible in comparison). On the right-hand side are the loss terms; the first loss term describes the effects of the molecular viscosity while the second term covers the losses due to turbulent exchange of momentum. Equation (8.1) represents a very general form of the conservation of momentum which comprises some technically significant special cases: 1. Equation (8.1) can be employed for investigating the flow in a rotating impeller in the relative system or, with ω = 0, for calculations of stationary components (diffusers, volutes, inlet casings). 444 8 Numerical flow calculations 2. By setting the right-hand side to zero, the Euler equation for non-viscous flow is obtained, which is solved in Euler-programs. 3. If the second term on the right-hand side is set to zero, the Navier-Stokes equations for laminar flows are obtained. The remaining term covers the shear stresses due to the molecular viscosity Ȟ. As long as constant temperature is assumed, Ȟ is a property of the fluid independent of the flow characteristics. Together with the continuity equation (8.2) the three component equations, exemplified by Eq. (8.1), yield a system of four partial differential equations for the four unknown functions p, wx, wy and wz. ∂w x ∂w y ∂w z + + =0 ∂x ∂y ∂z (8.2) If the velocities in Eq. (8.1) are considered as unsteady, the Navier-Stokes equations would be sufficient to calculate turbulent flows. This direct numerical simulation (“DNS”) of turbulent flows in pump components is, however, beyond present computing capabilities. In order to exactly describe all turbulent fluctuations, extremely fine calculation grids and excessive computing times are required. The number N of elements needed in a grid for DNS can be estimated from N ≈ Re9/4; in the order of 1013 elements would be required for Re = 106 [8.77]. Reynolds therefore replaced the unsteady velocities by w + w'(t) with w representing the time-averaged velocity and w'(t) its turbulent fluctuation. Equation (8.1) then gives the “Reynolds-averaged” Navier-Stokes-equations (referred to as “RANS”). Present-day Navier-Stokes-Programs work on this basis. The stresses caused by turbulent exchange of momentum (the “Reynolds-stresses”), are given by Eq. (8.3) (the over-bars signifying time-averaged quantities): σ'x = −ρ w 'x 2 τ'xy = −ρ w 'x w 'y τ'xz = −ρ w 'x w 'z (8.3) Since the fluctuating velocities wx', wy' and wz' are unknown, the system of the four equations cannot be solved (“closure problem”). Therefore, by means of a “turbulence model”, additional empirical equations must be devised for the velocity fluctuations. This can be done for example by establishing a relation between the Reynolds stresses and the average velocity components. 8.3.2 Turbulence models In general terms the turbulence model describes the distribution of the Reynolds stresses in the flow domain. All turbulence models in use are of an empirical nature. They contain constants and concepts which were selected so that CFD calculations agreed as well as possible with the test results in the considered particular geometry and flow regime (for example a one-sided plane diffuser). If a turbulence model contains five (or even more) empirical constants, a given set of test data gathered on a specific flow can be represented by different combinations of 8.3 Basics of Navier-Stokes calculations 445 said constants which are more or less equivalent for the specific case. It is however difficult to decide which of the combinations would be physically the most relevant and the best suited for extrapolation to different flow situations. A broader validation would require investigating a statistically relevant number of measurements, which appears scarcely practicable. The above discussion leads to the conclusion that there is no universally valid turbulence model which will yield optimum results for all applications. Instead it is necessary to select the turbulence model most suitable for the components to be calculated and to carefully validate it by comparing the CFD results to test data, Chap. 8.8.2. The choice of the turbulence model and the pertinent turbulence parameters thus harbors one of the main uncertainties of 3D-Navier-Stokes calculations of turbomachines. Turbulence models have therefore received considerable attention in a great number of publications, e.g. [8.8 to 8.10, 8.26, 8.33, 8.61 to 8.63, 8.79]. From the multitude of turbulence models only a few are briefly discussed below. Each of the models presented is available in different versions. This situation unmistakably testifies to the fact that the problem of turbulence modeling is not yet solved satisfactorily. Most turbulence models are based on the concept of the eddy viscosity νt (introduced in Chap. 1.5.1) which is determined solely by flow features (in contrast to the molecular viscosity ν which is a property of the fluid). The eddy viscosity is frequently expressed by a velocity scale Vt and a length scale Lt of the turbulence, through νt = Vt×Lt. Turbulence models can be classified by the number of transport equations employed. These are differential equations which attempt to describe the transport of the turbulence quantities in the flow domain. Zero-equation models (or “algebraic” turbulence models) assume that the eddy viscosity depends essentially on local flow quantities (such as velocity gradients) and on a given length scale for the energy carrying vortices. This dependence is expressed by algebraic equations. Zero-equation models are applicable to attached flows as well as to fluid jets and wake flows; their advantage is the comparatively low calculation time required. Single-equation turbulence models take into account (in addition to the local flow quantities) the “history”, i.e. the phenomena upstream of the fluid element considered. This is done by means of a transport equation for the velocity scale of the turbulence, for example in form of a transport equation for the kinetic energy k of the turbulence. Two-equation turbulence models use two transport equations for the turbulence parameters Vt and Lt (or k and ε). The currently popular k-ε-model belongs to this group. Usually it is combined with a wall function. The model is based on the specific kinetic energy k and the dissipation rate ε of the turbulent fluctuations. The production of turbulence is calculated from the local velocity gradients. The turbulence levels Tu and length scale of the vortices Lt are employed; these parameters are related by Eqs. (8.4) and (8.4a): 446 k= 8 Numerical flow calculations 1 § '2 3 '2 '2 · '2 ¨ wx + w y + wz ¸ ≈ wx 2© ¹ 2 Tu = w x '2 w Re f = 1 w Re f 2 k 3 ε= k3 / 2 Lt (8.4) The eddy viscosity is calculated from: ν t = cμ L t k = cμ k2 ε (8.4a) The “standard-k-ε-model” uses the constant cμ = 0.09. In addition, there are four more empirical constants in the transport equations for the turbulent kinetic energy k and the dissipation rate ε, which (among others) control the production of the kinetic energy of the turbulence. In two-dimensional attached boundary layers, the production and dissipation rates of the turbulent kinetic energy are in local equilibrium. In this case the relation of Eq. (8.4b) holds between the wall shear stress τw and the turbulence parameters (see also Table 8.1): wτ ≡ τw = cμ 1 / 4 k ρ (8.4b) The standard k-ε-model shows weaknesses when modeling the following types of flow: • • • • • • • Flows on curved paths Decelerated flows 3-dimensional boundary layers Rotating components, since the body forces influence the boundary layers Swirling flows Strong secondary flows Secondary flows which are solely induced by turbulence cannot be captured. Such secondary flows are encountered in channels with non-circular cross sections, see for example Fig. 1.6. • The production of turbulent kinetic energy is over-predicted in locations with strong velocity gradients (for example near a stagnation point). Therefore flow separations are not (or insufficiently) recognized. Almost all of the phenomena listed above are encountered in impellers, diffusers, volutes and inlet casings. As a consequence of the above deficiencies, the calculation of losses becomes unreliable and zones with flow separation are either predicted too small or they are not recognized at all, [8.32]. Furthermore, the velocity distributions in decelerated flow or curved streamlines are not calculated correctly, as was demonstrated with a diffuser in [8.48] and with bends in [8.27 and 8.49]. The k-ε-model failed equally when calculating the swirling flow in a pipe (as induced for example by a vortex rope): the decay of the swirl due to friction was over-predicted and the velocity distribution was not captured correctly, [8.49]. 8.3 Basics of Navier-Stokes calculations 447 It should be noted that the standard k-ε-model with the logarithmic wall function fails in apparently simple geometries such as diffusers or bends. It can therefore be concluded that the k-ε-model is in principle unsuitable for the calculation of pumps: after all, the flow through impellers, diffusers and volutes is decelerated, follows curved paths and is frequently subject to separation. The inability of the standard k-ε-model to capture these phenomena is not only due to the equations and parameters describing the turbulence. The situation is aggravated by the use of the logarithmic wall function which imposes rather than predicts the flow near the solid boundaries, Chap. 8.3.3.1 In spite of its severe shortcomings, the k-ε-model is widely used because convergence is better than with other turbulence models. Furthermore, as it is one of the first industrially applicable two-equation turbulence models, there is a broad validation basis. Some programs employ a modified k-ε-model which attempts to capture the influence of curved flow paths. Realizable k-ε-model: This modification of the standard k-ε-model prevents nonphysical solutions such as negative values of k and ε as well as negative normal turbulent stresses. Unrealistically high production of turbulent kinetic energy is thus avoided. k-ω-model: This model was developed specifically for flows against strong pressure gradients (such as encountered in a diffuser). It solves a transport equation for the frequency ω of “large” vortices (with Lt = k1/2/ω). The k-ω-model captures the flow near the walls more accurately than the k-ε-model, while the latter better describes the processes in the core flow. Therefore both models were combined providing a smooth transition from near-wall to core-flow, [8.61]. Close to flow separation there is a mismatch between turbulence production and dissipation. This situation is not well handled by the eddy viscosity concept employed in the k-ω-and k-ε-models. Shear stress transport model (SST): This model is designed to mitigate this shortcoming, [8.61]. It employs the k-ε-model in the core flow, the k-ω-model near solid surfaces and a modified relation for the eddy viscosity which limits the shear stresses due to pressure gradients. The k-ω-model uses five empirical constants. Additionally, the SST-model includes empirical functions for the transition between k-ω- and k-ε-model. Flows on curved paths are still not well simulated. Kato-Launder k-ε-model: This model is a modification of the standard k-εmodel. It avoids an over-production of turbulence in domains with high velocity gradients (e.g. near the stagnation point at the vane leading edge) by an alternative formulation of the production term, [8.81]. Since the model over-predicts turbulence production in rotating channels its application to pumps is not really attractive, [8.82]. Low-Reynolds k-ε-model: This model is able to better resolve the flow closer to solid surfaces. To this end, the constants of the standard k-ε-model are multiplied 1 In the future the standard k-ε-model will presumably lose significance for calculating pumps and compressors. 448 8 Numerical flow calculations by functions which depend on a Reynolds number defined with turbulence quantities. The maximum distance of the cells to the wall must be limited to y+ ≤ 2 (y+ is defined in Table 8.1). The first cells near a wall are to be located in the viscous sublayer. Since the sublayer in highly turbulent flow is very thin, a sufficient resolution of near-wall flow requires a high number of nodes. Two-layer model: This model equally attempts to more accurately describe nearwall flows including roughness effects. A comparison between calculations and test data has been reported in [8.62 and 8.63]. This procedure also relies heavily on empirical data. The limit between the near-wall layers and the core flow is often selected at νt/ν ≈ 20. Reynolds stress transport model: All components of the Reynolds stresses, Eq. (8.3), are calculated by transport equations without resorting to the eddy viscosity concept. The transport equations are derived by inserting w + w'(t) in the Navier-Stokes equation and subtracting the Reynolds-averaged equation. The resulting additional 7 coupled equations (6 Reynolds stress components and 1 length scale) require a considerable calculation effort and severely impair convergence. It is also necessary to define boundary conditions for the additional parameters. Therefore this model is not popular in the turbomachinery industry. Large eddy simulation (LES): Consider a local flow separation where vortices are generated in the shear layer between main flow and stalled fluid. At their origin the vortices are large. When traveling downstream of the stalled zone, the vortices decay through mixing. Large vortices have a high energy content and hence an impact on the flow distribution. When vortices break up into smaller sizes, energy is dissipated. In large eddy simulation the movements of the individual largescale vortices are followed by unsteady calculations, while the small eddies are handled by a statistical turbulence model, [8.36 and 8.37]. In LES calculations the wall treatment is difficult and wall functions are required to determine the unsteady wall shear stresses. LES calculations of a sharp 90°-bend with R/D = 1.0 in [10.52] and of a diffuser in [8.48] agreed well with test data in terms of flow separation, velocity distributions and friction coefficients. Also a diffuser opening to one side could well be calculated by LES, while parallel calculations with the k-εmodel yielded disappointing results because velocity profiles and friction losses deviated strongly from the measurements. However, because of contradictory investigations such results should not be considered as generic. LES calculations require very fine grids which should be as isotropic as possible. The number of elements needed increases with the square of the Reynolds number. Therefore LES applications in pumps are scarcely feasible at present because of limited computation power. Thus experience with pump calculations is lacking. Flows with Reynolds numbers of roughly 5'000, possibly up to 20'000, can be handled (depending on computer and geometry). In [8.48] a “v2-f” turbulence model (comparable with the SST model) was presented which (when used with a linear wall function) yielded results on the flow in a diffuser very similar to the LES calculations.