Shape optimisation of an integrated bulb-keel Thesis subtitle Master of Science Thesis KASPER LJUNGQVIST Department of Shipping and Marine Technology Division of Ship Design CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden, 2011 Report No. X-11/268 i A THESIS FOR THE DEGREE OF MASTER OF SCIENCE Shape optimisation of an integrated bulb-keel Thesis subtitle KASPER LJUNGQVIST Department of Shipping and Marine Technology CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2011 ii Shape optimisation of an integrated bulb-keel Thesis subtitle KASPER LJUNGQVIST © KASPER LJUNGQVIST, 2011 Report No. X-11/268 [1]Chalmers University of Technology SE-412 96 Gothenburg Sweden Telephone +46 (0)31-772 1000 Printed by Chalmers Reproservice Gothenburg, Sweden, 2011 iii Shape optimisation of an integrated bulb-keel Thesis subtitle KASPER LJUNGQVIST Department of Shipping and Marine Technology Division of Ship Design Chalmers University of Technology Abstract: The objective of the thesis was to shape optimise an integrated bulb-keel using a numerical optimisation method and to validate how well the used methods suit this kind of problem. The algorithm used for the optimisation was a Non-dominated Sorting Genetic Algorithm, NSGA- II. The optimisation was limited to the bulb part of the keel and 7 design variables were used together with 4 inequality constraints. The hydrodynamic performance of the keel was evaluated using a RANS viscous solver and the objective function used in the optimisation was the ratio of drag coefficient over lift coefficient CD/CL. The numerical CFD calculations were validated with experimental wind tunnel tests, which is a crucial part when investigating new types of numerical flow problems. The CFD calculations were also verified using a Least Square Root method in order to obtain the numerical error. The comparison between the wind tunnel results and the CFD results showed that the CL/CD-ratio was predicted very similarly for the numerical and experimental tests. The optimised keel design had a CD/CL-ratio improvement of around 4% compared to the initial keel design. Finally the keels were mounted on a hull of a 40-foot yacht and tested in a velocity prediction program (VPP) in order to get a more understandable performance indicator of the keels. Keywords: Keel, Bulb, NSGA-II, Optimisation, CFD, RANS, Integrated bulb-keel iv v Preface The thesis was carried out at Chalmers University of Technology in Gothenburg, Sweden, in cooperation with Aalto University, School of Science and Technology in Espoo, Finland. The thesis has therefore been published from both Universities. I would like to acknoledge some people who have been helping me throughout the thesis. Firstly I would like to like to express great gratitude to my examiner at Chalmers University of Technology Professor Lars Larsson for letting me carry out my thesis at Chalmers University of Technology even if I was not registered at the master’s program. I would also like to thank him for all his input in my work and his great knowledge in hydrodynamics and yacht design. Secondly I would like to thank my supervisors Master of Science (Tech.) Florian Vesting and Master of Science (Tech.) Michał Orych. Florian’s knowledge in the Friendship System and in optimisation and Michal’s knowledge in both CFD and yacht design have been have helped me a lot in throughout the thesis. Further I would like to thank supervisor at Aalto University Professor Jerzy Matusiak for being so flexible and letting me carry out the thesis at another university. He has also been very supportive throughout the project with both help concerning hydrodynamic issues but also with practical issues Finally I would like to thank my fellow students Björn Axfors and Hanna Tunander for the help in the wind tunnel and interchange of ideas, Friendship Framework’s software engineer Arne Bergmann for software support, Tommi Mikkola for improvement suggestions and Lotta Olsson for the help with all practical arrangements. Gothenburg, August, 2011 Kasper Ljungqvist vi vii Contents Nomenclature ------------------------------------------------------------------------------------- 1  1.  Introduction ----------------------------------------------------------------------------------------- 5  1.1  Background -------------------------------------------------------------------------------------- 5  1.2  Objective of the investigation ----------------------------------------------------------------- 6  1.3  Methodology ------------------------------------------------------------------------------------- 6  1.4  Limitations --------------------------------------------------------------------------------------- 7  2.  Theory ------------------------------------------------------------------------------------------------ 8  2.1  Introduction -------------------------------------------------------------------------------------- 8  2.1.1  Lift and drag ----------------------------------------------------------------------------- 9  2.1.2  Bulb -------------------------------------------------------------------------------------- 11  2.2  Computational hydrodynamics --------------------------------------------------------------- 12  2.2.1  Governing equations ------------------------------------------------------------------- 12  2.2.2  Turbulence models --------------------------------------------------------------------- 14  2.2.3  Boundary conditions ------------------------------------------------------------------- 15  2.2.4  Overlapping grids ---------------------------------------------------------------------- 17  2.2.5  Discretization --------------------------------------------------------------------------- 19  2.2.6  Verification and Validation ----------------------------------------------------------- 20  2.3  Optimisation ------------------------------------------------------------------------------------ 24  2.3.1  Genetic algorithm ---------------------------------------------------------------------- 24  2.3.2  Non-dominated Sorting Genetic Algorithm, NSGA-II --------------------------- 25  2.4  Velocity Prediction Program VPP ----------------------------------------------------------- 27  3.  Investigated cases ---------------------------------------------------------------------------------- 29  3.1  CFD ---------------------------------------------------------------------------------------------- 29  3.1.1  The initial keel model ----------------------------------------------------------------- 29  3.1.2  Grid generation ------------------------------------------------------------------------- 30  3.1.3  Boundary conditions ------------------------------------------------------------------- 32  3.1.4  CFD solver ------------------------------------------------------------------------------ 33  3.2  Optimisation ------------------------------------------------------------------------------------ 33  3.2.1  Work flow ------------------------------------------------------------------------------- 33  3.2.2  Design variables ------------------------------------------------------------------------ 34  viii 3.2.3  Constraints ------------------------------------------------------------------------------ 36  3.3  Experimental validation ----------------------------------------------------------------------- 37  3.3.1  Physical explanation of the wind tunnel and the keel models ------------------- 37  3.3.2  Alignment of the balance, the keel model and the flow -------------------------- 39  3.3.3  Correction and errors ------------------------------------------------------------------ 40  3.4  Velocity Prediction Program ----------------------------------------------------------------- 42  4.  Results ----------------------------------------------------------------------------------------------- 45  4.1  Grid variation results -------------------------------------------------------------------------- 45  4.1.1  Factor of Safety method --------------------------------------------------------------- 45  4.1.2  Least Square Root method ------------------------------------------------------------ 46  4.2  Results from the initial design ---------------------------------------------------------------- 47  4.2.1  Lift and drag ---------------------------------------------------------------------------- 50  4.3  Results from the optimisation ---------------------------------------------------------------- 51  4.3.1  Convergence ---------------------------------------------------------------------------- 52  4.3.2  Lift and drag ---------------------------------------------------------------------------- 53  4.3.3  Vortex development ------------------------------------------------------------------- 53  4.3.4  Pressure distribution ------------------------------------------------------------------- 54  4.4  VPP results -------------------------------------------------------------------------------------- 55  4.5  Wind tunnel results ---------------------------------------------------------------------------- 55  5.  Discussion ------------------------------------------------------------------------------------------- 57  5.1.1  Uncertainty analysis ------------------------------------------------------------------- 57  5.1.2  Optimisation ---------------------------------------------------------------------------- 57  5.1.3  Design comparison -------------------------------------------------------------------- 58  5.1.4  Problem area in the mesh ------------------------------------------------------------- 61  5.1.5  Comparison between the numerical and experimental results ------------------- 62  6.  Summary -------------------------------------------------------------------------------------------- 64  7.  Future work ---------------------------------------------------------------------------------------- 65  References --------------------------------------------------------------------------------------- 66  APPENDIX 1 -------------------------------------------------------------------------------------- I  APPENDIX 2 --------------------------------------------------------------------------------- VIII  APPENDIX 3 ----------------------------------------------------------------------------------- IX  APPENDIX 4 ----------------------------------------------------------------------------------- XI  ix APPENDIX 5 --------------------------------------------------------------------------------- XIII  APPENDIX 6 --------------------------------------------------------------------------------- XIV  x 1 Nomenclature  Symbols A Cross-section area of a control volume AR Aspect ratio ARe Effective aspect ratio B Centre of buoyancy BG Distance between B and G BM Metacentric radius C Mean chord length, Cross sectional area of wind tunnel C1 Root chord length C2 Tip chord length CD Drag coefficient CDi Induced drag coefficient CD0 Profile drag coefficient, Drag coefficient at zero degrees leeway Cf Local skin friction coefficient CF Total skin friction coefficient CL Lift coefficient D Drag force, Experimental solution Df Friction drag Di Induced drag D0 Form drag E Validation comparison error Fs Safety factor for numerical uncertainty Fx,y,z Body force in x-, y- and z-direction Fx Force in x-direction Fxp Pressure force in x-direction Fxf Frictional force in x-direction Fy Force in y-direction Fyp Pressure force in y-direction Fyf Frictional force in y-direction F(X) Objective function G Centre of gravity 2 g Acceleration of gravity GM Metacentric height GZ Righting arm h1 Typical cell size for the finest grid hi Typical cell size IT Transverse moment of inertia k Turbulent kinetic energy, blockage factor L Characteristic length, Lift force LCG Longitudinal centre of gravity LOA Length over all LWL Length water line M Metacentre m Mass ng Number of grids Ncells Number of grid cells P Distance metric p Pressure pa Atmospheric pressure pRE Observer order of accuracy pth Theoretical order of accuracy p∞ Undisturbed pressure R Convergence ratio Rij Reynolds stress tensor RM Righting moment Rn Reynolds number r Refinement ratio S Planform area SA Sail area Sij Rate of strain tensor SW Wetted surface area T Span Te Effective span TK Keel span  3 t Time TR Taper ratio u Time averaged flow velocity component in x-direction uτ Friction velocity U Free stream velocity Uϕ Numerical uncertainty UD Experimental uncertainty UG Grid discretization uncertainty Ui Instantaneous velocity component in i:th direction UI Iterative uncertainty Uinput Parameter uncertainty Uj Instantaneous velocity component in j:th direction Un Flow velocity’s normal component to a surface US Standard deviation of the numerical uncertainty’s curve fit V Ship velocity v Flow velocity v Time averaged flow velocity component in y-direction vn Normal velocity vparticle Particle velocity VCGkeel Vertical centre of gravity for the keel (from the root chord) w Time averaged flow velocity component in z-direction Wij Rotation rate X Design variable vector y+ Nondimensional wall distance Greek symbols α A constant for estimating the discretization error β* Modelling constant Γ Diffusion coefficient   ρ Density Δx Grid spacing δij  Kronecker delta  δRE  Discretization error 4 εsb Solid blockage correction εwb Wake blockage correction μ Dynamic viscosity μt Turbulent viscosity ν Kinematic viscosity ξB Parameter direction crossing the boundary ρ Density σij Stress tensor τw Wall shear stress φ Heel angle ϕ Dependent variable ΦB Distance function  Φexact Exact solution Φi Integral of local quantity Φo Estimated exact solution ω Specific rate of dissipation of turbulent energy ∇ Volumetric displacement Abbreviations CAD Computer-Aided Design CFD Computational Fluid Dynamics EASM Explicit Algebraic Stress Model FEM Finite Element Method FDM Finite Difference Method FS Factory of Safety FVM Finite Volume Method GA Genetic Algorithm ITTC International Towing Tank Conference LSR Least Square Root NACA National Advisory Committee for Aeronautics NSGA-II Non-dominated Sorting Genetic Algorithm II RANS Reynolds Average Navier-Stokes VPP Velocity Prediction Program 5 1. Introduction  1.1 Background   Yacht design has for long been more of an experience knowledge than an actual science. Little advanced research has been made compared to naval architecture for large ships. Keels, rudders and other parts of the yacht have often been selected based on experience and rules of thumb and the calculations have been more approximations than in depth analysis. It has mostly been racing yacht projects with large budgets that have been able to use advanced design tools in their analysis. Little of the results have been published for the public due to the rivalry between the racing teams. However, in the last years more research has been conducted in yacht design and the cruising yacht industry has become more involved in the academic research field. Computational Fluid Dynamics (CFD) has also become a more accepted method for hydrodynamic calculations mainly in the ship industry but it has slowly also become more accepted in the sailing yacht industry. Figure 1.1 The 4 different keels used in the bachelor thesis. The fin keel, the integrated bulb-keel, the L- bulb keel and the T-bulb keel [1] In the spring of 2010 Bergman et al. [1] carried out a bachelor thesis project at Chalmers University of Technology. In their thesis “Investigation of keel bulbs for sailing yachts” the students compared with CFD calculations the hydrodynamic performance of 4 different sailing yacht keels. The keels chosen in the investigation were 4 commonly used keel types (see Figure 1.1), a conventional fin keel, an integrated bulb-keel, an L-bulb keel and a T-bulb keel. The interested reader is directed to reference [1] for the complete description of all keels. The conclusion from the thesis was that the integrated bulb-keel had the best hydrodynamic performance of the 3 bulb keels; the fin keel had the best performance but it cannot provide the same righting moment as the bulb keels if all keels were made out of the same material. The keels were designed in a CAD program together with professional yacht designer Stefan Qviberg. Even if the keel types used in the investigation have been commonly used, little knowledge exists about their optimum shape. Hence the interest in optimising the keels evoked the idea of this thesis. It would be very interesting to optimise all the keels but due to time limitation and work load only one keel could be chosen. The integrated bulb-keel was chosen based on the results from the bachelor thesis and this keel is a very convenient choice even if the results have been questioned by Axfors & Tunander [2]. An integrated bulb- keel is practical because of its robustness, which makes the keel more impact resistant and less risk of getting items tangled up around it. Bergman et al. used the same CFD solver as the one that has been used in thesis, however the grid generator and the grid structure were different. 6 This thesis is divided into 7 chapters. First the introduction of the whole thesis will be presented together with the objective, methodology and limitation of the work. The theories behind the thesis are introduced in the second chapter followed by the investigated cases in chapter 3, where all the configurations for the calculations will be presented together with the setup of the wind tunnel. In chapter 4 the results from the calculations and from the experimental test of the wind tunnel are displayed and the discussion is carried out in chapter 5. Finally a summary of the thesis is presented in chapter 6 and the thesis is finished by some recommendations from the author for future study. 1.2 Objective of the investigation   The major aim of this thesis is to optimise the shape of a particular integrated bulb- keel. In this particular case shape optimisation means to improve the keels lifting force while decreasing its drag and maintaining a sufficient righting moment. Another objective of the thesis is to validate the numerical methods used in the process. The calculations are validated with experimental wind tunnel tests. Additionally an objective is to investigate whether the optimisation method used is sufficient for the problem of the thesis with the set limitations and constraints. 1.3 Methodology  The CFD solver that has been used is in this thesis is SHIPFLOW 4.4.04 and more specifically XCHAP, which is method for solving the Reynolds Average Navier- Stokes equations. XCHAP is using a finite volume method as discretization method and an EASM turbulence model has been implemented. Two methods have been used for the uncertainty analysis namely the Factor of Safety method and the Least Square Root method and the results have been compared and validated. A tool called Friendship Framework has been used for the optimisation. Friendship Framework integrates a number of useful computer programs needed in hydrodynamic optimisation. In this thesis, Friendship Framework’s CAD tool, generic integration tool and optimisation tool were used together with the integration of SHIPFLOW 4.4.04. The Non-dominated Sorting Genetic Algorithm-II (NSGA-II) was used for the optimisation algorithm. A velocity prediction program VPP has been used in order to get a more understandable performance indicator of the keels. The velocity prediction program takes also the stability aspect of the keels into account. WinDesign VPP developed by the Wolfson Unit was the chosen velocity prediction program for this thesis. Validation of results is a very important process when numerical calculations are carried out. Therefore wind tunnel tests were carried out in order to validate the CFD results calculated with SHIPFLOW. The wind tunnel tests were performed in the wind tunnel at the Department of Applied Mechanics at Chalmers University of Technology in Gothenburg, Sweden. 7 1.4 Limitations  In order to increase the quality of the thesis and to fulfil the requirements set for the thesis some limitations have been set. 1) The optimisation will be done for one angle of attack, since it would require too much computational effort to take more angles of attack into account. 2) The heeling angle was set to zero in both the numerical tests and in the experimental tests. 3) The fin part of the keel was not optimised in the optimisation process, since the optimisation was limited to the bulb part of the keel. The bulb was the part that was the focus in this thesis, although a global optimisation would be interesting. 4) Neither the strength aspect of the keels nor the manufacturing possibilities have been taken into account when designing the keels. 8 2. Theory  The theory of a yacht keel and the bulb will be presented in this chapter, together with optimisation theory and the computational hydrodynamics that has been used are explained. 2.1 Introduction  In this section the theory of the keel and the bulb will be presented. The hydrodynamic phenomena that occur around the keel will be introduced together with the stability aspects of the keel. The keel on a sailing yacht has two major tasks. The first one is to provide sufficient stability when the boat heels and the second one to counteract the aerodynamic force of the sail. The hydrodynamic force that counteracts the aerodynamic force is called lift L and it is developed by the pressure difference on the keel, the rudder and the hull. However the keel is the biggest contributor to the lift force in most cases. When a yacht sails upwind it will drift sideways due to the wind. Thus the yacht moves in a direction different from the longitudinal direction of the hull and the angle between the two directions is called leeway angle. The pressure difference on the keel is caused by the leeway angle. The keel will also develop a resistance component, which is called drag D [3]. The lift and drag forces will be explained in more detail in section 2.1.1. A sailing yacht has a significant difference compared to other seagoing vessels, since it heels considerably more. The reason for the large heel is the aerodynamic force acting on the sails of the boat. The aerodynamic force heels the boat over and this heeling moment has to be counteracted. Thus a righting moment RM is required. The righting moment is defined according to Rossell & Chapman in Principles of Naval Architecture [4] with equation (2.1), where m is the mass of the yacht, g the acceleration of gravity and GZ the righting arm. Equation (2.2) and (2.4) are also taken from Principles of Naval Architecture [4]. The righting arm is the distance between the centre of gravity of the boat G and the Z point in Figure 2.1. RM = m ∗ g ∗ GZ (2.1) The righting arm (2.2) is dependent on the heel angle φ and on the metacentric height GM. The metacentric height is the distance between the centre of gravity and the metacentre M, see Figure 2.1, and it is also a measurement of the initial stability of the boat. Equation (2.2) is only valid for small heel angles for floating bodies that do not have circular or spherical cross-sections [4]. However, the curve for the righting arm for most sailing yachts follows a sinusoidal path longer than the curve for large ships, because of the rounder transversal cross-section of the sailing yachts. Hence the equation is valid for larger angles of attack for most sailing yachts. GZ = GM ∗ sin φ (2.2) The metacentre is dependent on the shape of the hull and a high metacentre will give a high form stability. The exact formulation of the metacentric height can be seen in 9 equation (2.3), where BM is the distance between the metacentre and the centre of buoyancy and BG is the distance between the centre of gravity and centre of buoyancy B. BM is called the metacentric radius and is calculated with equation (2.4), where IT is the transverse moment of inertia of the water plane and ∇ is the volumetric displacement. The centre of buoyancy B will move to B’ when the boat heels and the buoyancy force will act at B’, see Figure 2.1. GM = BM− BG  (2.3) BM = !! ∇   (2.4) Figure 2.1 Transverse stability [3] The stability of the yacht can be increased by increasing the height of the metacentre or by lowering the centre of gravity. The metacentre can be difficult to move since it depends on the hull shape. Hence, a more convenient way of increasing the GM value is to lower the centre of gravity. The keel has a significant role of the yacht’s stability since it is has its centre of gravity very low. Therefore the keel is often made of material with high density in order to have a higher effect on the total centre of gravity. The total centre of gravity can further be lowered if a bulb is mounted on the keel [3]. 2.1.1 Lift and drag A keel can be seen as a wing and the fundamental aerodynamic theories can be applied. The lift is caused by the pressure difference on the keel as described above. Larsson & Eliasson [3] explain that when the flow meets the keel at an angle of attack (leeway angle) the flow becomes asymmetric over the keel. The point where the flow meets the keel is called the stagnation point and the velocity is zero at that point. The flow separates into two directions from that point, some of the flow flows around the upper side (suction side) of the wing (see Figure 2.2) and the rest of the flow follows the lower side (pressure side). The flow on the suction side is faster 10 than the flow on the pressure side and a higher velocity gives a lower pressure and a lower velocity gives a higher pressure respectively. Therefore the pressure difference is developed and the lift force L is created. The direction of the lift force is perpendicular to the flow direction. The lift force can be nondimensionalized with equation (2.5), where CL is the lift coefficient, U the velocity of the flow, ρ the density of the fluid and S is the planform area according to Abbot & von Doenhoff [5]. C! = ! !.!∗!∗!!∗! (2.5) The planform area is taken as mean chord length multiplied by the keel span. The resistance of the keel can be divided into three components, the induced drag Di, the form drag D0 and the frictional resistance DF. The induced drag is dependent on the pressure difference together with the aspect ratio. The form drag is an effect of the shape of the keel and the frictional resistance is caused by the friction between the keel and the water. According to Larsson & Eliasson [3] the pressure difference, which is caused by the leeway angle, makes the flow move from the pressure side to the suction side around the tip and this creates a slight upward flow on the suction side and a slight downward flow on the pressure side. The upward and downward flow are strongest close to the tip and weakest at the root. When the flow from the pressure and section side meet at the trailing edge they will have different directions due to the upward and downward flows and this will create vortices, which will be strongest close to the tip because the vertical motions are strongest in that area. All vortices merge into one big vortex, which has a lot of rotational energy. The rotational energy corresponds to the induced resistance for the keel. The drag D is nondimensionalized in the same way as the lift and the drag coefficient CD is given by equation (2.6) according to reference [5]. C! = ! !.!∗!∗!!∗! (2.6) The lift and the drag are dependent on the planform of the keel. The main parameters for a keel planform are the aspect ratio AR, the taper ratio TR and the sweep angle (the sweep angle for the 25% chord). The aspect ratio is defined by equation (2.7) according to Abbot & von Doenhoff [5] and the efficiency of the keel is very dependent on the aspect ratio. TK is the span of the keel and S is the planform area. AR = !!! ! (2.7) The taper ratio on the other hand is the ratio between the length of the root chord C1 and the tip chord C2 [3]. TR = !! !! (2.8) An elliptical force distribution over a wing is the optimal one (Abbot & von Doenhoff [5]). A keel with an elliptical force distribution over the span will have an effective aspect ratio ARe which is equal to the geometrical aspect ratio defined in 11 equation (2.7), any force distribution that differs from an elliptical will have a decreased effective aspect ratio. The location of the trailing vortex left behind the keel has an effect on the effective aspect ratio, since a higher vortex will decrease the effective aspect ratio. Thus the effective aspect ratio is ruled by the vortex and not the geometrical span of the keel. However the keel depth will have an effect on the vortex. The effective keel depth, which is dependent on the effective aspect ratio, can be calculated with equation (2.9) according to reference [6] T! = !"!∗! ! (2.9) The effective aspect ratio is further calculated with equation (2.10), where CDi is the induced drag. AR! = !! ! !∗!!" (2.10) The tip of the keel (on a fin keel) has also an effect on the efficiency of the keel. Different shapes of the tip have an effect on the location of the trailing vortex. The best tip shape is a square tip according to Larsson & Eliasson [3]. The square tip moves the vortex below the tip due to the sharp edge, thus the flow that is created around the tip separates quite early. Compared to a rounded tip that allows the flow to move around the tip and this will cause the flow to separate above the tip and decrease the effective aspect ratio. Figure 2.2 Concept picture of the pressure and suction side on a wing section There has been conducted much research on the sectional shapes of wings and there are recommendations which sections to use in different circumstances. The NACA- series are especially well known and their advantages and disadvantages are well documented. The NACA-series were developed for aircraft wings by the National Advisory Committee for Aeronautics [5]. 2.1.2 Bulb A bulb is attached to the fin keel if the keel itself does not have enough weight to give sufficient stability to the boat. A fin keel is the most typical keel type, since the keel looks like an airplane wing. The first keel in Figure 1.1 is a fin keel. The bulb will lower the vertical centre of gravity of a keel and thereby increase the stability. A bulb will in most cases have a negative effect on the hydrodynamic performance of the keel but the loss in hydrodynamic performance is compensated with the increase 12 in stability. There are many rules and theories (Larsson & Eliasson [3]) for yacht keels but most of the rules are for fin keels or for long keels. Rules and theories on the effect of the bulb has on the efficiency of a keel are rare. The frictional resistance is naturally increased if a bulb is mounted on a keel, due to the increased wetted surface. Another effect that bulb with a round sectional shape has is that it moves the trailing vortex further up and this will decrease the effective aspect ratio (Larsson & Eliasson [3]). 2.2 Computational hydrodynamics  In this section the governing equations used in computational hydrodynamics and more specifically in this thesis will be briefly explained. Further the turbulence models, boundary conditions, the overlapping grid concept and the discretization method will be explained. In the end of the section, theory of the verification and validation methods that have been used are described. 2.2.1 Governing equations One of the governing equations used is the continuity equation, which explains the conservation of mass. If the fluid is incompressible, which means that the density ρ is constant, the continuity equation can be written in accordance with reference [7]: ∇ • 𝐯 = 0, (2.11) where v is the flow velocity, in Cartesian coordinates the equation would be: !! !! + !! !! + !! !! = 0, (2.12) where u, v and w are the velocity components in x-, y- and z-direction. The motions of the fluid are explained with the Navier-Stokes equations, which illustrate the flow fields when solved. The interested reader is directed to reference [7] or [8] for the full derivation of the equations. Navier-Stokes equations in Cartesian coordinates, which can be found in most fluid dynamics literature, i.e. Matusiak [7], can for the different directions be written: !! !! + u !! !! + v !! !! + w !! !! = − ! ! !! !! + ν !!! !!! + !!! !!! + !!! !!! + ! ! F! (2.13) !! !! + u !! !! + v !! !! + w !! !! = − ! ! !! !! + ν !!! !!! + !!! !!! + !!! !!! + ! ! F! (2.14) !! !! + u !! !! + v !! !! + w !! !! = − ! ! !! !! + ν !!! !!! + !!! !!! + !!! !!! + ! ! F! (2.15) where the left-hand side is the so-called total acceleration in the given direction and the right-hand side is the pressure force, viscous force and the body force. All forces are per unit mass. The body forces are the external forces acting on the fluid. The pressure force is the first term on the right-hand side, where p denotes the pressure. The second term is the viscous force, where ν denotes the kinematic viscosity. The viscous force is due to the viscous stresses acting on the fluid elements [7]. 13 To solve the Navier-Stokes equations numerically would require tremendous amounts of computational time, and although the results would be very accurate it is not in reality feasible to use so much computational time. Hence some methods have been developed in order to lower the computational time. The Reynolds Average Navier-Stokes Equations (RANS) is one method for lowering the computational effort. The method is still very accurate and it can solve complex flow problems [8]. The RANS equations are time-averaged Navier-Stokes equations and the instantaneous values for the velocity ui, pressure p and stress tensor σij all consist of one fluctuating component and one average component. Thus the velocities would be denoted: u! = u! + u!′ (2.16) where u! is the average velocity components in x-, y- or z-direction and u!′ is the fluctuating component. The pressure and the stress sensor can be denoted respectively. The final form of the RANS equations can be seen in equation (2.17), (2.18) and (2.19) (Matusiak [7]). The complete derivation of the equation can be found in reference [7] and [8]. ρ !! !! + u !! !! + v !! !! + w !! !! = − !! !! + μ !!! !!! + !!! !!! + !!! !!! − ρ !!"! !! + !!"#" !! + !!"#" !! + F! (2.17) ρ !! !! + u !! !! + v !! !! + w !! !! = − !! !! + μ !!! !!! + !!! !!! + !!! !!! − ρ !!"#" !! + !!"! !! + !!"#" !! + F! (2.18) ρ !! !! + u !! !! + v !! !! + w !! !! = − !! !! + μ !!! !!! + !!! !!! + !!! !!! − ρ !!"#" !! + !!"#" !! + !!"! !! + F! (2.19) where the third term on the right-hand side of the equations is the turbulence contribution to the mean flow and μ is the dynamic viscosity. μ = νρ (2.20) The so-called Reynolds stress tensor Rij (2.21), which is included in the third term on the right-hand side, is unknown and it must be modelled [7]. R!" = − ρu!! ρv!u! ρw!u! ρu!v! ρv!! ρw!v! ρu′w′ ρv!w! ρw!! = −ρu!′v!′ (2.21) The Reynolds stress tensor is symmetric, thus ρu′v′ = ρv′u′, ρu′w′ = ρw′u′ and ρv′w′ = ρw′v′. Hence there are 6 unknowns but still only four equations (3 Navier-Stokes equations and the continuity equation), which mean that the unknowns cannot be solved; this is known as the closure problem [9]. In order to be able to solve the equations turbulence models are introduced. 14 2.2.2 Turbulence models Many turbulence models have been introduced and all of them are using empirical constants. The different models have different advantages and disadvantages for different flow problems. The turbulence models can be grouped in to five classes [8]: 1. Zero-equation models 2. One-equation models 3. Two-equation models 4. Algebraic stress model (ASM) 5. Reynolds stress model (RSM) In this thesis the Explicit Algebraic Stress Model (EASM) has been used, hence it is the only turbulence model that will be explained. The EASM is developed from the Boussinesq assumption but non-linear terms are added. The Boussinesq assumption assumes that the Reynolds stress tensor is dependent to the rate of strain tensor Sij in similar fashion as the viscous stresses (2.22), but instead of using the molecular viscosity µ the turbulent viscosity µt is used. σ!" = μS!" (2.22)   The turbulent viscosity is not constant over the flow field as the molecular viscosity; it will vary depending on the flow and the location according to Larsson & Raven [8]. According to the continuity equation (2.12) the rate of strain tensor would be 0 for Sii, this would give a zero value for the Reynolds stress tensor Rii as well. However if equation (2.21) and (2.23) are combined equation (2.24) can be derived. k = !!!!!! ! (2.23) R!! = −2ρk (2.24) where k is the turbulent kinematic energy. From the equations above the final Boussinesq assumption is derived (equation (2.25)) where δij is the Kronecker delta and it is 1 if i=j and 0 if i≠j which is presented by Larsson & Raven in reference [8]. R!" = μ!S!" − ! ! ρkδ!" (2.25) The linear Boussinesq assumption often fails to compute complex three-dimensional flows; therefore the EASM has been developed [10]. In the EASM the Reynolds stress tensor is extended with non-linear terms and it is calculated with equation (2.26) presented in XCHAP’s Theoretical Manual [9]. R!" = − ! ! ρkδ!" + μ! S!" + a!a! S!"W!" −W!"S!" − a!a! S!"S!" − ! ! S!"S!"δ!" (2.26) where 15 𝑎! = ! ! ! ! − 𝐶! (2.27) 𝑎! = ! ! 2− 𝐶! (2.28) 𝑎! = ! ! 2− 𝐶! (2.29) 𝑎! = ! !!!!!! !!/! !!!! (2.30) 𝛾! = !!! ! (2.31) γ! = !!! ! + !!"!!!" !!"!! (2.32) τ = ! ! (2.33) 𝜂! = S!"S!" (2.34) R! = !!"!!" !! (2.35) μ! = ν!ρ = max  −kα!, !.!!!"# !∗! (2.36) The specific rate of dissipation of turbulent energy is denoted ω and β* is a modelling constant. The rotation rate Wij is obtained with equation (2.37), where Ui and Uj are the instantaneous velocity components. W!" = ! ! !!! !!! − !!! !!! (2.37) α1 in equation (2.30) and (2.36) can be solved from equation (2.38) and the correct α1 is the solution with the lowest real part. !! ! ! − !! !!!!!! !! ! ! + !!!!!!!!!!!!! ! !! !!!!!!!!!!!!!!!!! !!!!!!! ! !! ! = − !!!! !!!!!!! ! (2.38) Cϵ1=1.44, Cϵ2=1.83, C2=0.36, C3=1.25, C4=0.4, C1 0=3.4 and C1 1=1.8 are all also modelling coefficients [9].   2.2.3 Boundary conditions The continuity equation and the Navier-Stokes equations are not enough to solve a flow field; boundary conditions have to be introduced in order to give boundary values for the equations. The boundary conditions can be grouped into three groups the solid boundary, the free surface and infinity as presented in reference [7]. Furthermore there will be boundary condition for the computational domain used in the calculations. These boundary conditions are modifications of the mathematical boundary conditions named above. The solid boundary conditions are the boundary conditions for the areas where the fluid meets solid surface. The surfaces could for instance be a hull or the seabed. The first condition is the no-leak condition, which states that there cannot be any flow through a solid surface. Therefore, the velocity’s normal component to the surface must be zero [7]. 𝐯! = 0 (2.39) 16 The velocities at the hull surface must be 0 if the coordinate system is moving along with the hull. Thus the tangential velocity at the surface is zero. 𝐯! = 0 (2.40) For surfaces that are not moving along with the surface the boundary condition can be replaced with u = 𝐕, v = w = 0 (2.41) where V is the ship velocity. The free surface boundary conditions are the boundary conditions for surface between two liquids, in most cases the water surface between air and water. It is usually assumed that there is a known atmospheric pressure in the air pa. The Reynolds stress and the surface tension are assumed to be zero. In reality they are very small and therefore they are neglected. Hence the atmospheric pressure is the stress and therefore the first free surface boundary condition, which is also known as the dynamic boundary condition, can be written [7]: p = p! (2.42) The second free surface boundary condition is known as the kinematic boundary condition and it states that a particle on the free surface has the same velocity in the normal direction as the free surface. v!"#$%&'( = v! (2.43) Where vn is the normal component of the free surface velocity and it is the only non- zero component because there is no shear stress at the surface. The waves that a ship or a boat generates are stationary if the coordinate system moves along with the vessel and in that case the free surface velocity is zero [7]. v!"#$%&'( = v! = 0 (2.44) The boundary conditions at infinity are used in order to have a place where the flow is completely undisturbed. In reality the place is of course not at infinity but a distance sufficiently far enough so that the flow is undisturbed and the following conditions are fulfilled: u = 𝐔 (2.45) v = w=0 (2.46) p = p! (2.47) where U is the free stream velocity and p∞ is the undisturbed pressure [8]. As mentioned above some of the boundary conditions have to be modified for computations because the computational domain might not be the same as reality. The computational boundaries may vary for different cases and they can be chosen arbitrary. 17 2.2.4 Overlapping grids In order to solve numerically the RANS equations with a turbulence model and the continuity equation one or more grids are required together with a discretization scheme (see section 2.2.5). The grids can be structured or unstructured and they must cover the whole computational domain. The structured grid is built up of cells in axial i, radial j and circumferential k direction. In a structured grid a mapping function can be used in order to transform the physical domain, which might be curve-linear, into the computational domain with Cartesian coordinates. In three-dimensional cases complete planes must be added when the cell number is increased, similarly complete lines must be added in two-dimensional cases. The cells in a three-dimensional structured grid are hexahedral. Structured grids are normally of H-, O- or C-type. An H-type grid is a grid where all the gridlines go from one edge to the opposite edge and an O-type grid is circular, see Figure 2.3. The O-grid can be a full circle or a half circle. Figure 2.3 H-grid (left) and O-grid (right). The topography is normally combined of two grid types and H-O-topography is especially common in hydrodynamics. The grid can be divided into several blocks or it can be a single block. The advantage of the structured grids is the faster computational time compared to the unstructured grid due to the regularity of the grid. The unstructured grid’s advantage is the flexibility; a single point can be added without having to create a complete plane. This will avoid areas with unnecessary cells, which often occur in structured grids due to the complete plane (or line) insert. This will however increase the computational effort [8]. Another method of increasing the flexibility is to use overlapping grids with two or more structured grids. By doing this the flexibility of the unstructured grid is almost obtained while keeping the most of the regularity and low computational effort of the structured grid. Overlapping grids are good in shape optimisation and for meshing difficult geometries. An example of an overlapping grid can be seen in Figure 2.4. 18 Figure 2.4 Overlapping grids Overlapping grid components can be both boundary fitted and non-boundary fitted. The non-boundary fitted components penetrate the boundaries and therefore special techniques are required to only use the boundary fitted grids close to a boundary. This calls for a robust prioritization of the grids. The boundary fitted grids are built up around the objects and this makes the discretization easier. Since overlapping grids are overlapping each other cells will also be overlapping each other. The cells that are from a grid with higher priority are normally prioritised, but this may cause problem at certain areas. The problem areas are the areas where a non-boundary fitted grid with higher priority is too close to a boundary of a boundary fitted grid with lower priority. Those areas are better represented by the grid that is fitted to the boundary, so therefore the priority rule cannot be used in those areas. In XCHAP this is solved with a distance function ϕB. This function calculates the shortest distance for a cell to the nearest boundary, thus at the physical boundaries the functions value is zero. If two cells overlap each other the cell with the lowest ϕB is chosen [11]. The grid close to the wall has to be very fine in order to resolve the boundary layer. In equation (2.48) y+ should not be larger than 1. y is the distance from the wall to the first discretization point and y+ is the nondimensional distance. Thus, the discretization point must be well inside the viscous sublayer, which is the inner layer of the turbulent boundary layer [8]. y! = !!! ! (2.48) uτ is the friction velocity and is derived from equation (2.49). 19 u! = !! ! (2.49) The wall shear stress τw can be expressed with: τ! = C! ! ! ρU! (2.50) where Cf is the local skin friction coefficient [8]. The skin friction coefficient is dependent on the Reynolds number Rn, equation (2.51), where L is the characteristic length of the analysed object. R! = !" ! (2.51) C! = !.!"#$ !! ! (2.52) Equations (2.48)- (2.52) are all presented in reference [8]. 2.2.5 Discretization There are three common discretization methods in numerical calculation methods, Finite Volume Method (FVM), Finite Difference Method (FDM) and Finite Element Method (FEM). FEM has hardly been used in computational hydrodynamics while the other two are more common. FVM is the most used within computational hydrodynamics and also the discretization method used in XCHAP and therefore it is the only method that will be explained in this thesis, it will be explained in accordance with reference [8]. Figure 2.5 Finite Volume Method 1D [8] In Figure 2.5 a conceptual picture of the FVM in one-dimension can be seen. The discretization is applied on node P and the capital letters W (west) and E (east) are the notations for the surrounding nodes in the control volume. The faces on each side of the node P are denoted with the small letters w (west) and e (east). The distances between the different nodes are also shown in the figure. The distance between W and P, δXWP and the distance between P and E, δXPE do not have to be equal. The objective is to link the dependent variable ϕ  (i.e.  ui,  k,  p,  ω) at P to the dependent variable ϕ at W and E and this is done by using Gauss Theorem. Thus the transport equation of the RANS equation, together with the turbulence models can be written in a steady state form: 20 ! !" (uϕ) = ! !" (Γ !! !" ) (2.53) The diffusion coefficient Γ is dependent on the dependent variable. Equation (2.53) applied on the one-dimensional case in Figure 2.5 with a cross-sectional area A of the control volume can be written: uAϕ ! − (uAϕ)! = (ΓA !! !" )! − (ΓA !! !" )! (2.54) The continuity equation (2.11) is also needed for the derivation, thus uA ! − (uA)! = 0 (2.55) (2.56) In a three-dimensional case 4 more nodes of the discretization point P have to be introduced n (north), s (south), t (top) and b (bottom) [8]. 2.2.6 Verification and Validation Since the numerical results from CFD calculations are just models which resemble real flow problems, the order of accuracy of the results is an important concern in order to know how close to the real solutions the model solutions are. In order to estimate the error a verification and validation method is used. P.J. Roache stated that verification is a process of “solving the equation right” while validation is a process of “solving the right equation” [12]. Thus verification is a method to estimate the numerical error and validation is a method of showing the model error. The validation is usually done through experimental testing such as wind tunnel tests and towing tank tests [13]. In this thesis two verification methods have been used and the reason for this is that they may lead to different results and the used method will be based on the outcome of the results. The two used methods are the Factor of Safety (FS) method and the Least Square Root (LSR) method. An assumption is made that the iterative uncertainty UI and grid discretization uncertainty UG are the only contributors to the numerical uncertainty Uϕ. The iterative error is dependent on the convergence of the calculated solution, thus even though the convergence criterion is met there is still a small error due to the lack of complete convergence. The discretization error is caused by the discretization of the equations used [8]. The round-off error is the last contributor to the numerical uncertainty but it is usually neglected in practical CFD problems [13]. Therefore the numerical uncertainty can be approximated with equation (2.57) as presented in reference [14]. 𝑈! = 𝑈!! + 𝑈!! (2.57) Both methods are based on grid convergence studies and compare all other grids to the finest grid. The purpose of the verification is to find a confidence interval where the exact solution lies within with a 95% probability, equation (2.58). ϕ− 𝑈! ≤ ϕ!"#$% ≤ ϕ+ 𝑈! (2.58) 21 Where ϕexact is the exact solution and ϕ is the numerical solution of the finest grid. The exact solution is the imagined solution for a grid with zero step size. [15]. Factor of Safety method The FS method is a method developed by Xing & Stern [15] where the uncertainty and error are estimated with a correction factor. This is similar to the method that the ITTC (International Towing Tank Conference) recommends. The method will be explained in accordance with reference [15]. The method requires the iterative uncertainty to be at least one order of magnitude smaller than the grid discretization error. Three grids (triplets), which have been systematically changed, are analysed in this method. The grids used are refined with a refinement ratio r, 𝑟 = ∆!! ∆!! = ∆!! ∆!! , (2.59) where Δx are the grid spacing for the three grids. Grid number 1 is the fine grid, grid number 2 is the medium grid and grid number 3 is the coarse grid. In order to decide the type of convergence a convergence ratio R is introduced. R is the ratio between the solution differences according to equation (2.60). 𝑅 = !!!!! !!!!! (2.60) The convergence is monotonic if 00 and oscillatory divergent if R<-1. The discretization error δRE, is calculated with Richardson Extrapolation equation (2.61), where ϕi is the local quantity, ϕo the estimated exact solution, α is a constant, hi the typical cell size and p is the order of accuracy. 𝜖 ≅ δ!" = ϕ! − ϕ! = αh! ! (2.61) Once the type of convergence has been determined the order of accuracy can be calculated. The theoretical order of accuracy pth =2 while the determined order of accuracy pRE is calculated with equation (2.62) followed by the error estimate equation (2.63). 𝑝!" = !"  !!!!!!!!!! !" (!) (2.62) 𝛿!" = !!!!! !!!!!! (2.63) The calculated order of accuracy should theoretically be the same as the theoretical order of accuracy pRE=pth. This means that the solution is within the asymptotic range. However this is rarely the case in reality, because pRE differs consistently in real cases from pth [15], and the ratio (distance metric) is denoted: 22 𝑃 = !!" !!! (2.64) The numerical uncertainty can generally be expressed as a factor of safety multiplied by the error, equation (2.65). 𝑈! = 𝐹!  𝛿!" (2.65) Xing et al. [15] introduces a factor of safety dependent on the distance metric and the final uncertainty is calculated according to equation (2.66). The safety factors are based on statistical analysis. 𝑈!" = 𝐹!  𝑃 𝛿!" = (2.45− 0.85𝑃) 𝛿!" , 0 < 𝑃 ≤ 1 16.4𝑃 − 14.8 𝛿!" ,𝑃 > 1 (2.66) Least Square Root method The LSR method presented below is developed by Eça et al. and described in reference [13] and [14]. The LSR method discards the iterative error and the condition for using the method is that the iterative error must be two orders of magnitude smaller than the discretization error. Thus if the iterative error is 100 times smaller than the discretization error the iterative error can be neglected. The LSR is also using Richardson Extrapolation equation (2.61). There are three unknowns ϕo, α and p in the equation, hence three grids are needed in order to solve the equations, but the method itself requires at least 4 grids [14] this will make the solution more certain. Thus one equation of form (2.61) is set up for every grid and the unknown’s ϕo, α and p are solved with a Least Square fit. By doing this a curve is fitted to the obtained results. The result is plotted against hi/h1, which can be calculated with equation (2.67). !! !! = !!"##$ ! !!"##$ ! ! (2.67) The convergence of the grid discretization is monotonic if p>0, oscillatory if nch≥INT(ng/3) where nch is the number of triplets with (ϕi+1‐ ϕi)(ϕi‐ ϕi‐1)<0 and the grid discretization is not behaving regularly in other cases. However equation (2.61) is only valid when the convergence is monotonic. Therefore the additional equations below are implemented, solved and curve fitted when the convergence criteria are not met. δ!" !" = ϕ! − ϕ! = α!"h! (2.68) δ!" !" = ϕ! − ϕ! =∝!! h+ α!"h! (2.69) 23 δ∆! = ∆! !!! !! !! (2.70) Where ΔM=max(|ϕi- ϕj|), 1≥i, j≥ng The discretization error δRE must be converted into a numerical uncertainty and this is done using the following formulas:   0.95 ≤ 𝑝 ≤ 2.05:  𝑈! = 1.25𝛿!" + 𝑈! (2.71)   0 < 𝑝 ≤ 0.95:  𝑈! = min (1.25𝛿!" + 𝑈!, 3𝛿!" !" + 𝑈!!") (2.72)   𝑝 ≥ 2.05:  𝑈! = max (1.25𝛿!" + 𝑈!, 3𝛿!" !" + 𝑈!!") (2.73) The numerical uncertainty is obtained from equation (2.74) if the convergence is oscillatory.   𝑈! = 3𝛿∆! (2.74) If the convergence is behaving randomly the numerical uncertainty is obtained with equation (2.75).   𝑈! = min (3𝛿∆! , 3𝛿!" !" + 𝑈!!") (2.75) The least squares fits standard deviations are denoted US, US 02 , and US 12 [14] and they are calculated with equation (2.76).   𝑈! = !!!(!!") ! !! !!! !!!! (2.76) Where δRE is replaced with equation (2.61) for US, with equation (2.68) for US 02 and with equation (2.69) for US 12 as presented by Eça et al. [16]. Validation The validation uncertainty is calculated, according to reference [14], with equation (2.77). UΦ is the numerical uncertainty and UD is the uncertainty obtained in experimental tests. Uinput is the uncertainty for the fluid properties and geometry and it is called the parameter uncertainty.   𝑈!"# = 𝑈!! + 𝑈!"#$%! + 𝑈!! (2.77) The validation uncertainty is compared to the comparison error E, which is taken as the numerical solution S subtracted by the experimental solution D.   𝐸 = 𝑆 − 𝐷 (2.78) 24 Eça et al. state [14] that the solution is validated if |E|>Uval. However, if |E|