Anna Lyhne Jensen

Henrik Sørensen

Lasse Rosendahl

S Y Mpo

O Tatin

Towards Simulation of Clogging Effects in Wastewater Pumps: Modelling of Fluid Forces on a Fiber of Bonded Particles using a Coupled CFD-DEM Approach

Keywords:

. The results show that fiber motion similar to the motion presented in literature is obtained using a coupling between EDEM and ANSYS Fluent.

INTRODUCTION

Clogging in wastewater pumps, due to objects such as textiles and plastic in wastewater, causes severe problems in wastewater pumping systems. Therefore, resilience to clogging is a key requirement in the design of wastewater pumps. The design process of wastewater pumps is commonly based on empirical knowledge and expensive and time consuming experiments on prototypes. Thus, a method for simulating clogging effects will significantly facilitate the design process.

Computational Fluid Dynamics (CFD) is widely used in the design process of pumps. However, simulation of clogging effects entails simulation of the objects in the wastewater causing clogging, such as plastic bags or textile. This fact significantly complicates the simulation, and therefore simulation of clogging effects in wastewater pumps has not previously been done. Different applicable approaches for simulation of the interaction between fluid and flexible materials exist. These include the Immersed Boundary Method [START_REF] Peskin | The immersed boundary method[END_REF], the Arbitrary Lagrangian-Eulerian method [START_REF] Sathe | Modelling of fluid-structure interactions with the space-time finite elements: Contact problems[END_REF], Smoothed Particle Hydrodynamics coupled with the Finite Element Method [START_REF] Attaway | Coupling of smooth particle hydrodynamics with the finite element method[END_REF] and coupled CFD-DEM, which has been discussed in [START_REF] Jensen | Towards simulation of clogging effects in wastewater pumps: A review of the state-of-the-art in cloth modelling and challenges in the simulation of clogging effects[END_REF]. Advantages of using the CFD-DEM approach are that selfcollision of the textile can be handled by the DEM soft-sphere collision model. Furthermore, it is advantageous that dynamic meshing is not required, since the DEM particles are not resolved by the CFD mesh. In the present study, a flexible fiber in shear flow is modeled using the DEM software EDEM from DEM Solutions in a two-way coupling to the CFD software ANSYS Fluent including fibre self-interactions. This model is the first step towards developing a CFD-DEM framework for simulation of flexible pieces of textile and the interaction of the textile with fluid applicable for simulation of clogging effects.

Bonded Discrete Elements DEM was introduced by Cundall and Strack (1979) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] for the simulation of solid flows with mechanical transport. The method computes the motion of a large number of particles of finite size using Newtons 2. law of motion as given in Eq. ( 1) and ( 2)

m i a i = F i,n + F i,t + F i,b , (1)
I i α i = r i,c × F i,t , (2)
where m i is the mass and a i is the acceleration of the i th particle. F i,n and F i,t are the normal and tangential particleparticle contact forces and F i,b is additional forces as gravity and fluid force. In Eq. 2, I i and α i are the moment of inertia and the angular acceleration of the i t h particle respectively, and r i,c is a vector from the center of the particle to the point where the force F i,t acts. When modeling a fiber as a chain of bonded spherical DEM particles, a second model for particle-particle interaction of the bonded particles must be introduced. Since the fiber is not completely flexible, normal and shear stiffness of the fiber are introduced in the model of the bonds between the particles. Thus, bonded particles resist normal and tangential displacements. Furthermore, the bond breaks if a specified normal or tangential limit for stress is exceeded. Existing bond models in DEM have primarily been developed for simulation of rock and concrete. This also applies to the build in bond model in the software EDEM, "Hertz-Mindllin with bonding", which is a modification of the bond model for rock proposed by [START_REF] Potyondy | A bonded-particle model for rock[END_REF] [START_REF] Potyondy | A bonded-particle model for rock[END_REF]. [START_REF] Guo | Validation and time step determination of discrete element modelling of flexible fibers[END_REF] [START_REF] Guo | Validation and time step determination of discrete element modelling of flexible fibers[END_REF] used DEM to model flexible fibers. The bond model used by [START_REF] Guo | Validation and time step determination of discrete element modelling of flexible fibers[END_REF] was also based on the bond model proposed [START_REF] Potyondy | A bonded-particle model for rock[END_REF] and the results showed good correspondence between theory and model. Further validation of the bond model in EDEM for the use for flexible fibers is necessary, but out of scope of this paper.

Fig. 1 illustrates a bond between two particles of a fiber. A bond is formed if the contact radii r c of two particles overlap. The radius r is the particle radius and is used for force calculation between particles.

bond Figure 1. The particle radius, r, and the contact radius r c is defined when modeling bonded particles. Bonds are formed if there is particle contact within the radius r c . Bond properties include normal and shear stiffness of the bond and critical shear and normal stress break criteria.

Flexible Fibers in Shear Flow

The DEM fiber model will be coupled to a CFD simulation of shear flow as sketched in Fig. 2. Previous work on fibers in shear flow includes [START_REF] Forgacs | Particle motions in sheared suspensions[END_REF] [START_REF] Forgacs | Particle motions in sheared suspensions[END_REF], who investigated the motion of flexible wood-pulp fibers in shear suspensions experimentally. More recently, the motion has been modeled by [START_REF] Yamamoto | A method for dynamic simulation of rigid and flexible fibers in a flow field[END_REF] [START_REF] Yamamoto | A method for dynamic simulation of rigid and flexible fibers in a flow field[END_REF] and [START_REF] Schmid | Simulations of fiber flocculation: Effects of fiber properties and interfiber friction[END_REF] [10] by dividing the fiber into spherical and cylindrical bonded segments respectively, and modelling the hydrodynamic force on these segments. [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF] [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF] also modeled fibers as chains of cylindrical fiber segments. In this model the flow was modeled by 3D incompressible Navier-Stokes equations using a two-way coupling between the fiber and the fluid. The model took into account both viscous and dynamic drag forces. An important conclusion of the work was that fiber motion cannot be described without taking two-way coupling into account. Based on this, a two-way coupling between the fluid and the fiber is implemented in the present work.

Different orbit types depending on the flexibility of the fiber was characterized by [START_REF] Forgacs | Particle motions in sheared suspensions[END_REF] [START_REF] Forgacs | Particle motions in sheared suspensions[END_REF]. These are rigid orbits, springy orbits, snake orbits and coiled orbits. These orbit types are sketched in Figure 3. The stiffness of the fiber increases from the coiled orbit type down to the rigid orbit type. Notice from the figure that rigid, springy and snake orbits are periodic. The present simulation of fibers in shear flow in EDEM and ANSYS Fluent is validated by qualitatively comparing the results to the different orbit regimes by [START_REF] Forgacs | Particle motions in sheared suspensions[END_REF], [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF]. [START_REF] Forgacs | Particle motions in sheared suspensions[END_REF], [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF].

The diameter of fibers in pulp and paper industry is in the order of 10 -5 m. The aim of this work is to simulate one yarn of a floor cloth which has a diameter in the order of 10 -3 m. Therefore the size scale differs significantly, but it is expected that the yarn will orbit in a way similar to that of the fibers. The fiber stiffness, density and the number of particles forming the fiber are changed, and the influence of these changes on the fiber motion is investigated to assess the applicability of the present bond model and coupling in EDEM and ANSYS Fluent, to simulation of flexible fibers or textiles in fluid flow.

METHODS

Navier-Stokes equations for an incompressible fluid with particles are given in Eq. (3) and Eq. ( 4) [START_REF] Kloss | Models, algorithms and validation for opensource dem and cfd-dem[END_REF]

∂α f ∂t + ∇ • α f u f = 0 , (3)
∂ α f u f ∂t + ∇ • α f u f u f = -α f ∇ p ρ f -R pf + ∇ • τ , (4)
where α f , u f and ρ f is the is the fluid volume fraction, fluid velocity and fluid density respectively. τ = v f ∇ u f is the fluid

Domain and Boundary Conditions

The domain and boundary conditions used in CFD are presented in Fig. 5. Periodic boundary conditions are used in the right and left part of the domain and symmetry is used on the front and back of the domain. The moving walls in the upper and lower part of the domain create a shear flow as sketched in Fig. 2. The velocity profile is linear and is given by Eq. 5

v x = a • y , (5)
where v x is the velocity in the x-direction and y is the position on the y-axis. A constant value of a = 2.5 is used in all the described simulations, and the role of shear rate will not be investigated. Furthermore, gravity is not included in the simulation.

The domain is meshed with hexahedral cells of equal size. The ratio between the particle volume and cell volume is approximately 1/15.

Fiber Model

The model of the mechanics of fibers include particle-particle interaction through bonds or collision, and drag force on the fiber. Morphological properties of the fiber such as length and diameter as well as the fiber stiffness and density are expected to have a major influence on resulting fiber orbit types. These parameters are varied and the influence on the fiber motion is evaluated. The diameter of the spherical particle segments forming the fiber is 2.5 mm, to resemble the size of the yarns of a woven floor cloth. The contact radius is set to be 20 % larger than the particle radius. Furthermore, the number of particle segments forming the fiber is 20 or 50, resulting in aspect ratios of 20 and 50 respectively. Fiber aspect ratios between 50 and 780 were investigated by [START_REF] Forgacs | Particle motions in sheared suspensions[END_REF], [START_REF] Schmid | Simulations of fiber flocculation: Effects of fiber properties and interfiber friction[END_REF] and [START_REF] Salinas | Bending and breaking of fibers in sheared suspensions[END_REF].

Drag Model

Using the CFD-DEM approach, the bonded particles are not resolved by the CFD mesh. Therefore a drag law must be ascribed to the particles. The coupling framework for EDEM and Fluent includes a freestream drag formulation as well as formulations of drag developed by [START_REF] Ergun | Fluid flow through packed columns[END_REF] [START_REF] Ergun | Fluid flow through packed columns[END_REF] and [START_REF] Wen | Mechanics of fluidization[END_REF] [START_REF] Wen | Mechanics of fluidization[END_REF] and di Felice (1994) [START_REF] Di Felice | The voidage function for fluid-particle interaction systems[END_REF] for fluid particle interaction in packed and fluidized beds. No drag formulation for the modeling of bonded particles is available. In the freestream drag formulation, the drag coefficient depends on the Reynolds number. The Reynolds number is based on the particle diameter, d p , and is calculated as given in Eq. 6.

Re = ρ f • v rel • d p µ f ( 6
)
Where v rel is the relative fluid velocity, and ρ f and µ f are the fluid density and viscosity respectively. The drag coefficient is calculated based on the classical formulation of drag on a sphere, using Eq. 7.

c d =           
24/Re, if Re ≤ 0.5; 24 • 1 + 0.15 • Re 0.687 /Re, if 0.5 < Re ≤ 1000; 0.44, otherwise.

(7) Using the drag coefficients determined for each particle by Eq. 7, the drag force on each particle of the fiber can be calculated by Eq. 8

F d = 0.5 • c d • ρ f • A • v rel , ( 8
)
where A is the particle area. The viscosity of the water is 0.001003 kg/m•s and the density is 998.2 kg/m 3 . The drag force on the particles is used to determine the particle-fluid momentum exchange.

RESULTS AND DISCUSSION

13 different fibers were simulated by varying the fiber density, bond shear and normal stiffness and number of fiber segments. The initial orientation of the fiber is vertical, perpendicular to the flow direction. The results are presented in Fig. 6 for a fiber formed by 50 segments and in Fig. 7 for a fiber formed by 20 segments. The figures show the initial state of the fiber and 10 following states. The time interval between the fibers in each case is 1 s. The properties of the fiber in each of the cases shown in Fig. 6 are listed in Tab. 1. Note that since the fiber consists of spherical particles, the aspect ratio of the fiber is equal to the number of particle segments.

The movement of the fibers in case 1-3 is very similar, especially for the first 4 s. Common for the three fibers is the density of 5000 kg/m 3 , which gives a density ratio between the fibers and the fluid close to 5. The difference between the three cases is the fiber stiffness. The fiber in the first case is the most flexible. This is evident from the results where the shear flow causes the fiber to continue the rotation in an orbit which resembles a cross between the snakelike motion and the coiled motion sketched in Fig. 3. The least flexible fiber (case 1), stops rotating in the area of low fluid velocity in the center of the domain.

A standard floor cloth has a dry density of around 150 kg/m 3 . When submersed in water the cloth will absorb water and the wet density of the cloth will be higher than the dry density. In case 4-6 it is assumed that the wet density of the cloth is around 500 kg/m 3 , and in the 7 th case the fiber and the fluid is assumed to have the same density. The results of case 4-7 are very similar, indicating that the difference in density from 500 kg/m 3 to 1000 kg/m 3 has a limited influence on the results. However, comparing to cases 1-3, with a fiber density of 5000 kg/m 3 and same range of fiber stiffness, the difference due to change in density is obvious, as the motion of the fibers in case 4-7 resembles the orbit of a rigid fiber until the fiber is parallel to the flow direction. After this, the fiber motion is similar to a coiled orbit with self entanglement. Case 4-6 differ in the fiber stiffness from 1e 7 to 1e 9 . But this change in fiber stiffness is not evident from the results.

The motion of the fibers formed by 20 fiber segments, illustrated in Fig. 7, is very similar to the fiber formed of 50 segments for the first 5 s. Both for a fiber density of 500 and 5000 kg/m 3 . This indicates that the number of segments forming the fiber is of small importance in this case. However, there are significant differences between fibers of 20 and 50 segments after 5 s of the simulation. This may be caused by a significant change in fiber aspect ratio, when reducing the number of fiber segments. The diameter of the fiber is kept constant at 2.5 mm in all simulations and consequently the aspect ratio of the fiber changes when the number of spherical segments forming the fiber is reduced from 50 to 20. A more direct comparison of the influence of fiber density and stiffness is presented in Fig. 8 and9. The figures illustrate the shape of the fiber after 4 s of the simulation. It is seen from both Fig. 8 and 9 that the fiber stiffness slightly changes the motion of the fiber. Furthermore, a significant difference is evident with a change of fiber density from 500 to 5000 kg/m 3 .

CONCLUSION

The software EDEM for DEM and the software ANSYS Fluent for CFD were used to simulate fibers in shear flow with a two-way coupling between fiber and fluid, to identify the importance of the fiber properties on the motion of the fiber in this flow condition.

A change in fiber density from 500 to 5000 kg/m 3 leads to a significant change in the fiber motion in shear flow. With a fiber density of 5000 kg/m 3 the motion of the fiber resembles the springy or snakelike orbit type reported in literature. In the case of a fiber density of 500 kg/m 3 the motion is more similar to rigid orbiting. However, in the case of 50 particle segments, this motion changes to coiled orbiting.

Changing the stiffness of the bonds between the particles from 1e 7 to 1e 9 only has a very limited influence on the results. A change in the number of particles forming the bond, from 50 to 20 gives similar results in the first half of the simulated time, but the difference between the results may be caused by the change in fiber aspect ratio.

This study has presented a qualitative comparison of fiber motion of fibers with different properties as well as comparison to fiber orbit types presented in literature. In order to determine if the built-in bond model and coupling in EDEM and ANSYS Fluent is applicable for clogging simulations, careful experimental validation of the motion of fibers with specific properties is needed. Furthermore, settling experiments or similar must be conducted for a model of a piece of textile. However, the presented results of fiber motion in shear flow, are similar to results presented in literature. This is promising for the future application of coupled CFD-DEM to model fiber motion and later simulation of clogging effects in wastewater pumps.

FUTURE WORK

The presented simulation of fibers in shear flow is the first step towards developing a framework for simulation of clogging effects in wastewater pumps. In order to accomplish this, experimental validation of both the coupling between CFD and DEM and the bond model in DEM for this application must be conducted. Based on the results of the validation, development and implementation of a new formulation of drag and DEM bonds may be necessary.

Figure 2 .

 2 Figure 2. Flexible fiber orbiting in shear flow.

Figure 3 .

 3 Figure 3. Sketch of four different orbit types of fibres in shear flow: rigid, springy, snakelike and coiled [8],[11].

Figure 4 .

 4 Figure 4. Block diagram of calculation process in a coupled CFD-DEM simulation. The blocks to the left represent the CFD solver and the blocks to the right represent the DEM solver.

Figure 5 .

 5 Figure 5. CFD domain. The upper and lower boundaries are moving walls with no slip. This creates a shear flow in the domain. The right and left part of the domain has periodic boundary conditions. The fiber of DEM particles orbits in the center of the domain.

Figure 6 .

 6 Figure 6. Development of the fiber shape with time in the different cases. The time interval between each image in a case is 1 s.

Figure 7 .

 7 Figure 7. Development of the fiber shape with time in the different cases. The time interval between each image in a case is 1 s.

Figure 8 .Figure 9 .

 89 Figure 8. Comparison of fibers formed by 50 segments after 4 s.

Table 1 .

 1 Properties of the simulated fibers.

			Fiber Simulation Cases
	Case Density	Stiffness Segments Diameter Length
		[kg/m 3 ] [N/m 3 ] [-]	[m]	[m]
	1	5000	1e 9	50	0.0025	0.125
	2	5000	1e 8	50	0.0025	0.125
	3	5000	1e 7	50	0.0025	0.125
	4	500	1e 9	50	0.0025	0.125
	5	500	1e 8	50	0.0025	0.125
	6	500	1e 7	50	0.0025	0.125
	7	1000	1e 9	50	0.0025	0.125
	8	5000	1e 9	20	0.0025	0.050
	9	5000	1e 8	20	0.0025	0.050
	10	5000	1e 7	20	0.0025	0.050
	11	500	1e 9	20	0.0025	0.050
	12	500	1e 8	20	0.0025	0.050
	13	500	1e 7	20	0.0025	0.050

ACKNOWLEDGMENTS

This work is funded by Grundfos in cooperation with the Department of Energy Technology in Aalborg University.


Marie Cecilie Pedersen
email: marie.ceciliepedersen@vattenfall.com

Henrik Sørensen

S Y Mpo

O Tatin

Towards a CFD Model for Prediction of Wind Turbine Power Losses due to Icing in Cold Climate

Keywords:

Icing induced power losses is an important issue when operating wind turbines in cold climate. This paper presents a concept of modelling ice accretion on wind turbines using Computational Fluid Dynamics (CFD). The modelling concept works towards unifying the processes of modelling ice accretion and the aerodynamic analysis of the iced object into one CFD-based icing model. Modelling of icing and obtaining ice shapes in combination with mesh update by surface boundary displacement was demonstrated in the paper. It has been done by expressing in-cloud icing in CFD by an Eulerian multiphase model, implementing an icing module into the CFD solution and finally by surface boundary displacement also included in the CFD solution. The model has been developed using ANSYS-FLUENT and user-defined functions. The naca profile, NACA64618, has been used to illustrate the functionality of the model. Running ice accretion for different meteorological boundary conditions has illustrated the capabilities of the model and the generated ice shapes showed agreement with the literature.

INTRODUCTION

Wind power in cold climate has become more common especially in Nordic countries, such as Sweden, Finland, Norway and Canada. The current total installed capacity of wind power in cold climate is approximately 70 GW, which is expected to increase up to 70% towards the end of 2017 [START_REF] Wallenius | Special Theme: Cold Climate Turbines[END_REF]. Cold climate areas are defined by experiencing temperatures below the operating limits of standard wind turbines and by the presence of icing events. Cold climate is attractive for wind energy extraction due to the low air density, remote location and most likely sparse population. Nevertheless, producing electricity in cold climate regions is challenging due to icing related risks, such as ice being shed from the turbine blade, blade fatigue due to additional ice load and finally production losses due to icing. In remote cold climate areas the most severe risk is the production losses, which can have tremendous impact on the profitability of a wind farm.

Production losses in cold climate happens during atmospheric icing, by the impact of moisture on a surface at temperatures below 0 o C resulting in ice accretion. The interest of this work is periods where the conditions are in favour of building up, which is called meteorological icing. Ice building up during meteorological icing is referred to as ice accretion. Opposite, periods where the accreted ice remains on the structure is called instrumental icing [START_REF]Athmospheric Icing on strutures[END_REF]. Atmospheric icing can be divided into 3 main categories; in-cloud icing, freezing rain and wet snow icing. During in-cloud icing super-cooled droplets are present in a cloud or fog at cold temperatures resulting in the formation of hard rime, soft rime or glaze ice depending on the temperature [START_REF]Athmospheric Icing on strutures[END_REF]. In this work in-cloud icing at cold temperatures and thereby the formation of rime ice has been modelled.

Modelling of icing has been studied for power lines, within aeronautics, for structures and objects and more recently also for wind turbines. Within wind power, the aeronautic icing models such as Messinger's thermodynamic model from 1953 [START_REF] Bernard | Equilibrium Temperature of an Unheated Icing Surface as a Function of Air Speed[END_REF] and the panel based code LEWICE [START_REF] Wright | User's Manual for LEWICE Version 3[END_REF] have been widely used and accepted. Also the CFD code FENSAP-ICE originally developed for the aircraft industry has been applied to wind turbine applications, as seen in [START_REF] Reid | FENSAP-ICE Simulation of Icing on Wind Turbine Blades, Part 1: Performance Degradation[END_REF] where performance degradation and power losses were studied. Other icing models highly used within wind power are the empirical tuned cylinder based model by Makkonen [START_REF] Makkonen | Models for the Growth of Rime, Glaze, Icicles and Wet Snow on Structures[END_REF] [2] and the panel based code TURBICE designed for de-icing applications [START_REF] Makkonen | Modelling and prevention of ice accretion on wind turbines[END_REF]. Recent studies by [START_REF] Davis | Icing Impacts on Wind Energy Production[END_REF] has managed to advancing Makkonen's cylinder based model to be more representative for a wind turbine. Common for LEWICE, TURBICE and Makkonen's cylinder based model is that the codes generate ice shapes. To study aerodynamic response of wind turbines during or after icing conditions the ice shape could be obtained of one of the latter mentioned methods followed by the use of CFD to study the aerodynamic changes. This was seen in [START_REF] Etemaddar | Wind turbine aerodynamic response under atmospheric icing conditions[END_REF] where LEWICE was used to generate icing shapes for a variation of meteorological conditions followed by a aerodynamic study using ANSYS-FLUENT. Thus, obtaining the ice shape was separated from the aerodynamic analysis, which was also seen in recent studies by [START_REF] Ali | In-Cloud Ice Accretion Modeling on Wind Turbine Blades Using Extended Messinger Model[END_REF]. In [START_REF] Villalpando | Numerical study of flow around iced wind turbine airfoil[END_REF] and [START_REF] Chi | CFD analysis of the aerodynamics of a business-jet airfoil with leading-edge ice accretion[END_REF] the ice shapes were obtained from experiments followed by CFD analysis of the aerodynamic response. The objectives of this work is to develop a concept, which unifies the ice accretion process and the aerodynamic analysis into one CFD-based ice accretion model.

METHODS

The aim of the ice accretion model is to mimic an icing event in the best possible way by employing a CFD-based frame and to make the model as generic as possible. The Eulerian multiphase flow model has been used to imitate icing conditions in various CFD-based icing models. This is seen by the commercial code FENSAP-ICE, originally developed for the aircraft industry [START_REF] Nti | User Manual FENSAP-ICE[END_REF], which has recently been applied for wind turbine applications [START_REF] Switchenko | FENSAP-ICE Simulation of Complex Wind Turbine Icing Events, and Comparison to Observed Performance Data[END_REF] and [START_REF] Reid | FENSAP-ICE Simulation of Icing on Wind Turbine Blades, Part 1: Performance Degradation[END_REF]. In other studies using ANSYS-FLUENT the Eulerian multiphase flow model has also been applied to imitate in-cloud icing conditions for ice accretion [START_REF] Wirogo | An Eulerian Method to Calculate the Collection Efficiency on Two and Three Dimensional Bodies[END_REF] and [START_REF] Mare | Preliminary modelling study of ice accretion on wind turbines[END_REF]. From the droplet solution of the Eulerian multiphase model, the variables of velocity, volume fractions, density and temperature for the phases will be available everywhere in the computational domain at any time.

The ice accretion model used in this study was inspired by Messiger's model [START_REF] Bernard | Equilibrium Temperature of an Unheated Icing Surface as a Function of Air Speed[END_REF], which has been applied for other CFD icing applications especially for aircraft purposes [START_REF] Nti | User Manual FENSAP-ICE[END_REF]. Nevertheless, the model in combination with an Eulerian multiphase flow model has also shown reasonable results and applications for wind turbines [START_REF] Ali | In-Cloud Ice Accretion Modeling on Wind Turbine Blades Using Extended Messinger Model[END_REF] and [START_REF] Switchenko | FENSAP-ICE Simulation of Complex Wind Turbine Icing Events, and Comparison to Observed Performance Data[END_REF]. The modelling approach is structured as illustrated by Figure 1, which consists of 3 main modules; 1) obtain multiphase flow solution, 2) ice accretion model and 3) surface boundary displacement. The ice accretion model and the surface boundary displacement is included in the ANSYS-FLUENT solver by using user-defined functions. Unique for the proposed modelling structure in Figure 1 is that the three modules are updated every time-step adding a high accuracy to the model. Updating the multiphase flow solution after every surface boundary displacement ensures the changes in geometry to be reflected in the flow solution. Other icing models, such as FENSAP-ICE [START_REF] Nti | User Manual FENSAP-ICE[END_REF] either uses an initial particle flow solution over the entire ice accretion time or uses a so-called quasi-steady approach, where the droplet solution are updated a given number of times during the ice accretion time [START_REF] Nti | User Manual FENSAP-ICE[END_REF]. Since the droplet solution and collection efficiency ( β) will change over time, see example in Figure 2, it is found necessary to update the Eulerian flow solution every time-step to ensure a reliable representation of the actual multiphase flow of air and particles at the given time of the ice accretion period. Figure 2 shows the change in collection efficiency ( β, Equation 4) over time and the corresponding shape of the airfoil exposed to icing over time.

Multiphase flow modelling

The Eulerian description of a multiphase flow involves the concept of phasic volume fractions, α q , which describes the space occupied by each phase. The laws of conservation of mass and momentum are satisfied by each phase individually [START_REF]ANSYS Fluent Theory Guide[END_REF]. The sum of all phases equals to one

n q=1 α q = 1 (1)
The volume fraction of each phase of the two-phase flow is solved by the continuity equation, shown here for phase q

∂ ∂t (α q ρ q ) + ∇ • (α q ρ q v q ) = n p=1 ṁpq -ṁqp + S q (2)
where v q is the velocity vector of phase q and m pq,q p characterize the mass transfer between the p th and q th phase. There is no mass exchange between the phases and no mass source are added for the phases, leaving the right hand side of Equation 2 to zero. The solution of Equation 2 for the secondary phase together with the conditions of Equation 1 allows for the calculation of the primary-phase volume fraction. The conservation of momentum for the fluid phase q is given as

∂ ∂t (α q ρ q v q ) + ∇ • (α q ρ q v q v q ) = -α q ∆p + ∆ • τq + α q ρ q g+ + n p=1 (K pq ( v p -v q )) + ( F q + F lift,q + F wl,q + F td,q ) ( 3
)
where g is the acceleration due to gravity, and K pq is the momentum exchange coefficient between the phases, ∆p, τq , F q , F lift,q , F wl,q , F td,q are the pressure gradient, stress-strain tensor, external body force, lift force, wall lubrication force and the turbulent dispersion force respectively. The energy conservation for the Eulerian multiphase model is solved by a separate enthalpy equation for each phase of the flow [START_REF]ANSYS Fluent Theory Guide[END_REF].

The k-ω SST turbulence is applied to the multiphase flow based on a preliminary study [START_REF] Mare | Preliminary modelling study of ice accretion on wind turbines[END_REF].

Ice accretion model For the simulation of ice accretion only rime ice has been considered. Figure 3 illustrates the mass balance of rime ice accretion on the surface boundary of an object, where m im p is the mass of droplets impinging on the object, m subl is sublimation of ice mass from solid to vapor and m ice is the accumulated ice mass. It is assumed that all mass which impinge on the body will freeze instantaneously and turn into ice, see Equation 5. In the model, the accreted ice mass becomes part of the object by applying a surface boundary displacement algorithm. Enclosing the accumulated mass of ice into the object geometry allows for studying the aerodynamic changes due to ice accretion. The mass of accumulated ice m ice , is found instantaneously over time based on the collection efficiency ( β)

β = α n (u d • n) U ∞ (4) m im p ≈ m ice = βU ∞ LWC (5)
where α n is the normalised droplet volume fraction

α d, i m p α ∞
, u d is the droplet impact velocity, n is the unit surface normal and U ∞ is the free-stream velocity. The normalised velocity u n is given as u d u ∞ . Thus the amount of accumulated ice mass is a function of the collection efficiency, the liquid water content in the cloud and the free-stream velocity.

Surface boundary displacement

The mass of instantaneously accumulated ice (m ice ) initiates the surface boundary displacement. The surface boundary displacement is constructed around the Dynamic-mesh package available in ANSYS-FLUENT [START_REF]ANSYS Fluent Theory Guide[END_REF], from where the DEFINE -GRID -MOTION macro has been used. The macro is transient and by an iterative process the node point of the surface can be updated. This means, that the mesh will be updated at every time step. The mesh is updated by a node displacement vector v n,i

v n, i = m ice ρ ice n (6)
where ρ ice is the density of ice. The node displacement vector gives a set of coordinates of the new location of the node. The process is obtained by a face-looping and node-looping approach as illustrated in of ice and the boundary displacement vector is calculated at each face (f) of the surface boundary as illustrated in Figure 5. By linear interpolation the contribution from the node neighbouring faces (f L , f R ) to the node displacement vector at node i (n, i) is enabled. The node positions are updated according to the following expression

p t+∆t n, i = p t n, i + v t n, i ∆t (7)
where p n, i is the current node positions, t is the current time and ∆t the time step.

n i n i+1 n i-1 f R f L vfL vfR vn,i

Model Set-up

A set of simulations have been carried out to illustrate the capabilities of the CFD icing model. The NACA64618 profile with a sharp trailing edge has been used for the simulation in a computational domain consisting of 4050 cells, in a C-grid shape, as seen in Figure 6. The NACA64618 profile was used, since it is part of the 5MW virtual NREL wind turbine also used for other icing studies [START_REF] Reid | FENSAP-ICE Simulation of Icing on Wind Turbine Blades, Part 1: Performance Degradation[END_REF] and [START_REF] Etemaddar | Wind turbine aerodynamic response under atmospheric icing conditions[END_REF]. The model inputs are presented in the following by specifying the conditions for the secondary phase (the super-cooled droplets) and for the primary phase (the airflow). Secondary phase inputs From the liquid water content and density the volume fraction of the secondary phase have been specified as

α d = LWC ρ d (8)
where ρ d is the density of the super-cooled droplets. Simulations of ice accretion have been carried out for the different cases of secondary phase settings presented in Table 1.

Meteorological inputs

The meteorological conditions of the rime ice events are inspired by [START_REF] Etemaddar | Wind turbine aerodynamic response under atmospheric icing conditions[END_REF] and presented in Table 2. The droplet size allows for a one-way coupling to exist between the airflow and the droplets.

Surface roughness

The surface roughness is a significant parameter to include in an ice model since it determines the transition location of the airflow from laminar to turbulent over the rough surface. The sand grain roughness depends

Meteorological boundary conditions

Phases: primary dispersed (airflow) (droplets)

U ∞ 25 25 m s T ∞ -10 -10 o C
on the actual icing condition seen by the iced object. A way to take roughness into account is by employing a sand grain roughness (K s ) [START_REF] Makkonen | Modelling and prevention of ice accretion on wind turbines[END_REF] and [START_REF] Shin | Prediction of ice shapes and their effect on airfoil drag[END_REF] based on empirical roughness parameters as During ice accretion the sand-grain roughness will change affecting the transitions area and flow field at the boundary of the iced object. In [START_REF] Makkonen | Modelling and prevention of ice accretion on wind turbines[END_REF] the time dependence of the surface roughness height is taken into account by multiplying by a factor depending on the freezing fraction ( f ), the fraction of particles which actually freeze at impact and the number of ice layers simulated. The freezing fraction is a number between 0-1, where 1 correspond to rime ice and 0 < f > 1 is glaze ice conditions, [START_REF] Battisti | Wind Turbines in Cold Climate[END_REF]. For f close to or equal to 1 the sand grain roughness is multiplied with 1 as soon as a reasonable rime ice layer is obtained. Thus, since only rime ice accretion is simulated in this study it is found reasonable to apply a constant sand-grain roughness based on Equation 9 from time equal to zero, even though this will enhance the ice accretion at the initial stage.

K s = k * s k * sb,LWC k * s k * sb,T ∞ k * s k * sb,U ∞ k * s k * sb,MVD k * sb c (

RESULTS AND DISCUSSION

The focus of this study has been on developing and testing the functionality of the presented CFD ice accretion model on the NACA64618 airfoil. The generated ice shapes have been compared to ice shapes presented in the work by Etemaddar et al. 2014 [START_REF] Etemaddar | Wind turbine aerodynamic response under atmospheric icing conditions[END_REF], where ice accretion has been simulated on the NACA64618 airfoil using the panel-based code LEWICE [START_REF] Wright | User's Manual for LEWICE Version 3[END_REF].

From [START_REF] Etemaddar | Wind turbine aerodynamic response under atmospheric icing conditions[END_REF] it was shown, that increasing LWC would increase ice accretion, especially in the direction of the chord. Likewise, increasing MVD would increase ice accretion but more severe along the thickness of the airfoil. In this study, ice accretion has been simulated for 20 minutes with a time-step size of 0.01 second for all cases of secondary phase settings presented in Table 1 with the meteorological boundary conditions given in Table 2. Furthermore a constant density of ice of 650 k g m 3 has been used an angle of attack of 0 o . The generated ice shapes is seen for the cases of LWC of 0.05 g m 3 , 0.10 g m 3 and 0.15 g m 3 in Figure 7 for cases with LWC of 0.20 g m 3 and 0.25 g m 3 in Figure 8. Figure 7 and 8 illustrate how the total amount of ice accretion increased with an increasing LWC. Furthermore, the figures illustrates that increasing MVD will result in more ice accretion, which was most severe for the high numbers of LWC. This can be seen in Figure 8. The results from the simulations was in agreement with the conclusions regarding LWC and MVD found in [START_REF] Etemaddar | Wind turbine aerodynamic response under atmospheric icing conditions[END_REF] and [START_REF] Matthew C Homola | Effect of atmospheric temperature and droplet size variation on ice accretion of wind turbine blades[END_REF].

Conclusions

From the study a promising concept towards a complete CFD model for prediction of wind turbine power losses due to icing in cold climate was established. The study has illustrated that it was possible to unify the ice accretion process with the aerodynamic analysis in one CFD ice accretion model. The concept has demonstrated the ability to model rime ice accretion and perform surface boundary displacement per time step for an icing period of 20 minutes. Furthermore, the simulated ice accretion cases showed agreement with the literature and the model was able to include the influence of droplet variation and variation in of the liquid water content.

Future work is to extend the icing time and include varying meteorological input conditions and finally to compare the model with site-measurements or data from a climatic wind tunnel, as seen by [START_REF] Peter M Blasco | Investigation of Wind Turbine Power Generation During Atmospheric Icing by Multi-Disciplinary Experimentation[END_REF]. Furthermore, changes in lift and draft forces for a variation of angle of attacks should be studied.

Figure 1 .

 1 Figure 1. Modelling structure.

Figure 2 .

 2 Figure 2. Top: variation in collection efficiency ( β) over 20 minutes of ice accretion. Bottom: corresponding ice shapes during 20 minutes of ice accretion.

Figure 3 .

 3 Figure 3. Mass balance of rime ice accretion on surface boundary.

Figure 4 .Figure 4 .

 44 Figure 4. Surface boundary displacement.

Figure 5 .

 5 Figure 5. Surface boundary displacement by surface displacement vectors.

Figure 6 .

 6 Figure 6. Computational domain.

9 )

 9 where c is the cord length in meters. The empirical roughness parameters are defined as follows k * s k * sb,LWC =(0.5714 + 0.2457(LWC) + 1.2571(LWC) 2 ) -0.0333(MVD), MVD > 20 k * sB =(0.001177)

Figure 7 .

 7 Figure 7. Generated ice shapes after 20 minutes of ice accretion using surface boundary displacement for U ∞ = 25 m s , T ∞ = -10 o C and LWC of 0.05 g m 3 , 0.10 g m 3 and 0.15 g m 3 .

Figure 8 .

 8 Figure 8. Generated ice shapes after 20 minutes of ice accretion using surface boundary displacement for U ∞ = 25 m s , T ∞ = -10 o C and LWC of 0.20 g m 3 and 0.25 g m 3 .

Table 1 .

 1 Table of CFD model inputs: secondary phase.

	Secondary phase conditions
	LWC g m 3	MVD µm	α d -
	0.05 8, 12, 16, 20 7.69 • 10 -8 0.1 8, 12, 16, 20 1.54 • 10 -7 0.15 8, 12, 16, 20 2.31 • 10 -7 0.20 8, 12, 16, 20 3.08 • 10 -7 0.25 8, 12, 16, 20 3.85 • 10 -7

Table 2 .

 2 Table of CFD model inputs: meteorological inputs.

ACKNOWLEDGMENTS

Marie Cecilie Pedersen would like to thank Vattenfall Vindkraft A/S for hosting and partly granting the Industrial Phd program and the Innovation Fund Denmark from the Danish Ministry of Higher education and Science for partly granting the program.


Alessandro Nocente
email: alessandro.nocente@ntnu.no

Tufan Arslan

Daniel Jasiński

Torbjørn K Nielsen

S Y Mpo

O Tatin

A Study of Flow inside a Centrifugal Pump: High Performance Numerical Simulations Using GPU cards

Keywords:

The present work reviews calculations of a steady three-dimensional (3D) flow past a centrifugal pump with its diffuser channels using different hardwares. The open source CFD software OpenFOAM has been ported to the GPU platform with double precision. The performance of a high-end GPU was reported. The considered computational domain is one of the three stages, since each has exactly the same design. In this work, simulations were carried out at the best efficiency point (BEP) for single phase flow with four different turbulence models. The results were compared with the performance report from the manufacturer. Only the realizable k-ε turbulence model has been taken into account for performance tests. The comparisons at BEP showed that integral quantity based results were not sensibly influenced by the turbulence model and the steady state simulations demonstrated to be a good approximation of the LDV data, always containing the error within an acceptable limit. The ability of solving numerical problems faster makes GPU hardware attractive to be used for further investigations which will involve two-phase flow studies on the same pump.

INTRODUCTION

Computational Fluid Dynamics (CFD) simulations are nowadays very accurate in predicting flow conditions in turbomachinery. The development and diffusion of CFD codes and the growth of available computer power have hugely increased the usage of CFD in turbomachinery [START_REF] Dick | Performance prediction of centrifugal pumps with cfd-tools[END_REF], [START_REF] Vande Voorde | Performance prediction of centrifugal pumps with steady and unsteady cfd-methods[END_REF]. It is less expensive and less time consuming than an experimental investigation, and allows testing new parts design before the actual realization of a prototype. Given a certain level of accuracy and confidence of the simulation results, the computational time has a great importance. Less time-consuming simulations are always preferable. GPU cards can reduce the computational time significantly. This reduction free resources which can be used to investigate deeper the internal flow by mean of a more accurate simulations (i.e. a finer mesh for the same domain). In this paper, we worked on a multistage centrifugal pump with closed impeller and a vaned diffuser specifically designed for offshore produced water treatment. It processes a dispersed phase flow which requires the least possible shear and a proper internal circulation [START_REF] Mark | Break-up of oil droplets in the production system[END_REF] in order to avoid the use of chemicals [START_REF] Sjöblom | Our current understanding of water-in-crude oil emulsions.: Recent characterization techniques and high pressure performance[END_REF]. Many researchers have used CFD for investigation of flow inside centrifugal pumps, most of them concentrating on pumps with volute casing. Shukla and Kshirsagar [START_REF] Shukla | Numerical experiments on a centrifugal pump[END_REF] calculated the internal flow using k-ε and k-ω turbulence models. Cheah et al. [START_REF] Cheah | Numerical flow simulation in a centrifugal pump at design and off-design conditions[END_REF] concentrated on the flow inside the impeller and Bacharoudis et al. [START_REF] Bacharoudis | Parametric study of a centrifugal pump impeller by varying the outlet blade angle[END_REF] realized a parametrical study on how performances are affected by different impellers. Shah et al. [START_REF] Shah | CFD based flow analysis of centrifugal pump[END_REF] carried over calculations on different operating points using both RNG k-ε and SST k-ω turbulence models. Among the researchers who worked specifically on vaned diffuser centrifugal pumps, Muggli et al. [START_REF] Felix A Muggli | Flow analysis in a pump diffuserpart 2: Validation and limitations of CFD for diffuser flows[END_REF] applied a Navier-Stokes code with a standard k-ε model to the study of the flow in a vaned diffuser of a high loaded centrifugal pump. Feng at al. [START_REF] Feng | Numerical investigation on pressure fluctuations for different configurations of vaned diffuser pumps[END_REF] investigated the flow and the pressure fluctuations on the impeller blade and the vaned surface of a diffuser. Zhou et al. [START_REF] Zhou | Numerical investigation of pressure fluctuation and rotor-stator interaction in a multistage centrifugal pump[END_REF] conducted a numerical investigation on unsteady flow field in a multistage centrifugal pump with focus on the pressure fluctuations while Pei et al. [START_REF] Pei | Numerical investigation of periodically unsteady pressure field in a high power centrifugal diffuser pump[END_REF] on a high power centrifugal diffuser pump for residual heat removal. Other researchers performed a comparison between the results obtained with CFD and the experimental data. Sinha and Katz [START_REF] Sinha | Quantitative visualization of the flow in a centrifugal pump with diffuser vanes -I: On flow structures and turbulence[END_REF] used PIV technique, Feng et al. [START_REF] Feng | Investigation of periodically unsteady flow in a radial pump by CFD simulations and LDV measurements[END_REF] used Laser Doppler Velocimetry (LDV) and Pedersen et al. [START_REF] Pedersen | Flow in a centrifugal pump impeller at design and off-design conditions -part I: Particle image velocimetry PIV and laser doppler velocimetry LDV measurements[END_REF] applied both experimental techniques.

The present study focuses on the characterization of the flow in the diffuser channels of a multistage centrifugal pump. One characterization of the flow passed the same multistage pump has already been achieved in previous works [START_REF] Nocente | Numerical simulation of flow inside centrifugal pump by two different solvers[END_REF], [START_REF] Nocente | Numerical prediction of a multistage centrifugal pump performance with stationary and moving mesh[END_REF]. The steady state and transient state simulation at the same working conditions have been carried on using commercial Ansys Fluent software. In this work, the open-source software RapidCFD [START_REF]RapidCFD[END_REF] has been used. The code is based on OpenFOAM [START_REF] Openfoam | [END_REF] solvers which are running on GPU cards rather than CPUs. The performance gain of GPUs over CPUs has been studied for the simulated case which is a steady three dimensional (3D) flow past the diffuser channels of a centrifugal pump.

NUMERICAL MODEL

The pump investigated in this work is a prototype realized for Typhonix AS. It is a three-stage end suction centrifugal pump for produced water. The computational domain consists of the inlet pipe, the first stage impeller and the first stage stator. After passing the stator, the fluid is addressed to the following stage by mean of return vanes. This part of the machine does not influence the performance significantly, therefore is not included in this work to save computational time. In Table 1, the main geometrical characteristics of the impeller and diffuser are listed along with head, flow rate and rotational speed at the operating point. Figure 1 shows the 3D computational domain. The computational grid is structured along the pump axis direction with hexahedral cells in the diffuser domain. The grid has no size refinement along the walls since we believe that the resolution were fine enough to take into account the wall effects on the flow. After the simulations, a check confirmed this assumption by showing that y + is in the range 30 < y + < 300. This allows us to use proper wall modelling in the calculation of the turbulence which will be mentioned in next section. The entire computational domain is shown on the left side of Fig. 2. The lower side of the figure shows a closer look to a small region of the diffuser mesh which is indicated as "A". Due to the complicated shape of the impeller vanes, the impeller grid is unstructured with tetrahedral cells. The total number of cells is about 10 millions (9,926,478) cells. The distribution of the cells in the different zones of the computational domain is summarized in Table 2. The mesh quality was checked on cell skewness which is defined as the difference between the shape of the cell and the shape of an equilateral cell of equivalent volume. The optimal value for skewness is 0. A maximum value greater than 0.95 may lead to convergence difficulties. The average and maximum values of mesh skewness for each mesh zone is shown in Table 2. The convergence behaviour will be mentioned in Section 4.

No-slip wall condition has been specified for all the wall boundaries such as pipes, hubs, shrouds, blades plus the upper and lower wall of both impeller and diffuser. At the inlet of the pipe, mass flow boundary condition is imposed. Zero pressure is imposed at the diffuser outlet in order to calculate the stage head. Turbulent intensity is considered moderate and it is imposed 5% at the inlet.

THE SOLVER RUNNING ON GPU

RapidCFD is a CUDA-ported version of a Finite Volume Method based code OpenFOAM, that runs all simulation stages on Graphic Processing Units (GPU). Running entire simulation on the GPU has a significant performance advantage over commonly used [START_REF] Malecha | Gpu-based simulation of 3d blood flow in abdominal aorta using openfoam[END_REF] GPU acceleration of solvers for linear equation systems. With the full acceleration, there is no need for constant memory copying between CPU and GPU memory space which is a huge bottleneck that significantly decreases overall simulation speed. 2. Mesh cells distribution (the parts of the computational domain can be see in Fig. 1) Porting OpenFOAM requires identification and modifications of some parts of the code that are responsible for most computationally intensive tasks, such as matrix assembly and solution. This requires creating data structures for storing all required information in GPU memory. A disadvantage of using GPUs is that all algorithms need to be parallelized. For this reason, the most efficient smoother algorithm -Gauss-Seidel -had to be replaced with a less efficient weighted Jacobi. For the same reason only relatively poor diagonal (Jacobi) preconditioner from the original implementation was ported to GPU. In order to improve performance of solvers based on Krylov subspace method, a new fully parallel preconditioner AINV [START_REF] Ilya | Algorithm for sparse approximate inverse preconditioners in the conjugate gradient method[END_REF] was implemented. This new method offers some performance improvement over the diagonal preconditioner, however it is still less effective than original Incomplete Cholesky and Incomplete LU algorithms. Despite the disadvantage of using less efficient algorithms, simulations on the GPU cards run significantly faster than on multi-core CPU. The advantage of GPU comes mostly from significantly higher memory bandwidth which is crucial for memory bound applications such as CFD.

The turbulence is simulated by four different Reynolds Averaged Navier-Stokes (RANS) models: Standard, Realizable, RNG k-ε and SST k-ω models. The pressurevelocity coupling is handled using SIMPLE scheme. Turbulent kinetic energy (k), turbulent dissipation ratio (ε) and specific dissipation ratio (ω) are discretized with a first order upwind scheme while the discretization scheme for the momentum equations is in second order. The discretized equation system is solved by Preconditioned Conjugate Gradient (PCG) methods [START_REF] Jacobs | Preconditioned conjugate gradient methods for solving systems of algebraic equations[END_REF] and Geometric Algebraic Multigrid Method (GAMG). The method which is used for steady state simulations is Multiple Reference Frame (MRF) since the rotational speed of the impeller is fixed and has a value of n = 1480 rpm. MRF, also known as Frozen Rotor (FR), considers the rotor as in a fixed position. The code is instructed on the impeller's motion and the calculation takes into account the forces that the motion imposes on the fluid due to centripetal and Coriolis acceleration.

RESULTS

In the first part of the work we focused on the turbulence modeling. The simulations are carried out at the best efficiency point (BEP) using the values shown in Table 1 and single-phase flow. Four different turbulence models have been taken into account: k-ε (standard, RNG and realizable) and SST k-ω. For each model, the flow The return vane between consecutive stages is not part of the calculation because of its computational demand, and the results demonstrate that they have negligible influence on head or efficiency. The return vanes have in fact no role in transferring or transforming energy.

In Table 3, head, torque and efficiency for each considered turbulence models are compared with the manufacturer's data. Data indicated as "Test" are to be referred to the commissioning test carried out from the manufacturer on the whole pump. The differences in the head is acceptable.

Concerning the internal flow in the diffuser channel, both the pressure and absolute velocity fields have been investigated. The overall ranges of velocities and pressures are the same, the flow structures are not affected by significant variation in the different turbulence models. In addition, there is not evident flow detachment. This means that the pump prototype is substantially well designed and explains the similarity of results in different turbulence models. According to these result, it is possible to state that for this simulation, the choice of turbulence model does not significantly affect the results. The comparison at this discharge level shows that the steady state model can approximate within 4-7% error the real flow characteristics as seen in Table 3. The same consideration can be done for the torque and the efficiency.

Because of the small influence of the turbulent model choice on this particular case, we decided to use only one of the turbulence models for the calculations. Although the standard k-ε model gave slightly more close results to the experimental data, the realizable k-ε model has been chosen because its robustness. In the meanwhile, an experimental campaign of Laser Doppler Velocimetry (LDV) measurements has been carried out in the Waterpower laboratory at NTNU, Trondheim. An experimental water loop has been designed and realized using a single stage of the pump. A special casing with plexyglass inspection window has been designed to guarantee an optical access for the laser probes. The velocity measurements have been repeated several times measuring on very fine grid of points. Hence it was possible to draw a quite accurate contour of the velocity along the mid plane of the diffuser. This contour is compared to a 2D contour obtained from the simulations in Fig. 3.

GPU PERFORMANCE

In the GPU cluster we have used, CPUs are Intel Xeon E5-2670 with 8 cores working at 2.60 GHz. The graphics cards are NVIDIA GeForce Titan X with 3072 CUDA cores and 12 GB GDDR5 memory. At each node, there are two cards installed via PCI Express connection and two CPUs. In this work, a single node is used. The performance tests were made by RapidCFD's turbulent steady state solver: simpleFoam. Both Geometric Algebraic Multigrid solver (GAMG) and Preconditioned Conjugate Gradient (PCG) solving for pressure equation were tested on both CPU and GPU in double precision.

Solving velocity and turbulent quantities is not having a significant role in computational effort therefore the solution method of pressure equation was focused more. As mentioned in the previous section, available preconditioner of PCG solver is AINV in RapidCFD (GPU) while standard OpenFoam solver (CPU) was using Diagonal Incomplete-Cholesky (DIC) which is the fastest one available. GAMG is using Gauss-Seidel smoother algorithm in OpenFOAM and Jacobi algorithm was used in RapidCFD. In the tests, simulations ran for 2 hours with several number of time steps (from 100 to 1000 based on the performance) and the first time step was discarded in the total wall time calculations. The convergence behaviour of the solver in the GPU simulations can be seen in Figure 4. As seen in Figure 5 the residuals are not fluctuating in CPU simulations. however the values for both simulations are smaller than 1e-5 for all equations.

The performance (calculated from wall clock time per iteration) for the GPU card (2 cards) over the 2 CPU unit (16 cores) is shown in Fig. 6. The GPU provided speed-up for PCG solver for two card is 3.4x which can be seen in Fig. 6. However, PCG is very slow compared to the GAMG solver for both GPU and CPU generally. In our case, GAMG shows superior performance over PCG as seen in Fig. 6. At the same time, GPU gives 2.3x speed up for GAMG solver. However, GAMG solver requires larger memory than PCG, therefore, our case could not be handled at a single GPU card when GAMG was used. The tests were performed with two cards by using MPI for data communication between the cards. The results were compared to the performance with 2 CPUs for both PCG and GAMG solvers.

CONCLUSIONS AND FUTURE WORK

This work aims to show that the new hardware platforms can be used for simulations at the industrial level and to demonstrate their importance for the future of CFD. The results demonstrates that OpenFOAM based RapidCFD shows promising speed up (around 2.6-3.5) with NVIDIA's GeForce Titan X graphics cards versus Intel Xeon E5-2670 CPU units for calculations in double precision. It brings significant advantage for simulating the flow because of the faster convergence in solving the pressure equation. However, the memory of the GPU As mentioned before in the results section, we can affirm that the simulation results are not influenced by the turbulence model used. The accuracy of the steady simulations has been evaluated in comparison with the test results provided by manufacturer. The difference in predictions of head tends to be small around the BEP, which is the design condition. As a future development, the steady state simulations will be used to simulate the flow in the pump at different load conditions and diffuser design. A dispersed phase flow will also be simulated in order study the behaviour of a dispersed phase in the diffuser for different load conditions and different dispersed phase density.

NOMENCLATURE

Figure 1 .

 1 Figure 1. Computational Domain[START_REF] Nocente | Numerical prediction of a multistage centrifugal pump performance with stationary and moving mesh[END_REF]

Figure 2 .

 2 Figure 2. Diffuser Mesh: overview (upper), closer look to the region near casing indicated as "A" (lower)

Figure 3 .

 3 Figure 3. Mean velocity magnitudes in mid-section of the diffuser: CFD (left), LDV (right) for BEP (60 m 3 /h) with Realizable k-ε model

Figure 4 .Figure 5 .

 45 Figure 4. Residuals in the simulations (GPU)

Figure 6 .

 6 Figure 6. Performance of GPU versus CPU with different matrix solvers

Table 1 .

 1 Pump Specification[START_REF] Nocente | Numerical prediction of a multistage centrifugal pump performance with stationary and moving mesh[END_REF]

	Impeller
	Number of blades	Z i	6
	Inlet diameter	D Ri	100 mm
	Outlet diameter	D Ro	264 mm
	Outlet height	h Ro	13 mm
	Diffuser
	Number of blades	Z d	10
	Inlet diameter	D Si	270 mm
	Outlet diameter	D So	478 mm
	Outlet height	h So	18 mm
	Design Operating Point
	Flow rate	Q BEP	60 m 3 /h
	Delivered head	H BEP	62 m
	Rotational speed	n BEP 1480 rpm

Table 3 .

 3 The cluster systems used in this work are Vilje [23] and Kongull Integral quantities values for the different turbulence models located in Trondheim, Norway and Abel[START_REF] Abel | [END_REF] which is located in Oslo, Norway.

		Head [m] Torque [Nm] Efficiency [%]
	Standard k-ε	21.8	32.7	70
	RNG k-ε	22.6	32.5	73
	Realizable k-ε	21.2	30.2	74
	SST k-ω	20.7	29.5	74
	TEST	21.0	32.6	66.7
	RATED	20.3	31.9	69.0

  Turbulence dissipation rate [m 2 /s 3 ]

			F R	Frozen Rotor
			GAM G Geometric Algebraic Multigrid solver
			H	pump Head [m]
			h Ro	Rotor outlet height [mm]
			h So	Stator outlet height [mm]
			k	Turbulence kinetic energy [m 2 /s 2 ]
			LDV Laser Doppler Velocimetry
			M RF Multiple Reference Frame
			n	Rotational speed [rpm]
			P CG Preconditioned Conjugate Gradient
			Q	Flow rate [m 3 /h]
	η	Efficiency [rpm]
	ω	Specific turbulence dissipation rate [m 2 /s 2 ]	T	Torque [Nm]
	AIN V Approximate Inverse	y+	Dimensionless wall distance
	BEP Best Efficiency Point	Z d	Number of diffuser blades
	D Ri	Rotor inlet diameter [mm]	Z

D Ro Rotor outlet diameter [mm] D Si Stator inlet diameter [mm] D So Stator outlet diameter [mm] DIC Diagonal Incomplete-Cholesky i Number of rotor blades

ACKNOWLEDGMENTS

The authors thank to NTNU via NOTUR (The Norwegian Metacenter for computational Science) for supporting the computing time.


C M Mazzoni
email: cosimo.mazzoni@eng.ox.ac.uk

R Ahlfeld

B Rosic

F Montomoli

Uncertainty quantification of leakages in a multistage simulation and comparison with experiments

Keywords:

The present paper presents a numerical study of the impact of tip gap uncertainties in a multistage turbine. It is well known that the rotor gap can change the gas turbine efficiency but the impact of the random variation of the clearance height has not been investigated before. In this paper the radial seals clearance of a datum shroud geometry, representative of steam turbine industrial practice, was systematically varied and numerically tested. By using a Non-Intrusive Uncertainty Quantification simulation based on a Sparse Arbitrary Moment Based Approach, it is possible to predict the radial distribution of uncertainty in stagnation pressure and yaw angle at the exit of the turbine blades. This work shows that the impact of gap uncertainties propagates radially from the tip towards the hub of the turbine and the complete span is affected by a variation of the rotor tip gap. This amplification of the uncertainty is mainly due to the low aspect ratio of the turbine and a similar behavior is expected in high pressure turbines.

.

INTRODUCTION

In order to accommodate rotor to casing relative movement during turbine operation, circumferential inlet and exit shroud cavities must be formed. These necessary alterations in the overall flow path geometry significantly modify the main blade flow field. All reported works (i.e. Wallis [START_REF] Wallis | Secondary and Leakage Flows in a Multistage Turbine[END_REF], Pfau et al. [START_REF] Pfau | Flow Interaction From the Exit Cavity of an Axial Turbine Blade Row Labyrinth Seal[END_REF], Rushton [START_REF] Rushton | High-Pressure Turbine Shroud Leakage[END_REF], Giboni et al. [START_REF] Giboni | Interaction of Labyrinth Seal Leakage Flow and Main Flow in an Axial Turbine[END_REF], Metzger and Rued [START_REF] Metzger | The Influence of Turbine Clearance Gap Leakage on Passage Velocity and Heat Transfer Near Blade Tips: Part I-Sink Flow Effects on Blade Pressure Side[END_REF][START_REF] Metzger | The Influence of Turbine Clearance Gap Leakage on Passage Velocity and Heat Transfer Near Blade Tips: Part II-Source Flow Effects on Blade Suction Sides[END_REF], Anto et al. [START_REF] Anto | Effects of Tip Clearance Gap and Exit Mach Number on Turbine Blade Tip and Near-Tip Heat Transfer[END_REF]) show the highly threedimensional nature of the flow in the cavities and the complex interaction between the leakage flow and the main passage flow. The cavities are fed by mainstream fluid and the flow inside the cavities is influenced by the blade-to-blade pressure field from the upstream and the downstream blade row and secondary flow structures from the upstream blade row (passage vortex, wakes, etc.). The presence of the shroud cavities, the mixing process when the leakage flow re-enters the mainstream, and the influence of leakage flow on the downstream blade row flow field all cause additional losses.

Shroud leakage loss mechanisms have been investigated by many authors. Traupel [START_REF] Traupel | Erbenisse von Turbinenversuchen[END_REF], Denton and Johnson [START_REF] Denton | The Tip Leakage Loss of Turbine Blades -A Review and a Theory for Shrouded Blades[END_REF], Denton [START_REF] Denton | Loss Mechanisms in Turbomachines[END_REF], Gier et al. [START_REF] Gier | Interaction of Shroud Leakage Flow and Main Flow in a Three-Stage LP Turbine[END_REF], Rosic and Denton [START_REF] Rosic | The Importance of Shroud Leakage Modelling in Multistage Turbine Flow Calculations[END_REF], Rosic et al. [START_REF] Rosic | The Control of Shroud Leakage Loss by Reducing Circumferential Mixing[END_REF], all tried to quantify these losses. However, the loss-generating mechanisms are highly dependent on the particular shroud and stage geometry, therefore it is difficult to develop universal design rules and to correlate different geometrical and operational turbine parameters to generic shroud design. A map of the possible turbine efficiency changes caused by different shroud modifications was obtained by Rosic et al. [START_REF] Rosic | The Influence of Shroud and Cavity Geometry on Turbine Performance -An Experimental and Computational Study, Part 1: Shroud Geometry[END_REF][START_REF] Rosic | The Influence of Shroud and Cavity Geometry on Turbine Performance -An Experimental and Computational Study, Part 2: Exit Cavity Geometry[END_REF] who investigated the influence of each geometric parameter on the mainstream aerodynamics. The intention was to summarize these effects and to highlight efficiency trends that can be used by turbine designers as guidelines for their particular shroud geometry rather than to create universal design rules.

Finally, blade tip clearance changes during operation due to rotor unbalances, dynamic deformations and blade vibrations may have a huge impact on the safety, the lifetime and the energy efficiency of turbo machines. Thus, the increasing need for reliability in gas turbine design, both for power generation and aero-propulsion (Fadlun et al. [START_REF] Fadlun | Measurement error influence on gas turbine operability for conditionbased maintenance and reliability/availability improvement[END_REF]), requires tools that are able to include the strong variability due to manufacturing tolerances, assembly process and in service operations in the design system [START_REF] Williams | The Use of Probabilistic Methods in Determining Turbine Disc Cyclic Life Uncertainty[END_REF]. The question arising in recent years is how to take into account these features with their associated aleatory distribution. Bunker [START_REF] Bunker | The Effect of Manufacturing Tolerances on Gas Turbine Cooling[END_REF] applied a Monte Carlo simulation to study the impact of geometrical variations on HP components. The author showed that there is a wide range of parameters that can vary and their relative impact on the component life is different. Manufacturing variations can reduce the nozzle life by 33%. A stochastic study on statistical variations introduced by in service degradation of the rotor tip was studied by [START_REF] D'ammaro | Uncertainty Quantification and Film Cooling[END_REF]. The authors proved that the tip heat transfer is strongly influenced by the stochastic distribution of fillet radius and tip gap. The study has shown that is crucial to use an uncertainty quantification approach in hot streak migration predictions. A summary of the current state of UQ for aircraft engines can be found in [START_REF] Montomoli | Uncertainty Quantification in Computational Fluid Dynamics and Aircraft Engines[END_REF].

Most of the uncertainty studies in literature have been obtained in simplified geometries. The study of Bunker [START_REF] Bunker | The Effect of Manufacturing Tolerances on Gas Turbine Cooling[END_REF], for example, was based on a flat plate assumption. More recently some studies have considered the impact or real geometries [START_REF] D'ammaro | Uncertainty Quantification and Film Cooling[END_REF][START_REF] Carnevale | Uncertainty Quantification: A Stochastic Method for Heat Transfer Prediction Using LES[END_REF][START_REF] Montomoli | Uncertainty Quantification and Conjugate Heat Transfer: a Stochastic Analysis[END_REF][START_REF] Montomoli | Geometrical Uncertainty in Turbomachinery[END_REF][START_REF] Büche | Uncertainty Analysis for large-scale industrial radial compressors[END_REF][START_REF] De Maesschalck | Performance Robustness of Turbine Squealer Tip Designs due to Manufacturing and Engine Operation[END_REF][START_REF] Panizza | Uncertainty Quantification of Hot Gas Ingestion for a Gas Turbine[END_REF][START_REF] Wunsch | Quantification of Combined Operational and Geometrical Uncertainties in Turbo-Machinery Design[END_REF]. However it is not clear how uncertainty is propagating in a flow. Carnevale et al. [START_REF] Carnevale | Uncertainty Quantification: A Stochastic Method for Heat Transfer Prediction Using LES[END_REF] have shown that the impact of uncertainty is modulated in a channel assuming a "transitional" like behavior. For this reason in this work we decided to simulate a multistage HP steam turbine to evaluate how the impact of uncertainty on seals clearance is affecting the downstream stages and is redistributing in the span-wise direction.

THE MODEL TURBINE AND EXPERIMENTAL METHODS

The experiments were carried out by Rosic [START_REF] Rosic | The Control of Shroud Leakage Flows in Low Aspect Ratio Multistage Turbines[END_REF] using a low speed multistage air turbine that was designed to represent the first few stages of a high pressure steam turbine (Fig. 1). Atmospheric air enters the inlet section radially and passes through a flow straightener before entering the turbine. Air is drawn through three turbine stages by a fan downstream of the turbine. Both stator and rotor blades were shrouded. The operating point, determined by the flow coefficient and rotational speed, is set by the fan power and dynamometer brake. The flow coefficient v x /U is calculated using the axial velocity at turbine inlet, which is measured indirectly using a mass flow measurement in the exhaust section and a density measurement at inlet. ). The key rig geometrical and operational parameters are presented in Tab.1. The blade parameters are taken at blade mid-height. The blading was designed to represent typical high pressure steam turbine conditions, and all parameters satisfy that condition. Only the Reynolds number, based on exit velocity and true chord, is significantly lower than that in a real turbine.

The rotor shroud cavities geometry and sealing arrangement under scrutiny is described in Figure 2. This configuration is representative of commercial turbine design as it allows for the axial movement of the rotor shaft relative to the casing. The shroud forms a radial clearance of 3mm with the casing. Using two approximately 0.75 mm clearance radial seals resulted in an over-shroud leakage flow of approximately 2%. A representation of the three stages experimental turbine is reported in Fig. 3.

To investigate the time-mean flow of the main passage, traverses were conducted using a five hole pneumatic probe. The properties were measured by performing full-span area (1) 2.53•10 5 2.97•10 5 (1) based on true chord and exit velocity Fig. 2. Shroud sealing arrangement traverses over one blade pitch at an axial distance of 6 mm and 17 mm downstream of the stator and rotor trailing edge tip respectively (Fig. 1). The traverse grid was 37 x 37 uniformly spaced points in both radial and pitch-wise directions. A detailed description of the experimental setup is reported in [START_REF] Rosic | The Control of Shroud Leakage Flows in Low Aspect Ratio Multistage Turbines[END_REF].

TBLOCK DESCRIPTION

TBLOCK is a multi-block structured grid solver developed by John Denton [START_REF] Denton | TBLOCK Manual[END_REF]. This is a steady or unsteady RANS solver and uses the finite volume method, with the explicit SCREE scheme [START_REF] Denton | The Calculation of Three Dimensional Viscous Flow Through Multistage Turbomachines[END_REF] for steady simulations and dual-time stepping method for unsteady calculations. Turbulence is modelled using a mixing length approach. Laminar to turbulent boundary layer transition can be modelled by assuming a laminar boundary layer up to a specified point on each blade and end-wall surface and a turbulent one downstream. The mixing plane approach allows reversed flow across the mixing plane, which may occur near the end-walls close to the leakage cavities. To accelerate convergence multigrid and spatially varied time steps are used. The code is aimed at predicting both main blade path and secondary gas path flows in turbomachinery although it can be used to calculate other types of flow. The solver overcomes most of the limitations of the simple shroud leakage model in calculating complex geometries, and was used for the calculations in this paper including the shroud leakage flow path. The flow field is divided into as many 'cuboid' blocks as are necessary to represent the geometry and each block is solved separately. Information passes between blocks using different types of interfaces treated as 'patches' (inlet and exit boundary, solid boundary, mixing plane, periodic boundary). The method is tolerant of extreme levels of grid distortion, and when necessary can interpolate the flow data from one patch onto another, which may lie on a non-contiguous face and may have a different number of grid points.

Four configurations with different radial seals clearance (0.25, 0.50, 0.75, 1.00 mm) on the shroud cavities for all the three stages were modelled using TLOCK. Simulations were completed modelling all three stages single passage domain, including rotor tip shroud and stator hub leakage paths. Figure 4 presents the flow domain decomposition.

The domain was divided into 72 blocks. In particular, the leakage flow paths were divided into the necessary number of rectangular blocks using structured H meshes. Two different levels of mesh refinement were analysed. For the finer mesh, the main blade passage was modelled using 124 grid points in the axial, and 114 grid points in both pitch-wise and radial directions (Fig. 5). The whole domain consists of 36.89 million grid points. The coarser mesh adopted 99 grid points in the axial direction and 91 points in pitch-wise and radial directions, for a total of 18.92 million computational nodes. Fig. 5. TBLOCK computational grid Fig. 6. TBLOCK -rotor and tip shroud grid structure

The simulation was run in parallel using 72 processors. A fully converged unsteady solution initialised from previous steady computations was obtained for 10 blade-to-vane passes. Averaging of the flow variables was then conducted for a further 10 blade-to-vane passes. The computational time required to complete one blade-to-vane pass for the finer mesh was about 7 hours.

The detailed grid and block structure used to model each rotor blade passage is presented in figure 6. The rotor tip shroud for all the configurations under scrutiny was modelled with different radial clearance gaps reproduced for every stage.

NUMERICAL RESULTS

Comparison with experiments

In comparing calculated and measured flow patterns it must be remembered that the real clearance is not known with great accuracy and that it varies significantly around the circumference. In the experimental configuration, the shroud forms a radial clearance of 3mm with the casing and the average radial gap for seals is assumed to be 0.75 mm.

Predicted distributions of total pressure coefficient downstream of the second and third stators and their comparison with experimental results are shown in Fig. 6. Both pitch-wise averaged profiles (left) and 2-D contour map (right) are presented. The measurement plane covers one pitch from 3% to 95% span. As it can be observed on the pitch-wise averaged profiles, numerical simulations on both coarse and fine mesh well reproduce the pressure levels and the main trend along the span. A comparison between 2-D contours show that TBLOCK calculations captured reasonably well the main flow features, such as the locations and intensity of the main loss core and of the hub end-wall secondary flow, and the migration towards the mid-span of the concentrated low momentum flow associated with casing end-wall boundary layers and leakage. Comparisons between experimental and numerical pitchwise averaged relative total pressure coefficients downstream the second and third rotor, are reported in Fig. 8. Numerical simulations are able to predict the radial locations of the mid-span and end-walls loss cores, although their levels were slightly overestimated.

Velocity field is analysed in Fig. 9 and 10, which show the span-wise distribution of pitch-wise averaged values of axial velocity (Fig. 9) and yaw angle (Fig. 10) at different downstream stator/rotor locations. Numerical simulations were in good agreement with experimental results in predicting magnitude and directions of the flow downstream both vanes and blades passages. In particular, in the region close to the casing and the shroud exit cavity of the third rotor, TBLOCK was able to reproduce the overturning trend in the main passage flow angle caused by mixing of the leakage flow.

Different seals clearance gaps

According to steam turbine industrial practice, due to manufacturing tolerances and assembly process, real shroud geometry presents some mismatches with their theoretical design and it is also very hard to have a reliable estimation on rotor shroud gap. Furthermore, during long-term turbine operation the sealing effectiveness deteriorates due to rubbing and worn sealing elements. Because of this, leakage fractions for shrouded blades can vary from 1% to 4 % within HP turbines [START_REF] Rosic | The Influence of Shroud and Cavity Geometry on Turbine Performance -An Experimental and Computational Study, Part 1: Shroud Geometry[END_REF]. Therefore it is important to investigate the influence of the different shroud radial gap on the main blade path. This was done in a numerical study using the same numerical flow domain as in the first part of this section varying the seal clearance over both radial fins simultaneously while the datum shroud geometry was kept unchanged.

The variations of pitch-wise averaged total pressure coefficient and yaw angle distributions with clearance gap are reported in Figures 11 and12, respectively. Absolute and relative values are shown at stators exits (the left sides) and downstream the rotors (right side), respectively.

It is worthwhile noting how the flow field in the exit rotor shroud cavity is greatly affected by the leakage jet generated in the upstream seal gap and re-entering the main passage. The leakage flow fraction is enhanced with increasing clearances, causing higher pressure losses near the casing (Fig. 9, right side). Also, as the clearance increases, the leakage jet increases its momentum (Fig. 11, right side), forming a region of higher relative yaw angle in the near casing of the main annulus (Fig. 12, right side), because of mechanism explained in [START_REF] Rosic | The Influence of Shroud and Cavity Geometry on Turbine Performance -An Experimental and Computational Study, Part 1: Shroud Geometry[END_REF].

The negative influence of the leakage flow on the flow field in the downstream blade row can be observed on the left side of Fig. 12, by analysing the yaw angle distributions downstream the stators exits. The leakage flow caused a strong underturning of the flow at the mid-span, and also the 'S' shape angle distribution close to the casing, associated with the casing secondary flows, was enhanced.

UNCERTAINTY QUANTIFICATION

A basic method to obtain the stochastic distributions of output variables like the stagnation pressure or the yaw angle if the clearance gap varies randomly is the Monte Carlo method. First, a probability distribution must be assumed for the clearance gap. Then, a suitable random number generator is used to create the individual input random samples to follow the prescribed distribution. Finally, the simulation only has to be run repeatedly for every random sample drawn from the PDF. Convergence is achieved by increasing the number of the simulations. Unfortunately, the convergence is slow: the sample size has to be quadrupled every time to half the error. For a cost intensive computational model, as it is the case here, it is therefore prohibitive to use a Monte Carlo method.

In the given simulation, only one random variable is present, the clearance gap, and it can be assumed to be Gaussian. This assumption is made based on the available literature on manufacturing uncertainties, which tend to display Gaussian behaviour [START_REF] Panizza | Uncertainty Quantification of Hot Gas Ingestion for a Gas Turbine[END_REF]. Both assumptions together can be used to achieve an immense cost reduction by using a Non-Intrusive Polynomial Chaos Method: the Probabilistic Collocation Method (PCM) as described in [START_REF] Hosder | Non-intrusive polynomial chaos methods for uncertainty quantification in fluid dynamics[END_REF].

The underlying idea of all Polynomial Chaos methods is that the functional form of the stochastic outputs of a model can be found with significantly lower computational effort, if they are approximated through a basis of optimal orthogonal polynomials defined through the input random variables. The additional concept of probabilistic collocation is that since the statistical moments are integrals of the output PDFs, they can be cheaply and accurately integrated by using optimal Gaussian quadrature on the found polynomial expansion. In consequence, the model only has to be run a few times using the Gaussian collocation points of the clearance gap to calculate accurate numerical approximations for all output distributions.

The number of the Gaussian collocation points is naturally connected to the polynomial order of the expansion. They are chosen as the roots of the next higher order polynomials. For example, 5 points are needed for a 4 th order expansion. Moreover, the number is also dependent on the number of input variables used. Since for multiple input variables the collocation points have to be combined through tensor products, Polynomial Chaos methods become quickly inefficient for increasing numbers of input random variables. For only one input variable, however, five collocation points can achieve an accuracy level comparable a million Monte Carlo runs [START_REF] Hosder | Non-intrusive polynomial chaos methods for uncertainty quantification in fluid dynamics[END_REF].

Probabilistic Collocation Method (PCM)

The basic concept of Polynomial Chaos expansions can be traced back to Wiener [START_REF] Wiener | The Homogeneous Chaos[END_REF], who found that the stochastic output y x ,ξ ( ) of a function can be decomposed into a series with deterministic coefficients α i x ( ) and stochastic orthogonal polynomials ψ i ξ ( ) , such that:

y x ,ξ ( ) ≈ α i x ( ) ψ i ξ ( ) i=0 N P -1 ∑ Eq.1
where N P is a truncation to limit the infinite expansion series.

For the one dimensional case N P can simply be chosen as the maximum polynomial order of the chosen expansion plus one. Moreover, ξ is a Gaussian random variable representing the clearance gap, and y is representative for either the stagnation pressure ot the yaw angle. For a Gaussian random variable, the optimal orthogonal polynomials ψ i ξ ( ) are the Hermite polynomials [33]. The choice of the polynomials is always a consequence of the probability distribution of the inputs variables. Thus, only the coefficients α i x ( ) are unknown. In order to find these, a linear system of N p equations can be solved in the optimal Gaussian collocation points ξ i : Eq.2

Here, the clearance gap is chosen as Gaussian random with distribution N 0.75,0.03356

(

) . This results in approximate gap widths collocation points 0.25,0.5,0.74,1.00 and 1.25 as used before and the PDF as shown in Fig. 14.

Fig. 14. Probability distribution used as input PDF including optimal Gaussian collocation points.

The mean µ y and variance σ y 2 of y can then be found using the coefficients by calculating

µ y = α 0 σ y 2 = α i 2 ψ i 2 i=1 N P -1 ∑ Eq.3
The PDF of y can be illustrated by Monte Carlo sampling the found expansion in eq. 1 using random samples with PDF N 0.75,0.03356 ( ) . Unlike direct sampling of the CFD model, this sampling process needs negligible computation time.

Convergence of the polynomial expansion can be found by increasing the order of the expansion. Here, orders from 3 to 5 were tested. Comparing 4 th and 5 th order a four decimal accuracy was found so that no higher orders were computed.

Results

Figure 15 shows the pressure coefficient distribution at the exit of each row. The dashed red line is the standard deviation. The uncertainty bars represent one standard deviation superimposed on the mean pressure distribution. In the first stage we can observe that the impact of uncertainty in the rotor tip is affecting only the higher part of the rotor. Only about 10% of the span (100-90%) is influenced by the rotor tip random variation. Similarly large random variation in this region can be observed for rotor 2 and rotor 3 as well. For the second stator, the peak of uncertainty is surprisingly at 50% and the region between 40% and 60% of the span is mostly affected. For the second rotor, the entire span is affected. In the last stage, the standard deviation for the stator peaks at 78%, with a region of increased variation higher between 65 and 80%. The last rotor is most strongly affected by the uncertainty with a strong variation from 20 to 80%.

There are different mechanisms involved. The changes in rotor tip gaps modifies the mass flow distribution and this alter the radial uncertainty. At the same time, the impact of the leakage flows modifies the work carried out by the rotor due to the variation of the inlet (and exit) angle. This is clearly shown in figure 16 where the exit angles indicates the variation in load associated to each row.

It is important to notice that the found regions with high standard deviation agree with the regions in which the CFD results differed largely from the experimental in the original work of Rosic et al [START_REF] Rosic | The Importance of Shroud Leakage Modelling in Multistage Turbine Flow Calculations[END_REF]. The performed uncertainty quantification study therefore gives suggest a reason of this disagreement between CFD and experiments: it can be associated to a small variation of the tip gap in the machine that is different from the CFD model.

Fig. 1 .

 1 Fig. 1. Schematic of model turbine Experiments were performed on 50% reaction blading at the design operation condition ( 384 . 0 = φ). The key rig geometrical and operational parameters are presented in Tab.1. The blade parameters are taken at blade mid-height. The blading was designed to represent typical high pressure steam turbine conditions, and all parameters satisfy that condition. Only the Reynolds number, based on exit velocity and true chord, is significantly lower than that in a real turbine.The rotor shroud cavities geometry and sealing arrangement under scrutiny is described in Figure2. This configuration is representative of commercial turbine design as it allows for the axial movement of the rotor shaft relative to the casing. The shroud forms a radial clearance of 3mm with the casing. Using two approximately 0.75 mm clearance radial seals resulted in an over-shroud leakage flow of approximately 2%. A representation of the three stages experimental turbine is reported in Fig.3.To investigate the time-mean flow of the main passage, traverses were conducted using a five hole pneumatic probe. The properties were measured by performing full-span area

Fig. 3 .

 3 Fig. 3. Representation of the experimental turbine

Fig. 7 . 3 Fig. 8 . 3 Fig. 9 . 3 Fig. 10 .

 73839310 Fig. 7. Measured and predicted total pressure coefficient downstream stator 2 & 3

Fig. 11 .

 11 Fig. 11. Predicted pitch-wise averaged total pressure coefficient profiles for different seals clearances

Fig. 12 .

 12 Fig. 12. Predicted pitch-wise averaged yaw angle profiles for different seals clearances

Fig. 13 .

 13 Fig. 13. Change in turbine efficiency with leakage fraction Finally, the influence of leakage flow on the mainstream flow in the three stage turbine analysed in this paper can be summarized in the diagram of Fig. 13, which shows the changes in turbine efficiency with leakage fraction. He resulting trend is in agreement with the classical shroud leakage theory which assumes that the loss in efficiency associated with leakage flow is directly proportional to the leakage flow rate

Fig. 15 .

 15 Fig. 15. Stagnation pressure coefficient distribution and uncertainty

Fig. 16 .

 16 Fig. 16. Yaw angle variation and standard deviation superimposed

Table 1 .

 1 Turbine Geometrical and Operational Parameters

	Design operating point data
	Flow Coefficient φ	0.384
	Stage Loading Ψ=∆h0/U 2	1.0
	Design Speed [rpm]	830
	Inlet Mach Number	0.045
	Inlet/Outlet Total Pressure [Pa]	100465.4	75545.5
	Inlet/Outlet Total Temperature [K]	300	279
	Number of Stages	3
	Inter Blade Row Spacing [mm]	25
	Blade Height h [mm]	75
	Mid-span Radius [mm]	462.5
	Hub to Tip Ratio	0.85
	Blade parameters	Stator	Rotor
	Blade Number	40	38
	True Chord l [mm]	89.73	105
	Aspect Ratio h/l	0.836	0.714
	Pitch to Chord Ratio	0.81	0.728
	Reynolds Number