Feasibility test of the simulation of a running turbine with openfoam

FEASIBILITY TEST OF THE SIMULATION OF A RUNNING TURBINE WITH OPENFOAM DANIELA FIGUEIREDO GOMES DE OLIVEIRA Dissertação submetida para satisfação parcial dos requisitos do grau de MESTRE EM ENGENHARIA CIVIL — ESPECIALIZAÇÃO EM HIDRÁULICA Orientador: Professor Doutor João Pedro Pêgo Coorientador: Elena-Maria Klopries, M.Sc. RWTH JULHO 2015 MESTRADO INTEGRADO EM ENGENHARIA CIVIL 2012/2013 DEPARTAMENTO DE ENGENHARIA CIVIL Tel. +351-22-508 1901 Fax +351-22-508 1446  [email protected] Editado por FACULDADE DE ENGENHARIA DA UNIVERSIDADE DO PORTO Rua Dr. Roberto Frias 4200-465 PORTO Portugal Tel. +351-22-508 1400 Fax +351-22-508 1440  [email protected].pt  http://www.fe.up.pt Reproduções parciais deste documento serão autorizadas na condição que seja mencionado o Autor e feita referência a Mestrado Integrado em Engenharia Civil - 2012/2013 - Departamento de Engenharia Civil, Faculdade de Engenharia da Universidade do Porto, Porto, Portugal, 2013. As opiniões e informações incluídas neste documento representam unicamente o ponto de vista do respetivo Autor, não podendo o Editor aceitar qualquer responsabilidade legal ou outra em relação a erros ou omissões que possam existir. Este documento foi produzido a partir de versão eletrónica fornecida pelo respetivo Autor. Feasibility test of the simulation of a running turbine with OpenFOAM To my Grandfather. Knowing is not enough, we must apply. Willing is not enough, we must do. Johann Wolfgang von Goethe Feasibility test of the simulation of a running turbine with OpenFOAM Feasibility test of the simulation of a running turbine with OpenFOAM i ACKNOWLEDGMENTS First of all, to my scientific mentor and the person who contributed the most to this work, Elena Klopries, M.Sc. RWTH. For all the never ending help, patience and teaching, a sincere thank you. It was for me a privilege to have learned from you. To Professor Dr. João Pedro Pêgo, for all the availability, since the first moment, to follow the development of my dissertation. For contributing with wisdom and knowledge when most necessary. To all the professors in the hydraulics department of the Faculty of Engineering of the University of Porto. In the end of the year, I was certain I had chosen the right specialization. To my friends and especially Marta Morais, for being always and unconditionally a source of inspiration and encouragement in good times, but, above all, in overcoming times. To Nuno Quelhas, my partner in this Erasmus experience. For his constant support and friendship a heartfelt thanks. To my grandparents, for their support, not only emotionally, but also financially. To my mother, for being the engine of our home, for all the effort made so that this path and especially this dissertation were possible to be performed under the conditions they were. Finally to my brothers, Ricardo and João. That the pride they can feel in me could be equivalent to the immensity they symbolize in my life. Feasibility test of the simulation of a running turbine with OpenFOAM ii Feasibility test of the simulation of a running turbine with OpenFOAM iii AGRADECIMENTOS Primeiramente à minha mentora científica e à pessoa que mais contribuiu para este trabalho, Elena Klopries, M.Sc. RWTH. Por toda a interminável ajuda, pela paciência e por todo o conhecimento transmitido, um sincero obrigado. Foi para mim um privilégio ter aprendido consigo. Ao Professor Dr. João Pedro Pêgo, por toda a disponibilidade desde o primeiro momento para acompanhar o desenvolvimento da minha tese. Por ter contribuído com sabedoria e conhecimento na altura certa. A todos os professores no departamento de hidráulica da Faculdade de Engenharia da Universidade do Porto. No final do ano, eu estava certa de que tinha escolhido a especialização certa. A todos os meus amigos, e em especial à Marta Morais, por serem sempre e incondicionalmente uma fonte de inspiração e de incentivo quer nos bons momentos, quer, acima de tudo, nos momentos de superação. Ao Nuno Quelhas, meu companheiro desta experiência de Erasmus. Pelo seu constante apoio e amizade, um agradecimento de coração. Aos meus avós, pelo apoio não só emocional, como financeiro. À minha mãe por ser o motor do nosso lar, por todo o esforço que fez para que este percurso e em especial esta dissertação fosse possível de ser realizada nas condições que foi. E por fim, aos meus irmãos, Ricardo e João. Que o orgulho que possam sentir em mim seja equiparável à imensidão que simbolizam na minha vida. Feasibility test of the simulation of a running turbine with OpenFOAM iv Feasibility testo f the simulation of a running turbine with OpenFOAM v ABSTRACT Field and laboratory studies reviewing the passage of fish through hydroelectric turbines have been object of research over the years. The biggest causes of fish mortality when passing a hydroelectric turbine are described as being the variation of pressure, cavitation, shear stress as well as strike and grinding. Most of the models used to calculate fish mortality are based on the probability of the injury to happen depending on parameters like flow, length of the fish or the angle of the turbine’s vanes. But most of these models are so-called black-box-models where the real, physical processes are not shown. This is where computational fluid dynamics plays an important part. It is possible to model and simulate different structures and in this way also the flow that may occur. Like this, it is possible to model what happens to a fish when passing a turbine, showing the real processes that happens during the passage. In this work, a four blade propeller turbine is modelled with the software OpenFOAM - Open Field Operation and Manipulation. The main focus of the work is on the simulation of the turbine’s movement. Two different approaches with regard to the simulation of the movement were made in order to better comprehend which one could be more suitable and more economical with regard to the given resources of the computers. The first one focused on an oscillating movement and the second one was done with a rotational movement. Both were solved resorting to turbulent resolution model SST k-ω. KEYWORDS: fish mortality, hydraulic turbines, Computational Fluid Dynamis, OpenFoam. Feasibility test of the simulation of a running turbine with OpenFOAM xii Fig. 33 – Values of pressure in laminar flow ...................................................................................... 55 Fig. 34 – Values of velocity in laminar flow ........................................................................................ 55 Feasibility test of the simulation of a running turbine with OpenFOAM xiii TABLE INDEX Table 1 – Fixed values for the Standard Model k-ε. ........................................................................... 28 Table 2 - Fixed values for the Standard Model k-ω ............................................................................ 30 Feasibility test of the simulation of a running turbine with OpenFOAM xiv Feasibility test of the simulation of a running turbine with OpenFOAM 1 1 INTRODUCTION With the improvement of electric and hydraulic technology, it’s almost impossible for fish to do their natural migration through a river without having to cross human made barriers such as dams and barrages. The studies regarding this fact are in constant development, not only in the field but also in laboratory. The fish are mechanically and hydraulically injured and because of this, the studies have special emphasis on tidal schemes in operation. Although the types of fish sufferer injuries are well documented and explained, their causes are not that clear. The specific hydraulic conditions that lead to these injuries and sometimes to mortality are yet not well known. So far, most of the known mortality models have been built on probability studies. This is where softwares like OpenFOAM can have a major difference. The possibility to model the most diverse structures as well as different flow models is the reasons why this software is adequate to give correct answers, regarding the purpose of this work. OpenFOAM or Open Field Operation and Manipulation is a software that, as the name indicates, is an open source software meaning that, besides the normal use of the libraries, it is also possible to change the codes, adapting them to the better form of the problem in question. The main purpose of this assignment is to test whether it is possible or not to build a numerical model of a running turbine resorting to this software. This work is divided in five chapters. The first one is merely introductory. The second one makes a detailed explanation about fish mortality, the studies made in this area and the influence of the turbines in the injuries and the mortality. The third chapter explains how computational fluid dynamics works in a way to better understand how all the boundaries were chosen to this work. The fourth chapter describes all the procedures regarding the modelling of the turbine. Finally, the fifth chapter exposes the conclusions about this model. Feasibility test of the simulation of a running turbine with OpenFOAM 2 Feasibility test of the simulation of a running turbine with OpenFOAM 3 2 THEORETICAL FRAMEWORK REGARDING FISH PASSAGE THROUGH RUNNING TURBINES 2.1. FISH MORTALITY ASSOCIATED WITH RUNNING TURBINES There are many different types of fish species that have migration as an important procedure on their basic needs, such as to get food and to breathe. This migration can differ on a scale of time and distance, varying from daily to annual and from meters to kilometres. During the migration process, fish usually need to go through some barriers. Some of these barriers are natural, like sandbars, landslides, waterfalls, and boulder cascades. Others are human made, such as dams and barrages, which usually are equipped with power schemes such as turbines. This barrier works, most of the times, as physical stressor that can have implication on fish mortality. If the studied species have economic relevance, this could be a serious problem and, therefore, the need for more efficient fish passage is a requirement. The field and laboratory studies on fish mortality are still being reviewed because it seems that little is yet known about which hydraulic conditions within the turbine directly affect fish. Although, the mechanical conditions responsible for the rates of mortality are already known, and they are:  The abrupt changes in pressure;  Water turbulence;  Shearing currents;  Mechanical contact with blades. On the other hand, the rates of mortality vary with fish size, turbine characteristics, such as head of water and runner diameter, and the operating conditions of the power scheme. The studies made so far on this subject are always related with the turbine passage, but it would be stimulating if the relationship between turbine characteristics, size and species of the fish could be studied. It is in fact a complex interaction with many variables. 2.2. CHARACTERISTICS OF RUNNING TURBINES The main causes of fish mortality are now subject of further explanation. The effect that pressure has on fish mortality is highly dependent on the turbine design and consequently on the head water. First, it is undoubtedly necessary to explain the operation in turbines and its consequences to the phenomena that causes injury and mortality on fish. Feasibility test of the simulation of a running turbine with OpenFOAM 4 Hydraulic turbines can be separated in two main categories: impulse and reaction turbines. The difference between these two types is mainly given by the action of the water and consequently the mechanical energy that is produced, later converted to electrical energy. Impulse turbines use the water movement and the associated kinetic energy of a high-velocity jet discharging at atmospheric pressure. (Turbak et al 1981). Fig. 1 – Operating system of an Impulse Turbine (Source: [1]) With this system, shown in figure 1, suction on the down side of the turbine does not occur because water flows through the bottom of the turbine after hitting the bucket. This type is usually used for high head and low flow power schemes. The PELTON turbine is one example of impulse turbines and the most used one. Pelton turbines are one of the most important turbines when it comes to the conversion of hydraulic energy into electricity, especially in mountain areas. As Dhakan P.K. and Chalil, A. B. P. (2013) refers this is due to the fact that this type of turbines can be easily used in places where the altitude difference between the water source and the location of the turbine is considerable. Pelton turbines can be used in large scale hydro installation for heads that can go from 20 meters to 150 meters. Usually these types of turbines are not chosen when dealing with lower water heads since the rotational speed becomes slower and the runner of the turbine required needs to be larger and it’s difficult to manage. In high water heads the flow rate tends to be lower going from 0,005 m3/s in the smaller systems till 1 𝑚3/s on larger systems. In consequence of these values the obtained power results can vary in a scale from a few kW up to hundreds of MW’s in those larger systems. Feasibility test of the simulation of a running turbine with OpenFOAM 5 In hydraulic systems with Pelton turbines, the reservoir is connected to a penstock head that is posteriorly linked to a penstock where the water flows and is directed by a nozzle against the buckets around the runner in a form of a thin jet. These turbines can be arranged in two forms: in a horizontal shaft or in a vertical shaft. The horizontal position allows the use of several runners leading to a higher specific speed and therefore a higher operational speed. The vertical position allows a multi jet construction leading to an improvement of the speed. In reaction turbines, the flow system is restricted to a closed conduit and at any point it gets in contact with air. This system is strictly closed from the headwater to the tail water. When water approaches the runner it is provided with pressure energy and kinetic energy. The first is due to the density of the water above this till the headwater surface and the second is due to the velocity. PROPELLER (Bulb Turbine, Straflo, Kaplan) and FRANCIS are some of the examples of reaction turbines. Fig. 2 – Operating system of a Reaction Turbine (Source [1]) As Turbak et al (1981) refers, fish mortality research has been almost entirely made with reaction turbines and for this reason those type of turbines will be object of a more extensive explanation. These turbines have a simple design that begins with an intake conduit that guides the water into the main casing where it gets in contact with the nozzles that are attached to the rotor. The acceleration caused by the water leaving the nozzles makes a reaction force on the pipes and consequently the rotor starts moving in the opposite direction of the water. The main shaft in the rotor can have different positions according to the different type of turbine. The Straflo uses only horizontal flow but Francis and Kaplan can use both horizontal and vertical flow. One of the most important and commonly used Propeller turbine is the Kaplan turbine. This turbine has a propeller with adjustable blades inside a tube. It is also known as an axial-flow turbine meaning that the flow direction doesn’t change while crossing the rotor. Feasibility test of the simulation of a running turbine with OpenFOAM 6 Theoretically, Kaplan turbines can work through an extensive spectrum of heads and flow rates, but, in practice, there are other types of more effective turbines for higher heads of water. Due to this, Kaplan turbines are usually chosen for lower water heads going from 1.5 to 20 meters with high flow rates from 3 m3/s to 30 m3/s. Hydroelectric power plants with Kaplan installations for these flows have output power values from 75 kW to 1MW. Fig. 3 – Basic layout of a Kaplan turbine (Source: [2]) These turbines have usually in their structure a nose cone, a rotor with adjustable blades and a vertical driveshaft. The inlet guide-vanes that lead the water to the turbine can regulate the flow rate in the turbine. This allows to fully stop the work of the turbine by closing completely the guide-vanes. It also allows, depending on the position of the guide-vanes, to control the flow and consequently guarantee that the oncoming flow can hit the rotor in the most efficient angle conducting to the highest efficiency. The blades can also be adjustable according to the dimension of the flow: a flat outline to very low flows and a heavily-pitched outline for high flows. Allowing the adjustment of the inlet guide-vanes leads to a big opening of the flow operating range and by adjusting the rotor blades it’s possible to optimize the turbine efficiency. Francis turbines are another type of really common reaction turbines and one of the most preferred since they are also the most reliable turbines used in hydroelectric power stations. Turbak et al (1981) states that the contribution of these turbines to the global hydropower capacity is around 60 per cent, basically because the efficiency is higher under a huge range of different operating conditions. Installations with Francis turbines can operate in heads from 30 to 300 meters with flows from 10 to 700 m3/s and the number of blades in the turbine can vary from 14 in lower heads to 20 for higher heads. In a simplified way, the Francis turbine is composed of a runner with complex shaped blades. The flow enters the turbine radially and leaves it axially as it is possible to notice on Fig. 4. Nevertheless there is a particularity about the blades of the Francis turbine: the cross-section is shaped like a thin airfoil. This means that when the inlet flow crosses the blades there is a low pressure in one side and a high pressure on the other side, leading to a lift force. Feasibility test of the simulation of a running turbine with OpenFOAM 7 It is possible to notice, as well, that the blade has a bucket kind of shape directed to the outlet. After the water flow hits the blades, it produces an impulse force to be able to leave the runner. The runner rotates due to both impulse and lift force. Fig. 4 – Basic layout of a Francis Turbine (Source: [3]) This may sound like a contradiction since Francis is a reaction turbine. The truth is that it’s not a pure reaction turbine and part of the force to move the runner is obtained by an impulse action. The runner is connected to a shat to posterior gain of electric energy production. As Turbak et al (1981) refers, fish mortality research has been almost entirely made with reaction turbines. These turbines have a simple design that begins with an intake conduit that guides the water into the main casing where it gets in contact with the nozzles that are attached to the rotor. The acceleration caused by the water leaving the nozzles makes a reaction force on the pipes and consequently the rotor starts moving in the opposite direction of the water. Feasibility test of the simulation of a running turbine with OpenFOAM 14 Feasibility test of the simulation of a running turbine with OpenFOAM 15 3 Computational Fluid Dynamics – CFD 3.1. INTRODUCTION TO COMPUTATIONAL FLUID DYNAMICS Computational Fluid Dynamics is the field in fluid mechanics that studies the resolution and analysis of the flow using algorithms and numerical methods. Fluids numerical modelling started, almost exclusively, with aeronautical study, but nowadays it is used in several areas of engineering and physics. It is employed in areas such as aircraft, turbo machinery, car and ship design and it can also be applied in areas such as meteorology, oceanography, astrophysics, oil recovery and also architecture. In its beginnings, it was a combination of physics, numerical mathematics and some computer science simulating fluid flows. With the advance in computer technology, it also occurred an advance in computational fluid dynamics. From simulating transonic flows based on the solution of non-linear potential equations, the CFD applications evolved to the first two-dimensional solutions and afterwards to three dimensional solutions. (Blazek. 2001) With the promptly increasing computers capacity and speed, it was possible to start developing some more robust simulations with inviscid flows passing complete aircrafts configurations or inside turbo machines. With this development, it came the need to start more demanding simulation with viscous flows governed by the Navier Stokes equations. As a result, turbulence models such as Reynoldsaverage Simulation (RAS) and the Large Eddy Simulation (LES) appeared. (Blazek. 2001) With the consequent advance of the complexity in flow simulations, the grid generation also had to follow this development. The progress started first with simple structured meshes constructed either by algebraic methods or by using partial differential equations. But with increasing geometrical complexity of the configurations, the grids had to be broken into a number of topologically simpler blocks and then the multiblocks approach appeared. Nevertheless, the generation of a complex structure using multiblocks grid would still take a long time in the order of weeks or even months. Due to this fact, the research on the field of unstructured grid generators and consequently its solver took a big growth by promising to reduce the setup times. (Blazek. 2001) With the use of the Navier-Stokes equations, the requirements for the meshes were more demanding. Prisms and hexahedra grids in viscous flows began to be used, improving the solution accuracy and saving the number of elements, faces and edges. The biggest problem on using fluids numerical modelling software is the greatest complexity and unpredictability on knowing how a fluid will behave since it has a non-linear behaviour. Despite this problem, with the appearance of more powerful computers, the properties of the fluids are more easily understood. Another problem is that the obtained solutions are only as correct as the physical models they are based on. Feasibility test of the simulation of a running turbine with OpenFOAM 16 The simulation with physical models involving fluids is very expensive. A physical model, depending on its complexity, needs a large set of material to better represent the prototype but it also needs measurement material and specialized staff. In spite of all these disadvantages, Bakker (2002) refers some important advantages, such as the cost reduction, since the computers became more advanced and faster in order to obtain results. One of the biggest advantages is the ability to simulate real conditions, such as flow and heat transfers and also to be able to control the physical process, isolating, if necessary, one specific phenomena. With the advance of technology, it is already possible to run some of these experiments on personal computers. Also with the growing of CFD simulations, there is already in the market some software that can be used without any cost to the user, such as OpenFOAM. Despite all the advantages Versteeg and Malalasekera (1995) refer that the numerical errors and the boundary conditions are two of the biggest problems when programming in CFD. Factors like the ambient temperature or the existence of air bubbles are not always included on the modelling process and that leads to calculation errors. Numerical errors are the most frequent and they usually occur because of the algorithms used in CFD software which use really complex mathematical equations. Such as the experiments and simulations in physical models, the CFD simulations need a critical review of the results in order to be validated. 3.2. SPECIAL FEATURES OF THE CFD TOOLS The base knowledge behind numerical modelling of fluids is that the numerical algorithms are able to simulate the fluid behaviour. Versteeg and Malalasekera (1995) refer three main elements: a preprocessor, a solver and a post-processor. 3.2.1 PRE-PROCESSOR This stage is where all the physical properties and information about the flow are defined. The data that is necessary is:  definition of the region in study;  definition of the geometry of the study case through meshes;  definition of both physical and chemical properties of the model in study;  definition of the fluid properties;  definition of boundary conditions. The definition of the mesh in a CFD problem is the first step to get the answer for a flow problem. The velocity, pressure and temperature, among other characteristics, are defined at the various nodes inside every cell of the mesh. As it is easy to understand, the precision of the result is largely dependent on the number of cells of the mesh. The larger the number of cells, the biggest is the accuracy of the solution. Efficient meshes should not be uniform: thinner in areas where big variations of characteristics happen and larger where this doesn’t happen. The evolution on CFD is to try to make the programs with the ability to create these different meshes by themselves but for now it is still up to the user to make the better judgement about the correct mesh to use. Feasibility test of the simulation of a running turbine with OpenFOAM 17 3.2.2. SOLVER It’s possible to obtain solutions by three different types of resolution: finite difference, finite element and spectral methods. As Versteeg and Malalasekera (1995) refer, the three are governed by the same assumptions: the ability to approximate the unknown flow variables with simple functions, to transform a continuous distribution in discrete units and then to approximate these into the governing flow equations and their respective mathematical manipulations. The finite difference method has the approximation to the Navier-Stokes equations in their simplified form and it’s of easy implementation. On the other side, it generates problems along the curved boundaries, such as the difficulty in achieving stability and reaching a convergence in the analysis. The mesh adaptation is also difficult. The finite volume method is an approximation of the Navier-Stokes equations as a system of conservation equations being the biggest advantage the fact that the special discretisation is carried out directly in the physical space. The disadvantages are that, like the previous one, it has difficult stability and convergence analysis and it’s difficult to adapt to unstructured meshes. The spectral method evaluates the spatial derivatives resorting to Fourier series or one of their generalizations. It is possible to obtain the biggest advantage from spectral methods if the functions used are periodic and the grid points uniformly spaced. Any change in geometry or boundary conditions requires a considerable change in the method, making these methods relatively inflexible. As it’s possible to conclude, the biggest difference between these three streams is the way the flow variables are approximated and consequently the discretisation process. These methods will be further explored in chapter 3.7. (Ferziger, J.H. 2002) 3.2.3. POST-PROCESSOR A big development has been done in the post-processing field. Almost all of CFD software is equipped with different visualization tools that can be divided into two forms: graphically or alphanumerically. In the graphical tools, there are vector plots, geometry and grid display contours, iso-surfaces, flowlines and animations. In the alphanumeric tools, there are some such as integral values, drag, lift and torque calculations, averages and standard deviations. These tools can facilitate the perception that the user has about work in progress. Fig. 8 – Flow around a cylinder: grid. (Source: Baker, 2002) Feasibility test of the simulation of a running turbine with OpenFOAM 18 For instance, for the mesh displayed in Fig. 8, representing the flow around a circular cylinder, it is possible to display for each one of the time-steps the vector plots and the surface variations regarding the velocity magnitude as shown in Fig. 9. The same is valid for the pressure distribution, which is shown in Fig. 10. Fig. 9 – Velocity vectors regarding the flow around a cylinder. (Source: Baker, 2002) Fig. 10 – Iso-surface of pressure regarding flow around a cylinder (Source: Baker, 2002) Feasibility test of the simulation of a running turbine with OpenFOAM 19 As it’s possible to acknowledge, the post-processor is as important as the rest of the procedure allowing the user to rightly obtain conclusions. 3.3. MESHES Analytical solutions are not usually a possible answer for the differential equations used to solve fluid flow and heat transfer problems due to the complexity and non-linearity of the governing equations. In order to better analyse the fluid flows, the domains are divided into smaller subdomains. These subdomains are characterized by geometric shapes like hexahedra and tetrahedral in 3 dimensions and quadrilaterals and triangles in 2 dimensions. The biggest help on doing this type of division is that the equations are now solved inside each one of these smaller domains. The methods used to solve this equation were previously mentioned in the chapter 3.2.2 - Solver. These subdomains or cells must have a logical continuity so the solution obtained across the different points can make sense when put all together. All these elements are called a mesh or a grid. Meshes can be rearranged in different forms and they can be structured or unstructured. Structured meshes (see Fig. 11) follow a i,j,k convention as opposite to unstructured meshes that don’t have any convention (see Fig. 12). It is also possible to have multiblocks that consist in a group of meshes each of which can be structured or unstructured. Fig. 11 – Structured mesh (Source: Yoganathan et al. 2004) Fig. 12 – Unstructured mesh (Source: Yoganathan et al. 2004) Feasibility test of the simulation of a running turbine with OpenFOAM 20 As it’s possible to notice, the biggest difference on the two types of meshes is that the first one has the topology of a square grid of parallelepipeds, but the same does not happen on the second. These two types of meshes can be unpredictably difficult to generate, especially because a mesh that can be used to evaluate fluid flows, for instance, needs to satisfy some limits that are hard to achieve. A mesh must represented as accurate as possible the object in study, but always with the restriction on size and shape of the elements of the mesh. Modelling objects with complex shapes, as is the case of a turbine, obligates many times the use of curved surfaces. And with these surfaces two types of boundaries appear: exterior and interior boundaries. As Shewchuk (1997) refers, exterior boundaries separate meshed and unmeshed portions of space, and are found on the outer surface and in the internal holes of a mesh. On a different point of view, internal boundaries appear within meshed portions of space, and enforce the constraint that elements may not pierce them. These boundaries are typically used to separate regions that have different physical properties. The size of the elements is an important control that should be also done while generating the mesh. This means that the variation of size of the elements should be done in a short distance. And why should this happen? The size of the elements has an influence on the finite element simulation. The smaller the elements are and more densely packed, the biggest is the accuracy of the results, but at the same time this is going to increase the required time for the computer to solve the problem. Of course that the physical phenomena in study has a highly influence on the size chosen for the elements. For instance, in a fluid flow simulation, in areas where turbulence occurs, the elements should be smaller than in areas of steadiness. It is possible to have the same size of elements through all the domain in study but that can lead to excessive computational demands. One way to solve this problem is to use a mesh generator that can make the gradation from small to large sizes of the elements inside the mesh. Mesh generator algorithms have the main goal to generate a mesh with as few elements as possible offering the option to refine some parts of the mesh where the elements are not as small as required. In conclusion, the last aim of generating a mesh is to try to make the elements as “round” as possible in shape because angles (as big or small they can be) degrade the quality of the solution. The quality of each element is still a relative subject, highly dependent on the model in study, the type of numerical method used and the polynomial degree of the piecewise functions used to interpolate the solution over the mesh. Feasibility test of the simulation of a running turbine with OpenFOAM 21 3.4. MAIN PRINCIPLES ABOUT FLUID DYNAMICS Before proceeding to the study of fluid dynamics, it is important to have some ground notions of how fluids work. Fluids are substances whose molecular structure is not affected by external shear forces. Common dynamic forces like pressure differences, gravity, shear, rotation and surface tension are the responsible for the fluid flow. The software that runs numerical modelling of fluids has its base ground on some fundamental conservation laws of fluids. Those are:  Mass conservation (Continuity law)  The rate of change of momentum equals the sum of the forces on a fluid particle (Second law of Newton)  Conservation of energy in a particle (First law of thermodynamics) The numerical formulation of these principles can be written in a differential way. When the study is done in an integral way, it’s necessary to consider the change of mass, movement or energy inside the volume where the inflow and the outflow of the fluid occur. In a differential study, it is used the Stokes equation where the conservative laws are applied to an infinitesimal volume. It’s in this last assumption that the majority of the fluid dynamics modelling is based on. The velocity of a flow is also an important characteristic and allied with the inertia of the fluid it defines if the flow is laminar or turbulent. As the velocity increases and it leads to instability and to a more random type of flow the flows goes from laminar to turbulent. 3.4.1. MASS CONSERVATION – LAW OF CONTINUITY The continuity law states that, in the mass conservation, the rate of increase of mass in the fluid element needs to be equal to the net rate of flow of mass into the same fluid element. This can be described by the following equation: Where ρ defines the fluid density, u 󰇍 , the absolute velocity of the fluid and n 󰇍 , the unit vector normal to the element with the area dA. The indices cv and cs have the meaning of control volume and control surface, referring to the fluid element. Defining a fluid element with sides such as δx, δy and δz , as its shown in the Fig. 13, can make the understanding of mass conservation easier. 𝜕 𝜕𝑡∫𝜌𝑑 𝑉 𝑐𝑣 +∫𝜌𝑢 󰇍 𝑐𝑠 ×𝑛 󰇍 𝑑𝐴 =0 (3.1) Feasibility test of the simulation of a running turbine with OpenFOAM 22 Fig. 13 – Fluid element for conservation laws (Source: Baker, 2002) Taking into account the mass flow across the face of the element in consideration, it’s possible to define the flow as it is shown in Fig. 14. Fig. 14 – Mass flows in and out of fluid element. (Source: Baker, 2002) Applying the Gauss Theorem to the equation 3.1, it’s possible to obtain: This equation defines now the unsteady and three-dimensional mass conservation at a certain point for a compressible fluid. For an incompressible fluid (for example, a liquid one) the second part of the equation which 𝜕𝜌 𝜕𝑡+𝑑𝑖𝑣(𝜌𝑢 󰇍 )=0 (3.2) Feasibility test of the simulation of a running turbine with OpenFOAM 23 describes the net flow of mass across the boundaries of the element and in a more specific way called convective term is null. This is explained by the fact that the density ρ of the fluid is constant making div(ρu 󰇍 )=0. Since this study is done with an incompressible fluid, the mass conservation equation is given by: 3.4.2. CONSERVATION OF THE MOMENTUM – SECOND LAW OF NEWTON Newton’s second law describes that the rate of change of momentum of a fluid particle is equal to the sum of the forces on the same particle. Following the same coordinates description, the rates of change of the momentum per unit volume are given in the three principal coordinates: Versteeg (1995) refers that there are two types of forces interacting with the particles of the fluid: Surface forces, such as pressure and viscous forces. Body forces, such as gravity, centrifugal, Coriolis and electromagnetic forces. With this knowledge and making use of the simplification of the newton second law: 𝐹 =𝑚×𝑎 it is possible to obtain the force vector as: Where τij represents the normal and tangential stress in the i and j vector direction respectively, ρ the density of the fluid, p the pressure and fi are the acting forces on the fluid in the direction yi (in this case we will only consider the gravitational action). Making use of the coordinates previously established, it’s possible to define the mass of the fluid as: Defining ui as the velocity in the direction xi and the time as t, the acceleration can be stated as: Combining the equations 3.5, 3.6 and 3.7 leads to the Navier-Stokes equation in its non-conservative form. According to Anderson (2009), it’s possible to obtain the conservation form of this equation by applying the derivation of the conservation law to the linear momentum, obtaining the following: 𝜕𝑢 𝜕𝑥+𝑑𝑣 𝑑𝑦+𝑑𝑤 𝑑𝑧=0 (3.3) 𝑥: 𝜌𝜕𝑢 𝜕𝑡 ; 𝑦: 𝜌𝜕𝑣 𝜕𝑡 ; 𝑧: 𝜌𝜕𝑤 𝜕𝑡 (3.4) 𝐹𝑖= (−𝜕𝑝 𝜕𝑥𝑖+∇𝜏𝑖𝑗)+𝜌𝑓𝑖 (3.5) 𝑚= 𝜌 𝑑𝑥 𝑑𝑦 𝑑𝑧 (3.6) 𝑎𝑖= 𝐷𝑢𝑖 𝐷𝑡 (3.7) 𝜕(𝜌𝑢𝑖) 𝜕𝑡 +∇×(𝜌𝑢𝑖)=(−𝜕𝜌 𝜕𝑥𝑖+∇𝜏𝑦)+𝜌𝑓𝑖 (3.8) Feasibility test of the simulation of a running turbine with OpenFOAM 30 Resorting to the previous equation, it is possible to write the two main solvers of the kinetic turbulent energy and the specific dissipation rate: Where σk, σω, β∗ and α assume the following value: Table 2 – Fixed values for the Standard Model k-ω. 𝜎𝑘 𝜎𝜔 β∗ β α 2 2 9100 ⁄ 340 ⁄ 59 ⁄ τij can be replaced by −ρ(umi′ 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 umj′ 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ) and by using the equation 3.16, it is possible to obtain the solution. The turbulent viscosity μT can be achieved through the resolution of the following: According to Versteeg (2007), the main advantage is the simple and exact way of how the behaviour of the flow near the walls is calculated, for low Reynolds numbers. For regions far away from the walls, the dissipation rate tends to zero and the viscosity turbulence starts tending to an infinite value. For these reasons, the values that are sometimes obtained by this model are not reliable and, to overcome this problem, Menter (1993) formulated a new model called Shear-Stress Transport (SST).  SHEAR-STRESS TRANSPORT MODEL k-ω. This model developed by Menter is the combination of the two previous models here explained: the kε and the k-ω model. With this junction, the main flaws from the previous models are eliminated. The SST k-ω model behaves like a k-ω model in the regions near the walls and then changes this behaviour to a k-ε model in the outer regions of the flow. Besides, the biggest difference from the previous model is that the SST k-ω model includes the transport effects in the calculation of the eddy viscosity: 𝜇𝑇=𝛼1𝑘 max (𝛼1𝜔𝑆𝐹2) (3.25) 𝜕 𝜕𝑡(𝜌𝑘)+𝜕 𝜕𝑥𝑗(𝜌𝑘𝑢𝑗)=𝜕 𝜕𝑥𝑗[(𝜇+𝜇𝑇 𝜎𝑘)𝜕𝑘 𝜕𝑥𝑗]+𝜏𝑖𝑗 𝜕𝑘 𝜕𝑥𝑗−𝛽∗𝜌𝑘𝜔 (3.22) 𝜕 𝜕𝑡(𝜌𝜔)+𝜕 𝜕𝑥𝑗(𝜌𝜔𝑢𝑗)=𝜕 𝜕𝑥𝑗[(𝜇+𝜇𝑇 𝜎𝜔)𝜕𝜔 𝜕𝑥𝑗]+𝛼𝜔 𝑘𝜏𝑖𝑗 𝜕𝑘 𝜕𝑥𝑗−𝛽𝜌𝑘𝜔 (3.23) 𝜇𝑇=ρ𝑘 𝜔 (3.24) Feasibility test of the simulation of a running turbine with OpenFOAM 31 By doing this, it is possible to predict the beginning and the size of the separation of the flow. The way it is possible to calculate the turbulent kinetic energy is the same as the previous model resorting to the equation number 3.23, but to calculate the rate of dissipation of specific energy it is necessary to apply the next equation: 𝜕 𝜕𝑡(𝜌𝜔)+𝜕 𝜕𝑥𝑗(𝜌𝜔𝑢𝑗)=𝜕 𝜕𝑥𝑗[(𝜇+ 𝜇𝑇 𝜎𝜔1)𝜕𝜔 𝜕𝑥𝑗]+𝛼𝑆2−𝛽𝜔2+2(1−𝐹1)𝜎𝜔21 𝜔𝜕𝑘 𝜕𝑥1𝜕𝜔 𝜕𝑥1 (3.26) The values of 𝐹1and 𝐹2 are calculated with the use of the equations presented below: 𝐹1=𝑡𝑎𝑛ℎ{{𝑚𝑖𝑛[𝑚𝑎𝑥( √𝑘 𝛽∗𝜔𝑦,500𝜇 𝑦2𝜔), 4𝜎𝜔2𝑘 𝐶𝐷𝑘𝜔𝑦2]}4} (3.27) 𝐹2=𝑡𝑎𝑛ℎ[[𝑚𝑎𝑥(2√𝑘 𝛽∗,500𝜇 𝑦2𝜔)]2] (3.28) Where 𝐶𝐷𝑘𝜔 is: 𝐶𝐷𝑘𝜔 =𝑚𝑎𝑥(2𝜌𝜎𝜔2 1 𝜔𝛼𝜕𝑘 𝜕𝑥1𝜕𝜔 𝜕𝑥1,10−10) (3.29) The values of 𝜎𝑘, 𝜎𝜔1, 𝜎𝜔2, 𝛽∗, 𝛽 and 𝛼1 are fixed values and are the following: Table 2 – Fixed values for the SST k-ω model 𝜎𝑘 𝜎𝜔1 𝜎𝜔2 β∗ β α1 1.0 2.0 1.17 9100 ⁄ 340 ⁄ 59 ⁄ Versteeg (2007) alerts to the fact that this precise model produces levels of turbulence extremely high in areas where the flow is static and also in areas with high acceleration. Despite this fact, all the other authors consulted present good behaviours and results when using this model. 3.6.2. LARGE EDDY SIMULATION Smagorinsky (1963) developed a method based on the Navier-Stokes equations for the simulation of big atmospheric scales. Later on, this model started to be used in engineering and nowadays it is used in the most diverse areas such as: combustion, acoustics, hydraulics, atmosphere simulations and so on. Large eddy simulations (LES) are turbulence models where the time-dependent flow equations are solved for the mean flow and the largest eddies and where the effects of the smaller eddies have an isotropic behaviour. For this reason, this last eddies can be neglected. Feasibility test of the simulation of a running turbine with OpenFOAM 32 The LES model uses a space filtration to separate the smaller vortices from the big ones. During this separation, the smaller ones get destroyed alloying the model to numerically obtain the behaviour of the flow resorting to the Navier-Stokes equations for the highest dimensions movements. 𝜕𝑢1    𝜕𝑡 +𝜕 𝜕𝑥𝑗(𝑢𝑖𝑢𝑗       )=1 𝜌𝜕𝑝 𝜕𝑥1+𝜕 𝜕𝑥1(𝜇𝜕𝑢1    𝜕𝑥𝑗)−𝜕𝜏𝑖𝑗 𝜕𝑥1 (3.30) Where 𝑢 and 𝑝 represent the velocity and pressure after filtration and 𝜏𝑖𝑗 is the residual stress it is formulated depending on the sub models of mesh used. According to Versteeg (2007), the LES model has a high precision regarding the calculation of the Reynolds stress and the terms of transportation. With this precision, there is the need of a bigger computational resolution and a much more fine mesh. The LES model needs the double of the processing comparatively with the RANS model. The need for a very fine mesh comes from the need that the simulations have to be addressed in such a way that the vortex shear stresses are more tension relieved. So that this can happen, filtering used cannot be larger than a small fraction of the fixed turbulent local scale. This scale decreased as the wall approaches and with it the width of the filter. When it comes to the resolution by the LES model, it should be noted that only a few years ago computers began to have enough processing power to simulate flow by this method. (Versteeg. 2007) 3.6.3. DIRECT NUMERICAL SIMULATION The direct numerical simulation (DNS) is based on the numerical modelling of flows through the Navier-Stokes equations without any help or with the approach of turbulence models. The full ranges of temporal and spatial turbulent scales have to be resolved, regardless of their size. This means that in the direct numerical simulation, the smallest vortex and the fastest fluctuations are settled, meaning that any movement, regardless of their irrelevance, is resolved. The advantages of direct numerical simulation can be summarized in three points:  Extremely precise details and parameters of movement of turbulence at any point of the flow;  Instant results, which are not possible to get through experimentation, can be generated in direct numerical simulation;  Modelling of turbulent flows that are impossible to happen in reality, for better physical perception of this type of flow (eg walls without tensions). Despite these advantages, direct numerical simulation is almost exclusively applied to research and at low Reynolds numbers. The fact that all the parameters are thoroughly modelled results in a very time-consuming and costly numerical modelling (Versteeg. 2007). Feasibility test of the simulation of a running turbine with OpenFOAM 33 3.7 PRINCIPLES OF SOLUTION OF THE GOVERNING EQUATIONS The discretization in space is done to better calculate the governing equations of movement. This discretization is done by transforming in smaller and less complex parts in order to facilitate the calculation. The continuity law states that matter is a continuous medium with no voids in its interior. This principle does not consider the molecular structure, but only their macroscopic form, making it possible to transfer the continuous governing equations for homologous discrete values. Spatial discretization methods are a tool which enables to approximate those governing equations for a system of algebraic equations, the variables of the problem. These variables will be obtained at discrete locations in space and time. The development for the different methodologies has gained a major improvement in the last years and is still an area in study. As it was mentioned on the chapter 3.3, there are two types of meshes (or grids): structured and unstructured. The methods used for the special discretisation schemes rely on these two types of grids and are the following:  Finite difference method  Finite volume method  Finite element method But before explaining each one of these methods, it is important to first explain the extent to which the different grids influence the used method. The highest advantage of using a structured grid relies on the fact that these meshes follow a logical sense represented by the indices i,j,k. These indices represent a linear address space that can also be called computational space and they relate to the way the flow variables are stored in the computer memory. Fig. 16 – Structured mesh (or grid) in two dimensions – physical space. (Source: Blazek, 2001) Feasibility test of the simulation of a running turbine with OpenFOAM 34 Fig. 17 – Structured mesh (or grid) in two dimensions – computational space. (𝜉 and η represent a curvilinear coordinate system). (Source: Blazek, 2001) Having this characteristic, it allows a quick and easy access to all the neighbour points in a grid by the simple act of adding or subtracting an integer number to the correspondent index (e.g. (i+1), (k-3)) as it’s shown on Fig. 16 and Fig. 17. This increases the facility of making the evaluation of gradients, fluxes and also to make the treatment of the boundary conditions. It also results in a more efficient memory usage since the neighbouring elements of a given cell are known a priori by the structure nature of the mesh. But using structured grids has also disadvantages and the biggest one is that it is difficult to perform structured grids to complex geometries. There is the possibility to simplify this procedure by dividing the space in smaller parts resorting to the multiblocks approach. In any methodology of this type, it is obvious to expect that having more blocks leads to an increase of the complexity of the flow solver, because it is essential that a different logic needs to be used to represent the exchange of physical quantities or fluxes between the blocks. (Blazek 2001) On the other hand, there are also the unstructured grids that offer a large flexibility in the treatment of complex geometries. The main advantage of this type of meshes is that the triangular for 2D meshes or tetrahedral, for 3D meshes, can be generated automatically regardless of the complexity of the domain. Of course it is still necessary to set some restrictions correctly in order to achieve a good quality mesh. There is still a big difference regarding the amount of time taken to construct an unstructured mesh comparatively to a structured mesh, in which the first one is much quicker. The disadvantage is that it is required to implement a sophisticated data structure inside the flow solver and this data structures, depending on the computer hardware, takes to a more or less reduced computational efficiency. Following these requirements, these structured schemes need a higher memory necessity. After presenting these two points of view about the grids it is comprehensible that most CFD software implements on their solvers the unstructured flows. After the generation of the grid, it is now possible to return to the beginning of this sub-chapter and present the methods to proceed to the discretisation of the governing equations briefly explained in the flowing subsections. Feasibility test of the simulation of a running turbine with OpenFOAM 35 3.7.1. FINITE DIFFERENCE METHOD This is the oldest and simpler method for the resolution of differential equations. This method is applied to the differential form of the governing equations. The main principle by which this method is led by is to apply a Taylor series expansion for the discretisation of the derivatives of the flow variables. The main advantage of this methodology is its simplicity and the possibility to easily obtain high-order approximations. By obtaining these approximations, it is possible to consequently accomplish high-order accuracy of the special discretisation. On the other hand, the disadvantage is that this method requires the application of a structured grid and this leads to a decrease on the application field. When faced with curvilinear coordinates, the finite difference method cannot be applied. Instead, the governing equations need to be transformed into a Cartesian coordinate system which means that they are going to go from the physical space to the computational space. This can be easily perceived by analysing the Fig. 16 and Fig. 17. This method is, nowadays, commonly used for the direct numerical simulation of turbulence – DNS but because of its harder application to complex grids it’s not of common use in industrial cases, for instance. 3.7.2. FINITE VOLUME METHOD This method uses directly the conservation laws – the integral formulation of the Navier-Stokes equations. The main principle that this method supports itself is in the division of the physical space into a number of arbitrary polyhedral control volumes. The control volume is the element in the mesh that has become the domain of the solution. The continuity law is applied in each one of these elements and the division of the mesh is done exactly to simplify the calculation. The precision of the special discretisation depends on the particular scheme with which one of these fluxes is evaluated. There are two possibilities to define the position and shape of the control volume regarding the mesh and they are as follows:  Cell centred scheme – the flow quantities are stored at the centroids of the grid cells and the control volumes are identical to the grid cells. Fig. 18 – Cell centred scheme. (Source: Kolditz, 2001) Feasibility test of the simulation of a running turbine with OpenFOAM 36  Cell-vertex scheme – the flow variables are stored at the grid points. Now the control volume can be one of two options: or the union of all cells sharing the grid point or some volume centred on the grid point. Fig. 19 – Cell-vertex scheme. (Source: Kolditz, 2001) The main advantage of the finite volume method is that the problem with any type of transformation between coordinate systems does not exist here because the special discretisation is carried out directly in the physical space. In comparison with other methods, like the finite difference method for instance, one important advantage is that this one is very flexible meaning that it is easily applied in both types of grids, structured or unstructured grids. And this feature is what turns this method in such an important tool on resolving discretisation problems because it can be used for the treatment of flows in complex geometries. Another big advantage of this method is that it directly discretises the conservation laws and by that, mass, momentum and energy are also conserved along the process. This is the most used method in CFD software and it also is used in this MS thesis. 3.7.3. FINITE ELEMENT METHOD On its basic grounds, the objective of this method was to analyse structural problems. With some evolution and research, the finite element method started to be used for the numerical solution of field equations in continuous media. As the previous method, this one can also be used in unstructured grids and its integral formulation makes is easier to be used in flows or around complex geometries. This method has a lot of mathematical resemblances to the finite volume method, although the numerical effort in this case is higher to resolve the same type of problems. By this reason, the two can sometimes be combined to resolve the boundaries and the discretisation of the viscous fluxes specially when encountered with unstructured grids. Feasibility test of the simulation of a running turbine with OpenFOAM 37 4 Computational Fluid Dynamics – Procedure 4.1. DESCRIPTION OF THE MODEL The main objective of this work was to validate the possibility of building a numerical model of a running turbine in the software OpenFOAM. Most of the studies, so far, were based on probabilities of occurrence of fish injury and achieving the result with this software could be a support for validation of the existing models. Most of the numerical studies made so far were accomplished resorting to Propeller turbines like Kaplan and Francis turbines and ,by that fact, the simulation developed in this project will be made with a four blade propeller turbine. An exact model of a complete turbine with all its components would be very extensive for the time available to perform this work and by this it was chosen to generate the part of the turbine that has movement: its rotor and the respective blades. It was also defined a box where the turbine was inserted that afterwards would have an inflow. This model allows a better comprehension of the values of velocity and pressure while the turbine is in movement. 4.2. PROCEDURE As any CFD analysis, the major point is to obtain a deeper understanding of the problem under consideration and to have valid results, the need to make appropriate choices when it comes to the solver applications is also required. In every computational fluid dynamics, the basic steps are usually common, regardless the software in question. In this work, making use of the OpenFOAM software, the mandatory steps followed were: 1. Problem definition 2. Pre-processing and mesh generation 3. Choice of turbulence model 4. Selection of boundary conditions 5. Selection of the solver 6. Execution of the numeric study 7. Discussion and verification of the results Feasibility test of the simulation of a running turbine with OpenFOAM 38 After defining the geometry in analysis the same geometry is divided into infinitesimal elements that all together represent the mesh. Once this is defined, the turbulence model is chosen among the ones discussed in chapter 3 and adjusted to the flow properties and to the boundary conditions. After having all this clearly distinct, it is necessary to define the solver that is going to run the problem. This process in question is a typical case of trial and error. The correct mesh generation and refinement and achieving all the correct parameters in all the dictionaries requires an elevated amount of time, especially for beginners in modelling software such as OpenFOAM. The computer used to run the program was an Intel Core 2 Quad Q9550 @ 2.83GHz with 4 Gb of RAM and the operating system used was Linux. This information is useful in order to better understand the required time to run every command including all the missed trials. 4.3. SOFTWARE PROCEDURE – OPENFOAM The software used in this work was OpenFOAM which means Open Source Field Operation and Manipulation, the version 2.2.0. OpenFOAM works with C++ libraries that are used to create executables that are named applications. These applications are divided into two categories: solvers and utilities. Solvers have the function to resolve a specific problem in continuum mechanics and utilities solve tasks that involve data manipulation. One of the most valuable features in OpenFOAM is that users can add new solvers and utilities as long as they have the ability of programming and understanding the underlying method. The main structure of OpenFOAM is shown in Fig. 20: Fig. 20 – Overview of OpenFOAM structure. (Source: OpenFOAM. The Open Source CFD Toolbox. User Guide. 2013) All the interfaces of pre and post processing are all utilities from OpenFOAM and this guaranties a consistent data handling. OpenFOAM uses the information kept in three main folders:  0  Constant  System Feasibility test of the simulation of a running turbine with OpenFOAM 39 On the 0 folder, that is a time folder, it is saved all the initial conditions defined by the user. On the constant folder, the information about the mesh generation as well as the physical properties of the flow are saved. All the information regarding the solvers is stored on the system folder. The main structure of this work folder was the following:  0  alpha1  k  nut  omega  p  U  constant  polyMesh  blockMeshDict  boundary  faces  neighbour  owner  points  trisurface  original.stl  dynamicMeshDict  g  mechanicalProperties  RASProperties  transportProperties  turbulenceProperties  system  controlDict  createPatchDict  decomposeParDict  fvSchemes  fvSolution  setFieldsDict  snappyHexMeshDict After running all the commands, more time folders will be created in the main folder. In these dictionaries there is all the necessary information to the development of the mesh and the flow. To create the simple meshes, it’s possible to introduce all the coordinates of the points in the dictionary blockMeshDict, but when it comes to meshes of higher complexity, such as a turbine, OpenFOAM provides a mesh generator called SnappyHexMesh. To use this generator, it is necessary to introduce its respective dictionary in the system folder and also a stereolithography file (.stl) in the folder triSurface in the constant folder. The available platform to visualize the results is called ParaView and it is accessible by running paraFoam in the terminal. Feasibility test of the simulation of a running turbine with OpenFOAM 46 In the library defined as libOscillatingVelocityPointPatchVectorField.C it’s possible to change the function that is already defined to move the patches. The one used in this work followed this mechanism: Fig. 26 – Mechanism of Oscillating Velocity. (source: Eslamdoost 2010) For each time step, it’s calculated a new angle for the oscillation arm after the user has defined the angular oscillation amplitude and frequency. This angle is calculated by: 𝑎𝑛𝑔𝑙𝑒=𝑎𝑛𝑔𝑙𝑒0+𝑎𝑚𝑝𝑙𝑖𝑡𝑢𝑑𝑒×sin (𝜔.𝑡) (4.1) As it’s possible to see in the Fig. 26, angle 0 is the initial angle of the oscillation arm and amplitude is the amplitude of the angular oscillation. Amplitude works as a weight factor, controlling the new angle in each time step. The frequency of the angular oscillation is represented by 𝜔 and t is the actual run time of the simulation. To use this library there needs to be a link on the controlDict dictionary that leads the program to use this application. On the controlDict dictionary, it was added the keyword libs (“libMyOscillatingVelocityPointPatchVectorField.so”). After altering the controlDict dictionary it is also important to add new keywords to the dictionaries fvsolution and fvschemes. In these two dictionaries, it needs to be specified which type of solver the program should use to simulate the movement. OpenFOAM simulates the movement resorting to the following solvers: PCG/PBiCG, smoothSolver and GAMG. The chosen one to solve the cellMotionU was the GAMG solver. This is a solver to generalised geometric-algebraic multi-grid and is much more inclusive that the others. It was also necessary to choose the respective preconditioner. The one used was DICGaussSeidel. Another important step before running the simulation is defining a dynamicMeshDict dictionary stored on the constant folder. Once again, there are two parameters of high importance: the solver and the diffusivity scheme. The solvers can be chosen from this three:  displacementLaplacian;  velocityLaplacian;  SBRStress. Feasibility test of the simulation of a running turbine with OpenFOAM 47 The chosen one was velocityLaplacian. The diffusivity models can be dived into two categories: quality-based methods and distance-based methods. Regarding the quality-based methods there is the possibility to choose from:  uniform;  directional;  motionDirectional;  inverseDistance. In the group of distance-based methods it’s possible to choose from:  linear;  quadratic;  exponencial. The selected one was the uniform. After defining each one of this points it is also necessary to define both pointMotionU and cellMotionU dictionaries. They are both stored in the 0 folder. The parameters defined in these two folders are the following:  axis – represents the axis of rotation, defined in meters.  origin – center of the rotation, defined in meters.  angle0 – angle which the oscillation occurs, defined in rad.  Amplitude – amplitude of the angular oscillation, defined in rad.  Omega – angular oscillation frequency, defined in rad/s. After running and adjusting several times this problem, the solution was never properly obtained. The solver couldn’t find convergence and that can be due to the mesh but more precisely to the movement. This was not a good strategy to impose movement to the turbine and another solution needed to be found. 4.8.2. SECOND APPROACH RESORTING TO ANGULARVELOCITY Although the first solution was not proper to the problem in question, the reasoning followed was not completely wrong. Like in every other open source software, the user has the possibility to change the code where the functions are defined. From the list of available boundary types, the one explained in the previous topic is an adequate starting point. After finding the directory for the angularOscillatingVelocity this was copied to the folder where the study was being made and it was renamed. In the .C and .H files, any mention of the word Oscillating was erased not only from the names of the files, but also inside them. The next point was defining the better solution to the problem in study. Inside the .C file it was defined the following equation: 𝑠𝑐𝑎𝑙𝑎𝑟 𝑎𝑛𝑔𝑙𝑒=𝑎𝑛𝑔𝑙𝑒0_+𝑎𝑚𝑝𝑙𝑖𝑡𝑢𝑑e_∗sin (omega_∗𝑡.value( )); (4.2) It is possible to infer from here how OpenFOAM is coded. Equations can be written in a much similar way from math notation. Feasibility test of the simulation of a running turbine with OpenFOAM 48 This is where modification should be done in order to obtain the desired movement. The new equation is: 𝑠𝑐𝑎𝑙𝑎𝑟 𝑎𝑛𝑔𝑙𝑒=𝑎𝑛𝑔𝑙𝑒0_+omega_∗𝑡.value( ); (4.3) The amplitude and the sin function were eliminated. In the two files .C and .H all the lines of code where the word amplitude appears were eliminated as well. After this was defined, the compilation of the parameters was determined. The files file inside the Make directory changed into: angularVelocity/angularVelocityPointPatchVectorField.C LIB = $(FOAM_USER_LIBBIN)/libAngularVelocityPointPatchVectorField The first line leads to the location of the .C file with the source code and the second to the directory where user libraries for OpenFOAM are stored followed by the name of the library that was created. Once all the files were correctly defined, it was necessary to compile the library by running: wmake libso After this new library was created, there is the need to define in the controlDic dictionary the path where the program needs to find the new code to solve the problem. So, in this dictionary, it is added: Libs (“libAngularVelocityPointPatchVectorField.so” “libOpenFOAM.so”) After this library is defined, the boundaries in the dictionaries pointMotionU needs to be corrected according to the new library created. The important part of this library is the part where the boundary of the turbine is defined and it’s the following: turbine_patch28539 { type libAngularVelocity; axis (1 0 0); origin (0 0 0); angle0 0; omega 20.94; //200rpm value uniform (1 0 0); } A propeller turbine can go from 79 rpm to 429 rpm. Since this is a numerical study, an intermediate value was chosen: 200rpm that correspond to 20.94 rad/s. Before running the problem, it is important to verify if the mesh is well constructed and so it’s necessary to run checkMesh. The log file that is acquired when running this is in the appendix E. Due to the high computational requisitions of the software, it was used the decomposePar tool that divides the main processor into 8 smaller processors and by that it is possible to achieve quicker results. Nevertheless, attention must be taken when doing this because a lot of boundary conditions may not be correctly handled and the user needs to correct them manually. Feasibility test of the simulation of a running turbine with OpenFOAM 49 Since the flow simulation requires a small time step, making the calculations longer in time it is possible to simply run moveMesh. This command will only calculate the movement of the mesh without taking into account the flow variables. After reconstructing the processors it is possible to visualize the results from Fig. 27 to Fig. 32. Fig. 27 - Time step 1 corresponding to 0.01s. (image generated by ParaView) Fig. 28 - Time step 15 corresponding to 0.15s. (image generated by ParaView) Feasibility test of the simulation of a running turbine with OpenFOAM 50 Fig. 29 - Time step 100 corresponding to 1s. (image generated by ParaView) Fig. 30 - Time step 200 corresponding to 2s. (image generated by ParaView) Feasibility test of the simulation of a running turbine with OpenFOAM 51 Fig. 31 - Time step 250 corresponding to 2.5s. (image generated by ParaView) Fig. 32 – Time step 300 corresponding to 3s. (image generated by ParaView) Feasibility test of the simulation of a running turbine with OpenFOAM 52 4.9. SOLUTION AND ALGORITHM CONTROL All the equation solvers, tolerances and other algorithm controls that the program has to follow to achieve a result are stored in the dictionary fvsolution stored in the system folder. In this dictionary there are some sub-dictionaries that include solvers, relaxationFactors, PISO and SIMPLE. The first sub-dictionary in fvSolution is the solver and the user needs to specify what the linear-solver for each discretised equation is. The main idea is that the main solver chosen to simulate the problem – in this case pimpleDyMFoam – solves each one of the equations specified here, like the velocity U and the pressure p. The dictionary for each one of these keywords also has the type of solver and its parameters. This last one includes tolerance, reTol, smoother, cacheAgglomeration, nCellsInCoarsestLevel, agglomerator and mergeLevels. The chosen solver for pressure (p) and the cell movement (cellMotionU) was GAMG – generalided geometric-algebraic multi-grid. Although in an early stage of the development of the work it was selected the PBiCG – Preconditioned bi-conjugate gradient for asymmetric matrices, this wasn’t allowing the simulation to even run till the end. This was caused by the fact that the solver could not make the correct convergence of the flow. Therefore, it was made the decision of changing the solvers of pressure and cell movement to GAMG and continuing to use the solver PBiCG in velocity and turbulent kinetic energy. Changing to the GAMG solver allowed that the simulation could run slightly faster since this solver acts by generating a solution resorting to only a small number of cells. The second step is mapping the solution onto a finer mesh assuming the results already obtained before. After choosing the solver, it is necessary to select the correspondent preconditioned conjugate gradient solver. The ones that OpenFOAM recognizes are:  DIC, Diagonal incomplete-Cholesky;  FDIC, faster diagonal incomplete-Cholesky;  DILU, diagonal incomplete-LU;  Diagonal;  GAMG;  None, when the choice is to have no preconditioning. The preconditioner chosen was the DILU for this work. Another keyword that needs definition is the smooth solvers that are:  GaussSeidel;  DIC;  DICGaussSeidel. GaussSeidel is the most reliable one and for that matter the one used with the GAMG solver, since this one requires the specification of a smoother. The tolerance defined started with the value of 1e-06 but due to the elevated time needed to run the simulation it was reformed to 1e-03. Finally, when it comes to algorithms OpenFOAM, it gives the possibility to used one of these:  PISO – Pressure implicit with splitting of operator;  SIMPLE – Semi implicit method for pressure linked equations;  PIMPLE – Pressure implicit method for pressure linked equations. Feasibility test of the simulation of a running turbine with OpenFOAM 53 Both PISO and SIMPLE are algorithms that are iterative procedures for solving equations of velocity and pressure. The first one is used to transient problems and the second for steady-state simulations. The common procedure is that both evaluate the initial conditions and then they corrected them. The difference is that SIMPLE only makes one correction and PISO executes more than one and usually less than four. The number of corrections can be defined in the fvsolution dictionary. The PIMPLE is a fusion of PISO and SIMPLE and it has two additional possibilities: to correct the loops and to slow down the iterations regarding exterior variables. The convergence is much more effective with this solver and by that matter was the chosen one. In a first stage, it was used the pimpleFoam to better authenticate all the choices made regarding the fvsolutions. At a later point, the solver was changed to pimpleDyMFoam that is a transient solver for incompressible flows of Newtonian fluids on a moving mesh. 4.10. DATA CONTROL The OpenFOAM solvers create some database solvers throughout the running process where the data generated is output in these database folders. The time conditions are all available for modification in the controlDict dictionary. The endTime was defined for 3 seconds and although this is a small number that can lead to not so good conclusions, the increase of this value leads to simulations in time not suitable for the extension of this work. The deltaT was set to 0.01 that supposedly would create folders for each one of this time steps, but that didn’t happen due to the high complexity of the mesh. The time-steps were being adapted during the simulation process automatically by the software. Another factor defined here is the Courant Number that can be calculated by: Where δt is the time step, |U| is the magnitude of the velocity through that cell and δx is the cell size in the direction of the velocity. Since the flow velocity varies across the domain, it’s important to ensure a Co < 1. Due to the complexity of the mesh, the chosen value was 0.7 for the Courant number. 𝐶𝑜=𝛿𝑡|𝑈| 𝛿𝑥 (4.2) Feasibility test of the simulation of a running turbine with OpenFOAM 54 Feasibility test of the simulation of a running turbine with OpenFOAM 55 5 CONCLUSION AND FUTURE WORKS 5.1. RESULTS OF THE STUDY The procedure done is explained in the previous chapter and it’s qualitatively summarized in the following paragraphs. The geometry of the casing where the turbine would be inserted was generated by resorting to a dictionary responsible to create a separate list of points for each boundary of the casing. When it comes to the boundaries, they can be easily changed and be adapted to any problem the user decides to solve. Working with a command line interface and not with a graphical user interface, hinders the process and makes it time consuming. The dictionary snappyHexMesh was written in order to obtain a mesh generated from a .stl file. The parameters have been adjusted several times in order to obtain not only a correct refinement but also a high mesh quality. A balance must be found, particularly when it comes to the refinement values. After the refinement, extrusion and snapping it was possible to visualize the generated mesh. While preparing the mesh, the user has to consider a lot of parameters in different dictionaries, which influence the performance of the whole model. Sometimes, these interdependencies are not apparent in the beginning and that can lead to performance problems during the simulation. After having tried several different approaches to improve the movement on the mesh, a new script, found in the Appendix D.1, has been written for this purpose. A new boundary condition and consequently a new library had to be implemented in OpenFOAM. The mesh was divided into eight parts that can be computed independently on their own processors in order to speed up the simulation. When running moveMesh, it is possible to correctly visualize the movement of the turbine without any calculations of the flow. The computing of the movement is the main objective of this work. The possibility of changing the number of rotations per minute and its consequent velocity is left to the user's discretion and, by that, this model can be used with different turbines and different running schemes. After the flow boundaries were correctly defined, pimpleDyMFoam was run in the terminal. This is where the problem was not well succeeded. The program could not solve the equations for the pressure when the mesh was in movement. To better understand the extension of the error, pimpleFoam was run in the terminal and in a first trial was not successfully solved. The flow model was then changed to laminar instead of turbulent and successfully solved. The problem could derive from the fact that the constant mesh deformation and the resolution of the equations to solve the pressure and velocity are not compatible together. This is still field of interest and development. Feasibility testo f the simulation of a running turbine with OpenFOAM 62 Wendt, J. E. (2008). Computational fluid dynamics: an introduction. Springer Science & Business Media. Williams, U.P., Scruton, D.A., Goosney, R.F., Bourgeois, C.E., Orr, D.C. and Ruggles, C.P. (1993). Proceedings of the Workshop on Fish Passage at Hydroelectric Developments. St. John’s, Newfoundland. Yoganathan, A. P., He, Z., & Casey Jones, S. (2004). Fluid mechanics of heart valves. Annu. Rev. Biomed. Eng., 6, 331-362. Websites: [1] http://www.daviddarling.info/encyclopedia/R/AE_reaction_turbine.html. April, 2015. [2] http://www.learnengineering.org/2013/08/kaplan-turbine-hodroelectric-power-gneration.html. April, 2015 [3] http://www.learnengineering.org/2014/01/how-does-francis-turbine-work.html. April, 2015 [4] http://www.fao.org/docrep/004/y2785e/y2785e03.htm. Março 2015. March, 2015. [5] http://www.ebah.pt/content/ABAAAATzIAE/tipos-turbinas-hidraulicas?part=2. May 2015 [6] http://www.bakker.org/. May 2015 [7] http://www.cfd-online.com/Wiki/Mesh_classification. May 2015 [8] http://www.pointwise.com/theconnector/March-2013/Structured-Grids-in-Pointwise.shtml. April 2015 [9] http://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture1.pdf. April 2015 [10] http://noc.ac.uk/f/content/science-technology/workshops/presentation_2_meshing.pdf. April, 2015. [11] http://www.renewablesfirst.co.uk/hydro-learning-centre/kaplan-turbines/ March, 2015. 63 Appendix A – Pre-processing A.1. Standart header of OpenFoam /*--------------------------------*- C++ -*--------------------------- -------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.1.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*-------------------------------------------------------------------- -------*/ FoamFile { version 2.0; format ascii; class dictionary; object blockMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // A.2 constant/polyMesh/blockMeshDict convertToMeters 1; vertices ( (-15 -15 15) (15 -15 15) (15 -15 -15) (-15 -15 -15) (-15 15 15) (15 15 15) (15 15 -15) (-15 15 -15) ); blocks ( hex (0 1 2 3 4 5 6 7) (100 80 80) simpleGrading (1 1 1) ); edges ( ); boundary ( inlet { type patch; faces ( (0 3 7 4) ); } outlet { type patch; faces ( (5 6 2 1) ); } down { type wall; faces ( (0 1 2 3) ); } up { type wall; faces ( (4 7 6 5) ); } front { type wall; faces ( (0 4 5 1) ); } back { type wall; faces ( (3 2 6 7) ); } ); mergePatchPairs ( ); A.3 system/snappyHexMeshDict // Which of the steps to run castellatedMesh true; snap true; addLayers true; // Geometry. Definition of all surfaces. All surfaces are of class // searchableSurface. // Surfaces are used // - to specify refinement for any mesh cell intersecting it // - to specify refinement for any mesh cell inside/outside/near // - to 'snap' the mesh boundary to the surface geometry { original.stl { type triSurfaceMesh; name turbine; } /* Verfeinerung { type searchableBox; min (-1 -10 -10); max (6 10 10); }*/ }; // Settings for the castellatedMesh generation. castellatedMeshControls { // Refinement parameters // ~~~~~~~~~~~~~~~~~~~~~ // If local number of cells is >= maxLocalCells on any processor // switches from from refinement followed by balancing // (current method) to (weighted) balancing before refinement. maxLocalCells 500000; // Overall cell limit (approximately). Refinement will stop immediately // upon reaching this number so a refinement level might not complete. // Note that this is the number of cells before removing the part which // is not 'visible' from the keepPoint. The final number of cells might // actually be a lot less. maxGlobalCells 2000000; // The surface refinement loop might spend lots of iterations refining just a // few cells. This setting will cause refinement to stop if <= minimumRefine // are selected for refinement. Note: it will at least do one iteration // (unless the number of cells to refine is 0) minRefinementCells 0; maxLoadUnbalance 0.10; // Number of buffer layers between different levels. // 1 means normal 2:1 refinement restriction, larger means slower // refinement. nCellsBetweenLevels 2; // Explicit feature edge refinement // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Specifies a level for any cell intersected by its edges. // This is a featureEdgeMesh, read from constant/triSurface for now. features ( /* { file "someLine.eMesh"; //level 2; levels ((0.0 2) (1.0 2)); }*/ ); // Surface based refinement // ~~~~~~~~~~~~~~~~~~~~~~~~ // Specifies two levels for every surface. The first is the minimum level, // every cell intersecting a surface gets refined up to the minimum level. // The second level is the maximum level. Cells that 'see' multiple // intersections where the intersections make an // angle > resolveFeatureAngle get refined up to the maximum level. refinementSurfaces { original.stl { level (1 2); } } resolveFeatureAngle 30; // Region-wise refinement // ~~~~~~~~~~~~~~~~~~~~~~ // Specifies refinement level for cells in relation to a surface. One of // three modes // - distance. 'levels' specifies per distance to the surface the // wanted refinement level. The distances need to be specified in // descending order. // - inside. 'levels' is only one entry and only the level is used. All // cells inside the surface get refined up to the level. The surface // needs to be closed for this to be possible. // - outside. Same but cells outside. refinementRegions { /*Verfeinerung { mode inside; levels ((1 4)); }*/ } // Mesh selection // ~~~~~~~~~~~~~~ // After refinement patches get added for all refinementSurfaces and // all cells intersecting the surfaces get put into these patches. The // section reachable from the locationInMesh is kept. // NOTE: This point should never be on a face, always inside a cell, even // after refinement. // This is an outside point locationInMesh (-0.033 -0.033 0.0033); locationInMesh (-10 -10 10); // Inside point // Whether any faceZones (as specified in the refinementSurfaces) // are only on the boundary of corresponding cellZones or also allow // free-standing zone faces. Not used if there are no faceZones. allowFreeStandingZoneFaces true; } // Settings for the snapping. snapControls { //- Number of patch smoothing iterations before finding correspondence // to surface nSmoothPatch 10; //- Relative distance for points to be attracted by surface feature point // or edge. True distance is this factor times local // maximum edge length. tolerance 1.2; //- Number of mesh displacement relaxation iterations. nSolveIter 10; //- Maximum number of snapping relaxation iterations. Should stop // before upon reaching a correct mesh. nRelaxIter 3; //- Highly experimental and wip: number of feature edge snapping // iterations. Leave out altogether to disable. nFeatureSnapIter 10; } // Settings for the layer addition. addLayersControls { // Are the thickness parameters below relative to the undistorted // size of the refined cell outside layer (true) or absolute sizes (false). relativeSizes true; // Per final patch (so not geometry!) the layer information layers { "turbine.*" { nSurfaceLayers 5; } } // Expansion factor for layer mesh expansionRatio 1.5; //- Wanted thickness of final added cell layer. If multiple layers // is the thickness of the layer furthest away from the wall. // See relativeSizes parameter. finalLayerThickness 0.8; //- Minimum thickness of cell layer. If for any reason layer // cannot be above minThickness do not add layer. // See relativeSizes parameter. minThickness 0.001; //- If points get not extruded do nGrow layers of connected faces that are // also not grown. This helps convergence of the layer addition process // close to features. nGrow 0; // Advanced settings //- When not to extrude surface. 0 is flat surface, 90 is when two faces // make straight angle. featureAngle 30; //- Maximum number of snapping relaxation iterations. Should stop // before upon reaching a correct mesh. nRelaxIter 5; // Number of smoothing iterations of surface normals nSmoothSurfaceNormals 1; // Number of smoothing iterations of interior mesh movement direction nSmoothNormals 3; // Smooth layer thickness over surface patches nSmoothThickness 10; // Stop layer growth on highly warped cells maxFaceThicknessRatio 0.5; // Reduce layer growth where ratio thickness to medial // distance is large maxThicknessToMedialRatio 0.4; // Angle used to pick up medial axis points minMedianAxisAngle 130; // Create buffer region for new layer terminations nBufferCellsNoExtrude 0; // Overall max number of layer addition iterations. The mesher will exit // if it reaches this number of iterations; possibly with an illegal // mesh. nLayerIter 50; // Max number of iterations after which relaxed meshQuality controls // get used. Up to nRelaxIter it uses the settings in meshQualityControls, // after nRelaxIter it uses the values in meshQualityControls::relaxed. nRelaxedIter 20; } // Generic mesh quality settings. At any undoable phase these determine // where to undo. meshQualityControls { //- Maximum non-orthogonality allowed. Set to 180 to disable. maxNonOrtho 60; //- Max skewness allowed. Set to <0 to disable. maxBoundarySkewness 20; maxInternalSkewness 4; //- Max concaveness allowed. Is angle (in degrees) below which concavity // is allowed. 0 is straight face, <0 would be convex face. // Set to 180 to disable. maxConcave 80; //- Minimum pyramid volume. Is absolute volume of cell pyramid. // Set to a sensible fraction of the smallest cell volume expected. // Set to very negative number (e.g. -1E30) to disable. minVol 1e-13; //- Minimum quality of the tet formed by the face-centre // and variable base point minimum decomposition triangles and // the cell centre. Set to very negative number (e.g. -1E30) to // disable. // <0 = inside out tet, // 0 = flat tet // 1 = regular tet minTetQuality 1e-30; //- Minimum face area. Set to <0 to disable. minArea -1; //- Minimum face twist. Set to <-1 to disable. dot product of face normal //- and face centre triangles normal minTwist 0.05; //- minimum normalised cell determinant //- 1 = hex, <= 0 = folded or flattened illegal cell minDeterminant 0.001; //- minFaceWeight (0 -> 0.5) minFaceWeight 0.05; //- minVolRatio (0 -> 1) minVolRatio 0.01; //must be >0 for Fluent compatibility minTriangleTwist -1; //- if >0 : preserve single cells with all points on the surface if the type processor; } procBoundary0to2 { type processor; } procBoundary0to3 { type processor; } procBoundary0to4 { type processor; } procBoundary0to5 { type processor; } procBoundary1to0 { type processor; } procBoundary1to2 { type processor; } procBoundary1to3 { type processor; } procBoundary1to4 { type processor; } procBoundary1to5 { type processor; } procBoundary2to0 { type processor; } procBoundary2to1 { type processor; } procBoundary2to3 { type processor; } procBoundary2to6 { type processor; } procBoundary2to7 { type processor; } procBoundary3to0 { type processor; } procBoundary3to1 { type processor; } procBoundary3to2 { type processor; } procBoundary3to6 { type processor; } procBoundary3to7 { type processor; } procBoundary4to0 { type processor; } procBoundary4to1 { type processor; } procBoundary4to5 { type processor; } procBoundary4to6 { type processor; } procBoundary4to7 { type processor; } procBoundary5to0 { type processor; } procBoundary5to1 { type processor; } procBoundary5to4 { type processor; } procBoundary5to6 { type processor; } procBoundary5to7 { type processor; } procBoundary6to2 { type processor; } procBoundary6to3 { type processor; } procBoundary6to4 { type processor; } procBoundary6to5 { type processor; } procBoundary6to7 { type processor; } procBoundary7to2 { type processor; } procBoundary7to3 { type processor; } procBoundary7to4 { type processor; } procBoundary7to5 { type processor; } procBoundary7to6 { type processor; } } A.6 0/k dimensions [ 0 2 -2 0 0 0 0 ]; internalField uniform 1; boundaryField { inlet { type turbulentIntensityKineticEnergyInlet; intensity 0.05; // 5% turbulent intensity value uniform 1; } outlet { type inletOutlet; inletValue uniform 1; } down { type kqRWallFunction; value uniform 0; } up { type kqRWallFunction; value uniform 0; } front { type kqRWallFunction; value uniform 0; } back { type kqRWallFunction; value uniform 0; } turbine_patch28539 { type kqRWallFunction; value uniform 0; } procBoundary0to1 { type processor; } procBoundary0to2 { type processor; } procBoundary0to3 { type processor; } procBoundary0to4 { type processor; } procBoundary0to5 { type processor; } procBoundary1to0 { type processor; } procBoundary1to2 { type processor; } procBoundary1to3 { type processor; } procBoundary1to4 { type processor; } procBoundary1to5 { type processor; } procBoundary2to0 { type processor; } procBoundary2to1 { type processor; } procBoundary2to3 { type processor; } procBoundary2to6 { type processor; } procBoundary2to7 { type processor; } procBoundary3to0 { type processor; } procBoundary3to1 { type processor; } procBoundary3to2 { type processor; } procBoundary3to6 { type processor; } procBoundary3to7 { type processor; } procBoundary4to0 { type processor; } procBoundary4to1 { type processor; } procBoundary4to5 { type processor; } procBoundary4to6 { type processor; } procBoundary4to7 { type processor; } procBoundary5to0 { type processor; } procBoundary5to1 { type processor; } procBoundary5to4 { type processor; } procBoundary5to6 { type processor; } procBoundary5to7 { type processor; } procBoundary6to2 { type processor; } procBoundary6to3 { type processor; } procBoundary6to4 { type processor; } procBoundary6to5 { type processor; } procBoundary6to7 { type processor; } procBoundary7to2 { type processor; } procBoundary7to3 { type processor; } procBoundary7to4 { type processor; } procBoundary7to5 { type processor; } procBoundary7to6 { type processor; } } A.7 0/omega dimensions [0 0 -1 0 0 0 0]; internalField uniform 0.001; boundaryField { inlet { type turbulentMixingLengthFrequencyInlet; mixingLength 0.01; // 1cm - half channel height k k; value uniform 1; } outlet { type inletOutlet; inletValue uniform 1; } down { type omegaWallFunction; value uniform 0; } up { type omegaWallFunction; value uniform 0; } front { type omegaWallFunction; value uniform 0; } back { type omegaWallFunction; value uniform 0; } turbine_patch28539 { type omegaWallFunction; value uniform 0; } procBoundary0to1 { type processor; } procBoundary0to2 { type processor; } procBoundary0to3 { type processor; } procBoundary0to4 { type processor; } procBoundary0to5 { type processor; } procBoundary1to0 { type processor; } procBoundary1to2 { type processor; } procBoundary1to3 { type processor; } procBoundary1to4 { type processor; } procBoundary1to5 { type processor; } procBoundary2to0 { type processor; } procBoundary2to1 { type processor; } procBoundary2to3 { type processor; } procBoundary2to6 { type processor; } procBoundary2to7 { type processor; } procBoundary3to6 { type processor; } procBoundary3to7 { type processor; } procBoundary4to0 { type processor; } procBoundary4to1 { type processor; } procBoundary4to5 { type processor; } procBoundary4to6 { type processor; } procBoundary4to7 { type processor; } procBoundary5to0 { type processor; } procBoundary5to1 { type processor; } procBoundary5to4 { type processor; } procBoundary5to7 { type processor; } procBoundary5to6 { type processor; } procBoundary6to2 { type processor; } procBoundary6to3 { type processor; } procBoundary6to4 { type processor; } procBoundary6to5 { type processor; } procBoundary6to7 { type processor; } procBoundary7to2 { type processor; } procBoundary7to3 { type processor; } procBoundary7to4 { type processor; } procBoundary7to5 { type processor; } procBoundary7to6 { type processor; } } B.3 0/cellMotionU dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { default { type empty; } inlet { type fixedValue; value uniform (0 0 0); } outlet { type fixedValue; value uniform (0 0 0); } down { type fixedValue; value uniform (0 0 0); } up { type fixedValue; value uniform (0 0 0); } front { type fixedValue; value uniform (0 0 0); } back { type fixedValue; value uniform (0 0 0); } turbine_patch28539 { type rotatingWallVelocity; axis (1 0 0); origin (0 0 0); angle0 0; omega 20.94 value uniform (1 0 0); } procBoundary0to1 { type processor; } procBoundary0to2 { type processor; } procBoundary0to3 { type processor; } procBoundary0to4 { type processor; } procBoundary0to5 { type processor; } procBoundary1to0 { type processor; } procBoundary1to2 { type processor; } procBoundary1to3 { type processor; } procBoundary1to4 { type processor; } procBoundary1to5 { type processor; } procBoundary2to0 { type processor; } procBoundary2to1 { type processor; } procBoundary2to3 { type processor; } procBoundary2to6 { type processor; } procBoundary2to7 { type processor; } procBoundary3to0 { type processor; } procBoundary3to1 { type processor; } procBoundary3to2 { type processor; } procBoundary3to6 { type processor; } procBoundary3to7 { type processor; } procBoundary4to0 { type processor; } procBoundary4to1 { type processor; } procBoundary4to5 { type processor; } procBoundary4to6 { type processor; } procBoundary4to7 { type processor; } procBoundary5to0 { type processor; } procBoundary5to1 { type processor; } procBoundary5to4 { type processor; } procBoundary5to7 { type processor; } procBoundary6to2 { type processor; } procBoundary6to3 { type processor; } procBoundary6to4 { type processor; } procBoundary6to7 { type processor; } procBoundary7to2 { type processor; } procBoundary7to3 { type processor; } procBoundary7to4 { type processor; } procBoundary7to5 { type processor; } procBoundary7to6 { type processor; C – Solving C.1 constant/turbulenceProperties FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object turbulenceProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // simulationType RASModel; C.2 constant/RASProperties FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object RASProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // RASModel kOmegaSST; turbulence on; printCoeffs on; C.3 constant/transportProperties FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object transportProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // transportModel Newtonian; nu nu [ 0 2 -1 0 0 0 0 ] 1e-06; rho rho [ 1 -3 0 0 0 0 0 ] 1000; CrossPowerLawCoeffs { nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; m m [ 0 0 1 0 0 0 0 ] 1; n n [ 0 0 0 0 0 0 0 ] 0; } BirdCarreauCoeffs { nu0 nu0 [ 0 2 -1 0 0 0 0 ] 0.0142515; nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; k k [ 0 0 1 0 0 0 0 ] 99.6; n n [ 0 0 0 0 0 0 0 ] 0.1003; } C.4 system/controlDict FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application pimpleDyMFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 3; deltaT 0.01; writeControl adjustableRunTime; writeInterval 1e-2; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable no; adjustTimeStep no; maxCo 0.7; libs ( "libAngularVelocityPointPatchVectorField.so" "libOpenFOAM.so" ); } C.4 system/decomposeParDict FoamFile { version 2.0; format ascii; class dictionary; location "system"; object decomposeParDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // numberOfSubdomains 8; method simple; simpleCoeffs { n (2 2 2); delta 0.001; } hierarchicalCoeffs { n (1 1 1); delta 0.001; order xyz; } manualCoeffs { dataFile ""; } distributed no; roots ( ); C.5 system/createPatchDict FoamFile { version 2.0; format ascii; class dictionary; object createPatchDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // This application/dictionary controls: // - optional: create new patches from boundary faces (either given as // a set of patches or as a faceSet) // - always: order faces on coupled patches such that they are opposite. This // is done for all coupled faces, not just for any patches created. // - optional: synchronise points on coupled patches. // - always: remove zero-sized (non-coupled) patches (that were not added) // 1. Create cyclic: // - specify where the faces should come from // - specify the type of cyclic. If a rotational specify the rotationAxis // and centre to make matching easier // - always create both halves in one invocation with correct 'neighbourPatch' // setting. // - optionally pointSync true to guarantee points to line up. // 2. Correct incorrect cyclic: // This will usually fail upon loading: // "face 0 area does not match neighbour 2 by 0.0100005%" // " -- possible face ordering problem." // - in polyMesh/boundary file: // - loosen matchTolerance of all cyclics to get case to load // - or change patch type from 'cyclic' to 'patch' // and regenerate cyclic as above // Do a synchronisation of coupled points after creation of any patches. // Note: this does not work with points that are on multiple coupled patches // with transformations (i.e. cyclics). pointSync false; // Patches to create. patches ( /* { // Name of new patch name turbine_patch; // Dictionary to construct new patch from patchInfo p0_ = p.localPoints(); } } libAngularVelocityPointPatchVectorField:: libAngularVelocityPointPatchVectorField ( const libAngularVelocityPointPatchVectorField& ptf, const pointPatch& p, const DimensionedField<vector, pointMesh>& iF, const pointPatchFieldMapper& mapper ) : fixedValuePointPatchField<vector>(ptf, p, iF, mapper), axis_(ptf.axis_), origin_(ptf.origin_), angle0_(ptf.angle0_), omega_(ptf.omega_), p0_(ptf.p0_) {} libAngularVelocityPointPatchVectorField:: libAngularVelocityPointPatchVectorField ( const libAngularVelocityPointPatchVectorField& ptf, const DimensionedField<vector, pointMesh>& iF ) : fixedValuePointPatchField<vector>(ptf, iF), axis_(ptf.axis_), origin_(ptf.origin_), angle0_(ptf.angle0_), omega_(ptf.omega_), p0_(ptf.p0_) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // void libAngularVelocityPointPatchVectorField::updateCoeffs() { if (this->updated()) { return; } const polyMesh& mesh = this->dimensionedInternalField().mesh()(); const Time& t = mesh.time(); const pointPatch& p = this->patch(); scalar angle = angle0_ + omega_*t.value(); vector axisHat = axis_/mag(axis_); vectorField p0Rel = p0_ - origin_; vectorField::operator= ( ( p0_ + p0Rel*(cos(angle) - 1) + (axisHat ^ p0Rel*sin(angle)) + (axisHat & p0Rel)*(1 - cos(angle))*axisHat - p.localPoints() )/t.deltaT().value() ); fixedValuePointPatchField<vector>::updateCoeffs(); } void libAngularVelocityPointPatchVectorField::write ( Ostream& os ) const { pointPatchField<vector>::write(os); os.writeKeyword("axis") << axis_ << token::END_STATEMENT << nl; os.writeKeyword("origin") << origin_ << token::END_STATEMENT << nl; os.writeKeyword("angle0") << angle0_ << token::END_STATEMENT << nl; os.writeKeyword("omega") << omega_ << token::END_STATEMENT << nl; p0_.writeEntry("p0", os); writeEntry("value", os); } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // makePointPatchTypeField ( pointPatchVectorField, libAngularVelocityPointPatchVectorField ); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam // ************************************************************************* // D.2 libAngularVelocityPointPatchVectorField.H Class Foam::libAngularVelocityPointPatchVectorField Description Foam::libAngularVelocityPointPatchVectorField SourceFiles libAngularVelocityPointPatchVectorField.C \*---------------------------------------------------------------------------*/ #ifndef libAngularVelocityPointPatchVectorField_H #define libAngularVelocityPointPatchVectorField_H #include "fixedValuePointPatchField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { /*---------------------------------------------------------------------------*\ Class libAngularVelocityPointPatchVectorField Declaration \*---------------------------------------------------------------------------*/ class libAngularVelocityPointPatchVectorField : public fixedValuePointPatchField<vector> { // Private data vector axis_; vector origin_; scalar angle0_; scalar omega_; pointField p0_; public: //- Runtime type information TypeName("libAngularVelocity"); // Constructors //- Construct from patch and internal field libAngularVelocityPointPatchVectorField ( const pointPatch&, const DimensionedField<vector, pointMesh>& ); //- Construct from patch, internal field and dictionary libAngularVelocityPointPatchVectorField ( const pointPatch&, const DimensionedField<vector, pointMesh>&, const dictionary& ); //- Construct by mapping given patchField<vector> onto a new patch libAngularVelocityPointPatchVectorField ( const libAngularVelocityPointPatchVectorField&, const pointPatch&, const DimensionedField<vector, pointMesh>&, const pointPatchFieldMapper& ); //- Construct and return a clone virtual autoPtr<pointPatchField<vector> > clone() const { return autoPtr<pointPatchField<vector> > ( new libAngularVelocityPointPatchVectorField ( *this ) ); } //- Construct as copy setting internal field reference libAngularVelocityPointPatchVectorField ( const libAngularVelocityPointPatchVectorField&, const DimensionedField<vector, pointMesh>& ); //- Construct and return a clone setting internal field reference virtual autoPtr<pointPatchField<vector> > clone ( const DimensionedField<vector, pointMesh>& iF ) const { return autoPtr<pointPatchField<vector> > ( new libAngularVelocityPointPatchVectorField ( *this, iF ) ); } // Member functions // Evaluation functions //- Update the coefficients associated with the patch field virtual void updateCoeffs(); //- Write virtual void write(Ostream&) const; }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif E – log file – checkMesh Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 275939 faces: 755576 internal faces: 721838 cells: 240822 faces per cell: 6.13488 boundary patches: 7 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 211485 prisms: 6552 wedges: 0 pyramids: 0 tet wedges: 42 tetrahedra: 0 polyhedra: 22743 Breakdown of polyhedra by number of faces: faces number of cells 4 2063 5 1418 6 6164 7 3358 8 2336 9 3378 10 68 11 18 12 3140 15 768 18 32 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology inlet 2500 2601 ok (non-closed singly connected) outlet 2500 2601 ok (non-closed singly connected) down 2500 2601 ok (non-closed singly connected) up 2500 2601 ok (non-closed singly connected) front 2500 2601 ok (non-closed singly connected) back 2500 2601 ok (non-closed singly connected) turbine_patch28539 18738 20574 ok (closed singly connected) Checking geometry... Overall domain bounding box (-15 -15 -15) (15 15 15) Mesh (non-empty, non-wedge) directions (1 1 1) Mesh (non-empty) directions (1 1 1) Boundary openness (-8.78672e-16 3.10042e-17 -1.2724e-16) OK. Max cell openness = 8.3837e-16 OK. Max aspect ratio = 13.9185 OK. Minimum face area = 0.000880194. Maximum face area = 0.411711. Face area magnitudes OK. Min volume = 8.75044e-05. Max volume = 0.236188. Total volume = 26890. Cell volumes OK. Mesh non-orthogonality Max: 59.9817 average: 9.02038 Non-orthogonality check OK. Face pyramids OK. Max skewness = 1.42542 OK. Coupled point location match (average 0) OK. Mesh OK. End