Valerius Boxberger
email: valerius.boxberger@tu-berlin.de

Roland Baar

Maike Gern

Rainer Zimmermann

Challenges of Turbocharging a Two-Cylinder Engine -A Computational Analysis of a Turbocharger Turbine by Transient CFD Methods

Keywords:

. As a result, the quasi-steady behavior, especially in two cylinder engines, does not seem to be valid and this has to be taken into account during the design process [2], [3]. This paper presents a numerical investigation of a pulse charged turbocharger turbine to gain an understanding of the aerodynamics at different pulse frequencies. A validated CFD-model is used to compute two and four-cylinder pulses. The time dependent boundary conditions for pressure and temperature are generated via 1D simulation. To make both pulse frequencies comparable, the turbine simulations were run with the same specific exhaust cycle enthalpy.

The analyses of the results shows a significant deterioration of the turbine operating behavior at low end torque due to the pulsating inflow of a two cylinder pulse compared to a four cylinder pulse. Higher engine speeds reduce the disadvantages. Furthermore it can be shown that the storage effect of the volute and the phase shift in temperature and pressure influences the throughput behavior of the turbine. For that the specific pulsating flow conditions must be considered for an optimization of the turbocharger of a two-cylinder engine.

INTRODUCTION

More stringent emissions standards and a reduction in fuel consumption represent the major challenges for automotive industry. Displacement reduction is considered to be a suitable way to achieve the legal requirements. Schernus et. al [START_REF] Weinowski | Zukünftiges Downsizing bei Ottomotoren-Potentiale und Grenzen von 2-und 3-Zylinder Konzepten[END_REF] shows that the lowest CO 2 emissions for vehicles with an inertia class less than 1130 kg occur at a displacement of 0.8 liter. Two-cylinder engines are at advantage because of their better single-cylinder volume and due to this light consumption benefits compared to three-cylinder engines. The radial turbines operate under strongly non-steady exhaust pulses. The boundary conditions of a turbocharger turbine therefore changes depending on the engine firing frequency [START_REF] Reuter | Erweiterung des Turbinenkennfeldes durch ein innovatives Messverfahren als Randbedingung für die Ladungswechselsimulation[END_REF]. Furthermore, it has been found that test bench measurements do not match the pulsed conditions. Deviations of cycle averaged turbine performance; the efficiency and the mass flow rate differ from stationary measured turbochargers [START_REF] Cao | Radial Turbine Rotor Response to Pulsating Inlet Flows[END_REF]. This makes an accurate the turbocharger modeling in the 1D simulation very complex. For two-cylinder engines, the pulsating exhaust gas flow has a particularly higher priority, since the bigger fluctuations in the pressure of the exhaust gas, the mass flow and the temperature have a far greater impact on the turbocharger behavior than for engines with more than two cylinders. Real engine conditions are difficult to realize on turbocharger test benches. The experimental measurement of high-frequency temperatures and mass flows is still a major problem [START_REF] Aymanns | A Revision of Quasi Steady Modelling in the Simulation of Pulse Charged Engines[END_REF]. This can be attributed to the thermal inertia during temperature measurements. A reliable measurement of the data can still not be carried out satisfactorily. Changes in the flow direction of the fluid further complicate the measurements. A detailed numerical analysis of the inlet conditions on a radial turbocharger turbine is performed. First a steady state validation with experimental test bench data is done. Furthermore transient real engine boundary conditions are generated for a two-cylinder engine in three different operating points via 1D simulation. Additionally, three artificial fourcylinder engine boundary conditions are generated which allow the usage as a reference to quantify the transient twocylinder pulse turbine performance.

EXPERIMENT

Numerical methods such as 1D and 3D simulations require a preliminary validation of the models. In this context, both the engine and the turbocharger used were stationary measured on test benches (engine and turbocharger).

Test engine

A 0.8 liter two-cylinder diesel engine with a performance power of 35 kW was used for the validation of the 1D simulation model. During the test bench measurements, the speed of the turbocharger, the stationary pressure and temperature were measured before and after the turbine for different engine operating points. The pressure signal was measured with high temporal resolution using pressure quartz at the same place as the steady-state values. The measurement result has been used for the validation of 1D simulation engine model. Figure 1 shows the pressure gradient over the crank angle in addition, the valve lift is shown. More engine details can be found in [START_REF] Zillmer | The Drivetrain for the New Volkswagen "1-Litre Car[END_REF].

Test turbocharger

The turbocharger used is currently the smallest available turbocharger on the market for diesel engines (type designation: GT0632SZ, manufactured by the company Honeywell Turbo Technologies). It is controlled via a waste gate. The diameter of the turbine impeller is 30 mm, the diameter of the compressor wheel is 32 mm. For the present study, a standard map on the hot-gas test bench for turbochargers was measured at the Department of Internal Combustion Engines at the TU Berlin with a turbine inlet temperature of T turbine inlet = 600°C. Furthermore a second adiabatic map was measured for the validation of the CFD model. Finally a sensitivity analysis of the pulse was carried out. In this, an artificial pressure pulse was generated by a cylinder head which was fed with compressed air. Different pulses frequencies could be modeled by changing the speed of the camshaft. The sensitivity analysis was used to examine the influence of the pulse frequency on the turbocharger speed.

NUMERICS METHODS AND VALIDATION

The detection of high-frequency flow phenomenon can currently be detected only on the numerical path. For this reason, a stationary model is presented here. This corresponds geometrically exact to the turbocharger test bench. Additionally a second transient model is presented. This uses the same mesh for the impeller and the housing, only the inlet and the outlet have been changed to save simulation time. The transient simulations run with 1D simulation generated boundary conditions and for this does not require the test bead measuring tubes (inlet and outlet).

Numerical model

The flow calculations of the turbocharger are based on a finite volume approach for solving steady and unsteady RANS/ URANS equations (Reynolds-Averaged-Navier-Stokes/ Unsteady-Reynolds-Averaged-Navier-Stokes). The simulations were carried out with the commercial code ANSYS CFX, version 15.0.7. The tool uses an implicit discrete approach. The simulations were run with the setting high resolution which uses the second order backward Euler scheme when it is required to maintain a bounded solution it switch to the first order backward Euler. The SST-turbulence model (SST = shear stress transport) was used in combination with the wall function. The transfer of information between the rotating and stationary domains was carried out by the frozen-rotor approach, which has been proven by [START_REF] Aymanns | A Revision of Quasi Steady Modelling in the Simulation of Pulse Charged Engines[END_REF] and [START_REF] Roclawski | Influence on Degree of Reaction on Turbine Performance for Pulsating Flow Conditions[END_REF] for steady and transient simulation. The walls were modelled as adiabatic and ideally smooth and the turbulence intensity was set to 5%. The break up criterion for the simulations has been set to the maximum residuals (R max = 1.0 • 10 -3 ). During the simulation the mass flow, pressure ratio and the temperatures at the inlet and outlet was monitored. Table 1 shows the mesh types which were used for the discretization. Stationary CFD model: The test bench structure is modelled as shown in Figure 2 (left). The time step is calculated as follows:

time step = 1 ω TC (1)
ω TC describes the turbochargers angular velocity. Measuring tubes and the manifold was integrated into the model to consider the right pressure loss in the system. In addition a mesh independence study was carried out. At the outlet, a constant steady pressure was defined. The prior evaluation of the engine tests showed that the rotational speed of the turbocharger during the 2C engine operation varies by less than 1%, see Figure 4. The variation during the 4C operating is estimated to be even lower. For that, the turbocharger speed was set as constant.

Figure 3: Boundary conditions transient simulation

The time steps were adjusted to the operating point. Since the simulation was run with the frozen rotor approach only the pulse frequency had to be resolved which led to a bigger time step. For the different operating points the time step has been set as follows: Δt OP1 =8.0 • 10 -5 𝑠, Δt OP2 =6.85 • 10 -5 𝑠, and Δt OP3 =5.33 • 10 -5 𝑠. The transient simulations were calculated with the 1.0 million mesh model. It has been used in order to decrease simulation time. Figure 2 (right) shows it. The initial conditions for the transient simulated were taken from a previous calculated steady state simulation.

Validation

High temperatures on the turbine side affect the oil in the bearing housing and also the compressed fresh air on the compressor side during the experimental measurement. It thus acts directly on the experimental determined efficiency. shows minor deviations, Figure 6. These were accepted to shorten the computation time.

TRANSIENT BOUNDARY CONDITIONS

To create a valuation basis for the operating behavior of the turbine during a 2C pulse, an additional artificial 4C pulse was modeled. A 4C pulse was selected as reference because of its high acceptance in the European market. Due to this a 4C pulse counts as the standard operating conditions for turbochargers. Two exhaust characteristics (2C and 4C) were calculated for three engine operating points.

The investigation of the thermodynamic behavior of the gas turbine requires highly resolved pressure and temperature gradients of the exhaust gas over time. These are set as boundary conditions for the numerical simulation. Because of the problematic detection of the experimental data, a 1D simulation model of the 2C engine was used. The model had previously been validated with the engine test bench data. Three different operating points of the 2C engine were calculated. 155.000 125.000 255.000 The artificial 4C pulse could be created by firstly adjusting the ignition angle. A change was made from 360 °CA (2C) to 180° CA (4C). Eq. 2 shows the approach for determining the specific exhaust gas enthalpy, equal to [START_REF] Vogt | Aufladesysteme für Ottomotoren im Vergleich[END_REF].

h exhaust,tt = h tt,3-4 = ∫ c p,3 • T 3t (φ) [1 -( p 4t (φ) p 3t (φ) ) κ-1 κ ] 720°0 °• dφ (2)
In addition the total inlet pressure and the total inlet temperature were integrated over a cycle and then scaled down to the level of the 2C pulse, see Eq. 3-7. As a result a 4C pulse with specific cycle exhaust gas enthalpy of a 2C engine was created.

p 3t,2 Cyl = ∫ p 3t (φ) 720°0 °• dφ (3)
T 3t,2 Cyl = ∫ T 3t (φ) 720°0 °• dφ (4)
For both engine pulses, an equal turbocharger speed was set for equal operating points. The mass flow of the turbine adjusts itself according to expansion ratio and the efficiency.

Dispersed engine operating points were selected to represent a wide operating range on the engine performance map. The selected points can be seen in Table 2. The first operating point (OP: 1) is a full load point at lower engine speed (1500 rpm), while OP: 2 is a partial load point (1750 rpm) and OP: 3 describes the point with the maximum torque (2250 rpm).

RESULTS AND DISCUSSION

Aerodynamics

Four-cylinder pulses show higher turbine mass flow rates; see Figure 8. The through flow difference is most evident in OP: 1; here the value is 28% higher compared to the 2C. The other operating points differ in the order of 16%. Schernus et. al [START_REF] Schernus | Turboaufladung für kleine Ottomotoren mit weniger als vier Zylindern, Motorprozesssimulation und Aufladung, III: Engine Process Simulation and Supercharging[END_REF] made comparable observations. They note that the turbine at 2C pulses seems to ''behave'' smaller than at pulses with a higher cylinder number. The reduced through flow can also be seen in Figure 9. Here the reduced turbine mass flow rate is plotted against the expansion ratio. Heat transfer influences on the housing, speed variations and temperature effects on the mass flow can be excluded because of the adiabatic boundary conditions and the reduced mass flow. Furthermore, the reduced turbine mass flow in all operating points is shown. The red hysteresis loop represents the system boundary in the sections 3-4 (complete stage) for the 2C pulse. Additionally the black hysteresis loop represents the system boundary in the sections 31-41 (impeller) for the 2C pulse. Figure 10 shows the diagrams for the 4C pulse. Comparable maximum through flows can be monitored between the 2C and the 4C pulse in every operating point. This behavior can be explained by the flow function. Here the reduced mass flow is a product of the flow function Ψ, the ideal gas constant 𝑅 and the effective turbine crosssection 𝐴 𝑇,𝑒𝑓𝑓 , see Eq. 8.

According to Eq. 9 the flow function is dependent on 𝜅 and the expansion ratio. It cannot be increased from a critical expansion ratio, see Eq. 10. The maximum through flow adjusts independent from the pulse frequency.

Ψ 31-4 = √ 𝜅 𝜅 -1 • [( 𝑝 4 𝑝 31 ) 2 𝜅 -( 𝑝 4 𝑝 31 ) 𝜅+1 𝜅 ] (9)
After reaching the critical value, a rise in the reduced mass flow can only be caused by changes of the effective turbine cross-section.

Ψ 𝑚𝑎𝑥 = ( 2 𝜅 + 1 ) 1 𝜅-1 • √ 𝜅 𝜅 + 1 (10)
It is remarkable that the expansion ratio range of a 2C pulse is much larger than that of the 4C. At very small expansion ratios, the 2C even reaches reduces mass flow values close to zero. At an expansion ratio Π 𝑇,𝑂 , Eq. 11, the turbine is running into a stagnation point, as shown in Figure 11. Because of the higher pressure level such conditions are usually not reached by 4C pulses. At the stagnation point, the gas impulse forces are no longer sufficient to the centrifugal forces caused by the rotation of the impeller, see Eq. 11, and the flow stagnates, see Figure 11. In this case, a constant speed has been specified for the impeller rotation; due to this negative mass flows can occur at high turbocharger speeds and very small expansion ratios. Over a cycle, the expansion ratio is fluctuating much more at a 2C pulse. At high expansion ratios the maximum through flow is reached and cannot be increased further. At small expansion ratios, the through flow stagnates. The 4C pulse makes the turbine run more often in a subcritical range. Thus the turbine has a better through flow.

Π 𝑇,𝑂 = [1 - 𝜋 2 • 𝑛 𝑇 2 2 • 𝑐 𝑝𝑇 • 𝑇 3𝑡 • (𝐷 3 ² -𝐷 4 ²)] 𝜅 𝑇 1-𝜅 𝑇 (11)
The emptying period shows the highest reduced mass flow. During a pressure pulse, different tip speed ratios can occur at the same expansion ratio 𝑝 4𝑡 𝑝 3𝑡

. The pressure propagates with the speed of sound. Due to the traveling time it comes to an offset of the expansion ratios in the sections 3 and 31. The boundary conditions at section 3 do not correspond to the boundary conditions at section 31. The turbine mass flow rate depends only on the conditions in sections 31. The offset in the reduced mass flow between 3 and 31 can be attributed to the storage effect in the volute. Furthermore, non-optimal flow conditions in front of the impeller influence the mass flow rate. The tip speed ratio is composed of the circumferential speed of the impeller and the isentropic nozzle velocity, see Eq. 12. An optimum tip speed ratio is in the range of approximately 0.6-0.7, depending on the design of the turbine. The tip speed ratio is an indicator for the optimum inflow angle. It affects both the turbine efficiency and the through flow performance.

𝑆 31 = π • D 3 • n 𝑇 √ 2 • c pT • T 31t • [1 -( p 41 p 31t ) κ-1 κ ] (12)
Here the through flow is dependent on the effective turbine cross-section. The effective turbine cross-section can be determined as follows:

ṁ𝑟 𝑒𝑑,31 = ṁ𝑇 • √𝑇 31𝑡 𝑝 31𝑡 = 𝐴 𝑇,𝑒𝑓𝑓 • √2 √𝑅 • Ψ ( 𝑝 41 𝑝 31𝑡 ) (13) 𝐴 𝑇,𝑒𝑓𝑓 = ṁ𝑟 𝑒𝑑,31 √2 √𝑅 • Ψ ( 𝑝 41 𝑝 31𝑡 ) ( 14
)
According to [START_REF] Pucher | Aufladung von Verbrennungsmotoren[END_REF] there is a relationship between the effective turbine cross-section and the tip speed ratio. At constant expansion ratios, the effective turbine cross-section decreases with increasing speed ratio. Thus high tip speed ratios reduce the turbine through flow. The tip speed ratio depends on c pT and κ and the total temperature. In accordance with the numerical boundary conditions, c pT and κ can be considered as constant. The temperature T 31t and the expansion ratio

p 41 p 31t
change during a pressure pulse. The expansion ratio has the greatest impact on the tip speed ratio.

Figure 12: Influence of the tip speed ratio (2C)

The comparison of the tip speed ratio in Figure 12 and Figure 13 show clearly the benefits of 4C pulse. The turbine is operating in a nearly constant tip speed ratio, close to the design point. Under 2C pulse the turbine is operating in poor regions, especially at low expansion ratios.

Specific power provision

Due to the 1st main equation of turbocharging there is a equilibrium of power between turbine and compressor. Both are linked by a common shaft. In a steady state point the equilibrium determines the turbocharger speed. In a transient simulation a deviation of the compressor or turbine power leads to a change in turbocharger speed. Furthermore the equation for the equilibrium has to be expanded by the inertia power which acts during acceleration or deceleration. The turbine power can be determined as follows:

P T = P C + P Fr (+P In ) (15)
The power output of the turbine is in charge for the intake mass flow, the specific heat capacities of the fluid and temperature. It can be calculated as: In the following the specific power output of the turbine will be analyzed. To evaluate the specific power provision of the radial turbine, the specific enthalpy over the entire turbine stage was measured, see Figure 14.

P T = 𝑚̇𝑇 • c p,exhaust • (T 3t -T 4t ) (16)
The overall highest values for the turbine enthalpy can be monitored at the 2C pulse. This is due to the partial higher expansion ratio and temperature.

For both pulse frequencies very small or negative power output could be determined. This can be traced back to small expansion ratios and partly reversed flow for the 2C pulse.

The equation for the specific exhaust gas enthalpy is shown in Eq. 17.

h T,is,tt = c p,exhaust • (T 3t -T 4t ) (17)
The consideration of enthalpy integral over the cycle shows that the biggest advantages of the 4C pulse are at low speeds and high loads. With increasing speed the specific enthalpy of both pulses is equal. The partially negative specific enthalpy of the 2C is compensated to the higher speeds and higher inlet pressures and temperatures.

Turbine efficiency

The determination of the isentropic efficiency under real pulsating conditions represents a challenge and has been a topic of research for many years; see [START_REF] Aymanns | A Revision of Quasi Steady Modelling in the Simulation of Pulse Charged Engines[END_REF], [START_REF] Cao | Radial Turbine Rotor Response to Pulsating Inlet Flows[END_REF], [START_REF] Lüddecke | Engine crank angle resolved turbocharger turbine performance measurements by contactless shaft torque detection[END_REF], [START_REF] Reuter | Erweiterung des Turbinenkennfeldes durch ein innovatives Messverfahren als Randbedingung für die Ladungswechselsimulation[END_REF]. In this chapter the problem of the determination of the temporal resolved efficiency is demonstrated based on the operating point BP: 3.

Turbine stage (section 4-3):

The standard determination of the pulsating efficiency leads to unphysical values as shown in Figure 15. Incorrect efficiency values are mainly due to storage effects. These occur in the components of the turbine (e.g. volute casing and turbine housing). Delayed gas exchange in the volute affects the total temperature after the impeller; see Figure 18 (top left).

η T,is,tt = 1 - T 4t T 3t 1 -( p 4t p 3t ) κ-1 κ (18)
While the temperature 𝑇 3𝑡 decreases in proportion to the falling pressure (120 ° CA) the temperature 𝑇 4𝑡 remains in this area until 𝑇 3𝑡 briefly falls below 𝑇 4𝑡 . The turbine housing is filled with old hot exhaust gas at this point. When the temperature of the gas is falling the mass flow is insufficient to push out the old gas. The fresh cold gas mixes with the old gas. A mixing temperature is measured at section 4, which does not allow conclusions to be drawn about the turbine efficiency. Figure 16 shows the moment when the space behind the impeller is filled with warm gas and fresh cold gas is exiting the impeller. The efficiency becomes negative at this point. Figure 17 shows the moment when the space behind the impeller is filled with old cold gas and fresh warm gas is entering the impeller. The old gas mixes with the fresh gas.

The resulting mixing temperature becomes lower than the real turbine outlet temperature. The consequences are isentropic efficiencies greater than one. The gas in front of the impeller starts recirculating and heats up. 𝑇 4𝑡,𝑖𝑠 increases in dependence to 𝑇 3𝑡 , and simultaneously the pressure level decreases over the turbine. 𝑇 31𝑡 remains temporarily in the stagnation region at an elevated temperature level, see Figure 19. This leads to isentropic efficiencies greater than one.

CONCLUSION

A numerical study of the turbine behavior of a turbocharger turbine at different engine pulse frequencies was carried out. The investigations were performed by transient 3D CFD simulations. A 2C engine pulse and a 4C pulse have been investigated. The artificial 4C pulse was created to give the analysis a quantitative reference.

2C engines have a doubled ignition angle compared to 4C engines. This results in a significantly increased time interval between the exhaust valve openings. The fluctuation of inlet pressure and temperature increases (more time for the impeller to reduce the pressure until the next valve opening). During this analysis the cycle averaged pressure and temperature and thus the specific exhaust enthalpy for both types of pulses were hold constant. As a result the pressure and temperature peaks were smaller at the 4C pulse (the pressure was spread over four instead of two cylinders.)

The analyses shows the biggest disadvantages of the 2C pulse in the area of low end torque. The variations in the expansion ratio are highest here. Small expansion ratios lead to stagnation or back flow, furthermore very high expansion ratios lead to choke. The influence on the cycle integrated reduced mass and specific power output is significant at high engine loads and small speeds. Operating points with higher engine speed reduce the 2C disadvantages.

The major finding in this paper is that the hysteresis loop is not only influenced by storage effects in the volute. Observing only the impeller in the system boundary does not make the hysteresis disappear. A second cause is the effect of phase shift in temperature and pressure in the volute. Upstream of the impeller the shift changes the blade speed ratio (high speed ratios lead to lower effective turbine cross-sections) and due to this reduces the swallowing capacity of the turbine.

The determination of the transient turbine efficiency is still a major problem. A reduced system boundary (only impeller) eliminates the storage problem and is simply done by the CFD methods. During the investigations it is shown that this method is only valid for medium and high expansion ratios. The calculation of the transient efficiency close to the stagnation point is not possible. Due to reverse flow and friction caused heatup of the gas upstream the turbine impeller.

Figure 1 :

 1 Figure 1: Pressure indication before and after turbine

  temperature were set at the volute inlet, see Figure3.

Figure 2 :

 2 Figure 2: 2.3 million model (left), 1.0 million model (right)

Figure 4 :

 4 Figure 4: TC speed reductionSince an accurate numerical simulation of heat transfer in the turbocharger is hard to achieve (detailed solid temperature, solid materials and y+=1-2 required) commonly adiabatic validation are carried out. In this case adiabatic conditions could be reached by aligning the temperature of fluid in the compressor, the turbine and the oil temperature to reduce heat flux (T T,inlet = T Oil = T 𝐶,𝑜𝑢𝑡𝑙𝑒𝑡 ). In this way, nearly adiabatic conditions are achieved for the turbocharger test bench measurements.

Figure 5 :

 5 Figure 5: Validation of the mass flowTo ensure a sufficient mesh density a grid independence study was carried out. Every domain of the turbocharger model was refined separately. The density was increased until the changes in the monitored values were negligible. The employed RANS code uses a wall-function for the boundary layer resolution. To ensure an optimum functionality of it, at least 10 nodes were places in the boundary layer. A preferably small model was needed for keeping the transient simulation feasible in a reasonable time. The validation shows the qualification of the 2.3-million nodes mesh for further investigations. To quantify the error of the 2.3-million-model a second model with a doubled mesh number has been generated. The two models were compared to each other. Only small deviations could be observed between the integral values (reduced mass flow rate and a deviation in the turbine efficiency of less than 2%, see Figure5. The detailed analysis of the internal velocity distribution

Figure 6 :

 6 Figure 6: Internal velocity distribution

p 3t, 2

 2 Cyl = p 3t,4 Cyl (5) T 3t,2 Cyl = T 3t,4 Cyl (6) h exhaust,tt,2 Cyl = h exhaust,tt,4 Cyl

Figure 7 :Figure 7 (

 77 Figure 7: specific exhaust enthalpy Figure 7 (left) shows the specific enthalpy gradients over the crank angle. The variations of enthalpy during a working cycle can be seen. The full load point shows the greatest variations.Figure 7 (right) shows the integrals of the specific exhaust gas enthalpy. The conformity of the integral values allows the use of these boundary conditions for further investigations.

Figure 8 :

 8 Figure 8: mass flow integral over cycle ṁ𝑟 𝑒𝑑 = 𝐴 𝑇,𝑒𝑓𝑓 • √2 √𝑅 • Ψ (8)

Figure 9 :

 9 Figure 9: reduced turbine mass flow 2C

Figure 10 :Figure 11 :

 1011 Figure 10: reduced turbine mass flow 4C

Figure 13 :

 13 Figure 13: Influence of the tip speed ratio (4C)

Figure 14 :

 14 Figure 14: Cycle integrated specific turbine power

Figure 15 :

 15 Figure 15: Isentropical turbine efficiency section (3-4) (left), section (31-41) (right)

Figure 16 :

 16 Figure 16: Efficiency values eta<0 (cold filling)

Figure 17 :Figure 18 :

 1718 Figure 17: Efficiency values eta>1 (hot filling) Turbine impeller (section 31-41): The influence of the storage effects can be avoided by reducing the system boundary to the sections 31-41. Figure 15 (right) and Figure 18 (right) show a corresponding attempt. The impeller can be treated as quasi-stationary, [1], [2].

Figure 19 :

 19 Figure 19: Temperature gradient at the stagnation point

Table 1 :

 1 Mesh overview

	stat. model	trans. model
	inlet (hexa)	volute (tetra)
	manifold (tetra)	impeller (hexa)
	volute (tetra)	body (tetra)
	impeller (hexa)	outlet (hexa)
	body (tetra)
	outlet (hexa)
	The operating points were set by specifying a constant total
	temperature, total pressure at the inlet and a static pressure at
	the outlet. The data for the boundary conditions was imported
	from the adiabatic turbocharger test bench measurements.
	Finally the turbocharger speed and the material values for
	exhaust gas were set.
	Transient CFD model: The boundary conditions for the
	transient simulation were derived from 1D simulation (crank
	angle resolved). They are modelled as a time dependent
	function. The temporally varying total pressure and the

Cyl. 1 Cyl. 2 Cyl. 1 Cyl. 2 Cyl. 1 Cyl. 2 OP:1 OP:2 OP:3

	CA [°]

Table 2 .

 2 Investigated engine operating points

	OP	rpm Engine	p me	rpm TC
	#.	[1/min]	[bar]	[1/min]
	1	1500	12.6
	2	1750	6.85
	3	2250	18.6


Zahir Hussain
email: z.hussain@mmu.ac.uk

Keywords: Rotating boundary-layer -crossflow instability -centrifugal instability -broad, slender rotating cone -co-rotating, counter-rotating spiral vortices

In this study, a new centrifugal instability mode, which dominates within the boundary-layer flow over a slender rotating cone, defined by half-angle ψ < 40 • , is used for the first time to model the problem when an enforced oncoming axial flow is introduced. The resulting similarity solution represents the basic flow more accurately than previous studies in the literature. This mean flow field is subsequently perturbed leading to disturbance equations that are solved via numerical and analytical approaches, importantly yielding favourable comparison with existing experiments. Meanwhile, a formulation consistent with the classic rotating-disk problem has been successful in predicting the stability characteristics of broad rotating cones, defined by half-angle ψ > 40 • , in axial flow.

INTRODUCTION

This article forms part of a series of studies which have used theoretical techniques to construct the correct models of governing instability for both broad and slender rotating cones. The current study represents a significant extension to the general problem in the slender cone case, introduced when enforcing an oncoming axial flow.

Physically, the problem represents an accurate model of axial airflow over rotating machinery components at the leading edge of a turbofan. In such applications, laminar-turbulent transition within the boundary layer can lead to significant increases in drag, resulting in negative implications for fuel efficiency, energy consumption and noise generation. Consequently, delaying transition to turbulent flow is seen as beneficial, and controlling the primary instability may be one route to achieving this. Ultimately, control of the input parameters of such a problem may lead to future design modifications and potential cost savings.

Our results are discussed in terms of existing experimental data and previous stability analyses on related bodies. Importantly, axial flow is seen to delay the onset of convective instability for both broad and slender rotating cones; the exact mechanism of interaction governing the transition process however is very different for both instabilities. Broad-angled rotating cones are susceptible to a crossflow instability visualised in terms of co-rotating spiral vortices, whereas slender rotating cones have transition characteristics governed by a centrifugal instability, which is visualised by the appearance of counter-rotating Görtler vortices. It is the relative com-petition of these two governing mechanisms that is explored in detail in this study, particularly with regard to the role of travelling modes in the breakdown process.

METHODS

We consider a cone of half-angle ψ rotating in a fluid of kinematic viscosity ν * with an angular velocity Ω * in an anticlockwise direction around the streamwise coordinate axis x * (where a * denotes a dimensional quantity in all that follows). We construct coordinate axes aligned along with and perpendicular to the spiral vortices ( x * and y * , respectively), as shown in Figure 1. Further details of the relationship between the coordinate systems, including a required Mangler transformation, are provided in [START_REF] Hussain | Stability and transition of three-dimensional rotating boundary layers[END_REF]. These are shifted from the conventional streamwise and azimuthal coordinates, x * and θ, which are based on cylindrical polar coordinates. In such a problem, there exists a boundary layer close to the rotating cone surface characterised by the distance along the cone l * and defined by the Reynolds number, R, such that:

R = Ω * l * 2 sin ψ ν * .
With the important distinction of the inclusion of the oncoming axial flow, the physical problem is subsequently altered such that there now exists a dimensional local slip velocity at the edge of the boundary layer, obtained via a wellknown potential-flow solution (see for example [START_REF] Rosenhead | Laminar Boundary Layers[END_REF]), given by U e = C * x * m , where C * is a constant. We subsequently compare this velocity to the rotational velocity of the cone surface, given by V w = Ω * x * sin ψ, to obtain an important ratio, which in part characterises the problem, known as the rotational-flow parameter, given by

s = V w U e 2 ,
and used in [START_REF] Kobayashi | Linear stability theory of boundary layer along a cone rotating in axial flow[END_REF][START_REF] Kobayashi | Boundarylayer transition on a rotating cone in axial flow[END_REF][START_REF] Garrett | Boundarylayer transition on broad cones rotating in an imposed axial flow[END_REF]. In this study, we will use s to facilitate comparison of our results with experiments, as well as make reference to physical cases where the cone is rotating 'quickly' and the axial flow is increased from a zero value (ie. s decreasing from ∞), for example during the take-off phase of an aeroplane, once the rotating turbofans have reached an optimum rotational velocity. Essentially, the boundary-layer flow undergoes competition between the streamwise flow component, due to the oncoming flow, and the rotational flow component, due to effect of the spinning cone surface, which can be described mathematically in terms of the control parameters, ψ and s. We present the results of convective instability analyses for the boundary-layer flow over broad and slender rotating cones in a variety of imposed axial flows, based on large Reynolds-number and short-wavelength asymptotics, as well as numerical solutions obtained via an Orr-Sommerfeld stability analysis.

Importantly, the basic flow quantities Ũ and Ṽ are expressed as projections along the shifted spiral coordinates. However, due to the introduction of an oncoming axial flow, these are now functions of both the non-dimensional streamwise and surface-normal coordinates, x and z, respectively (where x is the streamwise direction over the cone, scaled on l * ). Hence, obtaining the base flows now require the solution of a system of PDEs (see [START_REF] Garrett | Boundarylayer transition on broad cones rotating in an imposed axial flow[END_REF] and [START_REF] Hussain | Stability and transition of three-dimensional rotating boundary layers[END_REF] for full details and discussion) as opposed to the projected von Kármán solution of a system of ODEs presented in [START_REF] Hussain | The centrifugal instability of the boundary-layer flow over slender rotating cones[END_REF] for the still fluid problem. While there exist similar numerical basic flow for-mulations used by [START_REF] Kobayashi | Linear stability theory of boundary layer along a cone rotating in axial flow[END_REF] and [START_REF] Kobayashi | Boundarylayer transition on a rotating cone in axial flow[END_REF], experimental verification of the basic flow is currently planned by R. Lingwood (personal communication, 2015).

At this point, we outline the important link between the standard surface-normal coordinate η and the modified surfacenormal coordinate η 1 scaled on boundary-layer thickness according to the new velocity scales applied to the basic flow boundary-layer equations. The coordinate stretching yields

η 1 = η m + 3 2s 1 2 sin ψ 1 2 . ( 1
)
This relation enables the shifted velocity profiles Ũ (x, η) and Ṽ (x, η), expressed in terms of the standard boundary-layer coordinates, to be written in terms of f (s, η 1 ) and g(s, η 1 ), which depend on s and the modified boundary-layer coordinate. The shifted basic flow quantities are written correctly in the form:

Ũ (x, η) = U (x, η) cos φ + V (x, η) sin φ, Ṽ (x, η) = U (x, η) sin φ + V (x, η) cos φ.
Here, U (x, η) and V (x, η) can be expressed in terms of the solution functions f (s, η 1 ) and g(s, η 1 ) (where indicates ∂ ∂η 1 ) obtained in [START_REF] Hussain | Stability and transition of three-dimensional rotating boundary layers[END_REF] and presented in [START_REF] Garrett | Boundarylayer transition on broad cones rotating in an imposed axial flow[END_REF] for ψ = 50 • , 70 • . The numerical solutions for f (s, η 1 ) and g(s, η 1 ) are obtained via the D03PEF NAG routine using a Keller box scheme and the method of lines. However, in this study, we remain consistent with the formulation presented in [START_REF] Hussain | Stability and transition of three-dimensional rotating boundary layers[END_REF], pertaining to the shifted basic flow quantities, which are essential when considering the slender rotating cone problem for ψ < 40 • . Specifically, we may write

Ũ (x, η) = U e Ω * l * sin ψ f (s, η 1 ) cos φ + s 1 2 g(s, η 1 ) sin φ = s -1 2 Û (s, η 1 ), (2)
Ṽ (x, η) = U e Ω * l * sin ψ f (s, η 1 ) sin φ + s 1 2 g(s, η 1 ) cos φ = s -1 2 V (s, η 1 ), (3)
where Û and V are presented in Figure 2 for ψ = 15 • in a range of axial flows, increasing from s = 1.5 to s = 16 (corresponding to a 'quickly' rotating cone). We note that the x-component, Û exhibits a familiar inflexional nature, with its limiting value at the edge of the boundary layer increasing as s increases. However, for s ≥ 5, we observe from the results of [START_REF] Kobayashi | Boundarylayer transition on a rotating cone in axial flow[END_REF] that φ = 0 • , which is consistent with our basic flow solution where Û recovers the streamwise basic flow component, f , to within a factor of s -1 2 . In contrast, the y-component of velocity V exhibits a uniform shear and is consistently reduced as s is increased.

We assume that the spiral waves are periodic in the xdirection and introduce periodicity into the perturbation quantities of vortex x-wavenumber a and ȳ-wavenumber b. Scal-ing our perturbation quantities on the boundary-layer thickness, we introduce a combined flow of the form

ũ * = Ω * l * sin ψ[{ Ũ (x, η), Ṽ (x, η), R -1 2 W (x, η)} +R -1 2 { ũ(η), ṽ(η), w(η)}exp(ia x + ib ȳ)].
Similarly, the pressure perturbation term scales as

p * = ( ρ * Ω * 2 l * 2 sin 2 ψ)R -1 p(η)exp(ia x + ib ȳ). ( 4
)

RESULTS AND DISCUSSION

Asymptotic analysis

We proceed to solve the governing disturbance equations (identical to those given in [START_REF] Hussain | Stability and transition of three-dimensional rotating boundary layers[END_REF]) to determine leading-and next-order estimates of the scaled Taylor number for neutrallystable modes, which arise due to the scaling analysis and loosely follows [START_REF] Hall | Taylor-Görtler vortices in fully developed or boundary-layer flows: linear theory[END_REF] for the Taylor problem of flow between concentric rotating cylinders. Importantly, we are able to form comparisons with results in the literature expressed in terms of Reynolds numbers. The corresponding Taylor number, which characterizes the importance of centrifugal to viscous forces, is given by

T = 2 cot ψ cos φ sin 4 ψ . (5)
In the axial flow problem, for a fixed ψ, this quantity is an output of the analysis and represents a measure of how s affects the physical flow characteristics. Upon incorporating the basic flow expansions, we expand the perturbation quantities and pose a WKB solution for small values of , where a = -1 for the wavenumber a in the xdirection. As for the still fluid problem, the dominant terms in the governing disturbance equations balance if we scale T ∼ -4 and W/V ∼ O( -2 ), resulting in identical perturbation expansions to those presented in [START_REF] Hussain | The centrifugal instability of the boundary-layer flow over slender rotating cones[END_REF] (reproduced here for clarity when manipulating subsequent quantities):

ũ =E(u 0 (η) + u 1 (η) + 2 u 2 (η) + . . .), ṽ = 2 E(v 0 (η) + v 1 (η) + 2 v 2 (η) + . . .), w =E(w 0 (η) + w 1 (η) + 2 w 2 (η) + . . .), T = -4 (λ 0 + λ 1 + λ 2 2 + . . .)
,

where λ = λ 0 + λ 1 + λ 2 2 + . . ., E = exp( i ϕ K (τ)dτ), ϕ = sin ψ h1
η. and following [START_REF] Hall | Taylor-Görtler vortices in fully developed or boundary-layer flows: linear theory[END_REF], the growth K is a quantity to be determined.

After solving the leading-and first-order problems following [START_REF] Hussain | Stability and transition of three-dimensional rotating boundary layers[END_REF] and scaling out various dimensional quantities, we obtain an estimate for the scaled effective Taylor number, which is given by

T = -4 2 √ s m + 3 1 2 1 V (s, 0) + 2.3381 × 3 1 3 | V (s, 0)| 2 3 V (s, 0) + s -1 2 V (s, 0) 2 cos φ V (s, 0) 2 + . . . . (6)
(s) = 30.2 • , 22.5 • , 13.6 • , 6 • , 0 • , 0 • , 0 • .
Logarithmic plots of the scaled asymptotic Taylor number against vortex wavenumber, -1 = a, are shown in Figure 3 for ψ = 15 • and various values of s. We note that the asymptotic branches presented capture the effects of the leading-and first-order estimates for T. The unstable region is above the curves and the stable region below. In general, we observe that increasing s leads to a trend of reducing the asymptotic Taylor number branch. Physically, this can be interpreted as promoting the more dangerous centrifugal instability mode, and hence destabilising the flow, which leads to a larger unstable region above the neutral stability branch, as depicted in Figure 3.

Numerical analysis

In this section, we develop the corresponding numerical solution, outlining the major differences between the axial flow problem formulated in this study and the still fluid case presented in [START_REF] Hussain | The centrifugal instability of the boundary-layer flow over slender rotating cones[END_REF]. These arise due to the fact that the basic flow quantities Ũ and Ṽ are now functions of the logarithmic spiral coordinates x and ȳ, as well as η. We manipulate the governing disturbance equations and subsequently express the basic flow terms in terms of η 1 by making use of the coordinate stretching [START_REF] Rosenhead | Laminar Boundary Layers[END_REF]. The analysis involves neglecting Coriolis terms and viscous streamline-curvature effects. Importantly, we note that the centrifugal mode under investigation differs from the streamline-curvature mode for large half-angle cones, which arises due to viscous effects of the cone surface. In contrast, the centrifugal mode for small half-angle cones arises from the centrifugal forces present in the mean flow for small ψ, owing to the effects of surface-curvature. Such centrifugal curvature terms are not neglected in the analysis and contain the Taylor number as a factor. Proceeding in this fashion yields a modified Orr-Sommerfeld (OS) equation for stationary disturbances within the system, given by

i ∂ 2 ηη -k 2 2 + Re √ s α 1 Û + β 1 V ∂ 2 ηη -k 2 - Re 2s (m + 3) sin ψ α 1 Û + β 1 V w = 0, (7)
where

α 1 = a sin ψ Re , β 1 = b sin ψ, k = α 2 1 + β 2 1
represent the vortex wavenumbers in the x-, ȳ-and effective velocity-directions, respectively, and ∂ 2 ηη = ∂ 2 /∂η 2 . Furthermore, Re = x sin ψ is the local Reynolds number, interpreted as the local non-dimensional radius of the cone surface from the axis of rotation.

Importantly, we note from the basic flow profiles that for s ≥ 5, the vortex activity is located at the wall, with the minimum of V (s, η 1 ) existing at η 1 = 0. However, for s < 5, the curve has a minimum slightly departed from the wall, indicating the location of vortex activity will not be at η 1 = 0. This correlation results as a consequence of the requirement of obtaining valid real solutions for the growth rate K when solving the governing eigenvalue equation at leading order, which itself arises by following the study of [START_REF] Hall | Taylor-Görtler vortices in fully developed or boundary-layer flows: linear theory[END_REF] for the Taylor problem of flow between concentric rotating cylinders. For the case of s < 5, the solutions obtained are not the most dangerous modes available, but we include them as they provide useful information about non-zero wave angles (spiral waves) for a 15 • rotating cone in axial flow. Furthermore, an interesting observation pertains to the related study of [START_REF] Hussain | The convective instability of the boundary layer on a rotating disk in axial flow[END_REF] on the rotating disk in axial flow, where non-stationary travelling modes become more important as the strength of oncoming axial flow increases.

Indeed, the results of the numerical analysis are shown in figure 4, which presents the numerically calculated waveangle, φ = tan -1 ( β 1 /α 1 ), versus various s for a 15 • rotating cone. We compare results from the present study with the experiments of [START_REF] Kobayashi | Boundarylayer transition on a rotating cone in axial flow[END_REF], observing close agreement. Both studies observe that in the regime of a 'quickly' rotating cone (s > 1), increasing s leads to a reduction in φ, to the point where φ = 0 for s ≥ 5, physically corresponding to the transition from spiral waves to circular or 'Taylor' vortices. Interestingly, this appears consistent with our basic flow and asymptotic solutions where for s ≥ 5, the vortex activity remains located on the wall at η = 0. Conversely for s < 5, the stronger axial flow acts to sweep vorticity in the streamwise direction, which is again consistent with our asymptotic findings, namely that the vortex activity undergoes a slight departure from the wall in this regime. Specifically, for s < 5 in the current problem, it appears that the location of vortex activity departing slightly from the wall suggests that travelling modes may grow as s is reduced and in fact become the most unstable modes in this parameter regime. Indeed, physically, the departure of a vortex from the wall suggests that vorticity within the boundary layer is no longer fixed on the cone surface, but is instead propagating or travelling in the effective velocity x-direction. Ultimately, in order to confirm whether travelling instabilities may harbour the most unstable modes for the slender rotating-cone problem, a further investigation would be required, taking account of time-dependent terms within the governing disturbance equations.

In general, we observe close agreement between our asymptotic and OS numerical stability results, as well as with the numerical calculations of [START_REF] Kobayashi | Linear stability theory of boundary layer along a cone rotating in axial flow[END_REF] (see [START_REF] Hussain | The centrifugal instability of the boundary-layer flow over a slender rotating cone in an enforced axial free-stream[END_REF] for comparisons of various slender cones with ψ < 40 • ). Importantly, while we have used the asymptotic results to provide an envelope for the right-hand branch of the numerical neutral stability curve, they are unable to predict the effect of varying axial flow on the critical Reynolds numbers. Nevertheless, the asymptotic analysis has proved invaluable in this study, as it reveals the correct length-scalings on which to model the counter-rotating vortex pairs, which characterise the centrifugal mode. Furthermore, by expanding the shifted basic flows about the location η = 0, we were able to confirm that the vortex activity of the most dangerous modes is located at the wall. Subsequently as s was varied, we tracked the location of vortex activity, observing that it departs slightly from the wall for s < 5. As a result, we have posed the hypothesis that stationary modes could dominate in the region s ≥ 5, but below this non-stationary (or 'travelling') modes may begin to grow. Such an observation requires further investigation, but would not be possible through solely conducting a numerical analysis. Hence, the importance of an asymptotic analysis is clear in revealing the underlying physical mechanisms at work, along with how they might interact.

In contrast, the OS numerical stability results complement the asymptotics in confirming the existence of the neutral stability curve for the centrifugal mode. Furthermore, we observe a reduction in the critical Reynolds number Re c as well as an increase in the critical amplification rate α 1,c with increasing s (see [START_REF] Hussain | The centrifugal instability of the boundary-layer flow over a slender rotating cone in an enforced axial free-stream[END_REF] for various slender cones with ψ < 40 • ).

Hence, larger values of s are destabilising, suggesting that the centrifugal-instability mode is physically the most dangerous mechanism, despite alternatives being present, including the crossflow and Tollmien-Schlichting instabilities.

Ultimately, we propose a condition of 'optimal' stability existing around s = 1, where the competing effects of the rotational and streamwise flow components balance. For s < 1, the physical problem changes from a 'quickly' rotating cone (the parameter range considered in this study) to a 'slowly' rotating cone. In this regime, the physical effect of the oncoming axial flow strengthens, thereby promoting the streamwise Tollmien-Schlichting instability, which begins to dominate over the centrifugal mode. This optimal stability criterion of s = 1 appears to exist not only for the centrfiugal instability studied here, but also for the crossflow instability studied by [START_REF] Hussain | Stability and transition of three-dimensional rotating boundary layers[END_REF] for a cone with ψ = 70 • rotating in an oncoming axial flow (see figure 5). The flow for both inviscid type I and viscous type II modes stabilises for increasing s, with s = 1 corresponding to the most stable flow setup in the type II case. Hence, for both flow instabilities on the rotating cone the findings predict an optimal stability criterion, which appears to be close to s = 1, where the effects of the two competing instabilities balance. This conclusion has interesting implications for the design of spinning projectiles, for example in military and defence applications. Here, the streamwise component is often large due to the projected velocity of the missile. For example, projectile applications that involve high rotation rates, such as spinning bullets and spinning missiles, can spin up to 3000 • per second and higher. In such case, it is important to design a missile that spins at sufficiently high rotation rate in order to promote the centrifugal mode and obtain a suitable balance between the competing instabilities. In fact, the primary instability can break-down to a secondary instability, which has been observed, for example, by [START_REF] Kohama | Flow structures formed by axisymmetric spinning bodies[END_REF] in the formation of 'horseshoe-like' vortices. Essentially, the aim in such spinning body applications is to reduce the parameter-scope for transition-to-turbulence within the flow. Therefore, influencing the primary and, potentially, the secondary instability, over a longer streamwise distance along the spinning body may achieve a delay in turbulent-transition, which consequently leads to more accurate targeting and projectile control properties. It should be noted at this point that the current study neglects the effects of compressibility, which would play a significant role in accurately modelling such high-speed applications.

Figure 1 .

 1 Figure 1. Diagram of the spiral vortex instability on a rotating cone placed in an oncoming axial flow, showing coordinates in the x-and y-logarithmic spiral directions, as well as the corresponding vortex wavenumber, γ, and vortex waveangle, φ.

Figure 2 .

 2 Figure 2. Velocity profiles Û (s, η 1 ) and V (s, η 1 ) in the xand y-directions, respectively, at ψ = 15 • for s = 1.5, 2, 3, 4, 5, 10, 16 and φ(s) = 30.2 • , 22.5 • , 13.6 • , 6 • , 0 • , 0 • , 0 • (in the directions of the arrows).

Figure 3 .

 3 Figure 3. Asymptotic scaled Taylor number T as a function of non-dimensional vortex wavenumber -1 for ψ = 15 • , s = 1.5, 2, 3, 4, 5, 10, 16 and φ(s) = 30.2 • , 22.5 • , 13.6 • , 6 • , 0 • , 0 • , 0 • .

Figure 4 .

 4 Figure 4. A comparison between the experimental observations of [3] (dot-dashed line, ) and the current theoretical predictions (solid line, •) of the vortex orientation angle, φ, at the onset of instability. The diagram illustrates φ(s) reduces with increased rotational flow parameter, s, to a limiting value of φ = 0 • at s = 5.

Figure 5 .

 5 Figure 5. Plot of asymptotic neutral wavenumber predictions against rotational Reynolds number for inviscid type I and viscous type II modes for ψ = 70 • , s = 0.1 -0.6 (type I), s = 0.1 -10.0 (type II). Increasing s shifts the curves as shown.

ACKNOWLEDGMENTS

ZH wishes to acknowledge support from Manchester Metropolitan University.

Competing instabilities of rotating boundary-layer flows in an axial free-stream


B Alveroglu

A Segalini

S J Garrett

S Y Mpo

O Tatin

An energy analysis of convective instabilities of the Bödewadt and Ekman boundary layers over rough surfaces

Keywords:

An energy balance equation for the three-dimensional Bödewadt and Ekman layers of the so called "BEK family" of rotating boundary-layer flows is derived. A Chebyshev discretisation method is used to solve the equations and investigate the effect of surface roughness on the physical mechanisms of transition. All roughness types lead to a stabilization of the Type I (cross-flow) instability mode for both flows, with the exception of azimuthally-anisotropic roughness (radial grooves) within the Bödewadt layer which is destabilising. In the case of the viscous Type II instability mode, the results predict a destabilisation effect of radially-anisotropic roughness (concentric grooves) on both flows, whereas both azimuthally-anisotropic roughness and isotropic roughness have a stabilisation effect. The results presented here confirm the results of our prior linear stability analyses.

INTRODUCTION

Developing passive drag-reduction techniques has been a priority for many years. It was a classical belief that smooth surfaces are likely to generate less skin friction than rough surfaces until the recent studies have shown the opposite for some cases [START_REF] Sirovich | Turbulent drag reduction by passive mechanisms[END_REF][START_REF] Carpenter | The right sort of roughness[END_REF][START_REF] Choi | Fluid dynamics: The rough with the smooth[END_REF]. Indeed, the right sort of roughness on surfaces could lead to drag-reducing effects [START_REF] Carpenter | The right sort of roughness[END_REF][START_REF] Choi | Fluid dynamics: The rough with the smooth[END_REF][START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF].

A promising strategy for drag reduction is to delay the transition of the laminar boundary-layer flow with the help of surface roughness. A recent study by Cooper et al. [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF] suggests that this is possible in the case of the von Kármán [START_REF] Von Kármán | Über laminare und turbulente reibung[END_REF] boundary-layer flow established over a disk rotating in otherwise still fluid. Our previous study [START_REF] Alveroglu | The effect of surface roughness on the convective instability of the bek family of boundary-layer flows[END_REF] established the same drag-reducing effects for the broader BEK class of flows that includes the Bödewadt [START_REF] Bödewadt | Die drehströmung über festem grunde[END_REF], Ekman [START_REF] Ekman | On the influence of the Earth's rotation on ocean currents[END_REF] and von Kármán [START_REF] Von Kármán | Über laminare und turbulente reibung[END_REF] boundary-layer flows. In this current study we perform an energy analysis of the Ekman and the Bödewadt flows to verify those results.

The BEK class of flows occur in turbo-machinery and rotorstator devices, such as mixers, and their stability characteristics are therefore of practical importance in engineering design applications. Many studies have been performed to reveal the similarities in the stability characteristics across the entire BEK class (see, for example, the theoretical and experimental studies discussed in papers [START_REF] Savas | Stability of Bödewadt flow[END_REF][START_REF] Lopez | Flow between a stationary and a rotating disk shrouded by a co-rotating cylinder[END_REF][START_REF] Lingwood | Experimental study of the stability of the Bödewadt layer[END_REF][START_REF] Schouveiler | Instabilities of the flow between a rotating and a stationary disk[END_REF][START_REF] Cros | Transition to turbulence of the Batchelor flow in a rotor/stator device[END_REF][START_REF] Culverhouse | The hydrodynamic stability of crossflow vortices in the Bödewadt boundary layer[END_REF] for the Bödewadt layer; papers [START_REF] Tatro | Experiments on Ekman layer instability[END_REF][START_REF] Faller | An experimental study of the instability of the laminar Ekman boundary layer[END_REF] for the Ekman layer; and papers [START_REF] Lingwood | Absolute instability of the ekman layer and related rotating flows[END_REF][START_REF] Fedorov | Transitional flow conditions on a rotating disk[END_REF][START_REF] Corke | Stationary travelling crossflow mode interactions on a rotating disk[END_REF][START_REF] Lingwood | Instabilities of the von Kármán boundary layer[END_REF] for the von Kármán flow). These studies have shown the existence of two main flow instability mechanisms, commonly referred to as the Type I and the Type II modes. The Type I mode is inviscid in nature and results from an inflectional cross-flow component of the steady flow; the Type II mode is viscous in nature and results from streamwise curvature and Coriolis effects. The instability nature of a particular flow in the BEK system is determined by the dominant one of those instability modes.

Two distinct theoretical models for the steady boundarylayer flow over rough rotating disks have been established in the literature. They were developed by Miklavčič & Wang [START_REF] Miklavčič | The flow due to a rough rotating disk[END_REF] and Yoon et al. [START_REF] Yoon | Flow and heat transfer over a rotating disk with surface roughness[END_REF], henceforth referred to as the MW and YHP models, respectively. The MW model uses a partial-slip condition on the disk surface and is capable of modelling independent levels of roughness in both the radial and azimuthal directions. The MW model can therefore model all variations of anisotropic roughness (radial and azimuthal) and also isotropic roughness. In contrast, however, the YHP approach imposes a particular surface geometry on the lower disk and can model only radial anisotropic roughness. The limitations of the YHP model motivate the use of the MW model alone in this study and our previous linear stability analysis [START_REF] Alveroglu | The effect of surface roughness on the convective instability of the bek family of boundary-layer flows[END_REF]. This paper proceeds as follows: The mean-flow equations and calculations of the steady boundary-layer flows of the BEK system are summarized in Section 1. We then derive the relevant equations from the perturbed system and conduct energy analyses of the Ekman and the Bödewadt flows for particular disturbance eigenmodes in Section 2. Conclusions are then drawn in Section 3.

THE STEADY MEAN FLOWS

The steady-flow formulation of the BEK system for smooth surfaces is constructed in the related papers of Lingwood [START_REF] Lingwood | Absolute instability of the ekman layer and related rotating flows[END_REF] and Lingwood & Garrett [START_REF] Lingwood | The effects of surface mass flux on the instability of the bek system of rotating boundary-layer flows[END_REF]. Full details of the scalings can be found there. The formulation for rough surfaces under the MW model is detailed in our previous study [START_REF] Alveroglu | The effect of surface roughness on the convective instability of the bek family of boundary-layer flows[END_REF] and we discuss it only briefly here.

The Bödewadt and the Ekman boundary-layers are two distinct cases of a general class of fully three-dimensional boundary layers within the BEK family. Within this family flows are distinguished by a differential rotation rate between a lower disk and the fluid above generated by an additional rotating disk in the far-field. The system is assumed to rotate within an otherwise still, incompressible fluid. Both disks are assumed to be of infinite extent and rotate around a common axis with angular velocities Ω * D (lower disk) and Ω * F (far field). The Ekman flow is the case that Ω * D ≈ Ω * F and the Bödewadt flow has Ω * D = 0 and Ω * F 0. Although not considered here, the von Kármán flow has Ω * D 0 and Ω * F = 0. Using a rotating frame leads to the appearance of Coriolis terms and centrifugal effects in the governing equations. Therefore, the flows in the system can be classified by the Rossby number that is a constant of the flow relating to the ratio of Coriolis forces to inertial forces. It can be formulated as

Ro = Ω * F -Ω * D Ω * ,
where the system rotation rate Ω * defined as

Ω * = Ω * F 2 -Ro + Ω * D 2 -Ro = Ω * F + Ω * D 4 + Ω * F + Ω * D 4 2 + (Ω * F -Ω * D ) 2 2 1/2 .
In this study, we focus only on the Bödewadt (Ro = 1) and the Ekman (Ro = 0) boundary layers in the BEK system.

A cylindrical polar coordinate system that rotates with the lower disk is used in the formulation of the governing equations for the steady flow. The MW model is used to model the steady-boundary layer flow over a rough lower disk. The approach modifies the no-slip condition as stated below at the disk surface, z = 0, whereas the far-field boundary conditions at the upper edge of the boundary layer are identical to the smooth-disk formulation.

U (0) = λU (0), V (0) = ηV (0) and W (0) = 0, U = 0, V = 1, as z -→ ∞. (1)
Here, U, V, and W are the radial, azimuthal and vertical mean velocity profiles and the two parameters η and λ are empirical measures of roughness in radial and azimuthal directions, respectively.

The computed mean-flow profiles for different levels of surface roughness are presented in our previous study [START_REF] Alveroglu | The effect of surface roughness on the convective instability of the bek family of boundary-layer flows[END_REF].

A detailed discussion about the effects of roughness on the profiles can also be found in that paper. In this study, we proceed to perform an energy analysis on the computed meanflow profiles.

ENERGY ANALYSIS

A local linear stability analysis of the BEK system with surface roughness was conducted by Alveroglu et al. [START_REF] Alveroglu | The effect of surface roughness on the convective instability of the bek family of boundary-layer flows[END_REF]. Following the work of Cooper & Carpenter [START_REF] Cooper | The stability of rotatingdisc boundary-layer flow over a compliant wall. part 1. type I and II instabilities[END_REF] and Cooper et al. [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF], an integral energy equation for three-dimensional disturbances ( û, v, ŵ, p) to the undisturbed steady mean-flow profiles (U, V, W ) is derived in cylindrical polar coordinate system in rotating frame to analyse the underlying physical mechanisms behind Alveroglu et al.'s conclusions. The derivation of the governing energy equations begins by multiplying the linearized momentum equations by û, v and ŵ, respectively. Summing the resulting expressions leads to the kinetic-energy equation for the disturbances,

∂ ∂t + U ∂ ∂r + V Re ∂ ∂θ - RoW Re ∂ ∂z K = -û ŵ ∂U ∂z -v ŵ ∂V ∂z + Ro Re ŵ2 ∂W ∂z + RoU Re û2 (2)
+ RoU Re v2 - ∂( û p) ∂r + 1 Re ∂( v p) ∂θ + ∂( ŵ p) ∂z - Ro Re û p + ∂( ûj σ i j ) ∂ x i -σ i j ∂ ûj ∂ x i
Here K = 1 2 ( û2 + v2 + ŵ2 ) and σ i j are anti-symmetric viscous stress terms. Furthermore, Re is a local Reynolds number and it is written in terms of local radial position r a as Re = r a Ro. All O(1/r) viscous terms are now neglected to ensure consistency with the linearized governing stability equations.

We now average the perturbations over a single time period and azimuthal mode and integrate across the boundary layer. The steady, rotationally-symmetric nature of the energy is enforced by neglecting both t and θ, leading to

∞ 0          U ∂K ∂r a + ∂( û p) ∂r b - ∂ ∂r ( ûσ 11 + vσ 12 + ŵσ 13 ) c          dz = ∞ 0 -û ŵ ∂U ∂z + -v ŵ ∂V ∂z + ŵ2 Ro Re ∂W ∂z dz I - ∞ 0 σ i j ∂ ûj ∂ x i dz II + ∞ 0 Ro Re û p dz + ( ŵ p) w III (3)
-

ûσ 31 + vσ 32 + ŵσ 33 w IV + ∞ 0 Ro Re ∂K ∂z W dz + ∞ 0 Ro Re û2 Udz + ∞ 0 Ro Re v2 U dz V .
Overbars in the equations denote a period-averaged quantity, such that û v = û v * + û * v (where * indicates a complex conjugate). When Ro = -1, this equation is consistent with the energy equation for the von Kármán boundary layer [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF][START_REF] Cooper | The stability of rotatingdisc boundary-layer flow over a compliant wall. part 1. type I and II instabilities[END_REF]. We now proceed with a normal-mode analysis and impose perturbations of the form

( û, v, ŵ, p) = (u(z), v(z), w(z), p(z))e i (αr +β Reθ-ωt ) . (4)
Here the radial wavenumber α = α r + iα i ∈ C and the azimuthal wavenumber and frequency β, ω ∈ R. It is assumed that β O(1) and the number of spiral vortices around the disk surface is given by n = βRe. Furthermore, the orientation angle of the vortices with respect to a circle centred on the axis of rotation is defined as = tan -1 ( β/α r ). In all that follows we set the frequency ω = 0 as we are interested in stationary vortices rotating with the rough disk.

As discussed elsewhere [START_REF] Cooper | The stability of rotatingdisc boundary-layer flow over a compliant wall. part 1. type I and II instabilities[END_REF], the energy balance can be performed for any eigenmode of the perturbation equations and is derived after inserting (4) into (3) as

-2α i = (P 1 + P 2 + P 3 ) I + D II + (PW 1 + PW 2 ) III (S 1 + S 2 + S 3 ) IV + (G 1 + G 2 + G 3 ) V , (5)
where the terms can be interpreted physically as (I) the Reynolds stress energy production, (II) the viscous dissipation energy removal term, (III) pressure work terms, (IV) contributions from work done on the wall by viscous stresses, (V) terms arising from the streamline curvature effects and the three dimensionality of the mean flow.

The positive terms contribute to energy production and the negative ones remove energy from the system. The mode is amplified (α i < 0) when energy production outweighs the energy dissipation in the system. Using this formulation it is possible to discuss the effect of roughness on the instabilities of the Type I and Type II modes by calculating the total energy of the system which is the sum of all terms in Eq. ( 5). If the total energy increases as the roughness parameters are increased, then the modes become unstable. In contrast, a reduced total energy would indicate a stabilization effect. Underlying our approach to the calculations is a stability code that uses a Chebyshev polynomial discretization method in the wall-normal direction to obtain each term in Eq. ( 5) for various values of the roughness parameters. Consistent with the recent literature [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF][START_REF] Alveroglu | The effect of surface roughness on the convective instability of the bek family of boundary-layer flows[END_REF][START_REF] Miklavčič | The flow due to a rough rotating disk[END_REF], we use the parameters η and λ to denote empirical levels of roughness in the radial and azimuthal directions, respectively. The case η > 0 and λ = 0 corresponds to radial anisotropic roughness (concentric grooves); λ > 0 and η = 0 corresponds to azimuthal anisotropic roughness (radial grooves); and λ = η > 0 corresponds to isotropic roughness.

Results for the Ekman flow

The energy balance calculation for the Ekman boundary layer is carried out at the location of maximum amplifications of the both Type I and Type II modes at Re = R c + 25. Here R c is the critical Reynolds number for the onset of the Type I mode of instability for the particular rough surface being considered. Results for various levels of roughness are compared to smooth disk case in Figures 123.

Figure 1 shows the energy balance calculation for azimuthal anisotropic roughness (radial grooves). It reveals a stabilization effect on the Type I mode that arises mainly from the large reduction in energy production term P 2 . Note that the energy dissipation of the system, P 1 + D, is largely invariant, leading to a reduction in the total energy of the system with increased roughness of this type. Figure 2 shows a similar stabilizing effect of isotropic roughness on the Type I mode. Figure 3(a) demonstrates the major stabilizing effect of radial anisotropic roughness (concentric grooves) on the Type I mode. This is due to a reduction in the energy production term P 2 and largely invariant viscous stress dissipation D. As a result, the total energy of the system decreases. This confirms the stabilization effect of concentric grooves on the Type I mode. In contrast, however, Figure 3(b) presents a major destabilizing effect on the Type II mode. This arises from an increase in P 2 terms at the points of maximum amplification, together with a slight increase in the viscous dissipation term D. Collectively these bring about a increase in the total energy of the system that reveals the destabilization effect of concentric grooves on the Type II mode.

Supportive evidence for the above trends can be obtained from the form of the perturbation velocity profiles, that is, the eigenfunctions. The dominant eigenfunction is the axial perturbation velocity, w, which contributes to the dominant energy production term P 2 . Figure 4 shows the magnitude of the w-profile for a disk with radial anisotropic roughness. For the Type I mode shown in Figure 4(a), the form of the disturbance profile is largely unchanged as roughness is increased. The reduction in P 2 , however, results from the slight reduction in the amplitude of the profile as roughness is increased. Conversely, Figure 4(b) reveals that the increase of the P 2 term for the Type II mode is due to an increasing amplitude of this profile with roughness.

Results for the Bödewadt boundary layer

As predicted by Lingwood & Garrett [START_REF] Lingwood | The effects of surface mass flux on the instability of the bek system of rotating boundary-layer flows[END_REF], the very early onset of absolute instability in the Bödewadt flow leads to a "branch exchange" and it is impossible to find the location of maximum amplification of either convective mode. We therefore instead proceed with the energy balance calculations for the Bödewadt boundary layer at fixed Re = 400 with n = 85 for the Type I mode and n = 12 for the Type II mode. These values of n are chosen so as to be outside of the absolutely unstable region at Re = 400.

Figure 5 shows energy terms of Type I and Type II instabilities for n = 85 and n = 12, respectively, subject to increasing radial anisotropic roughness. Figure 5(a) reveals that the major stabilizing effect on the Type I mode is from reductions in the energy production by Reynolds stress and the effect of streamline curvature in the boundary layer. The total energy of the system decreases as roughness (η) is increased. This is interpreted as evidence of the stabilizing effect of concentric grooves on the Type I mode. Conversely, Figure 5(b) shows a destabilizing effect on the Type II mode that comes from a significant increase in the energy production term P 2 and subsequent increase in total energy. Figure 6 demonstrates the destabilizing effect of azimuthal isotropic roughness (radial grooves) on the Type I mode. In particular, increased roughness leads to growth in the energy production terms P 2 and G 1 . Even though the viscous dissipation is also increased, the net effect is an increase in total energy and a destabilization of the flow.

As shown in Figure 7, increasing the level of isotropic roughness promotes viscous dissipation from the system that causes a net reduction in total energy. This results in a stabilizing effect. In this case, the net energy production in the system is largely invariant. Furthermore, Figure 8 shows the magnitude of the wprofile for the case of radial anisotropic roughness. For the Type I mode in 8(a), a dramatic reduction in the amplitude of the profile as roughness is increased causes the large reduction in the P 2 energy production term. That results in a stabilization effect. In contrast, Figure 8(b) reveals that the destabilization of the Type II mode is due to an increase in the amplitude of this profile with roughness. This leads to an increased P 2 energy production term.

CONCLUSION

This study has been designed to establish the underlying physical mechanisms behind the effects of rough surfaces on the BEK system. The steady flow profiles over the rough disks are modelled using the partial-slip approach of Miklavčič & Wang [START_REF] Miklavčič | The flow due to a rough rotating disk[END_REF] without modification. The resulting energy balances are qualitatively consistent with the neutral curve calculations of Alveroglu et al. [START_REF] Alveroglu | The effect of surface roughness on the convective instability of the bek family of boundary-layer flows[END_REF].

Our study has revealed that radially-anisotropic roughness (concentric grooves) and isotropic roughness acts to reduce energy production of the Type I mode for both the Ekman and Bödewadt boundary layers. This is a stabilization effect on the Type I mode. Conversely, azimuthally-anisotropic roughness (radial grooves) is found to have the opposite destabilization effect on the the Type I mode. The energy balance calculation for the Type II mode demonstrates the destabilizing effect of radially-anisotropic roughness for each of the Ekman and Bödewadt boundary layers.

The conclusions arising from this study are easy to interpret: Surface roughness can be used as an effective passive flow-control mechanism for engineering flows of the rotorstator type modelled by the BEK system, and also other cross-flow dominant flows.

Figure 1 .

 1 Figure 1. Energy balance of the Ekman boundary layer at Re = R c + 25 for azimuthal anisotropic roughness (radial grooves). Type I instability mode.

Figure 3 .

 3 Figure 3. Energy balance of the Ekman boundary layer at Re = R c + 25 for radial anisotropic roughness (concentric grooves). (a) Type I instability mode, (b) Type II instability mode.

Figure 4 .

 4 Figure 4. Profiles of the axial perturbation velocity w of the Ekman layer at Re = R c + 25 with radial anisotropic roughness (concentric grooves). (a) Type I instability mode, (b) Type II instability mode.

1 Figure 2 .

 12 Figure 2. Energy balance of the Ekman boundary layer at Re = R c + 25 for isotropic roughness. Type I instability mode.

Figure 6 .

 6 Figure 6. Energy balance of the Bödewadt boundary layer at Re = 400 for azimuthal anisotropic roughness (radial grooves). Type I instability mode with n = 85.

Figure 5 .

 5 Figure 5. Energy balance of the Bödewadt boundary layer at Re = 400 for radial anisotropic roughness (concentric grooves). (a) Type I instability mode with n = 85, (b) Type II instability mode with n = 12.

Figure 8 .

 8 Figure 8. Profiles of the axial perturbation velocity w of the Bödewadt layer at Re = 400 with radial anisotropic roughness (concentrically grooved). (a) Type I instability mode, (b) Type II instability mode.

Figure 7 .

 7 Figure 7. Energy balance of the Bödewadt boundary layer at Re = 400 for isotropic roughness. Type I instability mode with n = 85.

ACKNOWLEDGMENTS

This research used the ALICE High Performance Computing Facility at the University of Leicester. BA wishes to acknowledge financial support from Ministry of National Education, Republic of Turkey. SJG is supported by a Senior Research Fellowship of the Royal Academy of Engineering, funded by the Leverhulme Trust.