Acta Polytechnica doi:10.14311/AP.2015.55.0373 Acta Polytechnica 55(6):373–378, 2015 © Czech Technical University in Prague, 2015 available online at http://ojs.cvut.cz/ojs/index.php/ap ANALYTICAL METHODS FOR CALCULATING FAN AERODYNAMICS Jan Dostál, Jan Kužel∗ VZLÚ, a.s., Prague ∗ corresponding author: kuzel@vzlu.cz Abstract. This paper presents results obtained between 2010 and 2014 in the field of fan aerodynamics at the Department of Composite Technology at the VZLÚ aerospace research and experimental institute in Prague – Letňany. The need for rapid and accurate methods for the preliminary design of blade machinery led to the creation of a mathematical model based on the basic laws of turbomachine aerodynamics. The mathematical model, the derivation of which is briefly described below, has been encoded in a computer programme, which enables the theoretical characteristics of a fan of the designed geometry to be determined rapidly. The validity of the mathematical model is assessed continuously by measuring model fans in the measuring unit, which was developed and manufactured specifically for this purpose. The paper also presents a comparison between measured characteristics and characteristics determined by the mathematical model as the basis for a discussion on possible causes of measured deviations and calculation deviations. Keywords: aerodynamics; fan; impeller; blade; total pressure; static pressure; flow tube. 1. Introduction The Department of Composite Technology at VZLÚ, a.s., formerly known as the Department of Aircraft Propellers, has developed algorithms in the field of fan aerodynamics and has designed a measuring unit. This development work is presented here. The first software was designed for calculating the aerodynam- ics of propellers. Other mathematical models and algorithms were also based on earlier models devel- oped for propellers, and these have gradually been replaced by models designed for fans. During the four-year development period, 22 versions of fan aero- dynamics computing programmes were developed and encoded. Each version has two branches: the first branch calculates separate modes with a more de- tailed output of the aerodynamic parameters along the blade, while the second branch calculates the aero- dynamic characteristics with the output of the integral parameters of the entire fan. Both branches share a common subroutine with coded algorithms. Each of the versions contains several so-called “constant” pa- rameters, i.e. unknown, optional and non-dimensional parameters that may depend on the blade geometry, the velocity profile in front of the impeller, etc. The non-dimensionality of the parameters should guaran- tee independence from the size of the computed fan. The closing part of Section 3 indicates the procedure for the design of the 14th constant parameter for the 23rd version of the computing programme, which was developed in 2015. The latest version of FanA22_7 [1] contains 13 parameters. Most of them are real numbers, while two are natural numbers. In order to determine them, optimization was proposed. This optimization requires, if possible, the verified aerody- namic characteristics of a suitable fan, together with its geometry. These are not easy to obtain, because the manufacturers of fans keep such information se- cret. The value function whose minimum is sought is defined by the sum of the absolute values of the weighted differences in the y-coordinates of the points of known and calculated characteristics. A measuring unit with the V936mod model fan has been designed and manufactured for these and other activities at VZLÚ. The measurement unit enables the required characteristics to be obtained and evaluated. 2. From propellers to fans The term “fan” is very general, and is used for simple devices with rotating blades used for ventilation and air conditioning, as well as for more complex machines for air transport. This paper deals with a single- stage fan without a stator and with a cylindrical hub and a case. This cylindrical case contributes to the significant difference between a propeller and a fan. Figs. 1 and 2 show models of the flow tubes used in the software for fans (from version 19 to 22) and for propellers. As each of the tubes of the fan model has the same volumetric flow rate, their thickness decreases away from the hub, by contrast with the propeller model with its cylindrical tubes of constant thickness. The proven whirl model for the components of induced velocity behind the propeller (axial, circum- ferential and radial) assumes that the radii of the flow helices are constant (Fig. 2), while in the actual flow behind the propeller with induced velocities close to the model velocities, we see contraction taking place, i.e. the helix radii are decreasing. Neither of the streamlines of the actual fan behaves exactly accord- ing to the model shown in Fig. 1, according to which 373 http://dx.doi.org/10.14311/AP.2015.55.0373 http://ojs.cvut.cz/ojs/index.php/ap Jan Dostál, Jan Kužel Acta Polytechnica Figure 1. Model of fan streamlines. the flow behind the fan should be stabilized in the computation plane and should follow the cylindrical flow area. In fact, the settling plane (perpendicular to the rotation axis) does not comply with the computa- tion plane as regards the constant axial distance, but moves away from it in the direction of the rotation axis; the difference increases with the extent to which the current mode differs from the computation mode. For the fans under consideration here, no flow con- traction occurs, as it is prevented by the cylindrical case behind the impeller. The flow section after the propeller and the static pressure are smaller than for the flow after a fan impeller of the same geometry. However, the flow speed after the propellers is higher. It should be noted in this context that there is a transition member between these two objects, i.e. a propeller in a ring, whose axial length determines which of the two objects is closer. The mathematical model of the propeller software is based on vortex theory for calculating the induced components of the velocity vectors after the propeller. These components, and also increases in the axial and circumferential components of the input speed, calculated for a propeller with the same geometry, were included in the input data for the computing programmes of the fan aerodynamics until the 18th version (inclusive). The subsequent versions replace the propeller-induced velocity by the solution of the equation of the radial balance between pressure and centrifugal forces. Figure 2. Model of propeller streamlines. 3. Basic assumptions and equations of the mathematical model of fan aerodynamics The assumptions on the model are the following: (1.) The flow is axially symmetric, and the boundaries of the flow tubes are rotating surfaces. (2.) Constant air density % (incompressible flow, with the exception of the input data of the lift and profile drag coefficients calculated for Ma and Re). (3.) The input absolute velocity is in the direction of the rotation axis (not necessarily constant along the radius). (4.) The impeller hub and the interior of the case have constant radii (cylindrical surface) [2]. (5.) A one-stage fan without stator vanes. (6.) Computation plane 1 before the impeller and computation plane 2 behind the impeller are located at a sufficient distance from the wheel, and so the radial velocity components are zero. The model is based on the followin three equations: Bernoulli Equation (energy conservation). pc1 + ∆pcT = pc2 + ∆pz, i.e., ∆pc = pc2 − pc1 = ∆pcT −∆pz, (1) 374 vol. 55 no. 6/2015 Analytical Methods for Calculating Fan Aerodynamics where pc1 [Pa] is the total pressure of the air mass before the impeller, pc1 = p1 + % 2c 2 1; p1 [Pa] is the static pressure of the air mass before the impeller; % [kg/m3] is the density of the flowing air mass; c1 [m/s] is the input speed of the air mass flow, which is equal to the axial component c1a; ∆pcT [Pa] is the theoretical increase in the total air mass caused by the energy supplied (disregarding losses); ∆pc [Pa] is the actual increase in the total air mass caused by the energy supplied (considering losses); pc2 [Pa] is the total pressure of the air mass after the impeller, pc2 = p2 + % 2c 2 2a (the ISO standard applies, but as the total pressure is considered variable, p̃c2 = p2 + % 2c 2 2a [3], where c2a [m/s] is the axial component of the output speed); p2 [Pa] is the static pressure of the air mass after the impeller; c2 [m/s] is the output air mass flow rate, c2 ≈ √ c2 2a + c2 2u applies and the radial component is neglected; c2a [m/s] is the axial component of the output speed; c2u [m/s] is the pe- ripheral component of the output speed; ∆pz [Pa] is the total loss of the pressure in the impeller plus the kinetic energy of the peripheral components of the absolute output speed (this applies to each flow tube: ∆pz = ∆pzp + ∆pzs + % 2c 2 2u; ∆pz [Pa] is the pressure loss in the entire fan is the sum of the pressure loss of all the flow tubes; ∆pzp [Pa] is the pressure loss is a profile loss; it depends on Re [4]; ∆pzs [Pa] is the pressure loss is a secondary loss [5]. Euler’s Equation (theorem of angular momen- tum). ∆pcT = %(u2c2u − u1c1u), (2) where u2 [m/s] is the transport velocity of the impeller rotation at the radius r2 of the intermediate plane of the flow tube under consideration after the impeller; c2u [m/s] is the circumferential component of the output velocity of the flowing air mass; u1 [m/s] is the transport velocity of the impeller rotation at the radius r1 of the intermediate plane of the flow tube under consideration before the impeller; c1u [m/s] is the circumferential component of the input velocity of the flowing air mass, according to assumption 3 we have c1u = 0. Equation of radial balance between pressure and centrifugal forces after the impeller. dp2 dr2 = %c2 2u r2 . (3) These three equations are taken from publications [6–10]; their derivation can be found in almost any turbomachine aerodynamics textbook. Equation (3) is differential; it is solved numerically in a mathemat- ical model, and the differentials along the radii r2 are replaced by differences equal to the thickness of the flow tubes in output computation plane 2. The number of flow tubes is one of the optimized constant parameters, and it is entered into a data file within the range of 3 to 20. The radial equilibrium equa- tion encoded in FanA22_7 is the simplest. It applies exactly to the cylindrical flow areas considered at a sufficient distance behind the impeller. It therefore disregards the inclination of the flow areas relative to the axis of rotation or their curvature. Determining the axial position of computation plane 2, i.e. the distance behind the impeller for which the assump- tion of cylindrical streamlines applies, should have been one of the tasks implemented in the measuring unit. However, it has not yet been determined, due to the complexity of the structural modifications of the existing solution. Equation (3) applies only to cylindrical flow areas considered at a sufficient distance behind the impeller. In the measurement of fan model V936mod carried out in 2014, the measuring probe was positioned relatively close to the trailing edges of the impeller, and so assumption 6 was not complied with. Another version of the fan aerodynamics computing programme for 2015, FanA23_7, must therefore also take into account radial components C2r of the output velocities in the general rotary flow areas: c2 = √ c2 2a + c2 2u + c2 2r. (4) The elementary centrifugal force acting against the mass element in the radial direction will therefore be dFO = dm (c2 2u r2 + c2 2m rc2 cos γ ) , (5) where dm = %r2 d%dr2 dl [kg] is the mass element from flow tube thickness dr2; r2 dϕ [m] is the rotation path element in radius r2; dl [m] is the path element in the axial direction; c2m [m/s] is the meridian velocity component, c2m = √ c2 2a + c2 2r; rc2 [m] is the radius of curvature in the meridian cross-section of the flow area; γ [rad] is the angle between the tangent and the meridian cross-section of computation plane 2, with the rotation axis. The elementary pressure force acting against the surface of face dS = r2 dϕdl (perpendicular to the radial direction) can be written as dFp2 = dS dp2. (6) where dp2 = p2(r2 + dr2/2) − p2(r2 − dr2) [Pa] is the elementary difference in pressure between the up- per and lower boundary rotation (axially symmetric) surface of the elementary flow tube. According to the radial balance of the centrifugal and pressure force at the flow tube element, the fol- lowing applies: dFO = dFp2. (7) By adding (5) and (6) in (7), we get dm (c2 2u r2 + c2 2m rc2 cos γ ) = dS dp2. (8) 375 Jan Dostál, Jan Kužel Acta Polytechnica Figure 3. Measuring unit. By another adding for dm and dS, %r2 dϕdr2 dl (c2 2u r2 + c2 2m rc2 cos γ ) = r2 d%dl dp2. (9) After reducing r2, dϕ, dl, we get the equation for the radial balance, expanded by the meridian element: dp2 dr2 = % (c2 2u r2 + c2 2m rc2 cos γ ) . (10) However, deriving equation (10) is not the end of our work. The second element in the brackets on the right side contains two unknown quantities: angle γ of the meridian velocity c2m and the meridian cross- section radius of curvature of rc2. Angles γ were measured on the model fan V936mod, at each point of the measured characteristics, with a sufficient density along the radius. What now needs to be done is to place computational plane 2 in the middle axial position of the measuring probe head, to evaluate graphically all measured curves and replace them with the analytical dependence γ(r2, Q, ϕ0.75), (11) where Q [m3/s] is the volumetric flow through the fan (see Figs. 6 and 5); ϕ0.75 [rad] is the angle of the impeller blades, at 75% of the blade length (see Figs. 6 and 5), between the rotation plane and the profile section chord. Finding the dependence (11) will not be an easy job. Paradoxically, it would be easier to find a better equation for radii of curvature rc2, for which we would need to carry out measurements on three successive planes, which is quite beyond realistic possibilities. However, if we look carefully at Fig. 1, we see that the meridian radius of curvature of the axisymmetric flow area could be inversely proportional to the difference of radii r2 − r1, in which the flow area intersects computation planes 1 and 2: rc2 = Crc r2 tr r2 − r1 , (12) where Crc [1] is the proportionality constant which can be included among the constant parameters of the up- coming 23rd version of the computing programme for Figure 4. Detail of the measuring unit at the rotor. the direct task of the aerodynamics of fan FanA23_7 ; rtr [m] is the mean radius of the flow tube between radii r1 and r2. Equation (12) can easily be verified for the extreme limit approach r2 → r1, occurring at the lower and upper limits of the flow path. The denominator of the fraction on the right side of (12) is in this case close to zero, which corresponds to the infinite radius of curvature of the straight lines of both cylindrical boundaries. 4. Measuring unit Figures 3 and 4 show the functional arrangement of the measuring unit, which has been completely de- signed and implemented at VZLÚ, a.s. It was created because established manufacturers are not willing, for obvious reasons, to share their “know-how” by pro- viding the aerodynamic characteristics of a particular rotor geometry, which are necessary for proper tuning of the “constant” parameters of calculation algorithms. Another reason for building our own measuring unit was to be able to verify the design calculations on real blades. The measuring unit, consisting of an actuator cov- ered with the outer and inner parts of the shaped input, a rotor, the shaped output, the casing, the attachment and throttle diaphragm, is equipped on its outer casing with a set of flanges allowing probes to be introduced into the space before and after the rotor. These probes attached to a traversing device allow the total pressure to be measured along the entire blade height. This allows the parameters for each flow tube to be verified. The measuring unit is provided with a cylindrical extension satisfying the condition of length L ≥ 8D necessary to measure steady flow using a ridge probe. The output of the unit is fitted with a movable throttle diaphragm, en- abling the aerodynamic loading of the rotor to be continuously controlled by simulating different levels of output losses. The whole measuring unit was de- signed as modular, so that it can be complemented, 376 vol. 55 no. 6/2015 Analytical Methods for Calculating Fan Aerodynamics Figure 5. V936mod – efficiency ηc (volumetric flow rate Q, angle of setting of blades). Figure 6. V936mod – pressure ∆pc (volumetric flow rate Q, angle of blade setting). without extensive interventions, with an output stator and a rotor with a different hub ratio. The mea- sured impeller can be removed from the measuring unit to allow the blades to be adjusted to the desired angle [11]. It was discovered during the measurements in 2014 that, due to the absence of a stator, the air mass exiting the fan has the character of a rotating vortex, adhering to the inner walls of the flow duct due to the effects of the centrifugal forces. No steady-state condition is achieved at any point of the attachment that would enable the use of a ridge probe for mea- suring the integral characteristics of the fan. For this reason, the steady-state condition measurement was replaced by measuring the pressure profile along the height of the blade in front of and behind the rotor, carried out using five-hole pressure probes installed on the crabbing device. These partial measurements were then used for determining the mean parameters of the flow in front of and behind the rotor, using integration. 5. Aerodynamic characteristics of the V936mod fan The aerodynamic characteristics of the measured V936mod fan are shown in Figs. 6 and 5. In our case (a single impeller without the stator), the total pressure and efficiency values are equal to the static values, since the circumferential velocity components cannot be used, and the axial components in the channel of a constant cross-section are the same in incompressible flow: ∆pc = pc2 − pc1 = p2 + 1 2%c 2 2a − p1 − 1 2%c 2 1a = p2 − p1 = ∆p, (13) where c2a and c1a [m/s] are axial velocity components in the measuring planes after and before the impeller. 377 Jan Dostál, Jan Kužel Acta Polytechnica In Fig. 5, the efficiency ηc = ∆pc ∆pc + ∆pz (14) is plotted. 6. Discussion on the correspondence between the calculations and the measurements The correspondence between the characteristics com- puted by the FanA22_7 program and the measured characteristics is not yet completely satisfactory, and the variations increase with the difference between the absolute values of the current and the calculation angle of the blade setting, which is 35°. There are several possible explanations for these variations: (1.) The constant parameters mentioned above are not yet firmly established, since the results of the measurements will be used to establish them. These parameters will be established by means of an op- timization process, which should be completed in 2015 after all the results have been evaluated. The evaluation will require the dependence on the radius along the blade, which will be processed graphically; (2.) In the FanA22_7 program, the module for cal- culating the dependence of the axial speed com- ponent on the radius was blocked. This blockage was created for reconciliation with the algorithm for calculating the circumferential component of the absolute speed by means of a simplified radial equi- librium (3), which does not take into consideration the radial component of the speed vector (see § 3). The resulting analytical model (consistent with the theory of the authors [2]) had a constant axial speed component along the radius, which corresponds to the cylindrical flow areas (similar to the propeller in Fig. 2). This model applies to the computation plane located at a sufficient distance behind the impeller, where the flow field is stabilized. How- ever, measurements close to the trailing edges of the blades showed non-zero radial speed components that do not conform to the analytical model. (3.) To achieve better correspondence of the calcu- lations with the measurements, the computation plane of the impeller will be placed in a position corresponding to the intermediate position of the inlet openings of the measuring probe. This will supplement the radial equilibrium equation with members containing the radial component of the flow speed and meridian radii of the curvature of the rotary flow areas (see the closing part of § 3). (4.) Neither the basic condition of axially symmetrical flow, nor the rotary limit area of the flow tubes (see article 3) is met at the present time [12, 13]. However, the mathematical model that is applied is not able to describe the necessary changes. 7. Conclusions The results of the comparison between the calculated aerodynamic characteristics and the measured aero- dynamic characteristics of the V936mod model fan in Figs. 6 and 5 are not yet completely satisfactory. However, they will serve as a good basis for ongoing research. Further development in 2015 of the comput- ing programme for the fan aerodynamics, FanA23_7, will focus on a numerical solution of the completed equation of the radial balance (10). This equation was derived in § 3, where the possible solution is also indi- cated, including the introduction of another constant parameter. It will also be necessary to complement the program with algorithms for calculating the complete stage, including the stator. References [1] Dostál, J.: Aerodynamický výpočet axiálního ventilátoru - manuál programu FanA22_7 ; Report of VZLÚ R-6093 [2] Cyrus, V.: Aerodynamic performance of rear axial compressor stage with annular diffuser and outlet chamber. Proceedings of ASME 98-GT-238; pp. 10 http://dx.doi.org/10.1115/96-gt-238 [3] Cyrus, V., Cyrus, J., Wurst, P., Panek, P.: Aerodynamic performance of advanced axial flow fan for power industry within its operational range, Proceedings of ASME GT2014-25339, Dusseldorf, Germany, pp.13 http://dx.doi.org/10.1115/gt2014-25339 [4] Cyrus, V.: The effect of the Reynolds number on the three-dimensional flow in a straight compressor cascade. Proceedings of ASME 88-GT-269; pp. 9 http://dx.doi.org/10.1115/88-gt-269 [5] Cyrus, V.: Prediction of secondary losses in axial compressors. Proceedings of ASME 92-GT-288; pp. 10 http://dx.doi.org/10.1115/92-gt-288 [6] Ušakov, K. A., Brusilovskij, I. V., Bušeľ, A. R.: Aerodynamika osových ventilátorů a jejich konstrukční prvky; Praha SNTL 1962; 444 pages (in Czech). [7] Čermák, J. a kolektiv: Ventilátory, Praha SNTL 1974; 412 pp. (in Czech). [8] Šmíd, V.: Ventilátory; Praha SNTL 1965, textbook, ČVUT FS, 100 s. (in Czech). [9] Nový, R.: Ventilátory; Praha Nakladatelství ČVUT 2007 (in Czech). [10] Kadrnožka, J.: Lopatkové stroje, Brno CERM 2003, 177 s. (in Czech) [11] Cyrus, V.: Aerodynamic performance of an axial compressor stage with variable rotor blades and variable guide vanes. Proceedings of ASME 98-GT-151, pp. 10 http://dx.doi.org/10.1115/98-gt-151 [12] Cyrus, V.: Experimental study of three-dimensional flow in an axial compressor stage. Proceedings of ASME 86-GT-118, pp. 9 http://dx.doi.org/10.1115/86-gt-118 [13] Beiler, M. G., Carolus, T. H.: Computation and Measurement of the Flow in Axial Flow Fans With Skewed Blades, Journal of Turbomachinery, vol. 121, Jan 01.1999, pp.59-66 378 Acta Polytechnica 55(6):373–378, 2015 1 Introduction 2 From propellers to fans 3 Basic assumptions and equations of the mathematical model of fan aerodynamics 4 Measuring unit 5 Aerodynamic characteristics of the V936mod fan 6 Discussion on the correspondence between the calculations and the measurements 7 Conclusions References