NUMERICAL SIMULATION OF TURBULENT FLOWS ...

67 downloads 0 Views 889KB Size Report
Jan 29, 1999 - von Karman Institute for Fluid Dynamics, .... visors - Professor Ulf H ll at Chalmers University of Technology and ...... The convective (Euler).
ISBN 91-7197-747-3 ISSN 0346-718X

Doctoral Thesis for the degree of Doctor of Philosophy

NUMERICAL SIMULATION OF TURBULENT FLOWS FOR TURBINE BLADE HEAT TRANSFER APPLICATIONS By

Jonas Larsson

Department of Thermo and Fluid Dynamics Chalmers University of Technology S41296 Gteborg, Sweden Gteborg, 1998

CHALMERS

NUMERICAL SIMULATION OF TURBULENT FLOWS FOR TURBINE BLADE HEAT TRANSFER APPLICATIONS By

Jonas Larsson Thesis submitted for the degree of Doctor of Philosophy at the School of Mechanical and Vehicular Engineering, Chalmers University of Technology, G teborg, Sweden. To be defended publicly on Friday 29 January, 1999, 10.15 a.m. in HA 2, H rsalsvgen 4, G teborg. Faculty opponent: Professor Tony Arts, von Karman Institute for Fluid Dynamics, Rhode-Saint-Genese, Belgium. Department of Thermo and Fluid Dynamics Chalmers University of Technology S41296 G teborg, Sweden G teborg, 1998

NUMERICAL SIMULATION OF TURBULENT FLOWS FOR TURBINE BLADE HEAT TRANSFER APPLICATIONS Jonas Larsson [email protected] http://www.cfd-online.com/Users/jola/ Department of Thermo and Fluid Dynamics Chalmers University of Technology S-412 96 Gothenburg SWEDEN

Abstract Turbine blade heat transfer is an important engineering problem characterized by complex ow elds and high turbulence levels. This thesis is focused on using a full Navier-Stokes solver with two-equation eddy-viscosity models to predict external heat-transfer in single-stage, linear, two-dimensional uncooled turbine cascades. The main application is supersonic space turbines, but most results presented are for subsonic and transonic cascades, for which there are measurements to compare with. The turbulence models evaluated include the algebraic Baldwin-Lomax model, three low Reynolds k ;  models (Chien, Launder-Sharma and Nagano-Tagawa) and two k ; ! models (Wilcox standard and transition). A new non-linear k ; ! model has also been developed, which improves the predictions in the leading edge region signicantly. Results are generally in good agreement with measurements. The main problem, which remains unsolved, is transition prediction. The numerical method is a block-structured explicit Runge-Kutta nite-volume scheme. A detailed description of this method and the governing equations is given in the thesis. The numerical quality of the simulations has been thoroughly investigated in order to ensure that the results are representative of the turbulence models and not the numerics. Grid and scheme independence has been veried, and a few general guidelines about the numerics are summarized.

Keywords: Turbine, external heat transfer, turbulence model, Navier-Stokes, numerical method, two-equation, k ; , k ; !, non-linear, low-Reynolds, explicit Runge-Kutta, nite volume, two-dimensional, stagnation ow, transition, supersonic, impulse blade i

Contents

Abstract Acknowledgments Preface Nomenclature 1 Introduction

1.1 Heat Transfer in Turbomachinery . . . . . . . . . . . 1.2 Turbine Blade Heat Transfer . . . . . . . . . . . . . . 1.2.1 Stagnation Region . . . . . . . . . . . . . . . 1.2.2 Laminar Regions . . . . . . . . . . . . . . . . 1.2.3 Transitional Region . . . . . . . . . . . . . . . 1.2.4 Turbulent Regions . . . . . . . . . . . . . . . 1.2.5 Separated Regions and Re-attachment Points 1.3 Heat Transfer Prediction . . . . . . . . . . . . . . . . 1.4 Supersonic Turbines . . . . . . . . . . . . . . . . . .

2 Governing Equations

2.1 Instantaneous Equations . . . . . . . 2.2 Favre Averaged Equations . . . . . . 2.2.1 Averaging . . . . . . . . . . . 2.2.2 Open Turbulent Equations . . 2.2.3 Approximations . . . . . . . . 2.2.4 Turbulence Energy Equation . 2.3 Turbulence Models . . . . . . . . . . 2.3.1 Baldwin-Lomax Model . . . . 2.3.2 k ;  Models . . . . . . . . . 2.3.3 k ; ! Models . . . . . . . . . 2.3.4 New Non-Linear k ; ! Model

3 Numerical Methods

3.1 Conservative Formulation . . . . . 3.2 Space Discretization . . . . . . . . 3.3 Flux Term . . . . . . . . . . . . . . 3.3.1 Convective Flux . . . . . . . 3.3.2 Diusive Flux . . . . . . . . 3.4 Source Terms . . . . . . . . . . . . 3.5 Time Integration . . . . . . . . . . 3.6 Turbulence Modications . . . . . . 3.6.1 TVD Limiter . . . . . . . . 3.7 Discussion and General Guidelines ii

. . . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

i iv v vi 1

1 1 3 4 5 6 6 6 7

11

11 13 13 13 14 18 19 21 23 25 27

29

29 30 32 33 35 38 38 40 42 43

4 Summary of Papers 4.1 4.2 4.3 4.4

Paper 1 Paper 2 Paper 3 Paper 4

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

5 Conclusions References Appendix A - Paper 1 Appendix B - Paper 2 Appendix C - Paper 3 Appendix D - Paper 4

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

44 44 45 47 47

48 49

iii

Acknowledgments This work was carried out at the Department of Thermo and Fluid Dynamics at Chalmers University of Technology. The project was funded by the Swedish National Space Board and Volvo Aero Corporation. It has been a part of a cooperation between Chalmers University of Technology (CTH), Volvo Aero Corporation (VAC), Ecole Central de Lyon (ECL) and SEP, division de snecma to develop new methods to predict external heat transfer in supersonic turbines. I would like to express my sincere gratitude and appreciation to my two supervisors - Professor Ulf Hll at Chalmers University of Technology and Professor Lars-Erik Eriksson at Volvo Aero Corporation. Ulf Hll's support and friendly advice has been very valuable. He always believed in me and he always helped me whenever I needed guidance in the world of turbomachinery. Lars-Erik Eriksson played an important role in the beginning of this project. His numerical expertise gave this work a head start and throughout the years he has been a constant source of good advice and friendly discussions. Thank you to you both! I would also like to take this opportunity to thank my mother, Stina Larsson. She introduced me to the world of mathematics and without her belief in me I would not be writing these words today. Thank you!

iv

Preface This thesis is based on the results presented in the following papers:

 Jonas Larsson, Lars-Erik Eriksson and Ulf Hll, External Heat Transfer

Predictions in Supersonic Turbines Using the Reynolds Averaged NavierStokes Equations, Proceedings from the 12th ISABE Conference, Melbourne, 1995, Copies distributed by AIAA.  Jonas Larsson, Turbine Blade Heat Transfer Calculations Using Two-Equation Turbulence Models, IMechE Journal A, Proc Instn Mech Engrs Vol 211 Part A, Special Turbomachinery Issue, pp 253-262, 1997, also in Proceedings from the 2:nd European Conference on Turbomachinery Fluid Dynamics and Thermodynamics, Antwerpen, 1997, pp 411-420.  Jonas Larsson, Two-Equation Turbulence Models for Turbine Blade Heat Transfer Simulations, Proceedings from the 13th ISABE Conference, Chattanoga, 1997, Vol 2, pp 1214-1222, Paper ISABE 97-7163, Copies distributed by AIAA.  Jonas Larsson with others, Simulations and Measurements on Impulse Blades for Heat Transfer Prediction in Supersonic Turbine Applications, ASME Paper 98-GT-154, Presented at the ASME Turbo Expo, Stockholm, 1998.

v

Nomenclature Cp Cv e e0 h k p Pr qi R Sij t T u ui   ij     ij !

Specic heat at constant pressure J=kgK] Specic heat at constant volume J=kgK] Internal energy J=kg] Total energy J=kg] Enthalpy J=kg] Turbulent kinetic energy J=kg] Static pressure Pa] Prandtl number ] Heat ux W=m2 ] Specic gas constant J=kgK] Trace-less viscous strain-rate tensor s;1 ] Time s] Static temperature K] Friction velocity m=s] Velocity m=s] Specic heat ratio = Cp=Cv ] Dissipation of k J=kgs] Kronecker's delta function ] Thermal conductivity W=mK] Dynamic viscosity Ns=m2 ] Kinematic viscosity m2 =s] Density kg=m3 ] Shear stress tensor N=m2 ] Specic dissipation s;1 ]

Subscript:

t

Turbulent property

Superscript:

conv diff lam tot turb

Convective part Di usive part Laminar part Laminar + turbulent part Turbulent part Reynolds uctuating part Favre uctuating part Density weighted (Favre) average Normal average (time or space)

0

00

e 



vi

1 Introduction The goal of this work is to contribute to the development of new numerical methods for the prediction of external heat-transfer in supersonic turbines. Most results presented are related to subsonic and transonic turbines, for which there are measurements available to compare with. However, the methods developed are suitable for supersonic cases. The next sections give an overview of turbine blade heat transfer and how it can be predicted. Then follows an introduction to the application, supersonic space turbines, and a detailed description of the physical models and the numerical techniques used. Finally the results from the papers are summarized and a few conclusions and thoughts about the future are given.

1.1 Heat Transfer in Turbomachinery Heat transfer in turbomachinery is an important engineering problem. By designing for a higher turbine inlet temperature it is possible to increase the eciency and performance of a gasturbine. A higher turbine inlet temperature usually gives a more ecient cycle and improves the power output and the fuel consumption of the engine. It can also allow for a lighter design. This is true for most gasturbines and especially for modern high-speed aircraft jet-engines. Consequently, there has been a general trend in the development of new jet-engines towards higher and higher turbine inlet temperatures. This has lead to a constant search for new materials, improved cooling techniques and a better understanding of the heat transfer processes involved. A modern military jet engine today can operate with a turbine inlet temperature above 1800K. The materials used in the blades can withstand metal temperatures up to a maximum of about 1300K. Internal cooling of the blades and external lm-cooling is therefore necessary.

1.2 Turbine Blade Heat Transfer The problem of predicting external heat transfer on turbine blades is essentially the same problem as to predict the boundary layer development on the blade surfaces. The thermal boundary layer is what determines the heat transfer between the blade and the gas. This is always good to keep in mind when analyzing and talking about turbine blade heat transfer. What happens outside, in the global ow eld, does not inuence the heat transfer unless it also inuences the 1

boundary layer development. However, as will be exemplied below, many global eects do have an important eect on the boundary layers and thereby also on the heat transfer. The external heat-transfer in a turbine is inuenced by many parameters, and detailed information about the ow is needed in order to obtain good results. The ow is characterized by regions with large continuous pressure gradients, causing strong acceleration and deceleration. Under certain operating conditions there might be separated regions. Some turbines operate under transonic conditions, and shocks, which interact with the boundary layers, can inuence the heattransfer. The mass-ow into the turbine comes from a combustor, and is strongly turbulent. This causes an early boundary-layer transition and an increased heattransfer in both the laminar 14] and the turbulent parts 49] 31]. Another important eect of the combustor is that it gives rise to high-temperature regions in the incoming ow and these hot-streaks can have an important inuence on the heat ux in the rst stages of a turbine. Multi-staging gives rise to wakes and shocks from the preceding stages, which sweep over the blades. As shown by many authors 5] 13] 15] 32] 41] this can increase the heat-transfer signicantly, although a high inlet turbulence level tends to decrease the relative importance of this eect. In a real turbine three-dimensional eects can also be important for the heat-transfer 8]. If the blades are lm-cooled this leads to even more complications, but this is not within the scope of the present work and will not be considered in this thesis. Summarizing, it is possible to isolate a number of dierent eects that are important for turbine blade heat transfer:

      

Unsteady eects due to multi-staging Hot streaks Pressure gradients Curvature Compressibility / shocks Free-Stream turbulence Secondary ows and 3D eects

Looking at this from another point of view the ow-path can be separated into the following dierent regions: 2

    

Stagnation region Laminar regions Transitional region Turbulent regions Separated regions and re-attachment points

1.2.1 Stagnation Region The stagnation region is a very critical region from the point of view of heat transfer - this is the region where the suction and pressure side boundary layers start to develop and eects occurring here are convected downstream and can inuence the heat transfer on large parts of the blade. The heat transfer in the stagnation region can become very high and this is often the region on the blade where the maximum heat ux occurs. The leading edge heat transfer is very dependent on the inlet turbulence. This is because the boundary layers here are thin and eddies from the inow can easily penetrate down to the blade surface and thereby inuence the heat transfer. Typically a high free-stream turbulence level (more than 10%) can increase the heat transfer in the stagnation point by 50% compared to the laminar level 40]. A simple physical explanation to this eect is that the turbulent eddies approaching a stagnation region are stretched by the diverging streamlines. This increases the vorticity and this increased vorticity is believed to be the cause of the amplied leading edge heat transfer. Note that the boundary layers can still remain in principal laminar. Not only the turbulence level but also the length scale is important for the stagnation region heat transfer. A smaller length scale, as shown by several authors 3] 19], usually gives a larger increase in heat transfer. This eect, however, is smaller than the inuence of the turbulence level, although it becomes more and more important at high turbulence levels. The simple explanation is that large scale eddies are seen by the blade as mean-ow variations since these eddies are too large to be stretched by the diverging streamlines in the stagnation region. Hence, large scale eddies are not able to inuence the leading edge heat transfer very much. Smaller eddies, which are stretched by the diverging streamlines, can, however, directly inuence the stagnation region heat transfer. Some authors 19] have also argued that for very small length scales viscous dissipation would destroy the turbulent eddies before they can interact with the boundary 3

layers, and that there must exist and optimum length scale for which the heat transfer reaches a maximum. However, this has not been veried experimentally - the heat transfer continues to increase with a decreased length scale even for length-scale over leading-edge-radius ratios as low as 0.05 19]. With this simple explanation of turbulent eddies stretched by the diverging streamlines it is easy to understand that the structure of the turbulence must also be important. Since the diverging streamlines mainly stretches the vorticity in the pitch-wise direction (normal to the stream-wise and span-wise directions) it can be assumed that pitch-wise vorticity has a larger inuence. This has also been show by authors like 3].

1.2.2 Laminar Regions The heat ux in a laminar boundary layer is lower than the heat ux in a corresponding turbulent boundary layer. The most critical aspect of the laminar parts of the boundary layers is therefore often when they become turbulent. For a laminar, incompressible boundary layer on a at plate the heat ux scales as Re1x=2 and Pr1=3. In laminar boundary layers there is a direct relationship between the velocity prole / skin-friction and the temperature prole / heat transfer. This means that in the simple laminar case when the coupling between temperature and velocity through density variations can be neglected the problem of predicting heat transfer is essentially the same problem as to predict skinfriction. For Pr = 1 the velocity boundary layer and the thermal boundary layer are exactly the same. On the pressure side the concave curvature and the centrifugal pressure gradient can sometimes give rise to G rtler vortices. These are stable laminar pairs of counter-rotating vortices that are aligned with the ow. The instabilities creating the G rtler vortices are similar to the instabilities creating Taylor vortices in Couette ows. G rtler vortices will increase the heat transfer. The occurrence of these vortices can be assessed by calculating the G rtler number, which is a dimensionless number relating viscosity, boundary layer thickness, curvature and free-stream velocity. For high turbulence levels the eect of these G rtler vortices can often be neglected. If a full Navier-Stokes solver is used to predict heat transfer G rtler vortices can be predicted automatically, provided that a full 3 dimensional simulation is made. The heat transfer in the laminar regions of the boundary layers is classically said to be independent of the free-stream turbulence. However, this is not true 4

at high turbulence levels. At high turbulence levels eddies from the free-stream interact with the laminar boundary layers and do increase the heat transfer also in these laminar regions. This eect is larger close to the leading edge of the blades, especially on the suction side. Further downstream on the suction side the convex curvature and the thicker boundary layers reduce the eect of the free-stream turbulence. On the suction side the most important eect of the freestream turbulence level on the laminar region is that it limits its extent. With a high inlet turbulence level diusion of turbulent energy into the laminar boundary layer will trigger an earlier by-pass transition and thus make the laminar region smaller. The inuence on the laminar heat transfer from a high free-stream turbulence level is usually larger on the pressure side of the blades than on the suction-side and the inuence also extends further downstream due to the concave curvature of the pressure side.

1.2.3 Transitional Region When talking about the transitional region in turbine blade heat transfer applications one usually refers to the transition occurring on the suction side. There is of course also a transition occurring on the pressure side, but this usually has a less dramatic eect on the heat transfer. The transitional process on the pressure side boundary layer is most often less obvious and occurs over a larger distance of the blade. Laminar instabilities like G rtler vortices etc. and the laminarizing acceleration close to the trailing edge tends to hide the transition on the pressure side. Transition is caused by a combination of dierent eects - boundary layer growth which eventually creates instabilities, diusion of turbulent energy from the freestream into the boundary layer, blade geometry triggering transition, surface roughness disturbing the laminar boundary layer, pressure gradients and shockboundary layer interaction and sometimes even laminar separation bubbles. Predicting transition in general is also very dicult. Suction-side transition on turbine blades usually leads to a signicant increase in the heat transfer - typically the heat ux can go up by 100% as the boundary layer becomes turbulent. The extent of the suction-side transitional region is usually quite small compared the extent of the other regions and the sharp increase in heat transfer over a short distance can therefore, for certain applications, create dangerous temperature gradients in the blades. For supersonic stages suction-side transition is often caused by the bow-shock impinging on the suction surface of the neighboring blade. 5

1.2.4 Turbulent Regions The turbulent region on the suction side is usually not very much inuenced by the free-stream turbulence. The convex curvature of the blade can sometimes be important though. However, this is for blades with very strong curvature, like, for example, the SFTB case used in this work, which on the suction side shoulder has a boundary layer thickness over radius of curvature ratio of about 0:10. A convex curvature has a stabilizing eect on the turbulence and reduces u0v0. This leads to a decreased heat transfer. For example, 50] and 34] showed that the heat transfer can be reduced by as much as 15% to 20% on a convex surface compared to the same at plate case. This eect also extends further downstream although the blade attens out. For example, 46] showed that in the recovery region following a curvature the Stanton number remains 15% to 20% below at plate values.

1.2.5 Separated Regions and Re-attachment Points Separated regions naturally aect the heat transfer signicantly. Separation can be caused by too strong adverse pressure gradients present under certain operating conditions. This most often occurs close to the leading edge on the pressure side or close to the trailing edge on the suction side. Another mechanism that sometimes can cause separation is shock - boundary layer interaction. This can for supersonic cascades occur either on the front part of the suction side, where the bow-shock interacts with the boundary layer, or on the rear part of the suction side, where the trailing-edge sh-tail shock/expansion-wave pattern interacts with the boundary layer. Separation on the pressure side does not normally occur due to shock-boundary layer interaction.

1.3 Heat Transfer Prediction Turbine heat-transfer can be calculated in many dierent ways. One way is to use a boundary-layer solver, with a prescribed global ow-eld obtained with a separate Euler or potential method. This is a fast and sometimes accurate method, which has been used by many researchers 12] 21] 48]. However, it usually breaks down when large separations occur, and it is dicult to give correct initial and boundary conditions. Thus, boundary-layer solvers are not good at predicting what happens at the leading edge and they often have diculties at the trailing edge due to separation. For supersonic impulse turbines these are 6

the regions where the heat-transfer is most critical because of the slender blade edges that are used. To handle eects like varying inlet conditions, high free stream turbulence levels, high Mach numbers, shock - boundary layer interaction, curvature, separation, transition and laminarization it is necessary to use a full Navier-Stokes method. The rapid development of faster computers now makes it possible to eciently use this kind of methods for engineering applications. The Navier-Stokes method needs a turbulence model to account for the small turbulent scales in the ow, which have a strong inuence on the heat-transfer. The main problem then becomes to nd or develop a turbulence model which is suitable for the cases of interest. Other researchers have used algebraic turbulence models 9], with good results for subsonic and transonic cases. However, the bowshock and the very high free-stream turbulence levels present in a supersonic turbine are dicult to take into account with an algebraic model. It is also dicult to predict separated regions correctly. With a transport model these problems can be solved, or at least reduced. Many researchers have also tried to use two-equation turbulence models to predict turbine blade heat transfer 1] 2] 11] 25] 26] 24] 27] 29] 36]. Using two-equation turbulence models in a turbine cascade simulation introduces several new problems. Many of these models behave badly in the stagnation region around the leading edge, and in order to obtain good results it is necessary to do something about this. Two-equation turbulence models can also potentially predict by-pass transition, but as will be shown this is very dicult in a turbine simulation. Most calculations done in this work are for stationary and two-dimensional singlestage turbine cascades. This greatly simplies the numerical analysis, and makes it possible to eciently study the inuence of other important factors like high free-stream turbulence levels, high Mach numbers and shock - boundary-layer interaction, blade curvature, transition and separation. In the future faster computers will make it possible to extend the methods to unsteady and three dimensional cases. The main principles used are the same but the amount of calculations necessary to obtain a solution is several orders of magnitudes greater.

1.4 Supersonic Turbines In rocket-engines using the gas-generator cycle, where turbines are used to power the fuel-pumps, the designs must be extremely light and it is usually not possible to cool the blades. This leads to problems with high temperature loads, especially at the thin leading and trailing edges of the blades. A method to predict the heatloads in the design phase would be very valuable. These space turbines are highly 7

Figure 1: Mach Number Contours, Supersonic Cascade specialized. Usually the blades are small and some stages are fully supersonic. The engine is started with a pyrotechnic powder starter, which gives a very quick transient. Further details about the ow in these turbines will be given below. Reference 52] contains an overview of a typical turbine used in the Vulcain rocket engine - the main engine in Ariane 5. The ow in these supersonic turbine is extremely complex, and it is dicult to predict the heat-loads on the blades. An example of the ow eld in a typical rst stage rotor can be seen in gure 1. The inlet Mach number for this cascade is 1.4, and there is a detached leading edge shock. This shock hits the suction side of the lower blade and interacts with the boundary layer. At the trailing edge a typical sh-tail shock/expansion-wave pattern is generated. This is reected in the suction side of the lower blade, where it in some cases might cause separation close to the trailing edge. The external heat-transfer is inuenced by many parameters, and detailed information about the ow is needed in order to obtain good results. The important characteristics of the ow can be summarized as follows:

8

 The blades are usually of impulse type, with a large deection angle (> 100), and there are large continuous pressure gradients in the cascade.     

Under certain operating conditions there might be separated regions in some parts. The strong acceleration in other parts has a laminarizing eect on the boundary-layers there. The mass-ow entering the turbine comes from a combustor, and is strongly turbulent. This gives rise to an early boundary-layer transition, and an increased heat-transfer in both the laminar 14] and the turbulent parts 49] 31] of the boundary-layers. The relative inlet velocity can be supersonic, and this gives rise to a bowshock in front of the blade. This shock interacts with the boundary layer on the suction side of the neighboring blade, where it might cause separation or transition. During startup and shutdown there are strong transients and the heattransfer varies rapidly over a short time. The aspect ratio of the blades is small (length=chord  2), and three dimensional eects are important 8]. Multi-staging gives rise to shocks and wakes from the preceding stages, which sweep across the blades. As shown by many authors 5] 13] 15] 32] 41] this can increase the heat-transfer signicantly, although a high inlet turbulence level tends to decrease the relative importance of this eect.

Clearly it is a very formidable task to predict heat-transfer in this kind of turbines. As can be seen from above several features of the global ow-eld (free-stream turbulence, separations, shocks, wakes, transients, 3D eects, ...) are important for the heat-transfer at the blade surface. Below follow an example of a prediction. However, note that no experimental data exists for this case and it is not possible to say if the predictions are good or not. Figure 2 shows heat transfer predictions side by side with mach contours for the typical supersonic SFTB cascade. The predictions in this case are with the Chien k ;  model. On the pressure side the heat transfer behaves fairly classically, with an early boundary layer transition and then a developing turbulent boundary layer further downstream on the blade leading to a declining heat-ux. More interesting things happen on the suction side though. First the strong acceleration keeps the boundary layer in a laminarizing state until the location where the bow-shock of the neighboring blade hits the suction side. The shock-boundary layer interaction triggers an immediate transition and a sharp increase in the heat transfer. Following after the shock there is a very strong acceleration and 9

Trailing Edge

Suction Side

Pressure Side Axial Distance

Leading Edge Q max

0

Heat Transfer

Figure 2: Heat Transfer Prediction, Supersonic Cascade the boundary layer nearly laminarizes again, but then close to the shoulder of the blade another shock-pattern interacts with the boundary layer and again causes transition. Then further downstream the sh-tail shock/expansion-wave pattern interacts with the boundary layer and causes an increased heat transfer by reducing the size of the boundary layer.

10

2 Governing Equations This section gives a detailed description of the equations used in this work. Included is also an overview of the underlying assumptions and their validity.

2.1 Instantaneous Equations The instantaneous continuity equation (1), momentum equation (2) and energy equation (3) for a compressible uid can be written as:

@ + @ u ] = 0 @t @xj j

(1)

@ (u ) + @ u u + p ; ] = 0 @t i @xj i j ij ji

(2)

@ (e ) + @ u e + u p + q ; u ] = 0 j i ij @t 0 @xj j 0 j

(3)

For a Newtonian uid, assuming Stokes Law for mono-atomic gases, the viscous stress is given by (4):

ij = 2 Sij

(4)

Where the trace-less viscous strain-rate is dened by (5): 

!

@ui + @uj ; 1 @uk  Sij = 21 @x 3 @xk ij j @xi 

def

(5)

The heat-ux, qj , is given by Fourier's law (6):

@T  ;C @T qj = ; @x p Pr @x j

j

11

(6)

Where the laminar Prandtl number Pr is dened by (7):

Pr def = C p

(7)

To close these equations it is also necessary to specify an equation of state. Assuming a calorically perfect gas the following relations are valid:

def = CCp p = RT e = Cv T Cp ; Cv = R v

(8)

Where , Cp, Cv and R are constant. The total energy e0 is dened by (9):

e0 def = e + uk2uk

(9)

Note that the corresponding expression 15 for Favre averaged turbulent ows contains an extra term related to the turbulent energy. Equations (1)-(9), supplemented with gas data for , Pr, and perhaps R, form a closed set of partial dierential equations, and need only be complemented with boundary conditions.

12

2.2 Favre Averaged Equations It is not possible to solve the instantaneous equations directly for the applications of interest here. At the Reynolds numbers typically present in a turbine these equations have very chaotic turbulent solutions, and it is necessary to model the inuence of the smallest scales. All turbulence models used in this work are based on one-point averaging of the instantaneous equations. The averaging procedure used is described in the next sections.

2.2.1 Averaging Let be any dependent variable. It is convenient to dene two dierent types of averaging of :

 Classical time average (Reynolds average): Z def 1

=

(t)dt T

T

0 def = ;

(10)

 Density weighted time average (Favre average):

e def =  

00 def = ; e

(11)

Note that with the above denitions 0 = 0, but 00 6= 0.

2.2.2 Open Turbulent Equations In order to obtain an averaged form of the governing equations, the instantaneous continuity equation (1), momentum equation (2) and energy equation (3) are time-averaged. Introducing a density weighted time average decomposition (11) of ui and e0, and a standard time average decomposition (10) of  and p gives the following exact open equations: 13

@ + @ uf] = 0 @t @xi i

(12)

@ (f @ hufuf + p + u00u00 ; i = 0 u ) + ij ji i j @t i @xj i j

(13)

@ (ef) + @ hufef + uf p + u00p + u00e00 + q ; u i = 0 j i ij j j 0 @t 0 @xj j 0 j

(14)

The density averaged total energy ef0 is given by (15):

ef0  e + uk2uk + k ff

(15)

Where the turbulent energy, k, is dened by (16): g 00 00 k def = uk2uk

(16)

Equation (12), (13) and (14) are referred to as the Favre averaged Navier-Stokes equations. , f ui and ef0 are the primary solution variables.

2.2.3 Approximations To analyze equation (12), (13) and (14) it is convenient to rewrite the unknown terms in the following way:

ji = fji + ji00

(17)

u00j p + u00j e000 = Cpu00j T + uiu00i u00j + uj 2ui ui

(18)

@T = ;C @ Te ; C @T 00 qj = ;Cp Pr p @xj Pr @xj p Pr @xj

(19)

00 00 00

14

ui ij = f ui fij + u00i ij + ufi ij00

(20)

Where the perfect gas relations (8) and Fourier's law (6) have been used. Note also that uctuations in the molecular viscosity, , have been neglected. Inserting (17)-(20) into (12), (13), and (14) gives:

@ + @ uf] = 0 @t @xi i

(21)

2 3 7 @ (f @ 66f u u =0 i) + iufj + pij + u00i u00j ; fji ; ji00 7 4 @t @xj | {z } |{z}5 (1 )

(2 )

2 @ (ef) + @ 666uf ef + uf p + C u00T + ufu00u00 + u00j u00i u00i @t 0 @xj 4 j 0 j | p {z j } | i {zi j} | {z2 } (3 ) (4 ) (5 ) 3 77 e 00 @ T

@T

00 00 ;Cp Pr @x ; Cp Pr @x ;ufi fij ; u| i{z ij} ; uf i ij 7 75 = 0 | {z } j | {z j}   (7 )

(6 )

(22)

(23)

(8 )

The terms marked with (1) ; (8 ) are unknown, and have to be modeled in some way. Term (1 ) and (4) can be modeled using an eddy-viscosity assumption for the Reynolds stresses, ijturb :

ijturb def = ;u00i u00j  2 tg Sij ; 23 kij

(24)

Where t is a turbulent viscosity, which is estimated with a turbulence model. The last term is included in order to ensure that the trace of the Reynolds stress tensor is equal to ;2k, as it should be. Term (2 ) and (8) can be neglected if: 15



j fij j >> ij00

(25)

This is true for virtually all ows, and has been used in all simulations here. Term (3), corresponding to turbulent transport of heat, can be modeled using a gradient approximation for the turbulent heat-ux:

t @ Te qjturb def = Cpu00j T  ;Cp Pr t @xj

(26)

Where Prt is a turbulent Prandtl number. In the present work a constant Prt equal to 0:9 has been used everywhere. Other researchers have reported no signicant improvement when using an algebraic expression for the variation of Prt 9]. Term (5) and (7 ), corresponding to turbulent transport and molecular diusion of turbulent energy, can be neglected if the turbulent energy is small compared to the enthalpy: k > @ 2 T 00 2 2

@xj

@xj

(29)

This is true for virtually all ows, and has been used in all simulations here. 16

To summarize, the governing equations (21)-(23), with assumptions (24), (25), (26), (27) and (29) can be written as in (30)-(39). These equations are valid for a perfect gas. Note also that all uctuations in the molecular viscosity have been neglected.

@ + @ uf] = 0 @t @xi i

(30)

i @ huf f @ (f g tot u ) + u + p ; @t i @xj j i ij ji = 0

(31)

i @ (ef) + @ huf f g g tot tot f f e + u p + q ; u i ij = 0 j @t 0 @xj j 0 j

(32)

g tot def g lam turb g ij = ij + ij

(33)

Where

!  fi @ ufj 2 @ ufk @ u def f g lam ij = ij =

+ ; ij

@xj

(34)

3 @xk

@xi

!  fi @ ufj 2 @ ufk @ u g turb 00 00 + ; ij ; 2 kij ij  ;ui uj  t

(35)

g tot def g lam turb qg j = qj + qj

(36)

@xj

@xi

3 @xk

" #

@ Te = ; lam def qg j = qej  ;Cp Pr @xj

@ p ; 1 Pr @xj 

" #

e turb  Cpu00 T  ;Cp t @ T = ; qg j j

Prt @xj 

t @ p ; 1 Prt @xj 

ff p = ( ; 1)  ef0 ; uk2uk ; k

17

3

!

(37) (38) (39)

If a separate turbulence model is used to calculate t, k and Prt, and gas data is given for , and Pr these equations form a closed set of partial dierential equations, which can be solved numerically.

2.2.4 Turbulence Energy Equation An equation for the turbulence energy, k, dened by (16), can be derived by multiplying the instantaneous momentum equation (2) with u00i and averaging. Rearranging the terms and making use of the instantaneous continuity equation (1) gives the following exact open equation for k: 2 @ (k) + @ 666uf k ; u00 | ji{z i} @t @xj 4 j

00 00 00 + uj 2ui ui | {z } Molecular Diffusion

;

0 u00 p|{z} j

Turbulent Transport PressureDiffusion

3 77 75 =

@ ufi ; @ui ; u00 @p + p0 @ui ;u00i u00j @x (40) ji i @x @x @x j j i i | {z } | {z } | {z } | {z } PressureWork PressureDilatation 00

00

Dissipation

Production

In turbine applications the maximum Mach numbers are well below the hypersonic range. This makes it possible to neglect the terms associated with the pressure - pressure diusion, pressure work and pressure dilatation. This has been done for all simulations presented here. The molecular diusion and turbulent transport can be estimated with a gradient approximation as in (28) and the production term follows directly from the eddy-viscosity assumption (24). The resulting modeled k equation becomes: "

#

@ (k) + @ uf k ; + t  @k = P ;  @t @xj j k @xj

(41)

Where the production, P , is dened by (42) and the dissipation, , by (43).

@ ufi  turb @ ufi P def = ;u00i u00j @x ij @x j

18

j

(42)

00 i  def = 1 ji @u @xj

(43)

Note that this equation is equivalent to the incompressible k equation. Also note that if this equation is used in a low-Re model extra terms and damping functions are often introduced to allow the model to be used also in the viscous sublayer close to the surface. For simplicity the Favre (e) and Reynolds () averaging notation will be dropped in all subsequent sections. Usually the meaning is apparent from the context, but if there is a risk of confusion it will be pointed out.

2.3 Turbulence Models None of the turbulence models used in this work have been specically developed for compressible ow calculations, and the only compressible modication is that the turbulence equations have been written in a conservative formulation. The major eort has been to test dierent models and compare with experimental data for subsonic and transonic turbine cascades. For these cases compressibility eects on the turbulence models are certainly negligible. If this is also the case for the supersonic cascades we are eventually interested in remains to be investigated. Good overviews of compressibility eects on turbulence can be found in references 51] and 30]. Wilcox's book 54] also contains a good chapter about more practical aspects of compressibility eects and turbulence modeling. Dierent types of turbulence models have been used. The simplest model is an algebraic eddy-viscosity model - the Baldwin-Lomax model 6]. This kind of model has been the standard turbulence model in turbomachinery Navier-Stokes calculations for some time. It usually behaves very well numerically, and seldom gives completely unrealistic results. For simple ows the results can even be quite good. However, the algebraic models can not account for any history eects on the turbulence. This makes it dicult to accurately and generally model the inuence of free-stream turbulence. This is a very important eect in heat-transfer calculations, and for that reason a number of two-equation transport models have also been tested. As seen in the following sections, going to a two-equation model introduces several new problems. However, the generality gained by using a transport model warrants this step. Going further to even more advanced dierential Reynolds stress models could potentially improve the predictions even more. However, the stability and convergence problems associated 19

with these models would be a major obstacle for using the methods for engineering purposes. Instead it seems better to concentrate on the well established two-equation models and try to improve them for our particular applications. The use of a non-linear stress-strain relation is one interesting improvement that is presented in this thesis. Algebraic Reynolds stress models might also be of interest. For heat-transfer calculations it is very important to model the eects close to the surface accurately. Thus, the turbulence models used should behave well also in the viscous sublayer. For k ;  models this means that extra low-Reynolds modications are necessary, whereas the k ; ! models can be used down to the surface without extra modications.

20

2.3.1 Baldwin-Lomax Model The Baldwin-Lomax model 6] is a two-layer algebraic model which gives t as a function of the local boundary layer velocity prole. The eddy-viscosity, t , is given by: 8 < t y  ycrossover

t = : inner

t outer y > ycrossover

(44)

Where ycrossover is the smallest distance from the surface where tinner is equal to

touter :

ycrossover = MIN (y) : tinner = touter

(45)

The inner region is given by the Prandtl - Van Driest formula:

t inner = l2 jj

(46)

Where

l = ky 1 ; e

;y +

A+



(47)

q

jj = 2ij ij 

@ui ; @uj ij = 12 @x @xi j

(48) !

(49)

The outer region is given by:

touter =  K CCP FWAKE FKLEB (y) Where 21

(50)



2 FWAKE = MIN yMAX FMAX  CWK yMAX FuDIF MAX

!

(51)

yMAX and FMAX are determined from the maximum of the function:

;y +

F (y) = y jj 1 ; e A+



(52)

FKLEB is the intermittency factor given by: 2 !6 3;1  y C KLEB 5 FKLEB (y) = 41 + 5:5

yMAX

(53)

uDIF is the dierence between maximum and minimum speed in the prole. For boundary layers the minimum is always set to zero.

uDIF = MAX (puiui) ; MIN (puiui) A+ CCP CKLEB CWK k

(54)

K

26 1.6 0.3 0.25 0.4 0.0168 Table 1: Model Constants, Baldwin-Lomax Model Table 1 gives the model constants present in the formulas above. Note that k is a constant, and not the turbulence energy, as in other sections. It should also be pointed out that when using the Baldwin-Lomax model the turbulence energy, k, present in the governing equations, is set to zero.

22

2.3.2 k ;  Models Good reviews of classical Low-Re k ;  models can be found in 39] and 43]. For this work three dierent low-Reynolds k ;  models have been used. These models were selected in order to give a good representation of classical low-Re k ; models. Future work will cover more advanced models. The models presented in this thesis are:

 Chien Model 10] - A very common model in turbomachinery applications. Has nice numerical properties.  Launder-Sharma Model 28] - A model which recently has attracted some attention for its ability to predict by-pass transition.  Nagano-Tagawa Model 37] - A model originally developed for heat-transfer applications.

Several other variants of low-Re models were also tested and was found to give very similar results. The three models chosen here have typical characteristics for this class of models. These models can be written in the general form of equation (55), (56), (57) and (58). C1 , C2 , C, k and  are model constants. The damping functions f , f1 and f2 and the extra source terms D and E are only active close to solid walls and makes it possible to solve k and  down to the viscous sublayer. Table 2 summarizes the constants, damping functions and boundary conditions for all k ;  models used. All of the models have been tested on simple at-plate benchmark problems to assure a correct implementation before using them in cascade simulations. "

#

@ (k) + @ ku ; + t  @k = P ;  ; D j @t @xj k @xj "

#

@ () + @ u ; + t  @ = (C f P ; C f )  + E j 1 1 2 2 @t @xj  @xj k

(55) (56)

2

t = C f k

(57)

@ui P = ijturb @x

(58)

j

23

Chien

Launder-Sharma

Nagano-Tagawa

c

0:09

0:09

0:09

k

1

1

1:4



1:3

1:3

1:3

D

2 yk2

2

E

; 2y2 e;0:5y+

2t

wall

0

0

C1

1:35

1:44

1:45

C2

1:8

1:92

1:9

f

;0:0115y+

1;e

f1

1

f2

Re 2 1 ; 0:22e ( 6 t )

1;

 @ pk 2 @y

0

@y2

0

 @ 2 u 2

2

1 ;

;3:4 e (1+Rt=50)2

;Re2t

1 ; 0:3e

Ret def = k2 , y+ def =

u y , 

1;e

;y +

26

 @ pk 2 @y

2

1 + Re3t =4

1





+ Re 2 1 ; 0:3e ( 6:5t ) 1 ; e ;y6 ;

kwall = 0

Table 2: Model Parameters, Low-Re k ;  Models

24

4:1

2

2.3.3 k ; ! Models The k ; ! models have remained in the shadow of the much more common k ;  models for a long time. However, recently there has been a renewed interest in these models, much thanks to the work by Wilcox 54]. The k ; ! models, unlike the k ;  models, do not need any extra damping functions close to the wall, and the boundary condition for ! can be a simple Dirichlet condition at the surface. This gives the k ; ! models nice numerical properties. Wilcox has also developed techniques to handle rough walls and mass injection by simply modifying the boundary condition for !. This makes the k ; ! models attractive. A problem often encountered with k ; ! models is that they are overly sensitive to the free-stream value of !. Menter 35], who in contrast to Wilcox, mainly uses full Navier-Stokes solvers, has developed a modied k ; ! model, which automatically switches to a k ;  model outside the wall regions. However, the sensitivity to the free-stream ! value can also be used to tune a simulation. For example, the inlet value of ! can be varied until the transition point on the suction side of a turbine blade is captured correctly. In lack of a good transition model this can be a way to estimate the heat-transfer if the transition point is known. For an example of this see the recent work by Chima 11]. Two dierent k ; ! models have been used in this work. The simplest model, herein referred to as the standard k ; ! model 53], contains no extra functions close to walls, and was developed by Wilcox in 1988. The other model, herein referred to as the Wilcox transition model 55], has a number of extra functions included to improve the results in transition calculations. Both of these models can be written in the general form of equation (59), (60), (61) and (62). Model constants and functions for the two models are given in table 3. "

#

@ (k) + @ u k ; ( +  ) @k = P ;  !k t @t @xj j @xj "

#

(59)

@ (!) + @ u ! ; ( +  ) @! =  ! P ; !2 t @t @xj j @xj k

(60)

t =  k !

(61)

@ui P = ijturb @x j

25

(62)

Standard Version

Transition Version



1

0 +(Ret =Rk ) 1+(Ret =Rk )



5 9



9 100



3 40

3 40



1 2

1 2



1 2

1 2

0 +(Ret =Rw )  1  1+( Ret =Rw )  5 +(Ret =R )4 9 18  1+(Ret=R )4 100

5 9

k ,  = ,  = 1 , R = 8, R = 6, R = 2:7 Ret def = ! 0

k ! 0 3 10

Table 3: Model Parameters, Wilcox k ; ! Models The theoretical asymptotic behavior of ! approaching a wall is that !  1=y2 as y ! 0. Usually a nite value of ! is used instead. This can be warranted by the fact that in the viscous sublayer the turbulent viscosity is negligible compared to the molecular viscosity, and it is not absolutely necessary that the turbulence has the correct limiting behavior. Using a nite value of !wall is good from a numerical point of view, and this is also what Wilcox recommends 54]. The surface roughness model gives a nite value for !wall as in equation (63) and (64). This model can and should be used also for smooth surfaces, and all simulations presented here were made with this model.

!wall = u SR 2

Where

8  > < 50=kR+ 2 kR+ < 25 SR = > : 100=kR+ kR+ 25

(63)  kR+ def = u kR

(64)

kR is the average height of the sand-grain roughness. Letting kR+ ! 0 corresponds

to the asymptotic behavior for a smooth surface. However, in practise it is sucient to set kR+ to a value below 5 in order to obtain a hydraulically smooth 26

surface. Tests by the author also veries that using an even lower value does not inuence the turbine heat-transfer results, provided a good grid is used. A too low value might cause numerical problems.

2.3.4 New Non-Linear k ; ! Model The nonlinear k ; ! model is based on an algebraic stress-strain relation derived from Dimensional Analysis, Rapid Distortion Theory (RDT) and Realizability. The formulation used here was rst proposed by Shih et al. 44, 45]. First dimensional analysis is used to obtain a truncated second order expression for the Reynolds stresses. This expression is then simplied using results from Rapid Distortion Theory 42, 33]. Realizability, dened as positiveness of all energy components and no violation of Schwarz' inequality between the uctuating velocities, then gives further constraints on the model constants (73). A detailed derivation of this model can be found in 45]. The nal stress-strain relation is given by equation (70), where the turbulent time scale T is dened by (69) and the trace-less strain and rotation parameters are dened in (71) and (72). The two remaining constants A0 and C0 are tuned using homogeneous shear ow and surface ow in the inertial sublayer. The values used here (74) are the same as those proposed in 44]. In this work the nonlinear stress-strain relation is recast in a k ; ! framework. This makes it possible to solve the equations down to the solid walls. The RDT technique used to derive the stress-strain relation is not strictly valid in the lowReynolds number region. However, in a k ; ! framework, the model performs reasonably well also in the inner parts of the boundary layers, and it has not been found necessary to use any further viscous corrections close to the walls. The same model constants as for the standard linear k ; ! model are used with the nonlinear stress-strain relation. "

#

@ (k) + @ u k ; ( +  ) @k = P ;  !k t @t @xj j @xj "

#

(65)

@ (!) + @ u ! ; ( +  ) @! =  ! P ; !2 t @t @xj j @xj k

(66)

@ui P def = ijturb @x

(67)

j

27

t def = C k T

(68)

= 1! = k

(69)

T

def

ijturb def = ;u00i u00j     ; 23 kij + C k 2 T Sij ; C2k 2 T2 ;Sik kj + ik Skj (70) 

!

@ui + @uj Sij = Sij ; 31 Skk ij Sij def = 21 @x @xi !  j @ui ; @uj ij def = ij ij def = 12 @x j @xi 

def

q

(71)

q

U  def = Sij Sij + ij ij S  def = Sij Sij    q    def  def Sij Sjk Ski  = ij ij W = (S )3 C = A + A1  U  T po s  As = 6 cos 1

C2 =

p1;(3C S T)2 

C0 +6S  T2

p

1 = 31 arccos 6W 

A0 = 6:5

C0 = 1:0

28

(72)

(73)

(74)

3 Numerical Methods All solvers used in this work are based on the G2DFLOW 18] series of codes developed by Lars-Erik Eriksson at Volvo Aero Corporation. A number of new turbulence models have been added to the original code, including low-Re k ;  models and k ; ! models. To be able to solve these models it is also necessary to modify and enhance the numerical methods. For simplicity all examples below will be for the k ;  model. The same basic techniques are also used for the other turbulence models. In short the numerical method is an explicit three-stage Runge-Kutta timemarching method using an optimized local time-step. The convective (Euler) parts are discretized with a third order accurate cell-centered nite volume scheme with upwind-biasing based on the characteristic variables and associated velocities. When solving turbulent quantities (k and ) a more stable rst order upwind scheme, or the third order scheme with a TVD limiter, making it second order accurate, is used. The viscous parts are discretized with a compact second order scheme. A detailed analysis of this scheme can be found in 17]. In addition to this a few extra tricks are necessary to avoid problems with k and  during the early stages of a simulation.

3.1 Conservative Formulation The governing equations can be written in a compact conservative form as follows:

@Q + @Fj = S @t @xj

(75)

Where 2 66 6 Q = 66 64

 ui e0 k 

3 77 77 77 5

2 uj 66 ui uj + pij ; ijtot 6 Fj = 66 (e0 + p)uj + qj(k;) ui ijtot 64 kuj + dj ()

uj + dj 29

3 77 77 77 5

2 66 6 S = 66 64

0 0 0

s s()

(k)

3 77 77 77 5

(76)

Q is the state vector to be solved for. Fj is the ux vector and S contains the source terms present in the turbulence model.

Note that the second term of each vector corresponds to the momentum equation, which in the general three-dimensional case consists of three components. Hence, in 3D equation (75) is a system of seven partial dierential equations. The diusion terms (d(jk) and d(j)) and the source terms (s(k) and s(k)) present in the equations for k and  are dened by: "  # @

t @k dj = ; @x +  @x j k j

(77)

"  #

@ t @ dj = ; @x +  @x j  j

(78)

@ui ;  ; D s(k) def = ijturb @x

(79)

(k)

()

def

def

j

def

s = ()



@ui C1 f1 ijturb @x j

!

; C2 f2  k + E

(80)

3.2 Space Discretization To obtain a nite-dimension equation-system in space the governing equations (75) are discretized on a structured non-orthogonal multi-block boundary tted mesh - see gure 3. Integrating equation (75) over a control volume  gives: Z @Q 

@t dV +

Z @Fj 

Z

@xj dV =  S dV

(81)

If the grid is stationary equation (81) can be rewritten with Gauss' theorem as: Z d V dt (Q) + @ Fj  dSj = V S

30

(82)



Boundary

Figure 3: A 2D structured non-orthogonal boundary tted mesh Where Q and S are cell averages of Q and S according to: Z Q = V1  Q dV

(83)

Z S = V1  S dV

(84)

V is the volume of the control volume  given by V = R dV . In words equation (82) says that the rate of change of the average state vector is equal to the sum of all uxes over the boundaries plus the source term. The ux term can be considered as the sum of the contributions from each face of the control volume : Z @

Fj  dSj =

XZ faces face

Fj  dSj

Assuming Fj to be constant on each face this can be rewritten as: 31

(85)

XZ

X

faces

faces

Fj  dSj = face

Fjface  Sjface

(86)

Where Sjface is a surface normal vector to the face with a length equal to the area of the face. To summarize, the governing equations (75) discretized on a structured mesh can be written as: X V dtd (Q) + Fjface  Sjface = V S

(87)

faces

The cell-averaged solution vector, Q, is stored at each cell, and the solution is advanced in time by evaluating the ux terms and the source terms in equation (87). The time-step is performed with an explicit 3-stage Runge-Kutta method. The problem then becomes to approximate the uxes and the source terms in terms of Q at each cell so that an accurate and stable method is obtained. For simplicity the over-line to denote cell-averaged values will be omitted in all subsequent sections.

3.3 Flux Term The ux term is split into a convective (Euler) part, and a diusive (viscous) part according to:

Fj = Fjconv + Fjdiff

(88)

Where 2 u 66 u u +j p ij 6 i j Fjconv = 66 (e0 + p)uj 64 kuj

uj

3 77 77 77 5

2 0 66 ; ijtot 6 Fjdiff = 66 qjtot ;(ku)i ijtot 64 dj

d(j)

32

3 77 77 77 5

(89)

Due to the dierent nature of the convective and diusive ux they are treated separately and in dierent ways.

3.3.1 Convective Flux A

B

C

D

Face Treated

Figure 4: Four Neighboring Cells of the Treated Face (2D) The convective ux across the cell wall in gure 4 is estimated with a four-node characteristic variable based scheme, which has a user-controllable up-winding based on the propagation direction of the characteristic variables at the cell wall. It is more convenient to work with the primitive variables, q: 2 66 6 q = 66 64

 ui p k 

3 77 77 77 5

(90)

Note that q in each cell can be calculated from Q in the same cell with second order accuracy. The inviscid parts of the governing equations are linearized around a point at the cell wall dened by (91). This gives a 1D system of hyperbolic equations.   qf = 21 qB + qC

(91)

By transforming to characteristic variables, r(n) , a decoupled set of equations can be obtained. The sign of the eigenvalues at the cell wall tells in which direction information about the corresponding characteristic variables travels. By choosing the direction of the interpolation stencil based on the sign of the characteristic velocities (eigenvalues) it is possible to construct an upwind-biased scheme. 33

The eigenvalues are given by:

1 2 3 4 5 6 7

ufj  Sj ufj  Sj ufj  Sj q ufj  Sj + cf qSj  Sj ufj  Sj ; cf Sj  Sj ufj  Sj ufj  Sj

= = = = = = =

(92)

Where, Sj = nj S , is the face area normal vector, and cf is the speed of sound at the cell face. The characteristic variables in the four neighboring cells are calculated with the following expressions:

r(1) =  ; (cfp)2 r(2) =

p

;S2 u1 +S1 u2 S12 +S22

r(3) =

;S3 (S1 u1 +S2 u2 )+(S12 +S22 )u3

r(4) =

f 2cf

pS 2+S 2 pS 2+S 2+S 2 1

2

1

2

3

3 u3 + pf 2  S1pu1S+S12+S2u222+S 2(c ) +S32

(93)

2 u2 +S3 u3 + pf 2 r(5) = ; 2cff  S1pu1S+S2+S 2 +S 2 2(c ) 1

2

3

r(6) = k r(7) =  The characteristic variables at the face are estimated from the characteristic variables in the four neighboring cells as in equation (94). Up-winding is obtained by making the direction of the interpolation stencil dependent on the sign of the eigenvalues.

34

(

( n) (n) (n) (n) C if n 0 1  rA + C2  rB + C3  rC + C4  rD rface = ( n) (n) (n) (n) C4  rA + C3  rB + C2  rC + C1  rD if n < 0 (n)

(94)

The interpolation coecients, Cn, can be set at runtime. The most common choices are:

 First order upwind: Cn = (0 1 0 0)  Third order upwind: Cn = (;1=6 5=6 2=6 0) The characteristic variables at the face can then be transformed back to primitive variables with the following expressions: (1) (4) (5) face = rface + rface + rface

(2) S2 rface S12 +S22

u1face = ; p

(3) S1 S3 rface S12 +S22 S12 +S22 +S32

p

;p

(2)

S1

 (4)

(5) rface;rface

+ cff  pS 2 +S 2+S 2

(3)

 (4)1



2

3  (5) S2 rface ;rface S12 +S22 +S32

cf u2face = + pS1Sr2face ; pS 2+SS22Sp3 rSface 2 +S 2 +S 2 + f  p +S 2 1

2

pS 2 +S 2

(3)

1

2

 (4)1

3  (5) S3 rface ;rface S12 +S22 +S32

u3face = rpfaceS 2+S 21+S 22 + cff  p 1

2



3

(4) (5) pface = (cf )2  rface + rface

2

(95)



(6) kface = rface (7) face = rface

3.3.2 Di usive Flux To calculate the diusive ux it is not only necessary to estimate the primitive variables at the faces, but also gradients of the primitive variables are needed.

 Primitive variables are estimated as the average of the two neighboring cells, as in equation (91).

35

 Gradients of primitive variables are calculated with a second order compact

scheme centered around the face. This scheme can either be interpreted as a nite dierence scheme or a nite volume scheme, but both interpretations are equivalent.

Finite Volume Interpretation (J)

Sj

D

C

Ω’ B

A

(I)

Sj

F

E

Face where the diffusive flux is calculated

Figure 5: Finite Volume Interpretation (2D) To calculate the gradient, @x@qi , at the face between cell A and B in gure 5, a new control volume 0 is dened by displacing the cell face along the centerline between the two neighboring cells. In 3D 0 is a parallelepiped, and in 2D it is a parallelogram. If @x@qi is assumed to be constant within 0 its value at the face can be calculated with Gauss' theorem as: 

@q @xi

!face

Z @q Z 1 1 = V 0 @x dV = V @0 q  dSi = i X 0 face 0 face = V1 q  Si 0 faces

(96)

Where V is the volume of 0 and Si face is the surface vector of each face. The 0

36

0 face value of q is taken as an average of the neighboring cells to each face. In 2D q0face is given by:

qface1 qface 2 qface3 qface4

= qB  = 14 qA + qB + qC + qD = qA  = 14 qA + qB + qE + qF

(97)

Finite Dierence Interpretation The gradient

@q @xi

at the face can be calculated with the chain rule as: 

@q @xi

!face



@q = @ j

!face 

@j  @x i

!face

(98)

Where the curve-linear coordinate system, i, is dened by the grid as in gure 6.

X2 ξ 2

Β Α

ξ1 X1 Figure 6: Physical and Computational Domain in 2D The two factors in the right hand side of equation (98) can be calculated as follows: 37



 @q face

is approximated directly using the closest cells.  @x face  @ j face ;1 i is calculated as the inverse of J = . Where J ;1  J = @x @ j i can be estimated using the nearest points in the grid, which is given by xi = xi(i). @ j

This leads to exactly the same second order compact scheme as the nite volume interpretation above.

3.4 Source Terms The source terms are evaluated as cell averages, using the state vector Q directly. To estimate the gradients present in the source terms a cell-centered nite volume scheme is used. This gives a second order accurate estimation of the gradients. Note that it is not good to use the gradient approximations at the cell faces. Instead a new cell-centered gradient approximation should be calculated, If the gradient @x@qi is calculated as the average over the control volume , it yields: 

@q @xi

!cell

Z @q Z X 1 1 = V  @x dV = V @ q dSi = V1 q  Si i faces

(99)

Where q on the faces is taken as the average of the two neighboring cells.

3.5 Time Integration The time step algorithm is a three-stage explicit Runge-Kutta scheme. In stationary simulations an optimized local time-step is used. The space-discretized governing equation (87) can be written as:

V  Qt = Sterm ; Fterm Sterm def = V S 38

(100)

Fterm def =

X faces

Fjface  Sjface

Where Sterm is the source terms and Fterm is the sum of all uxes. Qn+1 is calculated from Qn as follows:

 First step:

Q

 Second step:



= Qn + t 

1  (S ; F ) V term term Qn

(101)

1 1 1 1 n  Q = 2 Q + 2 Q + 2 t  V  (Sterm ; Fterm )  Q 

 Third step:

Qn+1 = 12 Qn + 21 Q + 12 t  V1  (Sterm ; Fterm)

Q

(102)

(103)

In stationary simulations the local time step, t, is determined from a local stability estimate in each cell. The scheme is stable for a CFL number up to about one, and t is given by (104). The eective spectral radius, SReff , is taken as the maximum of the convective and the diusive spectral radius (105). The convective spectral radius can be estimated with Fourier analysis of a local linearization of the Euler parts of the equations with locally frozen metrics 16]. The nal result is given in (106). An estimation of the diusive spectral radius is given in (107).

CFL t = SR

(104)

SReff = MAX (SRconv  SRdiff )

(105)

eff

SRconv =

(I ) (J ) (K ) p Si  ui + Si  ui + Si  ui + c M1 + 2M2

V

SRdiff = 2 + t 4M1 V+2 2M2 39

(106) (107)

Where the two metric terms are dened by:

M1 def = Si(I ) Si(I ) + Si(J )Si(J ) + Si(K )Si(K )





(108)

M2 def = Si(I ) Si(J ) + Si(I )Si(K ) + Si(J )Si(K )

(109)

Si(I ) , Si(J ) and Si(K ) are the surface normal vectors in each direction of the grid, V is the cell volume and c is the speed of sound.

3.6 Turbulence Modications It is a bit dicult to solve the turbulence equations (k and ) with an explicit time-stepping scheme. These equations are very source-term dominated, and thus they are very sensitive to small disturbances. At the beginning of a simulation it is necessary to take extra care to keep them from exploding in a few timestep. Many low-Reynolds k ;  models are also somewhat ill-behaved in that they can give negative values of k or  at some stages of a simulation. This is especially problematic when using a local time-step. Note that there are two dierent problems here! stability and too low k or  values, and that the latter often leads to the former. Under idealized conditions one can show that many low-Re k ;  models do not put any extra stability constraints on the local time-step. This has been done for the Chien model by Kunz et. al. 22] 23]. In reality it's not always that simple. Cases with low free-stream turbulence levels (< 1 %) can be dicult and demand extra care to converge. A large number of ways to improve stability have been tested and the experience can be summarized as follows:

 It's not possible to use the same third order scheme for k and  as for the

ow-equations. A TVD variant of the same scheme works for cases with high free-stream turbulence levels, but a few cases demand a rst order upwind scheme to be stable. The best convergence is often obtained by starting with a rst order upwind scheme and then switch to a high order TVD scheme as a physical ow eld has emerged.  During the early stages of a simulation it helps to set an upper limit on the turbulent viscosity. This limit can be relaxed in steps after a few thousand time-steps. It's not benecial to keep this limit too long. 40

 At startup, or when switching scheme or model, it helps to put an upper

limit on how much k and  are allowed to change in a time-step. This can be done by reducing the time-step locally if the magnitude of the total update is more than say 10 % of the cell value. It is better to treat k and  separately instead of reducing the time step for both equations if one of the variables change too much.  Point implicit treatment of parts of the source terms or the entire source terms does not improve stability signicantly, neither does various kinds of limitations on the production terms or =k.  The initial eld used to start a simulation is usually not very realistic, especially not the relation between the ow and the turbulent quantities. This, in combination with the use of a local time-step, causes k and  to reach unrealistically low levels at some stages of a simulation. The presence of shocks increases this problem, and cases with low free-stream turbulence levels are more problematic. If nothing is done this can cause a rapid fall in  and a following fall in k. Eventually this ends with a negative k or  value somewhere, if not instability has caused divergence before that. The stability problems are most pronounced in the  equation. The most eective way to prevent this from happening is to set a lower limit on k. A lower limit on  might seem more eective, but this only stops  from falling and does not help it to recover. Usually  gets stuck at the lower limit. The limit on k can be implemented by setting the update to zero if it is negative and the cell value is below a certain limit. The net eect of this is to inject some extra turbulent energy in problematic regions. This support of k has been sucient to prevent k and  from reaching too low levels in all cases considered. At startup the limit on k should be set to something like max(0.1% of free-stream k, Tu=0.001 %). Then it can be reduced in steps until it has no eect any more.  The k ; ! models have less numerical problems than the low-Re k ;  models. Setting the surface roughness, kR+, to a value below 1 can cause problems with some grids, but with a value above 1 convergence is good. Using these methods it has been possible to obtain converged solutions for all cases analyzed so far, including cascades with free-stream turbulence levels as low as 0.2 %.

41

3.6.1 TVD Limiter As mentioned above the third-order convective scheme is not stable enough to be used for the turbulent quantities k and . For most cases it is sucient to add an extra TVD (total variation diminishing) limiter to this scheme, which then becomes second-order accurate. The TVD limiter is applied only to k and , although it could be used for all the characteristic variables. Let r be a characteristic variable (k or  in this case). The parameters  and ' are then calculated as follows: 8 > < ' if = > 1 if > :

'