Analytical and numerical investigation of a small-scale radial-inflow steam turbine

Analytical and Numerical Investigation of a Small-scale Radial-Inflow Steam Turbine ÉCOLE POLYTECHNIQUE FÉDÉRALE DE LAUSANNE SCHOOL OF ENGINEERING INSTITUTE OF MECHANICAL ENGINEERING LABORATORY FOR APPLIED MECHANICAL DESIGN Author: Marcos Font Armenteras Supervisors: Professor J. Schiffmann P. H. Wagner Lausanne, Switzerland January 2019 1 2 ABSTRACT This thesis compares experiments to the analytical and numerical investigation of a small-scale radialinflow turbine (diameter of 15 mm). This turbine is part of the Fan-Turbine Unit (FTU). The turbine of this FTU propels the fan that recirculates the unused hydrogen and water vapor from the anode off-gas of a Solid Oxide Fuel Cell (SOFC) to its inlet. This recirculation improves the electrical efficiency of the SOFC system due to higher global fuel utilization and allows for an operation without external water supply for the steam reformer. Due to its high efficiencies, also at small-scale and the possibility for heat cogeneration, it is competitive compared with other energy generation systems. The main objective of this thesis is to improve the existing analytical and numerical models of the turbine and compare their respective results with experimental measurements. Good results are obtained with the analytical (38 Watt) and the numerical investigation (43 Watt), if the gas film bearing mechanical loss (18 Watt) is added to the experimental measurement (22 Watt). The turbine power at the design point is very low (40 Watt) compared to the bearing mechanical losses (18 Watt), so heat fluxes to the turbine impeller domain have a high impact. 3 CONTENTS ABSTRACT ............................................................................................................................... 2 LIST OF FIGURES ................................................................................................................... 5 LIST OF TABLES ..................................................................................................................... 7 NOMENCLATURE ................................................................................................................. 8 GREEK SYMBOLS .............................................................................................................................................. 8 ROMAN SYMBOLS ............................................................................................................................................. 8 1 INTRODUCTION .......................................................................................................... 10 1.1 MOTIVATIONS .................................................................................................................................... 10 1.2 CONTEXT AND OBJECTIVES ............................................................................................................ 10 2 THEORY .......................................................................................................................... 13 2.1 SOLID OXIDE FUEL CELL SYSTEM .................................................................................................. 13 2.2 TURBINE ............................................................................................................................................. 14 3 ANALYTICAL CALCULATIONS .................................................................................. 16 3.1 TURBINE POWER AND EFFICIENCY ............................................................................................... 16 3.2 PRESSURE LOSSES .............................................................................................................................. 20 3.2.1 Inlet pressure loss ........................................................................................................................... 20 3.2.2 Outlet pressure loss ......................................................................................................................... 24 3.3 VELOCITY ANGLE AT THE VOLUTE OUTLET ............................................................................... 25 4 EXPERIMENTAL MEASUREMENTS ......................................................................... 28 4.1 TURBINE EXPERIMENT AND MEASUREMENTS ............................................................................ 28 5 NUMERICAL CALCULATIONS ................................................................................... 30 5.1 TURBINE INLET AND VOLUTE SIMULATION ................................................................................ 30 5.1.1 Geometry ....................................................................................................................................... 30 5.1.2 Mesh .............................................................................................................................................. 30 5.1.3 Results for the design point (168 krpm) .......................................................................................... 31 5.1.4 Comparison between numerical and analytical results for the design point (168 krpm) ..................... 33 5.2 TURBINE SINGLE PASSAGE SIMULATION ...................................................................................... 34 5.2.1 Stationary simulation ..................................................................................................................... 34 5.2.2 Transient simulation ...................................................................................................................... 42 4 6 CONCLUSION ................................................................................................................ 45 6.1 SUMMARY ........................................................................................................................................... 45 A. MEAN ROUGHNESS EXPERIMENTAL MEASUREMENT (RA) ........................ 46 B. OTHER RESULTS CFD SIMULATIONS ................................................................. 48 C. MATLAB SCRIPTS ...................................................................................................... 50 C.1. VOLUTE INLET PRESSURE LOSS .......................................................................................................... 50 C.2. VOLUTE OUTLET VELOCITY ANGLE ................................................................................................. 51 C.3. VOLUTE OUTLET PRESSURE LOSS ...................................................................................................... 52 C.4. TURBINE POWER AND EFFICIENCY ................................................................................................... 54 BIBLIOGRAPHY ................................................................................................................... 56 5 LIST OF FIGURES Figure 1.1 - Meridional view of a schematic partial-admission radial-inflow turbine, as well as the nomenclature of the different sections and important components. Adapted from [1] ............... 11 Figure 2.1 – Stationary SOFC systems with steam reformer, schematic SOFC cell, burner, condenser, evaporator and anode off-gas recirculation propelled by a) a steam-driven ejector, b) a fuel-driven ejector, c) an electrically-driven fan, d) a patented thermally-driven fan, as well as e) direct steam supply for the steam reforming. [1] ........................................................................................................ 13 Figure 3.1 - Top view of a schematic turbine and the velocity triangles at the leading and trailing edge ...................................................................................................................................................................... 16 Figure 3.2 – Breakdown of the velocity vectors into the meridional (cm and wm) and transversal (cu and wu) components at the LE (section 7) and TE (section 8) of the rotor ........................................... 17 Figure 3.3 – The manufactured partial-admission turbine volute (without turbine stator) .................... 21 Figure 3.4 – Schematic turbine inlet with the nomenclature of the different sections. Adapted from [1]. ...................................................................................................................................................................... 21 Figure 3.5 – Tunnel-type volute section shape and nomenclature ............................................................. 22 Figure 3.6 – Graphic of conversion for the inlet pressure loss for the different number of discretization elements per part ....................................................................................................................................... 24 Figure 3.7 – Turbine volute with the removable turbine stator. [1] ........................................................... 26 Figure 3.8 – Schematic view of the turbine volute ....................................................................................... 27 Figure 4.1 – Overview of the Fan-Turbine Unit test rig for hot air at 200 ºC (fan inlet) without glass fiber insulation (left) and real implementation partly covered with glass fiber insulation without the oven cover. [1]..................................................................................................................................... 28 Figure 5.1 - Turbine volute and inducer geometry for the CFD simulation ............................................ 30 Figure 5.2 – Inlet volute and inducer domains with the generated mesh ................................................. 31 Figure 5.3 - Pressure loss graphic from the turbine inlet to the volute outlet (sections 1 to 4) of the numerical model compared to the analytical results (at 168 krpm) .................................................. 31 Figure 5.4 - Velocity angle at the volute outlet (section 4) graphic of the numerical model compared to the analytical result (at 168 krpm) ........................................................................................................... 32 Figure 5.5 – Velocity vectors at the volute outlet ......................................................................................... 32 Figure 5.6 – Schematic view of the turbine outlet with the velocity vectors calculated analytically (black) and numerically (red) ................................................................................................................................ 33 Figure 5.7 – Domain regions (right): inducer, stator, rotor and exducer domain from the top to the bottom and generated mesh (left) at the hub for each domain from the bottom view. Fluid-to- fluid interfaces are marked in green. [1] ................................................................................................ 35 6 Figure 5.8 – Optical microscopy with Hirox KH-8700 (left) of the turbine stator and rotor (digitally mirrored) and overview of turbine stator and rotor (upside-down) with turbine inducer (right). [1] ..................................................................................................................................................................... 37 Figure 5.9 - Turbine power graphic comparison between the experiments and the CFD simulations39 Figure 5.10 – Degree of reaction graphic comparison between the experiments and the CFD simulations ................................................................................................................................................. 40 Figure 5.11 – Isentropic efficiency graphic comparison between the experiments and the CFD simulations ................................................................................................................................................. 40 Figure 5.12 – Hole made between the rotor and the stator to measure the static pressure experimentally ..................................................................................................................................................................... 41 Figure 5.13 - Rotor-stator surface, simulating the real hole of the turbine, to calculate the static pressure at the CFD simulation (the span of the surface is higher for visualization, normally it is placed only at the shroud surface, span 0,99 to 1) ........................................................................................... 41 Figure 5.14 - Rotor-stator surface to calculate the static pressure at the CFD simulation with a higher span for visualization (for the calculation it is placed at the shroud surface, span 0,99 to 1) ...... 42 Figure 5.15 - Specific power graphic comparison between the experimental measurement, the stationary simulation and the transient simulation .............................................................................. 43 Figure 5.16 - Degree of reaction graphic comparison between the experimental measurement, the stationary simulation and the transient simulation (rotor-stator static pressure measured with the surface shown in Figure 5.13 with span 0.99 to 1).............................................................................. 43 Figure 5.17 - Degree of reaction graphic comparison between the experimental measurement, the stationary simulation and the transient simulation (rotor-stator static pressure measured as in Figure 5.15) ................................................................................................................................................ 44 Figure A.1 - Experimental measurement of the inlet pipe surface roughness with TESA Rugosurf 90G ..................................................................................................................................................................... 46 Figure A.2 – TESA Rugosurf 90G’s display after one of the measurements .......................................... 47 Figure B.1 – Boundary conditions for the different simulations ............................................................... 48 Figure B.2 – Static and total pressures for the different iterations ............................................................ 48 Figure B.3 – Power, degree of reaction and efficiency results for the different simulations ................ 48 7 LIST OF TABLES Table 3.1 - Geometrical parameters of the turbine needed for the analytical calculations .................... 18 Table 3.2 - Different initial conditions for the tested situations of the inlet analytical model .............. 23 Table 3.3 - Results of the analytical model for the different situations ..................................................... 23 Table 3.4 – Different initial conditions for the tested situations of the outlet analytical model ........... 25 Table 3.5 - Results of the analytical model for the different situations ..................................................... 25 Table 4.1 – Parameters measured during the turbine experiments by Wagner [1] .................................. 29 Table 5.1 - Analytical and numerical results comparison for the design point (168 krpm) ................... 33 Table 5.2 - CFD boundary conditions used for the nominal working point ............................................ 36 Table 5.3 – Total pressure calculated for the outlet boundary condition and total pressure at the outlet of the CFD simulation for each iteration at the design point (168 krpm) ....................................... 38 Table 5.4 - Experimental measurements and CFD simulation results for the different iterations for the degree of reaction and power (corrected 𝑝𝑠,9) at the design point (168 krpm) ............................ 39 8 NOMENCLATURE GREEK SYMBOLS ROMAN SYMBOLS  Absolute velocity angle rad  Relative velocity angle rad  Efficiency -  Absolute roughness mm 𝜖𝑎 Admission ratio - 𝑓𝐷 Friction coefficient -  Dynamic viscosity Pa·s  Degree of reaction -  Density kg/m3  Angular velocity rad/s A Area m2 b Channel width m c Absolute velocity m/s cs Sound velocity m/s Cu Transversal absolute velocity m/s Cm Meridional absolute velocity m/s Cp Heat capacity at constant pressure J/K C Courant number - D Diameter m d Distance between blades m h Blade height m h Specific enthalpy J/kg k Heat capacity ratio - 15 𝜂𝑖𝑠,𝑡𝑢𝑟𝑏𝑖𝑛𝑒=ℎ𝑡,8−ℎ𝑡,5 ℎ𝑡,𝑖𝑠,8−ℎ𝑡,5 (2.4) The stator LE is at plane 5, whereas the TE is at plane 6 (see Figure 1.1). The specific enthalpy difference can also be calculated with the temperature difference and the heat capacity at constant pressure (Cp) Δℎ𝑡𝑡=𝐶𝑝·(𝑇7−𝑇8) (2.5) Thus, the turbine power can also be expressed as 𝑃𝑡𝑢𝑟𝑏𝑖𝑛𝑒=𝑚󰇗8·𝐶𝑝·(𝑇7−𝑇8) (2.6) The degree of reaction is a non-dimensional number used for the turbine design. 𝛿ℎ=∆ℎ𝑟𝑜𝑡𝑜𝑟 ∆ℎ𝑠𝑡𝑎𝑔𝑒 (2.7) The degree of reaction is the fraction of the rotor specific enthalpy difference to the stage (rotor and stator) specific enthalpy difference. This reaction is always between 0 and 1, whereas turbines with a reaction of 0 or 1 are not achievable since there are always some frictional pressure losses within the rotor and the stator. Turbines with higher efficiencies are built at a degree of reaction of 0,5 [2]. Due to the high static pressure between the rotor and stator, a high reaction turbine is not capable for a partial-admission operation. The power density of a low reaction turbine is higher (about twice as much as a 0,5-reaction turbine with a similar mass flow rate). The 8 and 7 are more favorable, respectively the change of direction of the absolute velocity is higher. A higher change of velocity also involves higher aerodynamic losses; hence the low-reaction turbines have an inferior efficiency. Additionally, losses occur due to the increased velocities which are higher at the stator outlet, respectively at the rotor inlet compared to the 0,5 reaction turbines. 16 3 ANALYTICAL CALCULATIONS This chapter shows the analytical calculations made for the turbine power and efficiency, as well as the inlet and outlet pressure losses and the volute outlet velocity angle. 3.1 TURBINE POWER AND EFFICIENCY In this section, the turbine power is calculated with the two possible methods: (1) with the Euler equation (3.1) and the design parameters of the turbine and (2) with the first law of thermodynamics according to equation (2.6). The boundary conditions for the analytical calculations shown in this section are at the design point of the turbine (168 krpm). To calculate the turbine power with the Euler equation (3.1) it is necessary to calculate the velocity triangles of the turbine rotor 𝑃𝑡𝑢𝑟𝑏𝑖𝑛𝑒=𝑚󰇗8·Δℎ𝑡𝑡=𝑚󰇗8·(𝑐𝑢,8·𝑢8−𝑐𝑢,7·𝑢7) (3.1) The colour triangles (Figure 3.1) formed by the velocity vectors “u”, “c” and “w” explain how a turbine works and are used to calculate the turbine power with the Euler equation. The “u” vectors are the circumferential velocities of the blade, the “c” vectors are the absolute velocities of the fluid and the Figure 3.1 - Top view of a schematic turbine and the velocity triangles at the leading and trailing edge 17 “w” vectors are the relative velocities of the fluid with respect to the rotating blades. The absolute velocity is equal to the relative velocity plus the blade’s circumferential velocity. The inlet part of the rotor is the Leading Edge (7) and the outlet is the Trailing Edge (8). Figure 3.2 – Breakdown of the velocity vectors into the meridional (cm and wm) and transversal (cu and wu) components at the LE (section 7) and TE (section 8) of the rotor The boundary conditions needed to calculate this power analytically are the inlet and the outlet total pressures of the fluid, the inlet and the outlet total temperatures and the mass flow rate. Since the pressure ratio of the turbine (2,8) is higher than the critical pressure ratio (1,9), the Mach number in the critical section can be assumed equal to 1. This critical section is at the stator outlet, where the fluid velocity is equal to the sound velocity. The sound velocity (𝑐𝑠) is defined by the following equation (3.2) 𝑐𝑠=√𝐾·𝑅·𝑇 (3.2) Then, considering the inlet volute adiabatic and using the temperature at the inlet, the sound speed at the stator TE is 445 m/s. Thus, the fluid velocity at section 6 (𝑐6) is also 445 m/s. 18 Table 3.1 - Geometrical parameters of the turbine needed for the analytical calculations As shown in Figure 3.1  is the angle between the absolute velocity and the circumferential velocity of the blade and  is the angle between the relative velocity and the circumferential velocity of the blade. The components of the velocity triangle (Figure 3.1 and Figure 3.2) are calculated with the geometry parameters shown in Table 3.1. 𝑐𝑢,6=𝑐6·cos(𝛼6)=411,5 𝑚/𝑠 (3.3) 𝑐𝑢,7=𝑐𝑢,6·𝑟6 𝑟7=419,7 𝑚/𝑠 (3.4) 𝑢7=·𝑟7=132,1 𝑚/𝑠 (3.5) 𝑐𝑚,8=𝑚󰇗 3600·𝜌·𝐴8·10−6=162,9 𝑚/𝑠 (3.6) Parameter Nomenclature Value TE stator blade angle at section 6 𝛼6 0,38 𝑟𝑎𝑑 LE rotor blade angle at section 7 𝛽7 0,64 𝑟𝑎𝑑 Radius at section 7 𝑟7 7,5·10−3 𝑚 TE rotor blade angle at section 8 𝛽8 2,41 𝑟𝑎𝑑 Radius at section 8 𝑟8 6,5·10−3 𝑚 Distance between blades 𝑑8 0,695 𝑚𝑚 Height of the blades ℎ8 0,586 𝑚𝑚 Area between blades 𝐴8 5,3 𝑚𝑚2 Rotational speed of the rotor  17.614 𝑟𝑎𝑑/𝑠 19 𝑤8=𝑐𝑚,8 sin (𝜋 − 𝛽8)=244,4 𝑚/𝑠 (3.7) 𝑤𝑢,8=𝑤8·cos(𝜋−𝛽8)=182,2 𝑚/𝑠 (3.8) 𝑢8=𝑤·𝑟8=114,5 𝑚/𝑠 (3.9) 𝑐𝑢,8=𝑤𝑢,8−𝑢8=67,7 𝑚/𝑠 (3.10) With all the velocity components known it is possible to calculate the turbine power with the Euler equation 𝑃𝑡𝑢𝑟𝑏𝑖𝑛𝑒=𝑚󰇗·(𝑐𝑢,7·𝑢7−𝑐𝑢,8·𝑢8)=37,9 𝑊 (3.11) This turbine power can also be calculated with the heat capacity at constant pressure and the decrease of temperatures 𝑃𝑡𝑢𝑟𝑏𝑖𝑛𝑒=𝑚󰇗·𝐶𝑝·∆𝑇=21,9 𝑊 (3.12) This power is calculated with the measured temperature difference between the inlet (section 1 in Figure 1.1) and the outlet (section 12) of the turbine. The difference between both methods (16 W) is caused by the difference between the temperature at the outlet of the turbine (section 12) and at the outlet of the turbine rotor (section 8). Between the rotor outlet and the turbine outlet, a heat flux due to the bearing mechanical losses crosses the turbine fluid domain and thus the turbine outlet temperature is increased. This is the main reason of the difference between both methods for the power calculation. This heat flux coming from the bearings was measured and calculated by Wagner [1] and it is around 18 W. Then, assuming that this heat flux crosses the turbine, the adjusted turbine power is 39,9 W. With this corrected power the corrected total temperature at the rotor outlet is calculated as (3.13) 𝑇𝑜𝑢𝑡,𝑐𝑜𝑟𝑟𝑒𝑐𝑡𝑒𝑑=𝑇𝑖𝑛−𝑃𝑐𝑜𝑟𝑟𝑒𝑐𝑡𝑒𝑑 𝑚󰇗·𝐶𝑝 (3.13) 20 This corrected total temperature at the outlet is 171 ºC. The efficiency of the turbine is then calculated with this corrected temperature. 𝜂𝑖𝑠,𝑡𝑢𝑟𝑏𝑖𝑛𝑒=ℎ𝑡,8−ℎ𝑡,5 ℎ𝑡,𝑖𝑠,8−ℎ𝑡,5=(𝑇𝑡,𝑖𝑛−𝑇𝑡,𝑜𝑢𝑡,𝑐𝑜𝑟𝑟𝑒𝑐𝑡𝑒𝑑) 𝑇𝑡,𝑖𝑛·(1−(𝑃𝑠,𝑜𝑢𝑡 𝑃𝑡,𝑖𝑛)𝑘−1 𝑘) (3.14) The results of the isentropic efficiency with the corrected outlet temperature is 43,6 %, whereas with the non-corrected temperature it is 24 % (with the measured outlet temperature of 193,1 ºC). The Matlab scripts performed to automatize these analytical calculations are shown in Appendix C. 3.2 PRESSURE LOSSES The already existing numerical model for this turbine accounts for the inducer inlet (4) until the exducer outlet (9). In this model, it is assumed that there are no pressure losses between the turbine inlet (1) and the inducer inlet (4) and between the exducer outlet (9) and the turbine outlet (12). Although these losses are extremely low they should be considered to be as much accurate as possible. These losses are calculated in this section. The inputs needed for the existing numerical model of this turbine are the inducer inlet (section 4) total pressure and total temperature, the exducer outlet (section 9) static pressure, the rotational speed of the rotor and the inducer inlet (section 4) velocity angle which is the angle between the velocity and the volute outlet normal vector. These conditions were measured experimentally at the inlet (1) and at the outlet (12) of the turbine. For the case of the inlet velocity angle it was calculated analytically. The initial conditions for the numerical model are calculated accounting the pressure losses of the inlet and the outlet pipes of the turbine (from section 1 to 4 and from 9 to 12). 3.2.1 INLET PRESSURE LOSS For the inlet pressure loss calculation, the initial conditions needed are the mass flow rate, the inlet (section 1) total pressure and the inlet (1) total temperature (all of them measured experimentally). Figure 3.3 shows a photo of the manufactured prototype of the turbine volute. 21 Figure 3.3 – The manufactured partial-admission turbine volute (without turbine stator) For the analytical calculation of the turbine inlet pressure loss, the turbine inlet is divided into four sections (Figure 3.4). The first section (1-2) has a constant area and the other three sections (2 to 3, 3 to 3.1 and 3.1 to 3.2) are not constant, hence a discretization is necessary. The highest pressure-loss is in the volute (3.1 to 3.2) due to the higher change of area and the curve of the volute. Figure 3.4 – Schematic turbine inlet with the nomenclature of the different sections. Adapted from [1]. Inlet Outlet 22 The pressure loss of a pipe is calculated as ∆𝑝=12·𝑓𝐷·𝜌·𝑐2·𝐿 𝐷 (3.15) The friction factor (𝑓𝐷) used in equation (3.15) is obtained from the Moody Diagram, knowing first the Reynolds number (Re) and the relative roughness (𝜀 𝐷) of the pipe. The Reynolds number is a function of the fluid density, the dynamic viscosity, the flow velocity and the diameter of each part (equation (3.16). 𝑅𝑒=𝜌·𝑐·𝐷 𝜇 (3.16) The Moody Diagram is approximated by equation (3.17) in the straight pipes [3] 𝑓𝐷=0,0055·[1+(2·104· 𝐷+106 𝑅𝑒)1 3] (3.17) and by equation (3.18) in the curved pipes [4] 𝑓𝐷=0,314 0,95·𝑅𝑒0,25+0,0075·√ 𝐷 2·𝑟ℎ (3.18) where the Reynolds number is calculated with the equation (3.16) for each of the discretized elements and 𝑟ℎ is the curve radius. 𝐷 is the equivalent diameter of the pipe for each section and it is calculated as 𝐷=𝑅3.1+𝑏3.1 (3.19) Figure 3.5 – Tunnel-type volute section shape and nomenclature 23 To know the relative roughness of the material (𝜀 𝐷) it is necessary to know the absolute roughness () to be able to divide it by the diameter. The absolute roughness has been calculated with the following simplified expression (3.20) from Adams, Grant and Watson [5] with Ra (mean roughness). 𝜀=11,03·𝑅𝑎 (3.20) This mean roughness (Ra) has been measured experimentally with TESA Rugosurf 90G and it is shown and explained in appendix A. Once the absolute roughness is known, the friction factor can be calculated with the expressions shown before (3.17) (3.18). As it was stated before, the sections with non-constant area are discretized. Figure 3.6 shows the relation between the pressure loss at the turbine inlet (from section 1 to 4) and the number of discretization elements used per part. The result converges for 1000 elements per part. The pressure loss is also calculated for different initial conditions (mass flow rate and inlet total pressure) to verify that the model works for different situations. The different initial conditions tested are shown in Table 3.2, and the results are shown in Table 3.3. Table 3.2 - Different initial conditions for the tested situations of the inlet analytical model 1st situation 2nd situation 3rd situation Mass flow [kg/h] 2,86 1,87 1,14 Inlet pressure [bar] 2,74 1,94 1,40 Table 3.3 - Results of the analytical model for the different situations 1st situation 2nd situation 3rd situation Pressure loss section 1-4 [mbar] 8,9 5,7 3,2 24 Figure 3.6 – Graphic of conversion for the inlet pressure loss for the different number of discretization elements per part These results are used in section 5.2 to improve the inlet boundary conditions of the turbine numerical model. 3.2.2 OUTLET PRESSURE LOSS For the outlet part of the turbine the process has been similar as for the inlet. The outlet is divided into different sections to simplify the calculations depending on the shape of each section. The expressions used for the outlet pressure loss are the same as for the inlet (3.15), (3.16), (3.17) and (3.18), the only exception is the sudden expansion loss of the transition between the exducer and the diffuser (section 9) and from the outlet of the pipe to the ambient (section 12). This sudden expansion loss in the transition part has been calculated with the following expression (3.21) [6] ∆𝑝=𝜌·𝐴1 𝐴2·(1−𝐴1 𝐴2)·𝑐12 (3.21) where section 1 is the area before the transition and section 2 is the area after the transition. A1 is calculated with the factor of 33 % because it is considered that the fluid still has not expanded throughout the whole circumference. However, in section 2 the fluid is expanded to the whole tube. The expressions used to calculate these areas are 𝐴1=𝜋·(𝑟12−𝑟22)·0,33 (3.22) 𝐴2=𝜋·𝑟12 (3.23) 31 maximum element size and the same area, they will approximately have the same number of elements and they will match correctly. To ensure that both meshes match correctly the simulation has been solved for several meshes with different number of elements. 5.1.3 RESULTS FOR THE DESIGN POINT (168 KRPM) To verify the model and to ensure that the mesh is correct the simulation is solved for different number of mesh elements. The results are compared between the analytical model and the numerical simulation. Figure 5.3 - Pressure loss graphic from the turbine inlet to the volute outlet (sections 1 to 4) of the numerical model compared to the analytical results (at 168 krpm) 8 9 10 11 100 1.000 10.000 100.000 1.000.000 10.000.000 Pressure loss [mbar] Number of discretization elements Numerical model Analytical model Figure 5.2 – Inlet volute and inducer domains with the generated mesh 32 Figure 5.4 - Velocity angle at the volute outlet (section 4) graphic of the numerical model compared to the analytical result (at 168 krpm) Figure 5.3 and Figure 5.4 show respectively, the results of the pressure loss and the volute outlet velocity angle for different number of mesh elements. The simulation is solved for 28k, 35k, 125k, 350k, 900k and 1,8m mesh elements. The number of elements that optimized the accuracy of the results with a reasonable solving time is the one with 900k elements. The results of the simulation are 8,7 mbar (900k mesh elements) for the inlet pressure loss and 56,1º for the volute outlet velocity angle for the design situation (168 krpm). Figure 5.5 – Velocity vectors at the volute outlet 55 56 57 58 10.000 100.000 1.000.000 10.000.000 Volute outlet flow angle [º] Number of discretization elements Numerical model Analytical model 33 Figure 5.5 shows how the velocity vectors at the volute outlet are distributed equally. This would mean that the turbine volute has been designed correctly. 5.1.4 COMPARISON BETWEEN NUMERICAL AND ANALYTICAL RESULTS FOR THE DESIGN POINT (168 KRPM) After calculating the inlet pressure loss and the velocity angle with the analytical and the numerical models, they are compared to see the difference between them and to verify the models. Table 5.1 - Analytical and numerical results comparison for the design point (168 krpm) Analytical Numerical Difference Inlet pressure loss [mbar] 8,9 8,8 1,1 % Volute outlet velocity angle [º] 57,1 56,2 1,6 % As shown in Table 5.1 the results of both models are very similar, with an error of 1,1 % for the inlet pressure loss and 1,6 % for the velocity angle at the volute outlet (both based on the analytical value). The difference between the analytical and the numerical results are due to the assumptions made at the analytical calculation. These assumptions are explained in section 3.3. As in the numerical model the friction between the fluid and the walls of the pipe is considered, the angle is expected to be lower, Figure 5.6 – Schematic view of the turbine outlet with the velocity vectors calculated analytically (black) and numerically (red) 34 since 𝑟·𝐶𝑢 is not constant (see Figure 5.6). As the deviation between both models is very low, both numerical and analytical models are considered to be verified. Results are shown schematically in Figure 5.6, the analytical results are shown in black and the numerical results in red. The results used for the definition of the inlet boundary conditions of the numerical turbine model of the following sections will be the numerical ones, since these are considered as more accurate. The analytical model will be used to define the outlet boundary condition of the turbine. 5.2 TURBINE SINGLE PASSAGE SIMULATION The results of the inlet pressure loss, the outlet pressure loss, and the velocity angle at the volute outlet are used to solve the complete CFD simulation for the turbine domain. The boundary conditions at the inlet and at the outlet are adapted with the results of the previous chapters. The simulation only models one single passage including one stator blade and one rotor blade. This is explained in detail in section 5.2.1. 5.2.1 STATIONARY SIMULATION The simulation domain (see Figure 5.7) consists of the turbine inducer, the turbine stator, the turbine rotor and the turbine exducer. The stator has 12 blades and features an admission of 12+1 61 =0,21. Although the real turbine rotor has 59 blades, in the CFD simulation it has 61 blades to realize a frozen rotor interface between the stator and the rotor with no pitch change. The rotor pitch in the simulation is therefore 3,3% smaller than the pitch of the real turbine rotor. The inlet boundary conditions (section 4) are the total pressure, the total temperature and the absolute flow angle. The turbine inlet from the turbine section 1 to 4 is assumed adiabatic and thus the total temperature is constant from section 1 to 4, this total temperature was measured experimentally. The pressure loss between 1 and 4 (∆𝑝1−4) is calculated in section 5.1 and is 8,7 mbar (8,9 mbar analytically) for the nominal working point. So, the CFD inlet total pressure (𝑝𝑡,𝑖𝑛𝑙𝑒𝑡,𝐶𝐹𝐷) is the total pressure measured experimentally at section 1 (𝑝𝑡,𝑖𝑛𝑙𝑒𝑡) minus the pressure loss calculated (equation (5.1). 𝑝𝑡,𝑖𝑛𝑙𝑒𝑡,𝐶𝐹𝐷=𝑝𝑡,4=𝑝𝑡,𝑖𝑛𝑙𝑒𝑡−∆𝑝1−4 (5.1) The absolute flow angle at the inlet is calculated in section 3.3 analytically and in section 5.1 numerically. The result used for the turbine simulation is the numerical one and is 56,1º (angle between the absolute flow and the normal vector of the volute outlet surface). 35 Figure 5.7 – Domain regions (right): inducer, stator, rotor and exducer domain from the top to the bottom and generated mesh (left) at the hub for each domain from the bottom view. Fluid-to-fluid interfaces are marked in green. [1] The outlet boundary condition (turbine section 9) is the static pressure (𝑝𝑠,9), hence the mass flow rate through the passage is the result of the simulation. As for the inlet section, the outlet boundary condition takes the pressure loss that is calculated in section 3.2 into account (sudden expansion loss and frictional loss). This outlet pressure loss (∆𝑝9−12) is 16,2 mbar for the nominal point. So, the static outlet pressure is the outlet pressure measured experimentally (ambient pressure) plus the pressure loss from 9 to 12 calculated in the previous chapters minus the dynamic pressure (pressure caused by the velocity at section 9). For the dynamic pressure the velocity is calculated in the first place with the mass flow rate and the area at section 9. 𝑝𝑠,𝑜𝑢𝑡𝑙𝑒𝑡,𝐶𝐹𝐷=𝑝𝑠,9=𝑝𝑡,𝑜𝑢𝑡𝑙𝑒𝑡+∆𝑝9−12−12·𝜌·𝑐92 (5.2) 36 The stationary exducer domain hub assumes a rotating wall with the rotational speed of the rotor domain. The rotational speed at the nominal working point is 168,2 krpm. Table 5.2 - CFD boundary conditions used for the nominal working point Value Units Rotational speed 168,2 [krpm] Inlet total pressure (4) 2,726 [bar] Inlet total temperature (4) 220,3 [ºC] Inlet absolute flow angle (4) 56,1 [º] Outlet static pressure (9) 0,926 [bar] The measured mass flow in the experiment is 2,86 kg/h and the result of this CFD simulation is 13·0,264 = 3.43 kg/h which is equal to a difference of 0,57 kg/h. This could be caused by three different reasons: 1- The distance between the stator blades is lower due to the manufacturing limitation. This leads to a measured minimum distance between stator blades between 0,18 and 0,23 (see Figure 5.8), whereas it is 0,23 in the CFD simulation. Thus, there is an error between 0-22% because of the difference of passing area. 2- The difference in the geometry between the CFD inducer and the reality. In reality the inducer side walls have a displacement of the boundary layer of around 0,04 mm. This results in a partial blocking of the first and the last stator blade rows and thus in an overall blockage of all the stator blade rows. The CFD single passage simulation does not account for these two side walls effects. 3- The manufactured stator blades have a roughness based on the arithmetical mean deviation (Ra) of about 0,002 mm resulting in a 1,2% lower critical area. 37 These three facts could be the reason of the higher mass flow rate for the simulation and thus a higher power. Figure 5.8 – Optical microscopy with Hirox KH-8700 (left) of the turbine stator and rotor (digitally mirrored) and overview of turbine stator and rotor (upside-down) with turbine inducer (right). [1] In order to compare the CFD simulation and the experiment, the power is corrected with the mass flow relation between the simulation and the experiment. 𝑃𝐶𝐹𝐷,𝑐𝑜𝑟𝑟𝑒𝑐𝑡𝑒𝑑=𝑃𝐶𝐹𝐷 𝑚󰇗𝐶𝐹𝐷·𝑚󰇗𝑒𝑥𝑝 (5.3) In the previous case (for the nominal conditions) the power in the CFD for single passage is 3,71 W and the mass flow in the CFD single passage simulation is 0,26 kg/h. With the experimental mass flow of 2,86 kg/h, the corrected power result from the CFD is 40,2 W. The power that was measured experimentally by Wagner [1] is 21,9 W. If the gas film mechanical loss (18 W) is added to the measured power, the turbine power is 39,9 W. The degree of reaction is another measurement used to compare the results of the experimental measurements with the numerical simulations. The degree of reaction is defined in first approximation 38 as the ratio of static pressure drop in the rotor to the static pressure drop in the stage, since the static enthalpy difference in the rotor and the static enthalpy difference in the stage are not measured. 𝛿ℎ=∆𝑝𝑠,𝑟𝑜𝑡𝑜𝑟 ∆𝑝𝑡,𝑟𝑜𝑡𝑜𝑟+𝑠𝑡𝑎𝑡𝑜𝑟 (5.4) In the experiments performed by Wagner [1] the degree of reaction for the nominal working point is 0,15. For this calculation, Wagner [1] measured experimentally the static pressure between the rotor and the stator in a specific hole made between them. At the CFD simulation, the degree of reaction for the same boundary conditions was 0,22, a 47% higher. This difference is caused by the difference of the total pressure at the outlet (𝑝𝑡,9) between the CFD simulation and the real experiment, and this difference of total pressure is mainly caused by the difference of velocity. This velocity is higher in the CFD because the area where the fluid goes through is constant through all the single passage (21% of the circumference) and in reality, it expands to a higher proportion of the circumference. To correct this error the outlet static pressure (𝑝𝑠,9) is adjusted with the outlet total pressure (𝑝𝑡,9), and the velocity (Mach number) of the previous CFD simulation (equation (5.5). 𝑝𝑠,9=𝑝𝑡,9·(1+𝑘−1 2·𝑀𝑎2)−𝑘 𝑘−1 (5.5) where 𝑝𝑡,9 is the total pressure at the outlet which is calculated adding the outlet pressure loss from section 3.2.2 to the ambient pressure measured experimentally, and 𝑀𝑎 is the Mach Number which is calculated at the previous iteration by the CFD simulation. This iterative process is done two times to see the behavior of the degree of reaction, the power of the turbine and the efficiency. Table 5.4 shows the results of the two iterations for the reaction, the turbine power and the efficiency at the design point (168 krpm). Table 5.3 – Total pressure calculated for the outlet boundary condition and total pressure at the outlet of the CFD simulation for each iteration at the design point (168 krpm) 𝒑𝒕,𝒄𝒂𝒍𝒄𝒖𝒍𝒂𝒕𝒆𝒅 𝒑𝒕,𝒔𝒊𝒎𝒖𝒍𝒂𝒕𝒊𝒐𝒏 CFD simulation 0,98 1,18 CFD iteration 1 (corrected 𝑝𝑠,9) 0,98 1,09 CFD iteration 2 (corrected 𝑝𝑠,9) 0,98 1,05 39 Table 5.4 - Experimental measurements and CFD simulation results for the different iterations for the degree of reaction and power (corrected 𝑝𝑠,9) at the design point (168 krpm) Reaction Power [W] Efficiency Experimental measurement 0,15 39,9 0,39 CFD simulation 0,22 40,2 0,38 CFD iteration 1 (corrected 𝑝𝑠,9) 0,18 42,7 0,35 CFD iteration 2 (corrected 𝑝𝑠,9) 0,16 43,4 0,45 As it can be observed in Table 5.4, while the degree of reaction decreases to move closer to the experimental measurement, the power of the turbine increases and the efficiency decreases and both move further from the experimental measurement. This is caused by the change of the outlet total pressure. The single-passage simulation is also solved for other boundary conditions that were measured previously by Wagner [1] and the same iterative procedure has been carried out for some of them. Results are shown in the following figures. Figure 5.9 - Turbine power graphic comparison between the experiments and the CFD simulations 0 10 20 30 40 50 1,5 1,8 2,1 2,4 2,7 3,0 Power [W] Mass flow turbine inlet [kg/h] Experimental CFD CFD iteration 1 (corrected Ps9) CFD iteration 2 (corrected Ps9) CFD Transient 40 Figure 5.10 – Degree of reaction graphic comparison between the experiments and the CFD simulations Figure 5.11 – Isentropic efficiency graphic comparison between the experiments and the CFD simulations Figure 5.9 shows the turbine power for the different boundary conditions. As it can be observed, the turbine power is more similar between the CFD and the experiments as closer to the design point (168 krpm and 2,86 kg/h). The same behavior can be observed in Figure 5.10 with the degree of reaction and in Figure 5.11 with the isentropic efficiency. These three figures show the results shown in Table 5.4 for the rest of the measured points. There is also a deviation between the CFD simulation and the experimental degree of reaction which is caused by the difference in measuring between both models. In the experiment, the rotor stator static pressure is measured when the rotor is running at 168 krpm at one specific hole between the rotor and the stator. 0,0 0,1 0,2 0,3 0,4 0,5 1,5 1,8 2,1 2,4 2,7 3,0 Reaction Mass flow turbine inlet [kg/h] Experimental CFD CFD iteration 1 (corrected Ps9) CFD iteration 2 (corrected Ps9) CFD Transient 0,0 0,1 0,2 0,3 0,4 0,5 1,5 1,8 2,1 2,4 2,7 3,0 Efficiency Mass flow turbine inlet [kg/h] Experimental CFD CFD iteration 1 (corrected Ps9) CFD iteration 2 (corrected Ps9) CFD Transient 47 As it is explained in section 3.2.1 the friction coefficient fD (see equations (3.17 and (3.18) is calculated with a simplified expression of the Moody Diagram. These expressions calculate the friction coefficient from the Reynolds number (Re) and the relative roughness of the material (𝜀 𝐷). To calculate the absolute roughness () from the mean roughness (Ra) it has been used the simplified expression from explained in chapter 3.2 from Adams, Grant and Watson [5]. Figure A.2 – TESA Rugosurf 90G’s display after one of the measurements 48 B. OTHER RESULTS CFD SIMULATIONS As shown in section 5.2.1, the stationary single-passage CFD simulation was run for the different situations measured in the experiments. This section shows the boundary conditions and the results for all these different situations. Figure B.1 – Boundary conditions for the different simulations Figure B.2 – Static and total pressures for the different iterations Figure B.3 – Power, degree of reaction and efficiency results for the different simulations Situation 1 2 3 4 5 6 7 8 Rotational speed [krpm] 99,76 109,88 119,44 131,04 140,25 150,21 159,90 168,16 Inlet total temperature (1) [ºC] 221,4221,4221,1220,6220,8221,3221,4220,3 Measured inlet total pressure (1) [bar] 1,693 1,816 1,938 2,096 2,216 2,364 2,557 2,735 Inlet pressure loss numerical model (1-4) [bar] 0,0050 0,0055 0,0059 0,0065 0,0069 0,0074 0,0080 0,0087 Inlet pressure loss analytical model (1-4) [bar] 0,0047 0,0052 0,0058 0,0064 0,0069 0,0075 0,0082 0,0088 Corrected inlet (4) total pressure [bar] 1,688 1,810 1,932 2,089 2,209 2,357 2,549 2,726 Volute outlet velocity angle numerical model (4) [º]56,156,156,156,156,156,156,156,1 Volute outlet velocity angle analytical model (4) [º]57,157,157,157,157,157,157,157,1 Measured outlet total pressure (12) [bar] 0,961 0,962 0,963 0,964 0,964 0,964 0,964 0,965 Outlet pressure loss analytical model (9-12) [bar] 0,0048 0,0058 0,007 0,0085 0,0098 0,0115 0,0139 0,0162 Outlet total pressure (9) [bar] 0,966 0,968 0,970 0,972 0,974 0,976 0,978 0,981 Boundary conditions Situation 1 2 3 4 5 6 7 8 Corrected outlet static pressure (it. 0) [bar] (Ps9,0)0,950 0,948 0,946 0,943 0,941 0,937 0,931 0,926 Corrected outlet static pressure (it. 1) [bar] (Ps9,1)0,824 0,810 0,794 0,775 Corrected outlet static pressure (it. 2) [bar] (Ps9,2)0,780 0,764 0,735 0,703 CFD iteration 0 (Ps9,0)1,0508 1,0678 1,084 1,1031 1,115 1,1306 1,150 1,1793 CFD iteration 1 (Ps9,1)1,0311 1,0387 1,0629 1,0885 CFD iteration 2 (Ps9,2)1,0017 1,0055 1,0291 1,0543 CFD Transient (Ps9,2)1,0543 Total pressure [bar] Static pressure [bar] Situation 1 2 3 4 5 6 7 8 Experimental measurement 14,04 16,62 19,30 22,87 26,01 30,08 35,04 39,86 CFD iteration 0 (Ps9,0)10,28 13,15 16,20 20,46 24,11 28,55 34,58 40,18 CFD iteration 1 (Ps9,1)25,62 30,36 36,67 42,69 CFD iteration 2 (Ps9,2)26,14 31,07 37,68 43,40 CFD Transient (Ps9,2)42,38 Experimental measurement 0,116 0,121 0,125 0,131 0,136 0,139 0,146 0,151 CFD iteration 0 (Ps9,0)0,223 0,225 0,228 0,230 0,227 0,225 0,217 0,221 CFD iteration 1 (Ps9,1)0,173 0,170 0,173 0,176 CFD iteration 2 (Ps9,2)0,156 0,149 0,152 0,165 CFD Transient (Ps9,2)0,165 Experimental measurement 0,438 0,423 0,412 0,402 0,400 0,400 0,393 0,392 CFD iteration 0 (Ps9,0)0,313 0,325 0,335 0,349 0,360 0,368 0,376 0,381 CFD iteration 1 (Ps9,1)0,336 0,344 0,350 0,354 CFD iteration 2 (Ps9,2)0,327 0,336 0,340 0,338 CFD Transient (Ps9,2)0,458 Efficiency Power [W] Degree of reaction 49 Figure B.3 show the results of the power, the degree of reaction and the efficiency of all the CFD simulations run in this thesis. The CFD simulation with an assumed outlet velocity was run for the 8 different boundary conditions shown in Figure B.1. However, the iterative process to correct the outlet total pressure was only run for the last 4 situations, whereas the transient simulation just for the 8th simulation (design point: 168 krpm). Figure B.2 show the different static pressures (outlet boundary condition) for these simulations and the result of the total pressure at section 9, to see the difference with the calculated total pressure (in Figure B.1). 50 C. MATLAB SCRIPTS C.1. VOLUTE INLET PRESSURE LOSS %to calculate the pressure loss from the inlet of the turbine to the %inducer inlet (1-4) %inputs: T_inlet, P_inlet, m_dot, Ra (Roughness) m_dot=2.859; %input T_inlet=220.31; %input P_inlet=2.7347; %input Ra=3.172; %input epsilon=Ra*10^-3*11.03; %absolute roughness P_atm=0.9643; %atmospheric pressure Prel_inlet=P_inlet-P_atm; %inlet relative pressure %calculations R=286.9; %gas constant density=P_inlet*10^5/(R*(T_inlet+273.15)); %(kg/m^3) %section 1 to 2 r12=3; %(mm) A12=pi*r12^2; %(mm^2) c12=m_dot /(density*3600*A12*10^-6); %(m/s) D12=2*r12; %(mm) L12=34.25; %(mm) %moody Re12 = density*c12*D12*10^-3/(1.74*10^-5); %Reynolds number lambda12 = 0.0055*(1+(2*10^4*(epsilon/D12)+10^6/Re12)^(1/3)); %friction %coefficient from aproximation of the moody diagram AP12=lambda12*0.5*density*c12^2*L12*10^-5/D12; %pressure loss section 1- 2 APtotal_1=AP12; %total pressure loss Prel_2=Prel_inlet-AP12; %relative pressure at point 2 %section 2 to 3 %discretize in j sections APtotal_2=APtotal_1; AP23=0; j=1000000; for i=1:j r2i=r12-(i/j)*1; A2i=pi*r2i^2; c2i=m_dot/(density*3600*A2i*10^-6); D2i=r2i*2; L2i=21.866/j; Re2i = density*c2i*D2i*10^-3/(1.74*10^-5); lambda2i = 0.0055*(1+(2*10^4*(epsilon/D2i)+10^6/Re2i)^(1/3)); AP2i=lambda2i*0.5*density*c2i^2*L2i*10^-5/D2i; AP23=AP23+AP2i; APtotal_2=APtotal_2+AP2i; end 51 Prel_3=Prel_2-AP23; %section 3 to 3.1 %discretize in u sections APtotal_3=APtotal_2; AP34=0; u=1000000; for n=1:u r3n=2; A3n=((pi*r3n^2/2)+(2.2-((2.2-0.7)*(n/u)))*2*r3n); c3n=(m_dot)/(density*3600*A3n*10^-6); D3n=r3n+(2.2-((2.2-0.7)*(n/u))); L3n=10.575/u; Re3n = density*c3n*D3n*10^-3/(1.74*10^-5); lambda3n = 0.0055*(1+(2*10^4*(epsilon/D3n)+10^6/Re3n)^(1/3)); AP3n=lambda3n*0.5*density*c3n^2*L3n*10^-5/D3n; AP34=AP34+AP3n; APtotal_3=APtotal_3+AP3n; end Prel_4=Prel_3-AP34; %section 3.1 to 3.2 %discretize in w sections APtotal_4=APtotal_3; AP45=0; w=1000000; for t=1:w r4t=2-((2-0.56)*(t/w)); A4t=(pi*r4t^2/2)+(0.7*2*r4t); c4t=(m_dot-((t*m_dot)/(w+1)))/(density*3600*A4t*10^-6); D4t=r4t+0.7; L4t=12.98*2*pi*13/61/w; Re4t = density*c4t*D4t*10^-3/(1.74*10^-5); lambda4t = 0.3164/(0.95*Re4t^0.25)+0.0075*sqrt(D4t/2/(r4t)); AP4t=lambda4t*0.5*density*c4t^2*L4t*10^-5/D4t; AP45=AP45+AP4t; APtotal_4=APtotal_4+AP4t; end Prel_5=Prel_4-AP45; %total pressure loss from inlet to volute outlet (1 to 4) Total_pressure_loss_inlet=APtotal_4; C.2. VOLUTE OUTLET VELOCITY ANGLE %to calculate the angle between the velocity and the normal vector to the volute outlet surface epsilon = 13/61; %admission ratio r_1 = 12.775; %radius of the volute outlet surface r_2 = 2; %radius of the volute area at the volute inlet b = 0.7; %width of the volute area 52 r_0 = 12.975+2; %distace from turbine centre to volute inlet centre A_1 = epsilon * pi * 2 * r_1 * b; %volute outlet area A_0 = 0.5 * pi * r_2^2 + 2 * r_2 * b; %volute inlet area alpha = atan((r_0/r_1)*(A_1/A_0))*180/pi; %angle between velocity and the normal vector at the volute outlet C.3. VOLUTE OUTLET PRESSURE LOSS clear clc %inputs: m_dot, Ra (Roughness) m_dot=2.859; %input Ra=3.172; %input epsilon=Ra*10^-3*11.03; R=286.9; density=0.8; %(kg/m^3) %calculations %section 1 to 2 r12_1=4; %(mm) r12_2=1.135; %(mm) A12=pi*r12_1^2-pi*r12_2^2; %(mm^2) c12=m_dot/(density*3600*A12*10^-6); %(m/s) Dh12=2*r12_1-2*r12_2; %(mm) r12=Dh12/2; L12=4.7; %(mm) %moody Re12 = (density*c12*Dh12*10^-3)/(1.74*10^-5); lambda12 = 0.0055*(1+(2*10^4*(epsilon/Dh12)+10^6/Re12)^(1/3)); AP12=lambda12*0.5*density*c12^2*L12*10^-5/Dh12; APtotal_1=AP12; %bar %section 2 to 3 r23=4; %(mm) A23=pi*r23^2; %(mm^2) c23=m_dot/(density*3600*A23*10^-6); %(m/s) D23=2*r23; %(mm) L23=4.3; %(mm) %moody Re23 = density*c23*D23*10^-3/(1.74*10^-5); lambda23 = 0.0055*(1+(2*10^4*(epsilon/D23)+10^6/Re23)^(1/3)); AP23=lambda23*0.5*density*c23^2*L23*10^-5/D23; APtotal_2=APtotal_1+AP23; %bar %section 3 to 4 %discretize in t sections APtotal_3=APtotal_2; t=1000000; for i=1:t r3i=r23+(2*(i/t)); 53 A3i=pi*r3i^2; c3i=m_dot/(density*3600*A3i*10^-6); D3i=r3i*2; L3i=20.8/t; Re3i = density*c3i*D3i*10^-3/(1.74*10^-5); lambda3i = 0.0055*(1+(2*10^4*(epsilon/D3i)+10^6/Re3i)^(1/3)); AP3i=lambda3i*0.5*density*c3i^2*L3i*10^-5/D3i; APtotal_3=APtotal_3+AP3i; %bar end %section 4 to 5 r45=6; %(mm) A45=pi*r45^2; %(mm^2) c45=m_dot/(density*3600*A45*10^-6); %(m/s) D45=2*r45; %(mm) L45=126.2; %(mm) %moody Re45 = density*c45*D45*10^-3/(1.74*10^-5); lambda45 = 0.0055*(1+(2*10^4*(epsilon/D45)+10^6/Re45)^(1/3)); AP45=lambda45*0.5*density*c45^2*L45*10^-5/D45; APtotal_4=APtotal_3+AP45; %bar %section 5 to 6 %discretize in t pieces APtotal_5=APtotal_4; for n=1:t r5n=6-(n/t); A5n=pi*r5n^2; c5n=(m_dot)/(density*3600*A5n*10^-6); D5n=r5n*2; L5n=11.4/t; Re5n = density*c5n*D5n*10^-3/(1.74*10^-5); lambda5n = 0.0055*(1+(2*10^4*(epsilon/D5n)+10^6/Re5n)^(1/3)); AP5n=lambda5n*0.5*density*c5n^2*L5n*10^-5/D5n; APtotal_5=APtotal_5+AP5n; %bar end %section 6 to 7 r67=5; %(mm) A67=pi*r67^2; %(mm^2) c67=m_dot/(density*3600*A67*10^-6); %(m/s) D67=2*r67; %(mm) L67=188.6; %(mm) %moody Re67 = density*c67*D67*10^-3/(1.74*10^-5); lambda67 = 0.0055*(1+(2*10^4*(epsilon/D67)+10^6/Re67)^(1/3)); AP67=lambda67*0.5*density*c67^2*L67*10^-5/D67; APtotal_6=APtotal_5+AP67; %bar %Sudden expansion loss section 9 A_1=pi*(0.004^2-0.0028^2)*0.33; A_2=pi*0.004^2; c_1=m_dot/(3600*A_1*density); 54 sudden_loss_9=density*(A_1/A_2)*(1-(A_1/A_2))*c_1^2/100; %mbar %Sudden expansion loss section 12 A12=A67; c_12=m_dot/(density*3600*A12*10^-6); sudden_loss_12=(1/2)*density*c_12^2/100; Outlet_loss=sudden_loss_9+(APtotal_6*10^3)+sudden_loss_12; %mbar C.4. TURBINE POWER AND EFFICIENCY %Implement Refprop addpath('C:\Program Files (x86)\REFPROP\refprop-extras\MATLAB'); addpath(genpath('C:\Program Files (x86)\REFPROP\refpropextras\MATLAB\extras-matlab-functions')); m_dot_exp=2.859; %mass flow (kg/h) P_in=2.7347; %pressure inlet (bar) P_out=0.965+0.267465; %pressure outlet (bar) T_in=220.31+273.15; %temperature inlet (K) T_out=193.1+273.15; %temperature outlet density=refpropm('D','T',(T_in+T_out)/2,'P',(P_in+P_out)/2*100,'nitroge n','oxygen','argon',[0.755 0.231 0.014]); Cp=refpropm('C','T',(T_in+T_out)/2,'P',(P_in+P_out)/2*100,'nitrogen','o xygen','argon',[0.755 0.231 0.014]); R=287.13; %constant gases density_in=P_in*10^5/(R*T_in); %density inlet (kg/m^3) density_out=P_out*10^5/(R*T_out); %density outlet (kg/m^3) K=refpropm('K','T',(T_in+T_out)/2,'P',(P_in+P_out)/2*100,'nitrogen','ox ygen','argon',[0.755 0.231 0.014]); Cs=refpropm('A','T',(T_in),'P',(P_in)*100,'nitrogen','oxygen','argon',[ 0.755 0.231 0.014]); %geometry data c6=Cs; %(m/s) r6=7.65*10^-3; %(m) alpha6=22; %stator outlet blade angle (º) betha7=36.59; %rotor inlet blade angle (º) r7=7.5*10^-3; %(m) betha8=138.19; %rotor TE blade angle (º) r8=6.5*10^-3; %(m) d8=0.695; %distance between blades (measured with Catia) h8=0.586; %blade height A8=d8*h8*13; %area between the blades at the TE w=168200*pi/30; %(rad/s) %Mach = 1 in the critical section (stator outlet) %calculations Cu6=c6*cos(alpha6*pi/180); %(m/s) 55 Cu7=(Cu6*r6)/r7; %(m/s) u7=w*r7; %(m/s) Cm8=m_dot_exp/(3600*density_out*A8*10^-6); %perpendicular velocity to A8 W8=Cm8/sin((180-betha8)*pi/180); %relative velocity (m/s) Wu8=W8*cos((180-betha8)*pi/180); %(m/s) u8=w*r8; %(m/s) Cu8=Wu8-u8; %(m/s) P_exp=m_dot_exp*(Cu7*u7-Cu8*u8)/3600; %turbine power P_corrected=39.9; %(W) (adding the bearing heat losses) T_out_corrected=T_in-(P_corrected/((m_dot_exp/3600)*Cp)); P_temperatures=m_dot_exp/3600*Cp*(T_in-T_out); %turbine power experimental c8=sqrt(Cu8^2+Cm8^2); Static_pressure_out = P_out -(0.5*density_out*c8^2)*10^-5; Efficiency_corrected = (T_in-T_out_corrected)/(T_in*(1- ((Static_pressure_out/P_in)^((K-1)/K)))); Efficiency_non_corrected = (T_in-T_out)/(T_in*(1- ((Static_pressure_out/P_in)^((K-1)/K)))); 56 BIBLIOGRAPHY [1] P. H. Wagner, Integrated Design, Optimization, and Experimental Testing of the World's Smallest Steam-Driven Recirculation Fan for Solid Oxide Fuel Cell Systems, Neuchâtel: EPFL, 2018. [2] R. I. Lewis, Turbomachinery Performance Analysis, London, 1996. [3] L. F. Moody, An Appropriate Formula for Pipe Friction Factors, 1947. [4] P. M. a. S. N. Gupta, Momentum Transfer in Curved Pipes, 1979. [5] C. G. H. W. Thomas Adams, A Simple Algorithm to Relate Measured Surface Roughness to Equivalent Sand-grain Roughness, USA: Avestia Publishing, 2012. [6] H. Chanson, The Hydraulics of Open Channel Flow: An Introduction, 1999. [7] I. ANSYS, ANSYS ICEM CFD 11.0, USA, 2007. [8] I. ANSYS, Ansys CFX Tutorials, USA, 2010. [9] S. L. Dixon, Fluid Mechanics and Thermodynamics of Turbomachinery, UK: Elsevier Inc., 2014. [10] S. A. Korpela, Principles of Turbomachinery, USA: John Wiley and Sons, Inc., 2011.