Automated CFD Optimisation of a Small Hydro Turbine for Water Distribution Net- works Master’s thesis in Applied Mechanics Alexander Ohm and Hafþór Örn Pétursson Department of Applied Mechanics CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2017 Master’s thesis 2017:38 Automated CFD Optimisation of a Small Hydro Turbine for Water Distribution Networks ALEXANDER OHM AND HAFÞÓR ÖRN PÉTURSSON Department of Applied Mechanics Division of Fluid Mechanics Chalmers University of Technology Gothenburg, Sweden 2017 Automated CFD Optimisation of a Small Hydro Turbine for Water Distribution Networks ALEXANDER OHM AND HAFÞÓR ÖRN PÉTURSSON © ALEXANDER OHM AND HAFÞÓR ÖRN PÉTURSSON, 2017. Supervisor: Håkan Nilsson, Chalmers University, Department of Applied Mechanics Supervisor: Martin Holm, Undisclosed company Examiner: Håkan Nilsson, Chalmers University, Department of Applied Mechanics Master’s Thesis 2017:38 Department of Applied Mechanics Division of Fluid Mechanics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Pilot design and final design after optimisation, in a mean kinematic pressure field spanning 0 ≤ P/ρ ≤ 2.25. Typeset in LATEX Gothenburg, Sweden 2017 iv Automated CFD Optimisation of a Small Hydro Turbine for Water Distribution Networks ALEXANDER OHM and HAFÞÓR ÖRN PÉTURSSON Applied Mechanics Chalmers University of Technology Abstract Leakages in water distribution networks is a global problem, with leakage rates up to 50% of the flow [1]. One way to find leakages is to adopt more sensors in the network. Such sensors need electric power that is commonly not distributed with the network. An option is to generate the electricity where it is needed, using small hydro turbines that extract energy from the flowing water and produces just enough electric power for the sensor at each site. The present work provides an automated optimisation procedure for the design of such turbines, and employs it to optimise a turbine with respect to efficiency and kinematic pressure drop at a desired operating point. The automated optimisation loop handles the entire procedure of geometry creation, meshing, simulation, post-processing and optimisation using open source software. The geometry generation and solution procedure is based on a the work by Bergqvist [2], who implemented a Ruby code to generate a parameterised turbine geometry, generated the mesh with cfMesh, and used OpenFOAM to simulate the flow. The present work extends the method by using Dakota to create an automated loop for optimisation, and applies it for water network turbines. The optimisation is done by varying the pitch, number of rotor and stator blades, as well as the hub radius. A steady-state frozen-rotor approach with multiple ref- erence frames (MRF) is used in the simulations. The turbulence is modelled with the k− ω SST turbulence model. After convergence studies of the pilot design, the simulations are carried out on unstructured grids with 2 − 3 million cells for the varying geometries. To reduce the computational burden that direct optimisation with CFD simulations implies, a surrogate modelling technique is implemented using the Efficient Global Optimisation (EGO) algorithm [3] [4]. For this, Kringing modelling is used to create a response surface, and division of rectangles (DIRECT) to optimise it and find in- fill points for further evaluation. 51 sample points, generated with latin hypercube sampling, are simulated before the first abstraction of the objective function is cre- ated. After the learning period, 61 additional design iterations are simulated before the algorithm reaches the stopping criteria. In between each of the 61 optimisation designs, a surrogate optimisation is carried out to select the next infill point. The completed optimisation provides a new design that performs at 4.4 percent- age points higher efficiency and with 72% lower kinematic pressure drop. The new design features a bigger flow through area and a lower work output that is more v adapted to the application with a reduction of 70%. Keywords: Optimization, Turbomachinery, Hydro Turbine, OpenFOAM, Dakota, cfMesh, Ruby. vi Acknowledgements This project has been a valuable lesson and at certain times it proved to be tough. First and foremost, we would like to thank our supervisors Håkan and Martin for their input that made this project possible. We would also like to thank our better halves for their equally important support through all our ups and downs. Thank you also to Niklas and the senior staff at the company, for pushing us towards our goal and for showing interest in our progress. The very helpful PhD students at the Department of Applied Mechanics also deserve credit for their important tips and feedback and for putting up with our questions. Last but not least, we would like to thank Marcus, Karl and Fredrik for brightening up our days at the office. Thank you! Alexander Ohm, Hafþór Örn Pétursson Gothenburg, June 2017 viii x Nomenclature Abbreviation CFD Computational Fluid Dynamics EGO Efficient Global Optimisation EIF Expected Improvement Function MRF Multiple Reference Frames RANS Reynolds Averaged Navier Stokes URANS Unsteady Reynolds Averaged Navier Stokes Greek Symbols α Absolute flow angle β Relative flow angle εijk Permutation tensor η Efficiency γ Stagger angle ν Kinematic viscosity νt Turbulent kinematic viscosity Ω Angular velocity ω Specific dissipation φ Angle of attack ρ Density ε Dissipation Latin Symbols u Mean velocity hb Hub axial length from leading edge HE Hydraulic head Lb Blade chord length Nblades Number of blades p Instantaneous static pressure PH Hydraulic power Pt Power by turbine u Instantaneous velocity u′ Fluctuating velocity u′RMS root-mean-square of the turbulent velocity fluctuation U∞ Free stream velocity uAi Absolute velocity in stationary reference frame uRi Relative velocity in rotational reference frame w Weights wb Tangential chord length xi Z Zweifel number b Axial chord length c Absolute velocity g Gravitational acceleration hr Hub radius I Turbulence intensity k Turbulent kinetic energy l Turbulent length scale LB Lam-Bremhorst LS Launder-Sharma M Moment Q Volumetric flow rate r Radius rb Number of rotor blades Re Reynolds number rp Rotor pitch s Blades spacing sb Number of pre-swirl stator blades sp Pre-swirl stator pitch U Blade speed w Relative velocity Other Symbols i, j & k Coordinate directions in space x & θ Representing axial and tangential directions respectively xii Contents 1 Introduction 1 1.1 Aim of the project . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Method Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Theory 7 2.1 Turbine Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Blade Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1 Zweifel Number . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Turbulence Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5 Turbulence Quantities . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.6 Wall Treatment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.7 Optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.7.1 Efficient Global Optimisation Algorithm . . . . . . . . . . . . 16 2.8 Cavitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3 Methods 19 3.1 Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1.1 Geometrical Simplifications . . . . . . . . . . . . . . . . . . . 22 3.1.2 Geometrical Variables for Optimisation . . . . . . . . . . . . . 23 3.1.2.1 Hub . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.2.2 Pre-Swirl . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.2.3 Rotor blades . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Mesh Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 CFD setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.1 Flow Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3.2 Rotational Zone . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3.3 Boundaries and Wallfunctions . . . . . . . . . . . . . . . . . . 29 3.3.4 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 30 3.3.5 Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.3.6 Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4 Post Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.5 Optimisation Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 36 4 Method Validation 39 xiii Contents 4.1 Mesh Quality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2 Mesh Grid Independence . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2.1 Internal Mesh Comparison . . . . . . . . . . . . . . . . . . . . 41 4.2.2 Boundary Layer Resolution . . . . . . . . . . . . . . . . . . . 42 4.2.3 Domain Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.3 Rotor Position for the Pilot Design . . . . . . . . . . . . . . . . . . . 47 4.4 Turbulence Model Comparison . . . . . . . . . . . . . . . . . . . . . . 48 4.5 Flow Fields of the Pilot Design . . . . . . . . . . . . . . . . . . . . . 49 4.5.1 Casing Simplification Validation . . . . . . . . . . . . . . . . . 53 5 Results 55 5.1 Optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.1.1 All Designs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.1.2 Best Designs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.2 Evaluation of Final Design . . . . . . . . . . . . . . . . . . . . . . . . 61 5.2.1 Geometry Comparison . . . . . . . . . . . . . . . . . . . . . . 61 5.2.2 Rotor Position for the Final Design . . . . . . . . . . . . . . . 62 5.2.3 Flow Fields of the Final Design . . . . . . . . . . . . . . . . . 63 5.3 Propositions for and Discussion of Future Designs . . . . . . . . . . . 66 6 Conclusion 69 6.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 A Folder Structure I B Preliminary Simulation Test Using k − ε III C Additional Mesh Resolutions V D Flow Fields for Varying Rotor Position VII E Off Design IX xiv 1 Introduction Fresh water is part of people’s everyday life and in many parts of the world its availability in the tap is taken for granted. As global development continues, water demand increases. At the same time increasing water depletion is predicted [5]. This means that water should be used wisely; leakages in the water distribution sys- tems is however, a global problem with leakage rates between 3%-50% of the flow [1]. There are several methods to detect leakages in the distribution pipe networks. Eyuboglu et al. [6] divided the procedures of detecting leakages into three cate- gories; biological, hardware and software based methods. Biological methods detect leakages via on-checks, using human and animal senses, such as sight, hearing or sense of smell. Hardware methods uses tools, such as sensors, to measure any changes occurring in or around the leaking pipe, such as unwanted presence of water. Lastly, software methods analyse flow differences in the pipe system, such as mass equilib- rium, which may indicate leakages. One way to improve leakage detection is to introduce more sensors into the systems. These devices need energy and common solutions today are to provide energy with batteries or cables connected to the grid. Replacing batteries in remote locations or digging for cables in city environments is expensive and thereby limits the affordable number of sensors in the network. To make energy more available, a hydro turbine that provides the necessary energy directly from the water flow at the site is proposed by an undisclosed company. If many installations of such turbines are introduced into a connected network the downstream pressure influences from each unit adds up which could affect the function of the system. Contrary to most turbine applica- tions, it is therefore important to limit the static pressure drop over the turbine for this kind of application. Thus, many established design principles do not apply and the development of the turbine is therefore an exploratory process, characterised by trial and error. Such processes often benefit from replacing physical tests with simulations. Furthermore, optimisation can be added to the development process to aid with finding the right turbine geometry suitable for fresh water pipelines. Computational fluid dynamic (CFD) simulations has been the subject of much re- search to better understand turbine flows [7][8][9]. There are several different sim- ulation methods to simulate rotor stator interactions with CFD, each with its own drawbacks and advantages. Three popular approaches are evaluated here; the frozen rotor-, mixing plane - and sliding grid approach. 1 1. Introduction Frozen rotor. The frozen rotor approach is a steady state method where the rotor does not rotate physically. It is kept fixed and the rotation is applied to the fluid through a source term in the momentum equation. This is done by splitting the domain into multiple reference frames or MRF regions. A rotating reference frame for the rotating part of the fluid and a stationary reference frame is used for the rest of the domain. Interactions between rotating and stationary regions which cause unstable flows do not take place as no real rotation is present. Thus, as it is a steady state method, the results only show the flow at an instant of the rotation and do not accurately represent the real unsteady flow. However, this method is not computational expensive and gives reasonable estimates of the flow behaviour. These traits have made the frozen rotor approach popular in early design phases of hydro turbines [8]. Mixing Plane. The mixing plane method is another steady state method where the flow is divided into different sections. Averaging is utilised on the interfaces of these sections to update the boundary conditions between the stationary and rotat- ing sections. This method has the downside of not being able to predict the wakes generated by the blades due to the averaging. It is however more physical than the frozen rotor approach and can yield better results depending on the parameters of interest [8]. Sliding Grid. The sliding grid method is an unsteady method where the rotor moves in the domain. The real rotation of the flow is simulated by moving the mesh grid at the rotating region. The interaction which is not accounted for in the previous methods, is thus better predicted. As such it is more physical than the previously mentioned methods, but also computationally expensive [8]. Kaufmann and Bertram [7] compared the frozen rotor with a sliding grid approach and found good correlation for propeller efficiency when simulating surging forward motion of ships. Dubas [9] optimised a rim driven thruster and states that steady state based approaches show a good trade-off between computational cost and result accuracy. Because of its ability to handle varying geometries and its cheap compu- tational cost, the frozen rotor seems to be a viable option. In recent studies the so called evolutionary algorithms such as Genetic algorithms, non-dominant sorting genetic algorithm II and multi-objective genetic algorithm, have been used in combination with CFD [10][11][12][13][9]. All these methods explore the entire design space, to find a global optimum to the problems. Mohamed et al. [11] optimised an airfoil shape of a Wells turbine, using an evolutionary algorithm. His method managed to increase the work output by the blades by 11.3% as well as increase the efficiency by 1%. Using eleven variables describing the airfoil shape, it took 615 evaluations to reach an optimum. Öksüz and Akmandor [13] did an aerodynamic optimisation study on gas turbine blades. They used a modified multi-objective genetic algorithm combined with surrogate modelling to speed up the convergence towards an optimum. They concluded that the modified algorithm in association with surrogate model is a promising tool. Dubas [9] did an 2 1. Introduction optimisation study on a rim driven thruster, where he, like Öksüz and Akmandor, used a genetic algorithm in combination with a surrogate model to optimise the design. Generally, evolutionary algorithms are fast in finding prominent areas, their convergence to an optimum can however be slow [14]. The combination of a global optimisation tool with a surrogate model is therefore a favourable option when computational resources are scarce. 1.1 Aim of the project The present work is carried out at an undisclosed company where a pilot design of a turbine for the drinking water network had been developed, see pilot design in Fig- ure 1.1, where main components are shown. The goal is to develop an automated optimisation method that requires minimal manual input, and to apply it to im- prove the turbine designed by the undisclosed company. This is done by modifying the rotor blades, stator blades and the hub of the turbine, marked in Figure 1.1. Preferably the method should be cheap, both in terms of cost and computational resources. An optimisation approach based on open source software that can run and complete on a single workstation computer within a practical amount of time is therefore sought. The method should also be flexible such that it can be used in the future to carry out similar optimisation of new problems. Figure 1.1: The pilot design, with markings of important components. Artist: Martin Holm [15] 3 1. Introduction 1.2 Method Outline The optimisation method is divided into four steps, illustrated in in Figure 1.2. First, geometry is created, then a mesh grid is made using the geometry from pre- vious step. CFD simulations are carried out on the grid and lastly, the results are evaluated in an optimisation algorithm and a new design is proposed. At the Gothenburg Region OpenFOAM User Group Meeting in 2015, Bergqvist [2] shared a tutorial case that creates and simulates a turbine from a script. He created geometries with the code language Ruby and used the CfMesh and OpenFOAM software to mesh and simulate the flow around them. He suggested that an opti- misation software would be added to his loop to close the chain and form a fully automated loop. In the present work, the mathematical functions from Bergqvist’s [2] Ruby scripts are applied to create the geometry of the pilot turbine design from the company in a fully parameterised manner, with coded variables that can be used to modify the design during optimisation. The software CfMesh is then used to create a mesh to discretise the fluid volume. The OpenFOAM software is used to simulate the fluid flow with the frozen rotor approach because of its ability to handle varying geometries and its cheap computational cost. The software Dakota is used for the optimisation. With the goal of developing a computationally cheap method in mind, a surrogate optimisation technique is sought, similarly to the one used by [9]. The available derivative-free global method with surrogate modelling, the Efficient Global Optimisation (EGO) algorithm, is selected. It differs from the referenced papers in that it uses a sampling technique instead of an EA to optimise the response function of the surrogate model. Optimisation Geometry creation Spatial discretisation CFD Simulations Figure 1.2: Overview of the optimisation process 4 1. Introduction 1.3 Thesis Outline This thesis is divided into 6 chapters. After this introduction, chapter two describes the theoretical background relevant for the present work. In chapter three, the workings of the developed loop are described. Chapter four presents the test that are carried out to verify the loop. Chapter five then show the results that are generated from applying the method to improve the pilot turbine. Lastly, chapter six ends with conclusions and suggestions for the future. 5 1. Introduction 6 2 Theory 2.1 Turbine Performance Theory surrounding turbomachinery is a wide topic with many different kind of applications. The turbine considered in this present work is a hydro turbine which converts flow of water into usable electric energy. Turbines are often evaluated based on several parameters, the ones of concern in the present work are shown and dis- cussed in what follows. An important quantity for most turbomachines is hydraulic efficiency. It quantifies the percentage of the available energy stored in the flow that can be subtracted and later turned into electrical energy. The hydraulic efficiency of a turbine is more precisely defined as the ratio between the power output of the turbine,Pt, and the hydraulic power of the flow, PH , as [16], η = Pt PH . (2.1) The power of the turbine is expressed as Pt = ΩΣM , (2.2) and the hydrolic power of the flow as Ph = Qρg(H1 −H2), (2.3) where Ω is the angular velocity, ΣM is the sum of all moments acting on the turbine around the rotational axis, Q is the volumetric flow rate, ρ is the fluid density, g is the gravitational constant and H1 and H2 is the hydraulic head in the flow before and after the turbine, respectively. The specific hydraulic energy term, gH, is the available energy in the flow. Dixon [17] describes this as the sum of pressure, kinetic and potential energy and Stoessel [16] shows that it can be expressed as gH = p ρ + u2 2 + gz, (2.4) where p is the static pressure u is the velocity magnitude and z an elevation between a reservoir and the inlet of the turbine. 7 2. Theory Equation 2.4 can be rewritten as an expression of the total pressure p0 as [18] p0 = p+ 1 2ρu 2 + gzρ = gρH. (2.5) By combining equations 2.1, 2.3 and 2.5, the hydraulic efficiency can be written in terms of the total pressure drop before and after the turbine as η = ΩΣM Q∆p0 . (2.6) In the present work, the total pressure is calculated at the inlet and outlet of the pipe enclosing the system. Often, in turbomachinery analysis, it is desirable to rank or compare different tur- bines regardless of their sizes and configuration. This can be done by using dimen- sional analysis to create dimensionless quantities that describe the ratios between different relevant properties. One such number that is used in the present work is the specific speed, which can be used to evaluate what type of turbomachine is suitable for a given application. It is defined as [17] Ωs = ΩQ 1 2 (gH 3 4 ) . (2.7) The turbine in this task can be identified as a propeller turbine, which according to Dixon [17] typically operates in the specific speed range of Ωs ≈ 2 to 4. 2.2 Blade Design At an early design phase, so called velocity triangles, shown in Figure 2.1, are com- monly used to design the turbine blades. The analysis matches the flow angles and magnitudes to the geometry angles on the leading and trailing edges of all blades. Velocity triangles are drawn in the circumferential plane and cut through all blades at a constant radius. Analogue triangles can be drawn for every radius. By assuming a constant axial flow and that the flow occurs at a constant radius, the flow angles between the stations can be used to find the desired blade angles. The flow given by this analysis is not entirely representative since, in reality there is a radial flow and the flow does not perfectly follow the blade angles. The work output in relation to the flow can be expressed with the Euler equation as [17] Ẇt = ṁ∆(Ucθ), (2.8) where ṁ is the mass flow U is the blade speed calculated as U = rω, c is the absolute flow velocity and the index θ indicates the tangential direction. One way to increase the work output of a turbine is thus to increase the tangential redirection of the flow, ∆cθ. This can be accomplished by adding a pre-swirl in the flow by turning the first stator row. The stators themselves only introduce losses in the flow and do 8 2. Theory Figure 2.1: Velocity triangles for analysing the flow where c is the absolute velocity, w is the relative velocity, U is the blade speed, α is the absolute flow angle and β is the relative flow angle. Index 1 indicates the station upstream of the pre-swirl stator row, index 2 the station between the pre-swirl stators and the rotor, index 3 the station between the rotor and the third row, index 4 the station downstream of the third row and index x the axial direction. not produce any work, but combined with a proper rotor design the total output can be increased. Amin and Xiao [19] managed to increase the overall efficiency of an open water horizontal axis turbine by 13% by introducing pre-swirl. 2.2.1 Zweifel Number The most suitable number of rotor or stator blades for a certain application is not an obvious choice. The more blades there are, the better the alignment with the flow becomes since it is guided more precisely. However more blades also means more surface area which generates friction losses. The so called Zweifel criteraion is commonly used to evaluate this trade off and by assuming an incompressible flow and a constant axial velocity, can be expresseed as Z = 2s b cos2α2(tanα1 + tanα2), (2.9) where s is the spacing between the blades and b is the axial cord length. Through experiments it has been found that a value of Z ≈ 0.8 is optimal and hence an optimal value for s/b can be found given the flow angles [17]. 2.3 Governing Equations The famous Navier-Stokes equation, here given for an incompressible fluid, govern the momentum of viscous fluids as Dui Dt = −1 ρ ∂p ∂xi + ν ∂2ui ∂xj∂xj + fi, (2.10) where p stands for static pressure, u for velocity and f is the body force. It is derived from the conservation of mass, momentum and energy and together with the continuity equation 9 2. Theory ∂ui ∂xi = 0, (2.11) the equation set can describe the motions of fluid flows. By decomposing the velocity into a mean flow component and an instantaneous fluctuation, u = u + u′, and then applying a time average, the Reynolds averaged Navier-Stokes (RANS) equations are obtained as, derivation shown by Davidson [20], ∂uiuj ∂xj = −1 ρ ∂p ∂xi + ∂ ∂xj [ ν ∂ūi ∂xj − u′iu′j ] + fi. (2.12) In this process a new unknown term, the Reynolds stress tensor u′iu′j, appears. The Reynolds stress tensor, makes the equation set unsolvable without introducing further assumptions or modelling equations for closure. One way of modelling the Reynolds stress tensor is by applying the Boussinesq assumption [20], u′iu ′ j = −νt ( ∂ui ∂xj + ∂uj ∂xi ) + 1 3δiju ′ ku ′ k, (2.13) where νt is the turbulent viscosity and u′ku′k = 2k. In the present work, a turbine is solved using the frozen rotor approach, as intro- duced in Chapter [?]. The formulation of Petit et al. [21] showed the governing RANS equations implemented in the MRF approach. The equations do not have the Reynolds stress term, and thus only apply for steady flows as Continuity equation ∂uAi ∂xi = 0 (2.14) Momentum equation ∂uRiuAj ∂xi + ΩiuAjεijk = −1 ρ ∂p ∂xi + ν ∂2uAi ∂xj∂xj , (2.15) where the εijk is the permutation tensor and uAi is the absolute velocity in the stationary frame and uRi is the relative velocity, formulated as uRi = uAi − Ωirjεijk. (2.16) If the turbulence is included, using the Boussinesq assumption and using same derivation assumptions as used in the development of MRF [22], the momentum equation 2.17 can be formulated as ∂uRiuAj ∂xi + ΩiuAjεijk = −1 ρ ∂p ∂xi + ∂ ∂xj [ (ν + νt) ( ∂uAi ∂xj + ∂uAj ∂xi )] . (2.17) 10 2. Theory 2.4 Turbulence Models It is generally known that most flows in nature are turbulent. Turbulent flows are seemingly chaotic since they are full of disturbances in the form of swirling struc- tures. These structures are called turbulent eddies and exist at a wide variety of length scales. The biggest scales are of the same size as the biggest features of the flow. This is where energy is extracted from the mean flow by the turbulence. It is then transported through what is known as the cascade process to smaller and smaller scales, ultimately down to the smallest scale where the energy is dissipated as heat [23]. Turbulent flows are characterised by high Reynolds numbers, which are dimensionless numbers that describes the ratio between inertial and viscous forces in the flow. An example is that turbulence starts to build up in a pipe for Re ≥ 2000, and it undergoes a transition from laminar to turbulent until Re = 105, above that limit the flow is considered to be fully turbulent [24]. For the application in the present work, the Reynolds number in the pipe is about 3.5× 104. Flow structures smaller than the spatial discretisation, the mesh, cannot be resolved and thereby not simulated. Since the dissipative turbulent scales are extremely small compared to most applications relevant for engineers, properly simulating turbulence is often prohibitively difficult as it poses extreme requirements on the resolution. A different approach to deal with turbulence in simulations is to model it numerically. There are many differnt types of turbulence models [23], each of them are appropri- ate for different flow situations and come at different computational costs. The model most widely used in industry is the is the k − ε model. A reason for its popularity is its ability to model a wide range of flows [24]. It can however be in- accurate when simulating adverse pressure gradient flows [20] or flows which rotate, have curved boundaries or swirl [24]. The k − ε model is a high Reynolds number model and as such it is meant to be used with high Reynolds number wall functions, meaning that the near wall flow is modelled by typical behaviour. Modifications to the model has been made in order to further account for the effects of walls. These models are so called low Reynolds formulations for the k − ε model. The k − ε models are discussed by Davidson [23] and Moradnia [25]. Another popular model is the shear stress transport model, the k − ω SST model. This model implements the k − ε model in regions far from the wall and the k − ω model, discussed by Davidson [23], close to the wall. The k − ω model is known to be more accurate than the k − ε model in the near wall region [26] and is better at solving adverse pressure gradient flows [20], which occurs in a flow around wing profiles. The ε and ω equations are combined into one transport equation by the expression ε = β∗kω. The k−ω SST model shown by Menter et al. [26] is mainly used in the present work. Some of the coefficients are expressed with the same notation as in OpenFOAM, shown by Pirooz [25], but with a correction to the pressure term according to NASA [27]. The transport equations for k and ω are formulated as 11 2. Theory ∂k ∂t + ∂ūik ∂xi = min(Pk, 10β∗kω)− β∗kω + ∂ ∂xj [ (ν + αωνt) ∂k ∂xj ] (2.18) ∂ω ∂t +∂ūiω ∂xi = γ νt min(Pk, 10β∗kω)−βω2+ ∂ ∂xi [ (ν+αωνt) ∂ω ∂xi ] +2(1−F1)αω2 1 ω ∂k ∂xi ∂ω ∂xi , (2.19) where the production term Pk , which appears in both the k (2.18) and ω (2.19) equations, is calculated as Pk = νt ( ∂ūi ∂xj + ∂ūj ∂xi ) ∂ūi ∂xj , (2.20) where the turbulent viscosity νt is given by νt = a1k max(a1ω, SF2) . (2.21) F1 in Equation 2.19 and F2 in Equation 2.21 are blending functions written as F1 = tanh  { min [ max ( √ k β∗ωy , 500ν y2ω ) , 4ραω2k CDkωy2 ]}4  (2.22) F2 = tanh [max( 2 √ k β∗ωy , 500ν y2ω )]2 . (2.23) CDkw in Equation 2.22 is defined as CDkw = max ( 2ραw2 1 ω ∂k ∂xi ∂ω ∂xi , 10−10 ) , (2.24) note the similarities with the fourth term in the ω equation (2.19). The equation constants are given by αω1 = 0.5 αω2 = 0.856 αk1 = 0.85 αk2 = 1 β∗ = 0.09 β1 = 0.075 β2 = 0.0828 γ1 = 0.555̄ α2 = 0.44 Blending takes place between k − ω and k − ε through the quantity F1, which according to Menter et al. [26] is zero away from the wall where k− ε model is used and one in the boundary layer where k − ω is used. The constants are blended as ψ = ψ1F1 + ψ2(1− F1), (2.25) where ψ is an arbitrary constant. 12 2. Theory 2.5 Turbulence Quantities The turbulence intensity, I, of a flow is defined by Russo and Basse [28] as I = u ′ RMS U∞ , (2.26) where u′ RMS = √ 1 3u ′ iu ′ i is the root-mean-square of the turbulent velocity fluctuations and U∞ is the free stream velocity. The scaling of turbulence intensity with Reynolds number is investigated by Russo and Basse [28] for compressible and incompressible flows in smooth pipes using both measurements and simulations. The results are fitted to the power-law expression I = aReb, (2.27) where different values for the constants a and b are recommended for different flows. The turbulent kinetic energy is defined as [20] k = 1 2u ′ iu ′ i. (2.28) By combining equations 2.26 and 2.28 the turbulent kinetic energy can be expressed as k = 3 2(IU∞)2. (2.29) The dissipation, ε, can then be expressed as ε = k 3 2 l , (2.30) where l is the turbulent length scale [23]. The specific dissipation, ω, is then estimated from [20] ω = ε Cµk , (2.31) where the value of the constant Cµ = 0.09. The Boussinesq assumption allows the Reynolds stresses to be modelled with a turbulent viscosity. It has the same dimensions as the laminar viscosity, and from dimensional analysis the turbulent viscosity can be expressed as [20] νt = Cµ k2 ε . (2.32) 13 2. Theory 2.6 Wall Treatment The height of the first cell in the wall-normal direction needs to be considered when using turbulence models. The best way to capture the wall effect accurately would be to simulate the whole boundary using small cells close to the wall. This is not always practical with respect to computational time. Wall functions are thus in- troduced. They estimate the flow close to the wall by assuming that it behaves like a fully developed turbulent boundary layer. Using this formulation, a coarser grid close to the wall can be used, and as a result, less computationally expensive simulations are achieved [23]. A suitable wall mesh grid refinement is found by analysing the dimensionless quan- tity, y+. This quantity is used to describe the mesh resolution in the vicinity of the wall [29]. The boundary layer can be divided in three main regions; viscous region, buffer region and logarithmic region. According to Davidson [20], the viscous region spans y+ ≤ 5 and the logarithmic region y+ ≥ 30. In between those two regions is the so called buffer region (5 ≤ y+ ≤ 30). According to literature, when resolving the boundary, the most appropriate value is y+ ≈ 1, while when applying wall functions it should be around y+ ≈ 30 for high Reynolds models, [29], and below y+ = 5 for low Reynolds models [30]. The y+ of the cell closest to the wall is calculated as y+ = yuτ ν , (2.33) where y is the distance from the wall and uτ is the friction velocity [20]. The friction velocity can be expressed as uτ = ( τw ρ )1/2 , (2.34) where τw is the wall shear stress. The general expression for the wall shear stress is τw = µ ∂ū ∂y ∣∣∣∣∣ w , (2.35) where the subscript w denotes the wall. The y+ can only be evaluated in the presence of a flow field at a post processing stage. It is however good practice to have an initial estimate of the cell height y suitable for the flow. It is estimated by calculating the y+, using a formula for the wall shear stress. The wall shear stress is then expressed as [31] τw = 1 2CfρU 2 ∞, (2.36) where Cf is the skin friction. There are several ways of estimating the walls skin friction, one is to use the Schlichting relation for a turbulent boundary layer over a flat plate where [32] 14 2. Theory Cf = 0.455 (log10Re) . (2.37) As this is only an estimate, a good mesh grid convergence study is needed to find the appropriate cell height at the wall for each application. 2.7 Optimisation The goal with the deliberate design of an application is often to maximise or min- imise some relevant aspect of it, which can be called the objective. Alterations can be applied to certain aspects of the design, hereafter called variables, in order to change the performance or functionality of the application with regards to the objective. The objective can thus be described as a function of the variables, an objective function. Often, several aspects of the application are relevant and an ob- jective function that includes all the desired aspects then has to be formulated as a combination of multiple individual objectives. These objectives are often conflicting, like low cost and high performance, and sophisticated methods are often needed to evaluate the impact of any change in the variables on all the objectives. One such method that facilitates complicated design is to apply mathematical optimisation of the objective function. Additionally, there is seldom complete design freedom, as the variables and objectives are constrained by limits, e.g. maximum size, minimum performance and maximum cost. There are several ways of formulating multi-objective optimisation problems. One way is to use a weighted sum approach. When implementing this method, the objec- tive function is constructed from a combination of weighted individual objectives. The weight determines whether or not one of the objectives should be dominant during the optimisation [14]. If written on the negative-null form, which indicates that the constraints are fulfilled when their function value is negative, the problem can be formulated as minimize : f(x) = w1f1(x) + w2f2(x) Subject to : g1(x1, x2) = f2(x)− f2,max(x) ≤ 0 g2(x1) = x1 − x1,max ≤ 0 g3(x1) = −x1 + x1,min ≤ 0, (2.38) where f is the combined objective function of the individual objectives f1 and f2, w1 and w2 are the weights, x is the state vector of the variables x = [x1, x2, ..., xn] and g1 to g3 are the constraints. For many applications there is not only one optimum, there are often multiple local optimums complementing the global optimum. The maximum or the minimum ob- jective of all local optimum in a design space represents the global optimum, which is further explained in Papalambros and Wilde, [33]. In the present work the global optimum is the one of interest. 15 2. Theory As mentioned in Chatper 1, so called evolutionary algorithms are popular in CFD related optimisations. These type of algorithms mimic natural evolution by evalu- ating a population of points and then combines the most prominent points to create the next generation. They are good at finding truly global optimal values. The most optimal value of an objective function that is constructed from multiple ob- jectives can often be reached by several state configurations, there is a trade-off between the objectives. Since evolutionary algorithms works with entire genera- tions, and not one single point at a time, these algorithms are good for finding the trade-off between objectives. For the same reason, evolutionary algorithms require many function evaluations before converging to the optimal solution and are com- putationally expensive [34]. A less expensive alternative is to use a surrogate model for the objective function. The true objective function is then approximated with a model function, which is optimised instead. The Efficient Global Optimization (EGO) algorithm is such a method, and is used in the present work. 2.7.1 Efficient Global Optimisation Algorithm The EGO algorithm was first developed by Jones et. al. [3] and is a derivative-free, global, space filling algorithm that effectively evaluates the design space by refining the division of the variables in areas with high expected improvement and high uncertainty. First, the objective function is evaluated in a number of initial points. Kringing modelling is then applied to approximate the response function based on evaluated points, and its variance based on the spacing of the points. The objective function is modelled by the second degree polynomial f(x) ≈ G(x) = β0 + β1x + β2x2 + Z(x), (2.39) where the constants of the vectors βi are taken as the mean of the responses of previous function evaluations and Z is a Gaussian process that describes the model deviation from the true function [4]. An expected improvement function (EIF) is then formulated based on the current best design x∗ and the predicted values from the Kringing model as EIF (x) = max(G(x∗)−G(x), 0). (2.40) The point which maximises the value of the EIF is selected as the next infill point at which to evaluate the objective function. The Kringing model is then updated with this point and the process restarts. Once the EIF reaches a low enough value, the process stops [4]. 16 2. Theory 2.8 Cavitation Cavitation is a phenomenon that occurs when the static pressure in a liquid is reduced below the vaporisation pressure. The liquid then vaporises and transforms into gas bubbles. In hydraulic turbines this mainly occurs in vortices or close to the blades, in zones where the flow is accelerated. This can cause severe structural and performance problems [17]. When the gas bubble is transported out of the accelerated zone, it collapses due to the surrounding static pressure. This can cause liquid micro jets and shock waves which excert large forces on nearby solid surfaces resulting in pitting and erosive damage. The bubbles also constrict the liquid flow and thereby affect the performance. Turbines can be designed to take advantage of cavitation, so called super-cavitating turbines, but otherwise the phenomenon should be avoided [17]. In a highly pressurised system, like in the present work, the static pressure does not drop below the vaporisation pressure and cavitation does not occur. 17 2. Theory 18 3 Methods As previously explained in Chapter 1, the main goal of the project is to develop an automated optimisation procedure by coupling together several open source software and to use the method to improve the design of a turbine in fresh water pipelines. A study is carried out in the early stages of the project to determine what software are useful for the task. The method generated in the present work is based on a Ruby code to generate the geometry, CfMesh to generate the mesh grid, OpenFOAM to simulate the flow and Dakota to create a loop and perform optimisation. A schematic representation of the optimisation procedure is shown in Figure 3.1. The loop is initiated by Dakota in the first box where the structure for the rest of the procedure is defined. The algorithm then continues by creating the geometry of the first design with the Ruby scripts. The geometry is then sent to CfMesh that creates a mesh grid of the fluid domain around the current design. OpenFOAM then simulates the flow in the system and quantifies the performance of the turbine. If the optimisation is in the learning phase, this process repeats until the optimisation period is reached. If the learning process is completed the results from the simulations are used to update the response surface which models the system. The surrogate model of the response function is then optimised to find a new infill point. If the performance of the new point is predicted to be an improvement the process repeats, if not, the optimisation is complete. For further information about how the communication between the software is handled, see the folder structure in appendix A. Dakota:   Start Ruby:   Geometry  creation   CfMesh: Computational  grid  creation OpenFOAM: CFD  simulations Is  the  sampling complete?   Dakota:  Create/update   surrogate  model True False Is  there  an  expected   improvement? Yes No n=1 n=n+1 n≥k n