Martin Böhle

Annika Fleder

Matthias Mohr
email: matthias.mohr@mv.uni-kl.de

S Y Mpo

O Tatin

Study of the Losses in Fluid Machinery with the Help of Entropy

In the current design process of turbomachinery the use of CFD methods is state of the art. With CFD calculations we are able to predict the performance characteristics of the machines. In most of the cases the machines are rated based on their performance alone. A detailed evaluation of the losses involved often is not conducted, mainly because there is no established method to analyze the location and mechanisms of loss production. This leads to the fact that the loss mechanisms for some turbo machines are not known in detail, although they are very important to completely understand the operation principle of the machines and have a huge influence on their performance. If the location and quantity of losses were known, it would be possible to perform a selective optimization of the machine. The question is just how to calculate and localize the losses and evaluate them with CFD. In his thesis Kock introduces a method, to identify losses in CFD. His method is based on the entropy production rates in turbulent flow fields. In this paper the method is applied on hydraulic fluid machinery and two different machine types are discussed in detail. Those are a side channel pump and an inducer. The locations of the losses are detected and the size of losses will be evaluated by applying the method of Kock.

INTRODUCTION

In all fluid machinery there exist different kind of losses. We have to differentiate between internal losses and external losses. The external losses are, for example, losses in bearings, gears and shaft seals. Those losses are not included in the CFD calculation whereas the internal losses are considered by CFD. Among the latter are shock losses and friction losses. The different types of internal losses can be summarized in dissipation. The numerical simulation offers the opportunity to calculate the flow in the machines spatially and temporally dependent. Thus the flow in the machines can be observed three dimensionally at an arbitrary point in time. The calculated quantities are the fields of pressure, temperature and velocity. Derived quantities like pressure coefficient and efficiency can be calculated directly during post processing. To identify the different internal losses Herwig and Kock [START_REF] Kock | Local entropy production in turbulent shear flows: a high-reynolds number model with wall functions[END_REF][START_REF] Kock | Entropy production calculation for turbulent shear flows and their implementation in cfd codes[END_REF][START_REF] Herwig | Direct and indirect methods of calculating entropy generation rates in turbulent convective heat transfer problems[END_REF] proposed a direct method based on the evaluation of entropy production in turbulent flow fields. The advantage of the method is that it can be applied to a CFD simulation in a pure postprocessing step. The content of the present paper is the application of this method to two examples of hydraulic turbomachinery, a side channel pump and an inducer. For those two types of machines numerical investigations are conducted. To assess the quality of the numerical results they are compared to available experimental data. Eventually a detailed analysis of the losses based on the numerical data is performed and the influence of the turbulence model is examined.

METHODS

With the help of CFD it is possible to calculate the timedependent flow field in all kind of turbomachinery. The governing equations that are solved to describe the flow are the equations of conservation of mass, conservation of momentum and conservation of energy. The basis of the energy equation calculated in CFD is the first law of thermodynamics. Herwig and Kock use the production of entropy as a direct parameter for the efficiency of a technical system [START_REF] Herwig | Direct and indirect methods of calculating entropy generation rates in turbulent convective heat transfer problems[END_REF]. In the transport equation of entropy heat transfer and pressure drop are linked to one quantity. They describe two methods to calculate the entropy production in the context of a CFD simulation process. In this paper the so called direct method will be considered. One main feature of the direct method is that the entropy is treated as a postprocessing quantity. In the cases of hydraulic fluid machinery, the fluid is commonly regarded as incompressible and the heat transfer is negligibly. Therefore the energy equation is not calculated.

The specific entropy s is a state variable that, according to the second law of thermodynamics, grows in every real and thus irreversible technical process. Its transport or balance equation is given in the literature (e.g. [START_REF] Spurk | Fluid Mechanics[END_REF] or [START_REF] Schlichting | Grenzschichttheorie[END_REF]) and can be written as follows:

ρ δs δt + u • δs δx + v • δs δy + w • δs δz = -div q T + Φ T + Φ θ T 2 (1)
The first term on the right-hand side of equation ( 1) describes the reversible heat transfer div q T = divλ T grad T . This term is not important for the calculation of the losses. The function Φ

T describes the irreversible part produced by friction, and Φ θ T 2 = -q T 2 grad T the irreversible part produced by heat transfer. This means that in those two terms describing the irreversible part of entropy production all losses that occur in a machine by friction and heat transfer are included. As mentioned before, the heat transfer is negligible for the pumps under examination.

The terms Φ T andq T 2 grad T are always greater than zero because in an adiabatic system with q = 0 the entropy cannot decrease. In industrial applications of CFD analysis of turbomachinery usually the time-averaged Navier-Stokes (RANS) equations are solved. Denoting the time averaged quantities by an overbar, the time-averaged version of the transport equation for entropy, as shown in [START_REF] Kock | Local entropy production in turbulent shear flows: a high-reynolds number model with wall functions[END_REF], can be written in the following form:

ρ δ s δt + ū • δ s δx + v • δ s δy + w • δ s δz = -div q T -ρ δu s δx + δv s δy + δw s δz + Φ T + Φ θ T 2 (2)
Important for our studies is the term of entropy production caused by dissipation Φ T . This term is composed as follows:

Φ T = Ṡ D pr o + ṠD pr o (3)
The first term on the right-hand side of equation ( 3) stands for the entropy production by viscous or direct dissipation. It can be calculated directly in the CFD postprocessing procedure by using the following equation:

Ṡ D pr o = µ T 2 •       ∂ ū ∂ x 2 + ∂ v ∂ y 2 + ∂ w ∂z 2      + ∂ v ∂ x + ∂ ū ∂ y 2 + ∂ w ∂ y + ∂ v ∂z 2 + ∂ ū ∂z + ∂ w ∂ x 2 (4)
The second term on the right-hand side of equation ( 3) describes the entropy production by turbulent dissipation which is also called indirect dissipation. It can be expressed as follows:

ṠD pr o = µ T 2 •       ∂u ∂ x 2 + ∂v ∂ y 2 + ∂w ∂z 2      + ∂v ∂ x + ∂u ∂ y 2 + ∂w ∂ y + ∂v ∂z 2 + ∂u ∂z + ∂w ∂ x 2 (5)
This expression can't be computed directly, because the fluctuating quantities are not available in the CFD solution based on the RANS equations. To circumvent this problem the following model equation is proposed by Kock [START_REF] Kock | Local entropy production in turbulent shear flows: a high-reynolds number model with wall functions[END_REF][START_REF] Kock | Entropy production calculation for turbulent shear flows and their implementation in cfd codes[END_REF]:

ṠD pr o = ρ T (6)
It uses the dissipation rate of turbulent kinetic energy which is provided by the turbulence model and is readily available during postprocessing.

We have to consider that within this approach errors are introduced by the turbulence model itself which -as its name implies -is just a model and by the approximation involved in the formulation of model equation [START_REF] Kock | Bestimmung der lokalen Entropieproduktion in turbulenten Strömungen und deren Nutzung zur Bewertung konvektiver Transportprozesse[END_REF].

Following Kock [6] the entropy production rates caused by viscous and turbulent effects (equations ( 4)-( 6)) are connected to the viscous and turbulent dissipation respectively through the relations:

Φ D = T Ṡ D pr o (7) Φ D = T ṠD pr o (8)
Finally the dissipation portions Φ D and Φ D can be used to calculate the corresponding power losses by integrating over the volume of interest:

P D = V Φ D dV = V T Ṡ D pr o dV (9)
P D = V Φ D dV = V T ṠD pr o dV (10)
The calculation of the turbulent dissipation depends on the turbulence model. In the case of the side channel pump three different turbulence models, the k-model, the k-ω and the k-ω-SST model are used and their influence on the turbulent dissipation is shown. In the case of the inducer only the k-ω-SST model is utilized. For the remainder of the current paper we refer to the three applied models as k-, k-ω and SST respectively.

As will be seen later the viscous dissipation is small compared to the turbulent dissipation in the case of the examined machines. Hence to investigate the loss mechanisms and locations of loss production it is helpful to visualize the spatial distribution of turbulent dissipation inside of the machine. For this purpose we introduce the dimensionless loss coefficient ζ D . The loss coefficient is defined in the following equation:

ζ D = Φ D • d a ω 2 • ρ • V (11)
The loss coefficient depends on the volume flow rate V , the outer diameter of the impeller d a , the angular velocity ω and the turbulent dissipation Φ D . It can therefore be used to compare different operation points and machines.

APPLICATION 1: SIDE CHANNEL PUMP

Numerical set-up

In Figure 1 the simulation model of the side channel pump is depicted. It consists of three computational domains. Those are the inlet tube colored in grey, the impeller colored in red and the side channel including the outlet tube colored in grey.

The domains are discretized with hexahedral multiblock grids that were created using ANSYS ICEM CFD. The computations are performed with the commercial CFD code ANSYS CFX 14.5. Between the inlet tube and the impeller and between the impeller and the side channel transient rotor-stator interfaces are used (Figure 1). At the inflow the average total pressure is prescribed. The outlet boundary condition is realized as a velocity opening to prevent walls in the outlet caused by backflow of the fluid, which is generated by the strong swirl in the outlet tube of the pump. The default algorithm of velocitypressure coupling, implemented in CFX, is applied. We use a time step of 0.0001 s to simulate the pump at a speed of 1500 min -1 . For each operation point at least two full impeller turns are simulated to ensure periodic converged solutions. Pressure, momentum and velocities have to be time averaged over the pulsation length due to the characteristic pulsations of side channel pumps. The pulsation length is calculated as 360/z. To evaluate the influence of the turbulence model on the losses, the numerical calculations are carried out with three different turbulence models as already mentioned above. The influence of the turbulence model on the efficiency is larger as can be seen in figure 3. The best efficiency is calculated with the k-ω-SST model. Compared to the lowest efficiency that is calculated with the k-model there is a difference of about 2 %. In the measured efficiency mechanical losses of the bearings and the sealing are included. Because these external losses cannot be considered in the CFD simulations the measured efficiency curve lies below the numerically calculated points.

Dissipation Power and loss analysis

The occurring losses will be examined with the help of the dissipation power. As already discussed we have to distinguish between the power losses due to viscous dissipation (P D ) and the losses due to turbulent dissipation (P D ). Figure 4 shows the dissipation power loss, averaged over one pressure pulsation, for the side channel domain (index SK) and the impeller domain (index LR) measured in watt. The viscous dissipation power is a lot smaller in all operation points than the turbulent dissipation power. Their trend is the same. The analytically undetermined losses (inlet shock of the main and the circulating flow, friction and gap losses) decrease with increasing volume flow rate. Shock and gap losses are included in the dissipation power of the impeller, that decreases with increasing volume flow. The viscous and turbulent dissipation power in the side channel increase slightly. The losses in the side channel are mainly friction losses. Those increase with increasing volume flow rate and and the associated higher flow velocities in the side channel. The losses in the inlet tube are small and are neglected in the figure. The isosurface is especially developed in the outlet and in the inlet region (no. 3 in figure 5). In those regions the losses are large. They are induced by eddies, as can be seen from the vectors and streamlines on the section planes. In the left view the plane is located at an outer radius. At this radius the fluid exits from the impeller into the side channel with small values of ζ D . Only at the blade surfaces and between the impeller and side channel walls friction losses can be clearly recognized (no. 2 in figure 5). At the middle and lower radii illustrated in the two subsequent views, the fluid enters into the impeller with notable shock losses at the blade edges. These are indicated with no. 1 in figure 5. Inside the impeller regions with higher losses are situated. Due to the large difference between flow angle and blade angle flow separation occurs inside the impeller resulting in the development of vortices.

Those are responsible for the increase of losses. Figure 6 depicts the turbulent dissipation power in the impeller, figure 7 shows the turbulent dissipation power in the side channel. For the impeller as well as for the side channel the k-model predicts the highest dissipation power. The smallest losses are obtained when using the k-ω-SST model. This behavior is consistent to the efficiency characteristics of the pump obtained with the three different turbulence models and shown in figure 3. In general, the losses in the impeller decrease with increasing volume flow rate and the losses in the side channel increase with increasing volume flow rate. This is valid for all volume flows and turbulence models except for the simulation at overload with the k-model. While the losses in the impeller coincide quite well in case of the simulations using the k-ω and the k-model, the losses in the side channel agree better for the k-ω and k-ω-SST model. Using the kmodel results by far in the highest losses in this region. To Figure 10 shows the distribution of the loss coefficient on a plane in the middle of the gap between impeller and casing for all three turbulence models. The maximum loss coefficients found on this plane are substantially higher than in the middle of the impeller. Compared to figure 9 with a maximum loss coefficient of 30, values up to ζ D = 100 can be found here.

The distribution of the loss coefficient again is quite similar for all the considered cases. In the regions of high velocities and thus high friction the highest values of loss coefficient can be observed. Figure 8 depicts the loss coefficient in the 8). Higher losses occur in the annular gap near the interrupter. In this region higher velocities are caused by the large pressure gradient between suction-und discharge side of the pump (see no. 2 in figure 8).

From figure 11 it can be seen that independently of the used turbulence model there exist higher losses in the inlet and outlet regions while the overall distribution of the loss coefficient again looks roughly similar in all cases. In these regions the circulation of flow is not fully developed. Slightly increased losses can also be observed at the rotation axis of the circulating flow in the blade passages. The k-and k-ω model predict regions with somewhat higher losses compared to the result obtained using the k-ω-SST model.

At the wall of the side channel, shown in figure 12, we observe the largest differences between the turbulence models with respect to the calculated loss coefficient. Application of the k-model yields considerably higher losses compared to the other models. While for the k-ω and the k-ω-SST model large values of the loss coefficient are found only in the inlet region the k-model indicates high losses approximately along the first 45°of the side channel. In contrary to the other models the k-model predicts an additional region of higher losses near the outlet.

APPLICATION 2: INDUCER

Numerical set-up

As a second test case we consider a simple three-bladed helical inducer with axial hub shape under non-cavitating conditions. Steady state simulations of a periodical segment containing one blade passage of the inducer are performed using the commercial CFD code ANSYS CFX 15.0. Figure 13 shows the simulation setup consisting of three domains that represent the inlet pipe (blue), the inducer (grey) and the outlet pipe (blue). All domains are discretized with hexahedral elements only. To generate the grid for the inflow domain that includes the hub nose we used the software ANSYS ICEM The grids for the blade passage and the outlet domain were generated using ANSYS TurboGrid. To include the tip clearance losses in our numerical investigation we placed 25 cells in the gap between blade tip and shroud. The rotating domain containing the inducer is connected to the stationary inflow and outflow domains with frozen rotor interfaces. At the inlet boundary we prescribe the total pressure, at the outlet the mass flow is specified. We simulate five different volume flows at a speed of 3000 min -1 . As mentioned before we only use the k-ω-SST turbulence model for this case. Figures 14 and15 show the simulated data points in comparison to experimental results obtained at the test facilities of the TU Kaiserslautern. While a very good agreement exists for the pressure coefficient the numerically calculated efficiencies are throughout higher than the measured ones due to the mechanical losses that are not included in the CFD simulations. Especially under over-load conditions the efficiency is overpredicted by CFD. Figure 16 progression of the power losses as sums over all the simulation domains within the simulated operation range. While the power losses caused by viscous dissipation remain nearly constant a notable increase of the losses caused by turbulent dissipation occurs towards part-load. Like in the previous case of the side channel pump (see figure 4) the losses due to viscous dissipation are substantially smaller than the losses Figure 17 shows an isosurface for ζ D = 0.2 and contours of the loss coefficient on cross sections of the blade passage at three axial stations for the simulation at the BEP. It should be noted that the color levels of the contour plots are scaled logarithmically to simultaneously cover the full range of loss coefficient and emphasize the significant regions. The cut section on the left is located close to the leading edge, the cut section on the right at the trailing edge and the third approximately halfway between the the two of them. As expected the largest losses occur at the outer diameters of the inducer blade and in the tip clearance where flow velocities and friction reach the highest values. The greatest losses can be observed at the blade tip close to the leading edge where the typical backflow vortex (no. 1 in figure 17) is generated. As can be seen from the isosurface in figure 17 As already stated the loss coefficient increases with diameter and is especially high near the blade tip and along the shroud on the suction side of the blade. On the mid plane higher values of loss coefficient can be detected near the walls at the hub, blade and shroud that disturb the uniform distribution of loss coefficient in the core region. This indicates turbulent interaction of secondary flows with the core flow, causing an increase in losses. At the trailing edge position the core region containing small values of loss coefficient is further reduced and higher values can be found all over the hub and the blades.

DISCUSSION

In the first application, the side channel pump, we investigate the influence of the turbulence model on the predicted losses and their distributions. We can note that the occurrence of losses can be reasonably explained by all three turbulence models under examination. The higher losses in the impeller are due to the fact that the gap is included in this simulation domain and the complex flow pattern inside the impeller blade passages. Near walls the k-model predicts, in general, higher values of loss coefficient than the other two models. The spatial distribution of the losses and the tendency of dissipation power to increase and decrease with changed operation conditions can be physically justified and qualitatively measured. As already indicated the absolute values of the local losses depend highly on the assumptions made by the turbulence model used. Therefore no statements regarding the correctness of the quantitative losses can be made. However the method applied in this paper can be understood as a simple tool to understand loss production processes and locate the regions that are critical for the development and improvement of turbomachinery by using CFD simulations.

Our second application deals with the flow through an inducer. We could observe similar trends of the power losses within the investigated operation range. The power losses caused by viscous dissipation are comparable and equally small for the side channel pump as well as for the inducer whereas the losses due to turbulent dissipation are lower in the case of the latter. This smaller amount of turbulent dissipation explains the better (but still low) efficiency of the inducer compared to the efficiency of the side channel pump. By visualization of the losses with the help of the loss coefficient we could identify and localize the flow features that are responsible for a large portion of dissipation. Furthermore the development of the distribution of the loss coefficient inside the blade passages could be observed with this postprocessing method. The insights gained into the structure of the losses can help to better understand the mechanisms involved in loss production and improve the design of the machines.

Concluding we can state that a qualitative analysis of the losses in turbomachinery can be performed with the simple postprocessing procedure demonstrated in the current work. The ability to detect critical flow regions regarding the production of losses creates the opportunity to selectively optimize existing geometry. However it should be remarked that the quantitative balance of powers could not be shown with the current method for the cases under examination. This issue needs further investigation and is part of our current research work.

NOMENCLATURE Latin symbols

Symbol

Figure 1 . 2 . 2 Figure 2 .

 1222 Figure 1. Simulation set-up of the side channel pump

Figure 3 .

 3 Figure 3. Efficiency for the different turbulence models compared to experimental data

Figure 4 .

 4 Figure 4. Dissipation power for different volume flow rates for the impeller and the side channel

Figure 6 .Figure 7 .Figure 5 .

 675 Figure 6. Turbulent dissipation power in the impeller for the different turbulence models

Figure 8 .

 8 Figure 8. Loss coefficient and streamlines calculated with k-ω-SST turbulence model in the middle of the gap between impeller and side channel

Figure 9 .

 9 Figure 9. Loss coefficient in the middle of the impeller

Figure 10 .

 10 Figure 10. Loss coefficient in the gap between impeller and casing

Figure 11 .Figure 13 .

 1113 Figure 11. Loss coefficient in the plane of the exchange mass flow

Figure 14 .

 14 Figure 14. Measured and simulated pressure coefficient of the inducer

Figure 15 .Figure 16 .

 1516 Figure 15. Measured and simulated efficiency of the inducer

Figure 17 .

 17 Figure 17. Isosurface and contours of loss coefficient ζ D

ACKNOWLEDGMENTS

The simulations were conducted on the high performance cluster Elwetritsch at TU Kaiserslautern, which is part of the Alliance of High Performance Computing Rheinland-Pfalz (AHRP). We would like to thank for their kind support.


Alessio Abrassi

Alessandra Cuneo
email: alessandra.cuneo@edu.unige.it

Alberto Traverso

Off-design performance analysis of a micro-gas turbine under uncertainty

Keywords:

The analysis and design of complex energy systems is generally performed starting from a single operating condition and assuming a series of design parameters as fixed values. However, many of the variables on which the design is based are subject to uncertainty because they are not determinable with an adequate precision. For advanced energy system and other processes, the uncertainties associated with model input parameters can affect both the performance and the cost. Methods for system design under uncertainty thus become essential. As uncertainty is a broad concept, it is possible, and often useful, to approach it in several ways. One rather general approach, which is applied to a wide variety of problems, is to assign a probability distribution to the various uncertain input parameters of the model. Many studies have faced optimization problems, but almost no one has considered the uncertainty as a factor to be taken into account. An off-design static model of a micro-turbine has been studied: it is built on the configuration of a Turbec T100 actually installed at the University research laboratory. Stochastic analysis has been treated implementing the approximated method Response Sensitivity Analysis (RSA) based on Taylor series expansion. RSA methods have been used on this system model to estimate its main performance and economic parameters under the influence of uncertainties related to various operating parameters of the turbo-machinery.

INTRODUCTION

Almost all the existing design methods are deterministic since the presence of uncertainties in parameters and measurements are not taken into account. For realworld engineering systems, such simplification is generally impractical and this highly limits the applications of deterministic model updating methods. Therefore, stochastic updating methods involving uncertainty analysis are of great importance from a practical point of view. In a broad sense, uncertainties can be classified into two categories of aleatory (irreducible) and epistemic (reducible) uncertainties [START_REF] Fang | A stochastic model updating method for parameter variability quantification based on response surface models and Monte Carlo simulation[END_REF] and in many cases there is no strict distinction between these two categories. Aleatory uncertainty may result from geometric dimension variability due to manufacturing tolerances [START_REF] Kim | Stochastic Modeling and Uncertainty Analysis with Multi-objective Optimization Strategies for the Synthesis/Design and Operation/Control of a PEMFC Fuel Processing Subsystem[END_REF][START_REF] Kim | A hybrid multi-level optimization approach for the dynamic synthesis/design and operation/control under uncertainty of a fuel cell system[END_REF][START_REF] Heinze | A parametric model for probabilistic analysis of turbine blades considering real geometric effects[END_REF] or inherent variability of materials such as concrete, while epistemic uncertainty is caused by lack of knowledge (e.g. due to limited experimental data). The design of machinery, plants or any equipment connected to the world of energy, mechanical and civil engineering is often subjected to highly variable conditions; just think of how environmental condition can be aleatory for the working condition of energy plants, or how, during each design process, the uncertainties related to every parameter or chosen coefficients affect the final result. This is a very sensitive issue for the wide world of dynamic components, so does the need to refine methods able to take in account uncertainties of the studied system and to analyze their effects on monitored outputs. This work is part of a broader study that points to identify and develop one or more methods able to fulfill these tasks with more accuracy and speed as possible. In particular, this paper presents the study of a microturbine, the Turbec T100, at present installed at the research laboratory of the University of Genoa, located in Savona, under stochastic conditions, Figure 1.

METHODS

Due to different reasons, a system can have input data affected by uncertainty. This may be due to the inability to accurately estimate some of these parameters, or more simply it is related to a sudden change in operating conditions. The approach starts from the identification of the variables behavior, estimating their probability distribution. In most of cases, input data populations are organized in normal distributions or in distributions shaped very close to them. The next step was to find a method able to manage stochastic inputs. In open literature, there are many different possible approaches: Monte Carlo Simulation MCS, Response Sensitivity Analysis RSA, Fast Probability Integration FPI and Analysis of Variance ANOVA [START_REF] Kim | Stochastic Modeling and Uncertainty Analysis with Multi-objective Optimization Strategies for the Synthesis/Design and Operation/Control of a PEMFC Fuel Processing Subsystem[END_REF][START_REF] Kim | A hybrid multi-level optimization approach for the dynamic synthesis/design and operation/control under uncertainty of a fuel cell system[END_REF]. Both the sampling methods and the approximation panoramas have been explored: in fact, while the first one is strictly a sampling method, the others use an approximated approach or a mix of both as in case of ANOVA. The investigation led to the conclusion that MCS and RSA could be the two most appropriate methods for the scope of this work. The MCS starts sampling a single value from each initial distribution and uses it as actual variable for the model execution. At this point, the mGT model, working as a deterministic system, returns univocal results for each monitored output. So, in order to create an appropriate group of sampled individuals able to evaluate the uncertainty propagation on the outputs, the procedure has to be repeated many times leading to the construction of the stochastic distribution of outputs. MCS can be accurate but can require large computational efforts. RSA is an approximated method, which utilizes the Taylor series expansion to estimate all the probabilistic parameters for calculation of output distribution. In this work different implementation schemes about RSA have been tested. The RSA methodology was analyzed and validated against MCS by other researchers in [START_REF] Kim | Stochastic Modeling and Uncertainty Analysis with Multi-objective Optimization Strategies for the Synthesis/Design and Operation/Control of a PEMFC Fuel Processing Subsystem[END_REF][START_REF] Kim | A hybrid multi-level optimization approach for the dynamic synthesis/design and operation/control under uncertainty of a fuel cell system[END_REF]. In fact, RSA can be orders of magnitude faster than MCS, but risking in losing accuracy. Moreover, different approximations of the derivatives in the RSA method, for example central and forward finite difference schemes, have been used and compared. The results were evaluated in order to understand the potentiality of RSA method in the proposed case study.

RSA

Simulation of a deterministic model provides a set of outputs which may give an incomplete and frequently misleading representation of the system. To complete this information, the sensitivity (variability) and uncertainties in the results are needed. A deterministic single-point simulation gives no value of the range of sensitivity which may be expected in the system, and it also does not provide information about the uncertainties in the results. A deterministic model (i.e. one without uncertainty considerations) or nondeterministic/probabilistic model (i.e. one with uncertainty considerations) can be described in general terms by:

M f X , Y (1)
where the vector M represents a set of system output values, and vectors X and Y correspond to a set of synthesis/design and operation/control variables, respectively [START_REF] Kim | A hybrid multi-level optimization approach for the dynamic synthesis/design and operation/control under uncertainty of a fuel cell system[END_REF]. If the model is a nondeterministic or probabilistic one and the uncertainties on X and Y are known, their effects on system development and performance can be evaluated via several possible probabilistic design methods. Even though MCS produces exact solutions and is a powerful probabilistic design method for complex nonlinear energy systems, it is not very practical because of the very large computational effort or burden required. This computational difficulty can be overcome by approximated approaches. One of the practical uncertainty analysis methods is the sensitivity-based approximation approach, the RSA method, in which system outputs can be found by a Taylor series expansion [START_REF] Kim | A hybrid multi-level optimization approach for the dynamic synthesis/design and operation/control under uncertainty of a fuel cell system[END_REF]. The first order moment (mean value) and the second order moment (standard deviation) of the outputs can be estimated easily through the Taylor series expansion. If only the mean and variance of each system input (Zi) are known (the exact probability distribution functions are not known) and an implicit nonlinear functional relationship g Z between each system output M and the inputs Z is available, the approximated mean and variance of each system output M can be estimated by using a Taylor series expansion about the mean values μ of the inputs, i.e.,

M ≅ g μ Z μ ∂g Z ∂Z 1 2 Z μ Z μ ∂ g Z ∂Z ∂Z ⋯ (2)
The first and second order approximated means μ of each output, Mj, can be respectively expressed as shown in equation 3 and 4.

μ μ M ≅ g μ , μ , … , μ (3)
μ μ M ≅ g μ , μ , … , μ 1 2 ∂ g ∂Z ν Z (4)
While the first order approximated variance ν is as follow:

ν ν M ≅ ∂g ∂Z ν Z (5)
If there are no explicit functional relationships between the system responses and inputs, the partial derivatives of these functions (i.e., response sensitivities) with respect to the input variables cannot be determined analytically. However, they can be determined numerically using finite difference schemes. For example, for the i-th input variable, using the mean plus or minus one step size δ, the partial derivative in Equation ( 5) can be written as:

∂g ∂Z ≅ M M 2δ (6)
where:

M g μ , μ , … , μ δ , … (7)
M g μ , μ , … , μ δ , … (8)
The above first-order approximations may be successively improved by including higher-order terms in the Taylor series. In order to find the 2nd-order variance, the third and fourth moments of all the Z must be known, which is seldom the case. Thus, for practical purposes, the first-order variance Equation ( 5) and the second-order mean Equation ( 4) are generally used. In equations ( 2) to (11), the derivative term is called the system response sensitivity for M associated with Z .

Note that the perturbation step size δ should be small enough to minimize the truncation error but also large enough to avoid sensitivity to simulation error. If the system has highly nonlinear behaviors and the input values have non-normal distribution characteristics, the RSA method may have some errors in its results. However, it is nonetheless easily applicable to dynamic system simulation because it is a computationally inexpensive method compared to MCS. When system responses show piecewise linear characteristics, the RSA method provides high fidelity analysis results being very close to those provided by MCS. If the truncation error is large and/or step size is too small, it is helpful to use a 4th-order central difference scheme which has a truncation error term of the order (δ 4 ) and it is given by:

∂g ∂Z ≅ M 8M 8M M 12δ (9) ∂ g ∂Z ≅ 16M M 30M 16M M 12δ (10)
Where:

M g μ , μ , … , μ 2δ , … (11)
M g μ , μ , … , μ 2δ , … (12)

PROCEDURE

A first refinement of the model has been carried out by implementing the so called non-value added time costs.

Then, the deterministic model was translated into a stochastic one by inserting the frequency distributions related to the actual behavior of specific activities and by using the RSA method.

MODEL Figure 2. mGT model implemented in Matlab-Simulink environment

The model of the mGT, implemented in Matlab-Simulink environment, is based on the micro-gas turbine Turbec T100 installed at the research laboratory, feeding the electrical and thermal needs of the University campus [START_REF] Ferrari | Micro gas turbine recuperator: steady-state and transient experimental investigation[END_REF].In particular both compressor and turbine maps, based on experimental tests [START_REF] Caresana | Microturbogas cogeneration systems for distributed generation: Effects of ambient temperature on global performance and components' behavior[END_REF], were implemented in the model. The probability distributions of parameters and responses are assumed according to known information or actual measurements. For most engineering problems, a clear probability inference of parameters usually requires a large volume of experimental data, which is often impractical due to expense consideration or experimental limitations. Thus a normal (Gaussian) distribution is popularly adopted without losing the generality, which is, under such circumstance, more appropriate than other distributions. This is due to the facts that (a) normal distributions are often found in engineering problems when the data collection is adequate and (b) in many cases a detailed description of probability distributions is not required since only the means and variances are sought [START_REF] Friswell | Stochastic model updating in structural dynamics[END_REF].

Since the nominal performance, both in terms of electrical efficiency and power&heat outputs, is a critical issue in the performance analysis of a mGT, in this paper the nominal efficiency of compressor and turbine were taken into account as stochastic inputs (Table2).

All the random variables were assumed as independent. A scatter of ±5% was assumed for the efficiencies and 0.5% for rotational speed. Normal distribution was assumed for all the random variables taken into account. In addition, in order to consider the possible uncertainty that could affect the load demand, it was assumed that mGT works at three different load conditions along the year. All these uncertainties were implemented in the model in order to understand their impact in the performance (electrical power, fuel consumption and gross efficiency) and cost (pay-back period) of mGT. The assumptions for the calculation of the payback period are:

 an initial investment of € 130.000 for the purchase of mGT  revenue of € 0.18 for each kWh of electric power consumed by the University campus  revenue of € 0.08 for each kWh of thermal power consumed by the University campus  a cost of € 0.6 per kilogram of fuel consumed (pipeline natural gas)  1% of the total investment cost for the annual maintenance cost  an operating period of 8000 hours per year

C. F. C P C P C m (13) PBP F. C. 1 C C. F. (14)
In particular, the effect that each input has on the considered outputs was analyzed through the RSA method.

RESULTS AND DISCUSSION

The micro-turbine model gives the possibility to monitor various operating parameters, such as: electrical power, fuel mass flow and net efficiency of the system.

First of all, the model was implemented according to a deterministic approach in order to understand the behavior of the mGT in non-stochastic condition: 100% electrical load was assumed for 8000h, which is quite a non-realistic case. This first step is necessary in order to better compare and critically analyze the results obtained with the RSA approach. In Table 4 the output of the deterministic model at each percentage of load are reported: In a further step the probabilistic distribution of outputs were obtained using RSA analysis by introducing the effect of stochastic inputs described before. It is important to underline that all the distributions of parameters reported in Figure 3 were obtained maintaining the 100%±5% of electrical load for the whole year, in order to compare this case with the deterministic one summarized in Table4. The blue distributions are the results of stochastic analysis while each deterministic output is indicated with a squared marker on them. Since there is no difference between deterministic inputs and the means of stochastic input distributions used in RSA calculation, and because the boundary conditions for the two cases are very closed, as we could expect, the RSA approximates quite well the mean of the outputs in respect of the deterministic results, with differences not bigger than 3%. Moreover it returns more important information regarding the uncertainty of these values. Then, the model has been used to evaluate the performance also at part load. The calculation of PBP under stochastic conditions allows to estimate two different important aspects: on one hand, it allows to obtain the most probable period for the return of investment as the mean value of a normal distribution; on the other hand it is useful to evaluate the range of variation. The semi-amplitude of this interval is set here equal to the standard deviation.

The procedure done at 100%±5 load was repeated also at 50%±5 and 20%±5 of electrical load to evaluate the total PBP during the entire life of the mGT. In the end, the same economic analysis was developed, but this time taking in account all the three possible working conditions during the whole year. On the right hand side of Figure 4 the corresponding cash flow diagram is reported. In particular, for this mGT model, the combined effect of uncertainties of all the inputs seems to be more influent as the load is lower, and could explain the difference in estimation of standard deviation.

SENSITIVITY ANALYSIS

An important advantage of RSA method is its capacity to estimate the impact of each single input uncertainty on the monitored outputs at the same time of the probabilistic analysis. This is done through the parameter called sensitivity, defined as the first derivative of outputs with respect to the input variable, that can be non-dimensionalized as it follows:

Sensitivity ∂g ∂Z Z , g , (15)
In this way, it is possible to know how much influence each input can have on the monitored parameters; higher is this value higher in the influence. The series of pictures in Figure 5 shows, on one hand, what is the different effect of inputs at a fixed operating condition, and, on the other hand, how this effect evolves if the working condition is changed. It is clear for this system that reducing the load condition the influence of each input increases: this behavior is observable with different magnitudes but it appears as a common aspect for every parameter.

CONCLUSIONS

A stochastic analysis of a mGT was performed in order to understand how the uncertainty of the input parameters influence the outputs. In particular, the variability in load demand and in performance was taken into consideration to evaluate the impact on fuel consumption, net efficiency and cost (pay-back period). The analysis was performed using the RSA method. In general the RSA approximated method is able to return the most probable value of each monitored variable and the range, described in terms of standard deviation of the outlet probabilistic distribution, where it could fall. The advantage of RSA is that it provides the possibility of treating problems with high level of uncertainty linked to their internal variables, in particular:

 it shows an important saving in computational time if compared to traditional stochastic analysis like Monte Carlo simulation, keeping good accuracy.  it can be used to investigate the direct effect of a single uncertainty on final results, in order to know what are the most influent inputs and what are the negligible ones.

From the economic point of view, in this specific case study, the value of 2.34 years was found to be the most probable Pay Back Period, with a range of ±0.36 years of variation where it could fall with 68% level of confidence.

Figure 1 .

 1 Figure 1. Energy Hub at laboratory of University of Genoa.

Figure 3 .

 3 Figure 3. Outputs distributions: electrical power, net efficiency and fuel consumption.

Figure 4 .

 4 Figure 4. (a) Pay Back Period [yr] distribution and (b) cash flow diagram.

Figure 5 .

 5 Figure 5. Sensitivity of the outputs to each input variable at (a) 100%, (b) 50% and (c) 20% of electrical load condition

Table 1 .

 1 Compressor and turbine data at 100% of load condition

	Variable	Compressor Turbine
	[kg/s]	0.9003	0.9038
	[1/min]	71000	71000
	T [K]	288	1123
	R [J/kgK]	288.2	294.4
	P [Pa]	101300	430100
	;	4.3	4.1
	;	0.78	0.83

Table 2 .

 2 Uncertainty in input variables

	rpm		,		,
	μ	σ	μ	σ	μ	σ
	71000 355	0.78	0.039	0.83	0.0415

Table 3 .

 3 Operating condition per year

	% of the load Operating hours
	100%±5%	1/3* O.P.
	50%±5%	1/3* O.P.
	20%±5%	1/3* O.P.

Table 4 .

 4 Deterministic results

	OUTPUTS	100% load 50%load 20% load
	Pel [kW]	97	49	21
	ηg [%]	0,34	0,22	0,11
	[nl/min] 0,000414 0,000329 0,000278
	PBP [year]	1,17	2,04	3,46

Table 5 .

 5 Mean and standard deviation of outputs at different load conditionsThe total PBP is calculated considering the combination of the three cases reported in Table5: so the mean

		Pel [kW]	/
		μ	σ	μ	σ
	100%	97.52	20.19 4.01E-04 2.80E-05
	50%	48.12	18.29 3.18E-04 2.65E-05
	20%	20.40	16.59 2.76E-04 2.32E-05
		ηg		PBP
		μ	σ	μ	σ
	100%	0.35	0.05	1.18	0.16
	50%	0.22	0.07	2.10	0.40
	20%	0.10	0.08	3.74	1.00


P T Griffiths
email: paul.griffiths@le.ac.uk

O Tatin

Quantifying non-Newtonian effects in the context of transitional rotating boundary layers

Keywords:

The stability of the boundary-layer on a rotating disk is considered for fluids that adhere to a non-Newtonian governing viscosity relationship. For fluids with shear-rate dependent viscosity the base flow is no longer an exact solution of the Navier-Stokes equations, however, in the limit of large Reynolds number the flow inside the three-dimensional boundary-layer can be determined via a similarity solution.

The convective instabilities associated with flows of this nature are described both asymptotically and numerically via separate linear stability analyses. Akin to previous Newtonian studies it is found that there exists two primary modes of instability; the upper-branch type I modes, and the lower-branch type II modes. Results show that both these modes can be stabilised or destabilised depending on the choice of non-Newtonian viscosity model. A number of comments are made regarding the suitability of some of the more well-known non-Newtonian constitutive relationships within the context of the rotating disk model. Such a study is presented with a view to suggesting potential control mechanisms for flows that are practically relevant to the turbo-machinery industry.

INTRODUCTION

The stability and transition of the boundary-layer on a rotating disk is a classical fluid mechanics problem that has attracted a great deal of attention from numerous authors over many decades. In the Newtonian limit there exists a vast wealth of material concerning the solutions of the rotating disk equations.

The pioneering study of Gregory et al. [START_REF] Gregory | On the stability of three-dimensional boundary layers with applications to the flow due to a rotating disk[END_REF] contains the first observation of the stationary cross-flow vortices on a rotating disk. The flow is convectively unstable, within certain regions, to disturbances stationary in the frame rotating with the disk. These disturbances are excited by roughnesses on the surface and, because these roughnesses are fixed in time in the rotating frame, the stationary disturbances are consistently excited and reinforced such that they are evident in flow-visualisation experiments.

Theoretical studies of convective modes within the rotatingdisk system have shown that the flow is susceptible to two distinct modes of instability. The type I (upper-branch) mode due to the cross-flow instability and the type II (lower-branch) mode attributed to external streamline curvature. Because of this the neutral curve has a characteristic two-lobed structure, as noted by Malik [START_REF] Malik | The neutral curve for stationary disturbances in rotating-disk flow[END_REF] and Lingwood [START_REF] Lingwood | Absolute instability of the boundary layer on a rotating disk[END_REF], amongst others. Hall [START_REF] Hall | An asymptotic investigation of the stationary modes of instability of the boundary layer on a rotating disc[END_REF] demonstrated that in the limit of large Reynolds number the upper and lower branches of the neutral curve can be described (with excellent agreement) using an asymptotic theory. More recently, the aforementioned pioneering studies have been extended to include rotating spheres, see Garrett & Peake [START_REF] Garrett | The stability and transition of the boundary layer on a rotating sphere[END_REF] and rotating cones, see Garrett et al. [START_REF] Garrett | The crossflow instability of the boundary layer on a rotating cone[END_REF]. In both cases the flow can be stabilised (or destabilised) with control of the governing parameter, that being the angle of latitude, in the case of the sphere, and the half-angle, in the case of the cone. Indeed, research on the von Kármán boundary-layer remains a topic of active study, the interested reader is referred to the recent review article by Lingwood & Alfredsson [START_REF] Lingwood | Instabilities of the von Kármán Boundary Layer[END_REF].

Rotating disk flows have practical relevance to turbomachinery where non-Newtonian fluids are commonplace. However, far less attention has been given to the corresponding non-Newtonian rotating disk problem. Very recently Griffiths [START_REF] Griffiths | Flow of a generalised newtonian fluid due to a rotating disk[END_REF] determined the base flow profiles for numerous generalised Newtonian fluid models. Using these results we investigate the convective instability of the boundary-layer on a rotating disk for shear-thinning power-law and Carreau fluids. Results show that accurately modelling the variation of viscosity within the boundary-layer is of paramount importance.

FORMULATION

We consider the flow of a steady incompressible non-Newtonian fluid due to an infinite rotating plane located at z * = 0 (here the superscript * denotes dimensional quantities). The plane rotates about the z * -axis with angular velocity Ω * . The motion of the fluid is in the positive z * direction, the fluid is infinite in extent and the only boundary is located at z * = 0. In a rotating frame of reference the continuity and Navier-Stokes equations are expressed as

∇ * • u * = 0, ( 1a
)
∂ u * ∂t * + u * • ∇ * u * + Ω * × (Ω * × r * ) + 2Ω * × u * = - 1 ρ * ∇ * p * + 1 ρ * ∇ * • τ * . (1b)
Here u * = ( Ũ * , Ṽ * , W * ) are the velocity components in cylindrical polar coordinates (r * , θ, z * ), t * is time, Ω * = (0, 0, Ω * ) and r * = (r * , 0, z * ). The fluid density is ρ * and p * is the fluid pressure. For generalised Newtonian models, such as the power-law model, the stress tensor is given by τ * = µ * ( γ * ) γ * , where γ * = ∇ * u * + (∇ * u * ) T is the rate of strain tensor and µ * ( γ * ) is the non-Newtonian viscosity. The magnitude of the rate of strain tensor is γ * = ( γ * : γ * )/2. The governing relationships for µ * ( γ * ) that will be considered herein are:

µ * = m * ( γ * ) n-1 , (2a)
µ * = µ * ∞ + (µ * 0 -µ * ∞ )[1 + (λ * γ * ) 2 ] (n-1)/2 , (2b)
these are the power-law and Carreau models, respectively. Here m * is the consistency coefficient and n is the fluid index, for n > 1 the fluid is said to be shear-thickening, whilst for n < 1 the fluid is said to be shear-thinning. The Newtonian viscosity relationship is recovered when n = 1 and µ * 0 = µ * ∞ , respectively. The infinite-shear-rate viscosity is µ * ∞ , the zeroshear-rate viscosity is µ * 0 and λ * is the characteristic time constant, often referred to as the 'relaxation time'.

In the Newtonian limit an exact solution of (1) exists, as was first determined by von Kármán [START_REF] Von Kármán | Über laminare und turbulente Reibung[END_REF]. However, no such solution exists for flows with n 1. It is only in the large Reynolds number limit that the leading order boundary-layer equations admit a similarity solution analogous to the exact Newtonian solution. The governing boundary-layer equations are given in Griffiths [START_REF] Griffiths | Flow of a generalised newtonian fluid due to a rotating disk[END_REF], for brevity we exclude these here.

We introduce the generalisation of the classic Newtonian similarity solution in order to solve for the steady mean flow relative to the disk. The dimensionless similarity variables are defined by

U (η) = Ũ * 0 r * Ω * , V (η) = Ṽ * 0 r * Ω * , W (η) = W * 0 χ * , P(η) = P * 1 ρ * ( χ * ) 2 , (3) where χ * = [ν * /(r * 1-q Ω * 1-2q )] 1/(q+1) . Here ( Ũ * 0 , Ṽ * 0 , W * 0 )
are the leading order velocity components, P * 1 is the leading order fluid pressure term and ν * = σ * /ρ * is the kinematic viscosity. The dimensionless similarity coordinate is η = r (1-q)/(q+1) z, where here r and z have been made dimensionless with respect to L * = (ν * /Ω * 2-q ) 1/2 . For power-law fluids q = n and σ * = m * , whereas for Carreau fluids q = 1 and σ * = µ * ∞ . Thus the laminar-flow profiles are determined from the following set of non-linear ordinary differential equations:

2U + ηU + W = 0, (4a)
U 2 -(V + 1) 2 + (W + ηU)U -(µU ) = 0, (4b) 2U (V + 1) + (W + ηU)V -(µV ) = 0, (4c)
where η = η(1q)/(q + 1) and the primes denote differentiation with respect to η. The power-law viscosity function is given by

µ = (U 2 + V 2 ) (n-1)/2 , (4d)
whilst the Carreau viscosity function takes the form

µ = 1 + c 0 [1 + k 2 (U 2 + V 2 )] (n-1)/2 . ( 4e
)
Here c 0 = (σ * -µ * 0 )/σ * is the viscosity ratio, and k = r * λ * Ω * ( ρ * Ω * /σ * ) 1/2 is the dimensionless equivalent of λ * . Throughout this study we fix the values of c 0 and k such that c 0 = 1 and k = 100, this is consistent with previous investigations.

Using a fourth-order Runge-Kutta quadrature routine twinned with a Newton iteration scheme to determine the values of the unknowns U (0) and V (0) the set of ordinary differential equations are solved subject to

U (0) = V (0) = W (0) = 0, U (η → ∞) → 0, V (η → ∞) → -1. (5)
Although the formulation here is different these results are an exact reproduction of those presented by Griffiths [START_REF] Griffiths | Flow of a generalised newtonian fluid due to a rotating disk[END_REF]. The laminar-flow profiles are presented in figures 1 (power-law) and 2 (Carreau) for a range of shear-thinning values of n.

The differences between the base flow profiles, and their significance in relation to the resulting linear stability calculations, will be discussed in section 4.

STABILITY ANALYSES

Type I asymptotic analysis

This analysis is based on the assumption that the Reynolds number is large and that the disturbances have wavelengths scaled on the boundary-layer thickness, given by δ = Re -1/(q+1) . Here Re is the asymptotic representation of the Reynolds number. When considering the type I modes we observe the existence of three distinct layers. An inviscid layer, or zone, a wall layer and a critical layer. The inviscid zone encompasses the entirety of the boundary-layer, the wall layer is needed to ensure the no-slip condition at the wall is satisfied and the critical layer exists so that the singularities that arise within the inviscid zone are smoothed out.

By solving for zero-order perturbations in the linear disturbance equations we arrive at a modified form of Rayleigh's equation that is dependent on the parameter n. At the next order we determine a nonhomogeneous modified form of Rayleigh's equation that is required to match our leading order solution in the wall layer. This wall layer solution is given in terms of the decaying Airy function. Despite the appearance of additional viscous terms in the leading order governing equations we find that analytic solutions are obtainable. Matching these solutions between the two layers enables us to plot solutions for the asymptotic wavenumber and wave angle of the disturbances as functions of R = r 2/(q+1) Re 1/(q+1) . Here R is the numerical representation of the Reynolds number. Doing so enables us to make direct comparisons between the asymptotic and numerical solutions. In this study we present only the type I asymptotic analysis with the viewpoint of validating our numerical results when the Reynolds number is large. Indeed initial investigations have begun on the corresponding type II analysis, we hope to report on this investigation in due course. For more details regarding the type I analysis see Griffiths et al. [START_REF] Griffiths | Stability of the boundary layer on a rotating disk for power law fluids[END_REF].

Numerical analysis

The stability analysis, applied at a specific radius, involves imposing infinitesimally small disturbances on the steady mean flow, in the form of scaled normal-mode quantities

(u, v, w, p) = ( û, v, ŵ, p)(η; α, β, ω; R, n)e i(αr +βθ-ωt ) . (6)
The frequency of the disturbance in the rotating frame is ω (taken to be zero in this stationary study), the complex radial wavenumber is α = α r + α i and β is the real azimuthal wavenumber.

After making an approximation akin to the parallel-flow assumption and a viscous assumption stemming from the asymptotic results, the stability equations may be written as sets of six first-order ODEs using transformed perturbing variables φ i (η) with i = 1, 2 . . . 6. The governing power-law equations can be found in Griffiths et al. [START_REF] Griffiths | The neutral curve for stationary disturbances in rotating disk flow for power-law fluids[END_REF], the equivalent Carreau fluid equations can be easily inferred from there. In all that follows the eigenvalue problem defined by the stability equations is solved with the homogeneous boundary conditions

φ i = 0 at η = 0, (7a) φ i → 0 as η → ∞, (7b)
for all i. This eigenvalue problem is solved for certain combinations of values of α, β and ω at each Reynolds number, R, and for the specified value of n. From these we form the dispersion relation, D(α, β, ω; R, n) = 0, at each n, with the aim of studying the convective instabilities. The step size in η was reduced and the value of infinity increased until there were no discernible differences in the numerical results. The values taken were such that the boundary layer was approximated by 2000 equally spaced data points between η = 0 and η = 20. This discretization is known to be consistent with Lingwood [START_REF] Lingwood | Absolute instability of the boundary layer on a rotating disk[END_REF] and Garrett & Peake [START_REF] Garrett | The stability of the boundary layer on a rotating sphere in a uniform axial flow[END_REF], for example, and represents an appropriate balance between accuracy and computational effort for each n.

In order to investigate the structure of the spatial branches at each n, we solve the dispersion relation for α whilst marching through values of β at fixed R. For each n in the particular range of interest two spatial branches determine the convective instability characteristics of the system. Neutral curves, defined by α i = 0, have been calculated for a variety of shearthinning values of n using both the power-law and Carreau fluid models.

DISCUSSION & CONCLUSIONS

Our results, both asymptotic and numerical, predict differing behaviour for shear-thinning fluids described by the power-law model when compared to those described by the Carreau model. We observe from figure 3 that the power-law model predicts a strong stabilising effect as the fluid index is decreased, whereas the Carreau model suggests that shear-thinning destabilises the boundary-layer flow. Under the power-law regime the value of the critical Reynolds number significantly decreases with n, as does the area encompassed by the neutral curves. Contrary to this, using the Carreau fluid model, we predict that the critical Reynolds number will be marginally increased as the shear-thinning parameter is increased, and moreover the area encompassed by the neutral curve will also significantly increase.

Figure 4 outlines the excellent agreement obtained between the results of the asymptotic analysis when compared to our numerical solutions. The large Reynolds number asymptotic analysis does indeed predict boundary-layer stabilisation in the power-law limit and destabilisation in the Carreau limit. This is wholly consistent with the numerical predictions at much lower Reynolds numbers.

The opposing results for our two shear-thinning models are indeed quite striking. These surprising results are attributed to the inability of the power-law model to describe shearthinning flows for vanishing shear-rates, i.e. far from the disk surface. It is well known that for shear-thinning power-law fluids µ * ( γ * → 0) → ∞, which is exemplified by our base flow solutions where we see that the viscosity function grows without bound as we move away from the surface of the disk, or equivalently move towards a region of low shear-rate. This essentially has the effect of increasing the predicted boundarylayer thickness. Contrastingly, we have that µ * ( γ * → 0) → µ * 0 for shear-thinning Carreau fluids. Thus, far from the disk, in regions of low shear-rate, we find that µ → 1 + c 0 , and in this case the boundary-layer thickness is reduced. It is this failing of the power-law model in regions of low shear rates that dramatically affects the predicted velocity profiles and therefore our linear stability results.

It transpires that the inability of the power-law model to accurately describe shear-thinning flows in the limit as γ * → 0 has such a significant effect on the base flow profiles that it in turn effects the predicted stability characteristics. Therefore it is with relative confidence that we can say that the results owing from the Carreau model do indeed provide a much better description of the observed cross-flow instabilities.

These results clearly show the importance of accurately modelling the variation of viscosity within the boundarylayer. The power-law model may be useful for describing experimental results in regions of moderate shear-rate, however, it fails within this theoretical framework. As such, we conclude that the Carreau model provides a better physical representation of the boundary-layer flow and hence that the introduction of shear-thinning fluids will have a destabilising effect. This work clearly motivates the need for detailed experimental results with which to compare our theoretical analysis.

To the best of the author's knowledge no such experiments have yet taken place. Certainly in the Newtonian regime a wealth of literature exists, and this is currently a topic of particular interest, see Imayama et al. [START_REF] Imayama | A new way to describe the transition characteristics of a rotatingdisk boundary-layer flow[END_REF][START_REF] Imayama | An experimental study of edge effect on rotating-disk transition[END_REF][START_REF] Imayama | On the laminar-turbulent transition of the rotating-disk flow: the role of absolute instability[END_REF], for example. Personal communication with these authors has revealed the difficulty of obtaining consistently accurate experimental results; therefore we can only envisage the introduction of non-Newtonian fluids would serve to significantly complicate any experimental procedure.

It is acknowledged that the parallel-flow approximation utilised within the numerical stability analysis means that the perturbation equations solved here are not rigorous at O(R -1 ). Although it is clear that the approximation will lead to inaccuracies at the predicted critical Reynolds numbers, it is our opinion that these will be small. The excellent agreement obtained between our numerical neutral curves and our asymptotic predictions shows that the effects of this approximation are negligible at high Reynolds number.

Figure 1 .

 1 Figure1. Plots of U, V , W and µ versus η for power-law fluids with n = 1, 0.9, 0.8, 0.7, 0.6. The η-axis has been truncated at η = 20.

Figure 2 .

 2 Figure 2. Plots of U, V , W and µ versus η for Carreau fluids with n = 1, 0.95, 0.75, 0.5, 0.25. The η-axis has been truncated at η = 20.

Figure 3 .

 3 Figure 3. Wavennumber neutral stability curves for shear-thinning power-law fluids and Carreau fluids.

Figure 4 .

 4 Figure 4. Neutral stability curves for large R with corresponding asymptotic predictions. Excellent agreement is observed in both cases for every n in the range of interest.

ACKNOWLEDGMENTS

The author gratefully acknowledges support from the Engineering and Physical Sciences Research Council UK.


S J Garrett

A J Cooper

M Özkan

P J Thomas

S Y Mpo

O Tatin

Towards roughness-based drag reduction in cross-flow dominated flows

Keywords:

 and the surface-geometry approach of Yoon et. al (2007). An energy analysis reveals that for both instability modes, the main contributors to the energy balance are the energy production by the Reynolds stresses and conventional viscous dissipation. For the Type I mode, energy dissipation increases and the Reynolds-stress energy production decreases with roughness under both models. This suggests a clear stabilising effect of the anisotropic roughness on the Type I mode. For the Type II mode, the Reynolds-stress energy production increases with roughness under both models. However, the energy dissipation of the Type II mode decreases with the roughness under the surface-geometry model and increases under the partial-slip model. This sensitivity to the precise form of the anisotropic roughness suggests that maximising dissipation by an appropriately designed roughness can theoretically lead to an overall beneficial stabilisation of both the Type I and Type II modes. This is a potential route to overall boundary-layer-transition delay and drag reduction in cross-flow dominated flows.

INTRODUCTION

The von Kármán boundary layer represents a typical, generic example of a general class of fully three-dimensional boundary layers that share the common, characteristic feature of a cross-flow velocity component [START_REF] Kármán | Ueber laminare und turbulente Reibung[END_REF][START_REF] Lingwood | Instabilities of the von Kármán boundary layer[END_REF][START_REF] Saric | Stability and transition of three-dimensional boundary layers[END_REF]. Similar boundary layers are encountered in many applied contexts such as, for instance, on the blades of wind turbines or over the highly swept wings of aircraft.

Here we study effects of distributed surface roughness on the laminar-turbulent transition process of the von Kármán boundary-layer flow. All boundary layers with a cross-flow component display similar laminar-turbulent transition characteristics due to the existence of an inflection point on the cross-flow velocity profile [START_REF] Lingwood | Instabilities of the von Kármán boundary layer[END_REF][START_REF] Saric | Stability and transition of three-dimensional boundary layers[END_REF]. Hence, the results presented here are of general, direct practical relevance to all applied flow configurations where a boundary layer with a cross-flow component is established. Our motivating, long-term aim is to develop theoretical methods to enable the energeticallyoptimal design of surface-roughness that can be exploited, for boundary layers with a cross-flow component, in the context of the development of new, passive drag-reduction techniques.

It has now been firmly established that, contrary to the classic belief, the interaction of boundary-layer flow with the right sort of roughness on surfaces can result in energetically beneficial, drag-reducing effects [START_REF] Sirovich | Turbulent drag reduction by passive mechanisms[END_REF][START_REF] Carpenter | The right sort of roughness[END_REF][START_REF] Choi | The rough with the smooth[END_REF]. The challenge that remains, however, is to identify what represents the right sort of roughness that leads to such drag-reducing effects in any particular application. One fundamental, general strategy known to result in reduced drag is to control the laminar flow and delay its transition to turbulence. This method of the stabilisation of the boundary-layer flow exploits the fact that laminar flows are subject to smaller dissipative energy losses than turbulent flows. The theoretical results to be discussed here will reveal clear evidence of stabilising effects of the investigated surface structures on von Kármán boundary-layers on rotating disks.

Two distinct models for the steady boundary-layer flow over rough rotating disks exist in the literature. They were developed by Yoon et al. [START_REF] Yoon | Flow and heat transfer over a rotating disk with surface roughness[END_REF] and Miklavčič & Wang [START_REF] Miklavčič | The flow due to a rough rotating disk[END_REF] and are henceforth referred to as the YHP and MW models, respectively. Both models show how successively increasing roughness levels lead to deviations from the classic similarity solution for the flow over a smooth disk due to von Kármán [START_REF] Kármán | Ueber laminare und turbulente Reibung[END_REF]. The main goal of this continuing study is to demonstrate that theoretical predictions for the stability characteristics of the von Kármán boundary layer, based on the two fundamentally different modelling approaches, lead to overall consistent, and energetically beneficial, results.

The YHP approach models roughness by imposing a particular surface distribution as a function of radial position and assumes a rotational symmetry. The YHP approach therefore models roughness in the radial direction only, i.e. the flow over disks with concentric grooves of a particular crosssectional shape profile. Alternatively, rather than imposing a particular mathematical form of the surface roughness, the MW approach models roughness empirically by replacing the usual no-slip boundary conditions with partial-slip conditions at the disk surface. The MW model is therefore capable of modelling independent levels of roughness in the radial and azimuthal directions by separately modifying the boundary conditions in these directions.

The case of uniform levels of roughness in every direction over a surface is typically referred to as isotropic roughness, and non-uniform levels of roughness as anisotropic roughness. The MW model is clearly capable of modelling both anisotropic and isotropic roughness; the YHP model, however, is only capable of modelling the roughness in the radial direction, i.e. a particular case of anisotropic roughness.

It is important to note that the analytical results of the YHP and MW models have not yet been corroborated by experiment and so the two models should be considered as competing models of an actual flow. For this reason, it is instructive to consider a comparative study of the stability characteristics of the steady flows arising from both models. Comparisons will be made as appropriate between the stability characteristics of the YHP flows and those of 'equivalent' levels of radial anisotropic roughness under the MW flow.

In our previous publication [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF] we have considered all classes of flows arising from the full MW model. In particular, our previous study considered isotropic roughness and two forms of anisotropic roughness: concentric grooves (radial anisotropic, as will be considered here) and radial grooves (azimuthal anisotropic). The results obtained reveal that both anisotropic and isotropic roughness results in a stabilisation of the Type I instability mode. In addition, in the case of radial anisotropic roughness, a significant destabilisation of the Type II mode is found. The emergence of the Type II mode with increasing roughness should therefore be understood as halting the beneficial roughness effects on the Type I mode. The interested reader is also referred to that paper for an extensive literature review of previous rotating-disk and surface-roughness studies.

This current paper should be considered a summary of our other paper [START_REF] Garrett | On the stability of von Kármán rotating-disk boundary layers with radial anisotropic surface roughness[END_REF] and full details can, of course, be found there.

THE STEADY FLOWS

The surface-geometry model due to YHP [7]

Full details of this formulation can be found in papers [START_REF] Yoon | Flow and heat transfer over a rotating disk with surface roughness[END_REF][START_REF] Garrett | On the stability of von Kármán rotating-disk boundary layers with radial anisotropic surface roughness[END_REF]. Here it is sufficient to understand that, under the YHP approach, the lower surface is described by a particular functional form expressed within the dimensional cylindrical polar coordinate system (r * , θ, z * ). The disk is considered to be rotating about its axis of symmetry at a constant rotation rate Ω * and we formulate the analysis in the rotating frame. All dimensional quantities are scaled on a characteristic length-scale given by the boundary-layer thickness, d * = √ ν * /Ω * , where ν * is the kinematic viscosity, and a velocity scale given by r * Ω * . This leads to the Reynolds number Re = r * Ω * d * /ν * = r.

The particular surface function used here scales to

s(r) = δ cos 2πr γ (1)
which gives two non-dimensional control parameters: δ, the height of the roughness, and γ, the pitch of the roughness. It is useful to define the aspect ratio a = δ/γ which we henceforth refer to as the roughness parameter within the YHP model. It is necessary to transform out the surface distribution before attempting to solve the governing equations. To this end we use a new coordinate system, (r, θ, ζ ), a modified form of the non-dimensional cylindrical polar co-ordinates, defined by the transformation ζ = zs(r). These lead to transformed, rotationally symmetric flow components U (r, ζ ), V (r, ζ ) and W (r, ζ ) in the radial, azimuthal and normal directions, respectively. The governing boundary-layer equations for the steady flow are then obtained after introducing variables closely related to the von Kármán similarity variables,

f (r, ζ ) = 1 r U (r, ζ ), g(r, ζ ) = 1 r V (r, ζ ), h(r, ζ ) = W (r, ζ ) (2)
and given by

2 f + r ∂ f ∂r + ∂h ∂ζ = 0, (3)
r f ∂ f ∂r + h ∂ f ∂ζ + 1 + r s s 1 + s 2 f 2 = (4) (1 + s 2 ) ∂ 2 f ∂ζ 2 + (1 + g) 2 1 + s 2 , (5)
r f ∂g ∂r + h ∂g ∂ζ = (1 + s 2 ) ∂ 2 g ∂ζ 2 -2 f (1 + g), (6)
These PDEs are subject to the usual no-slip and quiescent boundary conditions in the rotating frame and can be solved using the commercially available NAG routine D03PEF. The transformed flow field across (r, ζ ) is found to vary at two distinct spatial scales in the radial direction. At the scale characterised by γ, we have a response dependent on where r is within the oscillatory cycle of the surface function [START_REF] Kármán | Ueber laminare und turbulente Reibung[END_REF]. In addition to this oscillatory behaviour, we see a similaritytype solution scaling with r (as per von Kármán) at the larger spatial scale. For γ < O(10 -1 ), as is envisaged here, we argue that the small-scale response of the viscous flow will not occur in practice and it is a reasonable approximation to take an ensemble average of the flow field over any complete cycle in r. This approach leaves only the similarity-solution variation within the averaged flow field, ( f (ζ ), ḡ(ζ ), h(ζ )). Note that overbars have been introduced to denote averaged quantities. The ensemble average acts to 'average away' the oscillatory surface distribution, that is s(r) = 0 and so ζ → z. Under our approach the surface roughness is therefore seen to lead to a modified von Kármán flow, denoted ( f (z), ḡ(z), h(z)).

Note that throughout this study we compute all ensemble averaged quantities at 100 regularly spaced locations over one wavelength and the results have been confirmed to be independent of the starting radial position. Our results also show that the aspect ratio a = δ/γ determines the flow response and so, despite having two control parameters, we can work in terms of the single roughness parameter, a.

The partial-slip model due to MW [8]

Rather than imposing a particular mathematical form for the surface roughness, the MW approach assumes that roughness can be modelled by a modification of the no-slip conditions at the disk surface. In particular, the model assumes partial slip at the disk surface but is otherwise identical to the von Kármán formulation [START_REF] Kármán | Ueber laminare und turbulente Reibung[END_REF]; full details can be found elsewhere [START_REF] Miklavčič | The flow due to a rough rotating disk[END_REF][START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF]. The full MW model has two parameters η and λ (giving empirical measures of the roughness in the radial and azimuthal directions, respectively). These appear in the surface boundary conditions of the von Kármán ODEs such that

f (0) = λ f (0), ḡ(0) = η ḡ (0),
where a prime denotes differentiation with respect to the normal spatial variable. Here we are concerned with the particular case of radial anisotropic roughnesses, consistent with the capabilities of the YHP model, and so set λ = 0.

Comparison of steady flows

Despite both models predicting a reduced radial jet and a thickened boundary layer, the areas enclosed by the radial profiles are found to increase with roughness under the YHP model and decrease in the MW model. This area is a measure of the volume of fluid transported outwards in the radial direction and accounts for the different behaviour of the axial entrainment between the two models. Given the different physical predictions arising from the two models and the empirical definition of roughness in the MW model, a direct quantitative comparison between 'equivalent' levels of roughness is not possible. Instead we proceed with a qualitative comparison of the effects of increasing roughness under both models. The particular values of η used here are therefore reasonably arbitrary and we have opted to use the maximum value of the radial jet as a matching parameter. That is, for each value of a in the YHP model, the value of η in the MW model is chosen such that the maximum values of the radial wall jet, max( f ), agree. Note that the azimuthal and wall-normal components can never be matched between the two models and we emphasise again that direct quantitative comparisons should not be made.

The paired parameter values are found to be a = 0.1 ∼ η = 0.14, a = 0.2 ∼ η = 0.57 and a = 0.3 ∼ η = 1.18 and comparisons between the resulting flow profiles can be seen in Figures 1(a), 2(a) and 3(a). We observe good agreement between the two models for small values of the roughness and an increasing discrepancy as the roughness increases.

CONVECTIVE INSTABILITY ANALYSIS

The two approaches used to calculate the steady flows have both resulted in similarity solutions in the scaled physical space (r, θ, z). The resulting flows are therefore related to the von Kármán flow, and, importantly, the stability analyses of the rotating-disk flow presented elsewhere [START_REF] Lingwood | Absolute instability of the boundary layer on a rotating disk[END_REF][START_REF] Lingwood | The effects of surface mass flux on the instability of the BEK system of rotating boundary layer flows[END_REF] are directly applicable in this current study. Full details of the governing perturbation equations can be found in those references. Here it is sufficient to understand that we conduct a normal-mode analysis with perturbations of the form ( û, v, ŵ, p) = (u(z), v(z), w(z), p(z))e i(αr +β Reθ-ωt ) .

The wavenumber in the radial direction, α = α r + iα i , is complex, as required by the spatial convective analysis to be conducted; the frequency, ω, and circumferential wavenumber, β, are real. It is assumed that β is O(1) and the integer number of complete cycles of the disturbance around the azimuth is n = βRe. We identify n with the number of spiral vortices around the disk surface. Furthermore, the orientation angle of the vortices with respect to a circle centred on the axis of rotation is = arctan( β/α). The quantities n and can be compared directly to experimental observations. In what follows we are concerned with stationary vortices that rotate with the rough surface and so set ω = 0.

The governing perturbation equations are solved using a Chebyshev polynomial discretisation method in the wallnormal direction to obtain solutions of the dispersion relation D(α, β; Re, [a, η]) = 0 with the aim of studying the occurrence of convective instabilities for various values of the roughness parameters. The results of our code in the smooth case (a = 0 = η) have been compared against those in the literature and the predictions for the critical parameters of the Type I mode are found to be entirely consistent with other published results [START_REF] Malik | The neutral curve for stationary disturbances in rotating-disk flow[END_REF][START_REF] Cooper | The stability of rotating-disc boundary-layer flow over a compliant wall. Part 1. Type I and II instabilities[END_REF][START_REF] Appelquist | Direct numerical simulations of the rotating-disk boundary-layer flow[END_REF][START_REF] Turkilmazoglu | Convective and absolute instability in the incompressible boundary layer on a rotating-disk[END_REF][START_REF] Garrett | The crossflow instability of the boundary layer on a rotating cone[END_REF].

As with existing analyses of smooth rotating disks, modes of Type I and II are found to determine the convective instability properties of the disturbance modes over rough rotating disks. The Type I mode, appearing as the upper lobe in Re-α r neutral curves, is known to arise from the inflectional nature of the steady-flow profiles, and the Type II mode, appearing as the lower lobe, is known to arise from streamline curvature and Coriolis effects. The neutral curves arising from the analysis of both models at comparable roughness parameters are shown in Figures 1(b), 2(b) and 3(b). The neutral curves display the same qualitative behaviour: the Type I lobe is diminished (both in terms of critical Re and width) with increased roughness, and the Type II mode exaggerated. This is entirely consistent with the results of our previous study of radial anisotropic roughness [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF]. The results of the YHP model appear much more sensitive to increased roughness, although this conclusion depends somewhat on the equivalence of the parameter values used for the comparisons. Critical parameters at the onset of unstable Type I and Type II modes are given in Table 1.

ENERGY ANALYSIS

Following previous work [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF][START_REF] Garrett | On the stability of von Kármán rotating-disk boundary layers with radial anisotropic surface roughness[END_REF][START_REF] Cooper | The stability of rotating-disc boundary-layer flow over a compliant wall. Part 1. Type I and II instabilities[END_REF] an integral energy equation for three-dimensional disturbances ( û, v, ŵ) to the undisturbed three-dimensional boundary-layer flow (U,V,W ) is derived in order to extract possible underlying physical mechanisms behind the effects of roughness on the stability of rotating disk boundary-layer flow. Essentially, the energy-balance approach enables one to assess the relative influences of the various energy transfer mechanisms affecting the destabilisation of fluid disturbances. The method was used in extensively for the full MW model in our previous publication [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF] and full details are presented there. As demonstrated elsewhere [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 rotating-disc boundary-layer flow over a compliant wall. Part 1. Type I and II instabilities[END_REF], the energy equation that applies to a particular eigenmode is given by

-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 , (7)
where, physically, the terms on the right-hand side are identified as follows:

(I) the Reynolds stress energy production term, obtained from {P i },

(II) the viscous dissipation energy removal term, obtained Calculations have been carried out for both roughness models, and for both the Type I and II modes at Re = 400. The corresponding growth rates are shown in Figure 4. This emphasises the stabilising effect of roughness on the Type I mode, the destabilising effect on the Type II mode and the stronger effect on both of these modes for the YHP model. For the MW model the amplification of the Type II mode is more modest, even though the Type II lobe of the neutral curve shows similar augmentation to the YHP case, and the stabilising effect on the Type I mode is not quite so strong.

By calculating all terms in the energy equation [START_REF] Yoon | Flow and heat transfer over a rotating disk with surface roughness[END_REF] it is possible to identify where the effects of roughness are the greatest. Given the boundary conditions for the YHP model some terms in the energy equation are identically zero (PW 2 , S 1 , S 2 , S 3 ). The results of the energy balance for the three roughness values a = 0.1, 0.2 and 0.3 are compared to those for a smooth disk (a = 0) in Figure 5. For both modes the main contributors are energy production by the Reynolds stress (P 2 ) and conventional viscous dissipation (D). Terms P 1 , P 3 , PW 1 and G 2 are found to be negligible and the geometric terms G 1 and G 3 remove energy from the system. The strongly stabilising effect of roughness on the Type I mode is manifested in a striking reduction in P 2 and a slight increase in viscous dissipation. Conversely, the growth with roughness of the Type II mode arises from a net increase in energy production through increased Reynolds stress alongside a reduction in viscous dissipation.

Results of the energy balance calculation for the MW model are shown in Figure 6. Again both modes are dominated by contributions from P 2 and D. The main difference from the YHP model in the case of the Type II mode is that, although the Reynolds stress energy production term increases as before, the MW model also shows an increase in viscous As detailed in publication [START_REF] Garrett | On the stability of von Kármán rotating-disk boundary layers with radial anisotropic surface roughness[END_REF] where a presentation of the disturbance eigenmodes is given, the Type II disturbances generally extend further into the boundary layer than the Type I disturbances. The further stretching of the disturbance profile of the YHP model as roughness is increased, together with the thickening of the boundary layer with roughness, would appear to contribute to the augmentation of the Type II instability mode. The MW model results in a similar stretching of the Type II eigenmodes to the YHP case. However, since there is no significant increase in the boundary-layer thickness for the MW steady flow, the fluid disturbances do not extend quite as far into boundary layer as for the YHP case. These affects appear to be related to the subtly different responses of the energy balance of the Type II mode across each roughness model. The Type I disturbance profiles are very similar to those for the YHP model.

CONCLUSION

We have summarised our theoretical study investigating the effects of radial anisotropic surface roughness resulting from concentric grooves on the stability of the von Kármán boundarylayer flow over a rotating disk. Our theoretical analysis was based on the two alternative approaches suggested by Yoon et al. [START_REF] Yoon | Flow and heat transfer over a rotating disk with surface roughness[END_REF] and by Miklavčič & Wang [START_REF] Miklavčič | The flow due to a rough rotating disk[END_REF] for modelling the steady boundary-layer flow over a rough rotating disk. Our results have shown that both these fundamentally different approaches yield qualitatively consistent results for the two dominant steady-flow velocity components (radial and azimuthal) of the rotating-disk boundary layer for increasing roughness levels. Similarly, the subsequent linear stability analysis based on the steady-flow profiles obtained from the two modelling approaches revealed that in both cases the type of roughness investigated results in a stabilisation of the Type I instability mode of the rotating-disk boundary layer and in a destabilisation of the Type II mode. However, it was found that, overall, the effects arising when following the YHP model are somewhat more pronounced than when using the MW approach.

The energy analysis has revealed that for both the YHP and the MW approach, and for both the Type I and Type II instability modes, the main contributors to the energy balance are the energy production by the Reynolds stresses and conventional viscous dissipation. For the Type I mode dissipation increases with the roughness level and the increase is more pronounced for the MW model than for the YHP model. For the Type I mode Reynolds-stress energy production decreases with the roughness level and the decrease is more pronounced for the YHP model than for the MW model. Thus, in summary, increased energy dissipation and decreased energy production by Reynolds stresses implies a stabilisation of the Type I mode by increasing roughness levels.

For the Type II mode Reynolds-stress energy production increases with the roughness level and the increase is slightly less pronounced for the YHP model than for the MW model. The main qualitative difference is observed for the energy dissipation of the Type II mode. For the Type II mode energy dissipation decreases with the roughness level for the YHP model whereas there is a slight increase under the MW model. Yet, the overall increased Reynolds-stress energy production and decreased energy dissipation result in a destabilisation of the Type II mode for both the YHP and the MW model.

Our study suggests that the beneficial stabilisation of the Type I mode by concentric roughnesses becomes suppressed when the Type II mode is destabilised as it moves upstream and eventually becomes the critical mode at the lowest Reynolds number. This is consistent with our previous result [START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF]. The results of the energy analysis imply that the dissipation of the Type II mode is sensitive to the precise form of the steady-flow base profile. Consequently, maximising dissipation by an appropriately designed surfaceroughness pattern, that leads to the energetically optimal base profile, can theoretically lead to an overall beneficial stabilisation of the Type II mode. Provided that the Type I mode is not adversely affected, this could result in a boundary-layertransition delay and drag reduction. This points to a possible way forward for exploiting the beneficial stabilising effects of radial anisotropic roughness on the Type I mode in future drag reduction techniques relying on transition delay for boundary layers with a cross-flow component. Full details of this study into radial anisotropic roughness can be found in publication [START_REF] Garrett | On the stability of von Kármán rotating-disk boundary layers with radial anisotropic surface roughness[END_REF].

Figure 1 .

 1 Figure 1. Model comparison at a = 0.1 and η = 0.14.

Figure 2 .

 2 Figure 2. Model comparison at a = 0.2 and η = 0.57.

Figure 3 .

 3 Figure 3. Model comparison at a = 0.3 and η = 1.18.

Figure 4 .

 4 Figure 4. Type I and II growth rate curves at Re = 400.

  Type II mode with n = 12.

Figure 5 .

 5 Figure 5. Results of energy analysis for YHP model at Re = 400.

  Type II mode with n = 8.

Figure 6 .

 6 Figure 6. Results of energy analysis for MW model at Re = 400.

Table 1 .

 1 Critical values of measurable parameters at the onset of instability under both models. Type I and (Type II). Bold text indicates the most dangerous mode in terms of critical Reynolds number. dissipation (opposite to YHP) which would account for the more modest growth observed in Figure4. The Type I mode shows a less pronounced decrease in P 2 , but more viscous dissipation compared to the YHP model.

			YHP model
		Re	n
	a = 0	286.1 (461.5) 22.2 (21.3) 11.4 (19.2)
	a = 0.1	311.5 (394.4) 20.7 (16.7) 11.1 (19.5)
	a = 0.2	426.8 (283.3) 20.6 (11.1) 10.8 (19.0)
	a = 0.3	593.9 (220.6) 21.8 (8.5) 11.1 (18.5)
			MW model
		Re	n
	η = 0	286.1 (461.5) 22.2 (21.3) 11.4 (19.2)
	η = 0.14 300.6 (390.3) 19.6 (17.5) 9.8 (16.9)
	η = 0.57 343.7 (311.5) 15.4 (9.2)	7.3 (12.3)
	η = 1.18 393.8 (284.9) 12.4 (6.3)	5.6 (9.2)

ACKNOWLEDGMENTS

SJG is supported by a Senior Research Fellowship of the Royal Academy of Engineering, funded by the Leverhulme Trust. SJG also wishes to acknowledge the hospitality of the Faculty of Science, School of Mathematics at the University of Sydney, Australia where the final version of this paper was completed. MÖ wishes to acknowledge financial support from Republic of Turkey Ministry of National Education.