Prediction of High-Speed Planing Hull Resistance and Running Attitude A Numerical Study Using Computational Fluid Dynamics Thesis for the Degree of Master of Science DAVID FRISK LINDA TEGEHALL Department of Shipping and Marine Technology CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2015 Master’s thesis X–15/320 Prediction of High-Speed Planing Hull Resistance and Running Attitude A Numerical Study Using Computational Fluid Dynamics DAVID FRISK LINDA TEGEHALL Department of Shipping and Marine Technology Division of Marine Technology Chalmers University of Technology Gothenburg, Sweden 2015 Prediction of High-Speed Planing Hull Resistance and Running Attitude A Numerical Study Using Computational Fluid Dynamics DAVID FRISK, LINDA TEGEHALL © DAVID FRISK, LINDA TEGEHALL, 2015. Supervisors: Fredrik Carlsson, FS Dynamics Sweden AB Arash Eslamdoost, Department of Shipping and Marine Technology Examiner: Rickard Bensow, Department of Shipping and Marine Technology Master’s Thesis X–15/320 Department of Shipping and Marine Technology Division of Marine Technology Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Free surface around a planing hull at 37.5 knots. Typeset in LATEX Printed by Chalmers Reproservice Gothenburg, Sweden 2015 iv Prediction of High-Speed Planing Hull Resistance and Running Attitude A Numerical Study Using Computational Fluid Dynamics DAVID FRISK, LINDA TEGEHALL Department of Shipping and Marine Technology Chalmers University of Technology Abstract Accurate predictions of the resistance and running attitude are key steps in the process of hull design and manufacturing. The predictions have traditionally relied on model testing, but this technique is both expensive and time consuming. In this study, the performance of CFD simulations of planing hulls is evaluated us- ing two commercial software: ANSYS FLUENT, developed by ANSYS, Inc., and STAR-CCM+, developed by CD-adapco. This was done by predicting the steady resistance, sinkage and trim angle of one semi-planing and one planing hull in calm, unrestricted water. The Reynolds averaged Navier-Stokes equations with the SST k-ω turbulence model was used along with the volume of fluid method to describe the two-phase flow of water and air around the hull. Furthermore, a two degrees of freedom solver was used together with dynamic mesh techniques to describe the fluid-structure interaction. The simulations were performed with both fixed and free sinkage and trim to make careful comparisons of the software and with experimental data. The results from the fixed sinkage and trim simulations of the planing hull in FLUENT and STAR-CCM+ show a good consistency. However, there is a sig- nificant difference in the pressure resistance obtained from the two codes that could not be explained. The free sinkage and trim simulations were mainly conducted in STAR-CCM+ due to problems with obtaining a stable solution in FLUENT. Froude numbers between 0.447 and 1.79 were simulated and the results follow the same trends as what is seen in the experimental data. The calculated resistance, sinkage and trim angle show good correspondence to experimental data in the planing region, where the errors of the predicted values are below 10%. Keywords: Computational fluid dynamics (CFD), planing hull, resistance, sink- age, trim, fluid-structure interaction (FSI), ANSYS FLUENT, CD-adapco STAR- CCM+, Reynolds averaged Navier-Stokes equations (RANS), Volume of fluid (VOF). v Acknowledgements This thesis is a result of our Master’s thesis project which has been carried out at FS Dynamics in the spring 2015. We would like to thank our supervisors at FS Dynamics, Fredrik Carlsson, Julia Claesson, Ulf Engdar and Mattias Wångblad, for their continuous encouragement, support and feedback throughout the project. Also, we thank our supervisor at the Department of Shipping and Marine Technol- ogy of Chalmers, Arash Eslamdoost, for the very rewarding discussions of the details of hull simulations and for his valuable feedback. Many thanks to all the people at FS Dynamics for a warm welcome and for showing interest in our project. Finally, we want to express our gratitude to Björn Henriksson at Swede Ship Marine for providing us with data and showing interest in the results. David Frisk and Linda Tegehall Gothenburg, June 2015 vii Contents List of figures xiii List of tables xv Acronyms xvii Nomenclature xix 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Purpose and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Demarcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.5 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.6 Report structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Theory 5 2.1 Planing hull theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 Vessel resistance . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2 Hull induced waves . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.3 Planing hull characteristics . . . . . . . . . . . . . . . . . . . . 7 2.1.4 Towing tank tests . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.4.1 Model test scaling . . . . . . . . . . . . . . . . . . . 8 2.2 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.1 Turbulent flow . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.2 Turbulence modelling . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2.1 The standard k-ε model . . . . . . . . . . . . . . . . 11 2.2.2.2 The k-ω model . . . . . . . . . . . . . . . . . . . . . 11 2.2.2.3 The SST k-ω model . . . . . . . . . . . . . . . . . . 12 2.2.3 Boundary layers . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Free water surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.1 The volume of fluid method . . . . . . . . . . . . . . . . . . . 14 2.4 Fluid-structure interaction . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.1 Rigid body motion . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.2 Dynamic hull simulations . . . . . . . . . . . . . . . . . . . . 15 2.5 Non-dimensional coefficients . . . . . . . . . . . . . . . . . . . . . . . 16 ix Contents 3 Numerical methods 17 3.1 The finite volume method . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 Spatial discretization schemes . . . . . . . . . . . . . . . . . . . . . . 18 3.2.1 VOF discretization schemes . . . . . . . . . . . . . . . . . . . 18 3.2.1.1 HRIC scheme . . . . . . . . . . . . . . . . . . . . . . 18 3.2.1.2 Compressive scheme . . . . . . . . . . . . . . . . . . 19 3.3 Temporal discretization schemes . . . . . . . . . . . . . . . . . . . . . 19 3.4 Pressure-velocity coupling . . . . . . . . . . . . . . . . . . . . . . . . 19 3.5 Dynamic meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.5.1 Diffusion-based smoothing . . . . . . . . . . . . . . . . . . . . 20 3.5.2 Overset method . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.6 Convergence criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.7 Grid dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.7.1 The grid convergence index procedure . . . . . . . . . . . . . . 23 4 CFD simulations 25 4.1 Computational domain definition . . . . . . . . . . . . . . . . . . . . 25 4.2 Mesh generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.3 Model definitions and properties . . . . . . . . . . . . . . . . . . . . . 27 4.3.1 Mathematical models . . . . . . . . . . . . . . . . . . . . . . . 27 4.3.1.1 Two-phase flow models . . . . . . . . . . . . . . . . . 28 4.3.1.2 Hull motion . . . . . . . . . . . . . . . . . . . . . . . 28 4.3.2 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . 28 4.4 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.5 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.6 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.7 Simulated operating conditions . . . . . . . . . . . . . . . . . . . . . 30 4.7.1 Athena hull . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.7.2 Swede Ship hull . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5 Results and discussion 35 5.1 Athena hull . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.1.1 Grid dependence study . . . . . . . . . . . . . . . . . . . . . . 35 5.1.2 Calculated properties . . . . . . . . . . . . . . . . . . . . . . . 36 5.2 Swede Ship hull, fixed sinkage and trim . . . . . . . . . . . . . . . . . 38 5.2.1 Grid dependence study . . . . . . . . . . . . . . . . . . . . . . 38 5.2.2 Volume fraction of air at the hull . . . . . . . . . . . . . . . . 39 5.2.3 Pressure coefficient . . . . . . . . . . . . . . . . . . . . . . . . 40 5.2.4 Skin friction coefficient . . . . . . . . . . . . . . . . . . . . . . 41 5.2.5 Dimensionless wall distance . . . . . . . . . . . . . . . . . . . 42 5.2.6 Free surface wave pattern . . . . . . . . . . . . . . . . . . . . 43 5.2.7 Calculated properties . . . . . . . . . . . . . . . . . . . . . . . 44 5.3 Swede Ship hull, free sinkage and trim . . . . . . . . . . . . . . . . . 45 5.3.1 Grid dependence study . . . . . . . . . . . . . . . . . . . . . . 45 5.3.2 Centre of gravity sensitivity analysis . . . . . . . . . . . . . . 47 5.3.3 Calculated properties . . . . . . . . . . . . . . . . . . . . . . . 47 5.3.4 Free surface wave pattern . . . . . . . . . . . . . . . . . . . . 51 x Contents 5.3.5 FLUENT simulations of free sinkage and trim . . . . . . . . . 52 6 Conclusion 53 7 Future work 55 References 57 A Data analysis I B Athena hull data III C Swede Ship hull data V D Grid dependence studies VII xi Contents xii List of figures 1.1 Work flow of the CFD simulations. . . . . . . . . . . . . . . . . . . . 4 2.1 Characteristic parts of a hull. . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Characteristic hull measurements. . . . . . . . . . . . . . . . . . . . . 6 2.3 Kelvin wake pattern behind a moving object. . . . . . . . . . . . . . . 7 2.4 A planing hull in front view. . . . . . . . . . . . . . . . . . . . . . . . 8 2.5 Schematic illustration of a boundary layer at a flat plate. . . . . . . . 12 2.6 Coordinate system showing the 6 degrees of freedom of a rigid body. . 15 2.7 Iterative procedure of a hull simulation used to describe the fluid- structure interaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1 Schematic illustrations of dynamic meshes around a tilted square. . . 20 4.1 Dimensions of the computational domain. . . . . . . . . . . . . . . . 25 4.2 Schematic illustrations of the FLUENT mesh structure. . . . . . . . . 26 4.3 Schematic illustrations of the STAR-CCM+ mesh structure. . . . . . 27 4.4 Boundaries of the computational domain. . . . . . . . . . . . . . . . . 29 4.5 Sketch of the Athena hull showing the two different centres of gravity, CG1 and CG2, used in the simulations. . . . . . . . . . . . . . . . . . 31 4.6 Sketch of the Swede Ship hull showing the centre of gravity (CG) and the centre of buoyancy (CB) at zero speed. . . . . . . . . . . . . . . . 32 5.1 Convergence of the total resistance coefficient with grid refinement for the free sinkage and trim simulations of the Athena hull performed in STAR-CCM+ at Fn = 0.545. . . . . . . . . . . . . . . . . . . . . 36 5.2 Volume fraction of air at the symmetry plane of the Athena hull at Fn = 0.545. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.3 Cross sections at different y-coordinates used in the illustrations of the results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.4 Convergence of the total resistance coefficient with grid refinement for the fixed sinkage and trim simulations of the Swede Ship hull performed in FLUENT at Fn = 1.68. . . . . . . . . . . . . . . . . . . 39 5.5 Contour plots of the volume fraction of air on the hull for the fixed sinkage and trim simulations at Fn = 1.68. . . . . . . . . . . . . . . . 39 5.6 Contour plots of the pressure coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. . . . . . . . . . . . . . . . 40 5.7 Longitudinal cross sections of the pressure coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. . . . . . . . . . 40 5.8 Transverse cross sections of the pressure coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. . . . . . . . . . . . . 41 xiii List of figures 5.9 Contour plots of the skin friction coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. . . . . . . . . . . . . . . . 41 5.10 Longitudinal cross sections of the skin friction coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. . . . . . . . 42 5.11 Contour plots of the dimensionless wall distance on the hull for the fixed sinkage and trim simulations at Fn = 1.68. . . . . . . . . . . . . 42 5.12 Contour plots of the dimensionless free surface height for the fixed sinkage and trim simulations at Fn = 1.68. . . . . . . . . . . . . . . . 43 5.13 Cross sections of the dimensionless free surface height at Fn = 1.68. . 44 5.14 Convergence of the total resistance coefficient with grid refinement for the free sinkage and trim simulations of the Swede Ship hull performed in STAR-CCM+ at Fn = 1.68. . . . . . . . . . . . . . . . . . . . . . 45 5.15 Convergence of the sinkage with grid refinement for the free sinkage and trim simulations of the Swede Ship hull performed in STAR- CCM+ at Fn = 1.68. . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.16 Convergence of the trim angle with grid refinement for the free sinkage and trim simulations of the Swede Ship hull performed in STAR- CCM+ at Fn = 1.68. . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.17 Computational and experimental results for the Froude number de- pendency of the total resistance for the Swede Ship hull with free sinkage and trim. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.18 Computational and experimental results for the Froude number de- pendency of the sinkage for the Swede Ship hull with free sinkage and trim. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.19 Computational and experimental results for the Froude number de- pendency of the trim angle for the Swede Ship hull with free sinkage and trim. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.20 Computational and experimental results for the Froude number de- pendency of the wetted area for the Swede Ship hull with free sinkage and trim. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.21 Computational and experimental results for the Froude number de- pendency of the wetted length for the Swede Ship hull with free sink- age and trim. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.22 Contour plots of the dimensionless free surface height for the free sinkage and trim simulations obtained in STAR-CCM+ at Fn = 0.447, Fn = 0.894 and Fn = 1.68. . . . . . . . . . . . . . . . . . . . . 51 5.23 Dimensionless free surface height around the free hull at Fn = 1.68 showing the unphysical dislocation of the free surface when the hull moves downwards and upwards. . . . . . . . . . . . . . . . . . . . . . 52 xiv List of tables 4.1 Hull dimensions of R/V Athena. . . . . . . . . . . . . . . . . . . . . . 31 4.2 Fluid properties used in the simulations of the Athena hull. . . . . . . 31 4.3 Properties and dimensions of the Swede Ship hull at a water density of 1025.87 kg/m3 and a shell plating thickness of 8.00mm. . . . . . . 32 4.4 Fluid properties used in the simulations of the Swede Ship hull. . . . 32 4.5 Simulated cases for the Swede Ship hull with fixed sinkage and trim. . 33 4.6 Simulated cases for the Swede Ship hull with free sinkage and trim. . 33 5.1 Experimental and calculated properties of the Athena hull with free sinkage and trim at Fn = 0.545. . . . . . . . . . . . . . . . . . . . . . 36 5.2 Calculated properties of the Swede Ship hull with fixed sinkage at LPP/2 and trim angle of 0.685m and 3.19◦, respectively, at Fn = 1.68. 44 5.3 Calculated properties of the Swede Ship hull with free sinkage and trim at Fn = 1.68. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 B.1 Hull dimensions of R/V Athena. . . . . . . . . . . . . . . . . . . . . . III B.2 Fluid properties used in the simulations for the Athena hull. . . . . . III C.1 Dimensions of the Swede Ship hull. . . . . . . . . . . . . . . . . . . . V C.2 Properties of the Swede Ship hull at a water density of 1025.87 kg/m3 and a shell plating thickness of 8.00mm. . . . . . . . . . . . . . . . . V C.3 Fluid properties used in the simulations of the Swede Ship hull. . . . VI D.1 Number of cells in the meshes used in the grid dependence study of the Athena hull with free sinkage and trim at Fn = 0.545. . . . . . . VII D.2 Number of cells in the meshes used in the grid dependence study of the Swede Ship hull with fixed sinkage and trim at Fn = 1.68. . . . . VII D.3 Number of cells in the meshes used in the grid dependence study of the Swede Ship hull with free sinkage and trim at Fn = 1.68. . . . . . VII xv List of tables xvi Acronyms In the following table, expansions of the acronyms used in the report are presented. Acronym Expansion AP Aft perpendicular CAD Computer-aided design CFD Computational fluid dynamics DNS Direct numerical simulation DOF Degrees of freedom EFD Experimental fluid dynamics FP Forward perpendicular FSI Fluid-structure interaction FVM Finite volume method GCI Grid convergence index HRIC High resolution interface capturing scheme PDE Partial differential equation RANS Reynolds averaged Navier-Stokes SIMPLE Semi-implicit method for pressure linked equations SST Shear stress transport UDF User defined function VOF Volume of fluid xvii Acronyms xviii Nomenclature In the following tables, the nomenclature used in the report is explained. For the Greek and Roman symbols, the physical dimensions are expressed in the base di- mensions mass (M), length (L) and time (T). When a variable is not connected to a certain physical quantity, the dimension field is marked with a star (?). Accents Accent Description Time average −→ Vector Greek symbols Symbol Description Dimensions α Diffusion parameter 1 δij Kronecker delta 1 δRE Discretization error estimation ? ε Dissipation of turbulent kinetic energy L2T−3 Γ Diffusivity L2T−1 γ Colour function used in the VOF method 1 Γm Mesh diffusivity L2T−1 ∇ Volume displacement L3 ν Kinematic viscosity L2T−1 νt Turbulent viscosity L2T−1 Ω Angular velocity T−1 ω Specific dissipation of turbulent kinetic energy T−1 φ General flow variable ? φnb Value of general flow variable in neigh- bouring cells ? xix Nomenclature φP Value of general flow variable in present cell ? ρ Density ML−3 σ Surface tension MT−2 τ Torque ML2T−2 τw Wall shear stress ML−1T−2 θ Trim angle 1 Roman symbols Symbol Description Dimensions anb Discretization coefficients associated to neighbouring cells ? aP Discretization coefficient associated to present cell ? AW Wetted area L2 Cε1, Cε2 Cµ, σk, σε Constants of the k-ε turbulence model 1 Cω1, Cω2 β∗, σωk , σω Constants of the k-ω turbulence model 1 d Normalized boundary distance 1 ED Discretization uncertainty ? EI Iterative uncertainty ? EN Numerical uncertainty ? F Force MLT−2 FS Factor of safety used in GCI 1 g Gravitational acceleration LT−2 hi Characteristic cell length of grid i 1 k Turbulent kinetic energy L2T−2 L Length L LOA Length overall L LPP Length between perpendiculars L LW Wetted length L M Moments of inertia (tensor) ML2 m Mass M xx Nomenclature N Number of cells 1 P Pressure ML−1T−2 p Fluctuating pressure ML−1T−2 P∞ Pressure of free stream ML−1T−2 Pk Production of turbulent kinetic energy L2T−3 q Richardson extrapolation constant ? r Richardson extrapolation observed order of accuracy 1 Rnorm Residual normalization value ? Rφ Residual of flow variable φ ? Rφ,s Residual of flow variable φ, globally scaled ? Rpres Present value of residual ? rr Grid refinement ratio 1 Rrms Root mean square residual ? Rtot Total resistance MLT−2 S Source term in transport equation ? s Sinkage L S0 Extrapolated solution for an infinite grid density ? Si Solution on grid i ? SP Variable depending part of source term ? SU Constant part of source term ? T Draught L t Time T U Velocity LT−1 u Fluctuating velocity LT−1 Uhull Hull velocity LT−1 u∗ Friction velocity LT−1 xi Spatial coordinate in dimension i L y+ Dimensionless wall distance 1 xxi Nomenclature Dimensionless numbers Symbol Name Definition Cf Skin friction coefficient Cf = τw 1 2ρU 2 hull CP Pressure coefficient CP = P − P∞ 1 2ρU 2 hull CT Total resistance coefficient CT = Rtot 1 2ρU 2 hullAW0 Fn Froude number Fn = Uhull√ Lg Rn Reynolds number Rn = UhullL ν Wn Weber number Wn = ρU2 hullL σ Subscripts Subscript Description AP Aft perpendicular FP Forward perpendicular xxii 1 Introduction Hull design is a field of engineering that has been developed for hundreds of years. Today, the main focus is to meet the contract speed requirements with minimum fuel consumption. It is important to make an accurate prediction of the hull resistance to be able to choose an appropriate engine power of the propulsion unit. It is also of interest to investigate the running attitude of the hull to get desirable seakeeping properties. Traditionally, the hull resistance and running attitude have been determined by performing experiments on down scaled model hulls in towing tanks. These experi- ments have proved to predict the behaviour of the full scale hull very well, but the method is time consuming, expensive and only applicable to the tested model. A more universal method for hull performance predictions is to use computational fluid dynamics (CFD) which is a branch of computer simulations where the behaviour of a system involving fluid flow is analysed using numerical methods. There are well established methods for CFD simulations of displacement hulls and the results predict the behaviour of the hull with high accuracy. Resistance pre- dictions of planing hulls are more difficult, and a lot of today’s research is focused on improving CFD methods for planing hulls. In contrast to displacement hulls, which are supported by hydrostatic forces acting on the hull, planing hulls utilize hydrodynamic forces from the water to reduce the resistance and thereby be able to reach high speeds with relatively low fuel consumption. This study is done in cooperation with Swede Ship Marine AB which is a company that designs and manufactures high speed planing vessels such as rescue vessels, pilot boats, coast guard vessels and navy vessels. If Swede Ship could make use of CFD in their hull development, they could improve their vessel design process and reduce the production time and cost. 1.1 Background Hull development has traditionally relied on model scale experiments in towing tanks. Instead of doing these experiments, empirical models based on regression analysis of experimental data have been developed to predict the behaviour based on the characteristics of the hull shape. Two widely used empirical models are that of Holtrop and Mennen [1], used for displacement hulls, and that of Savitsky [2], 1 1. Introduction used for planing hulls. The main drawback with these empirical models is that they are restricted to simple hull shapes which are similar to the hulls used to develop the models. During the last decades, several numerical studies on planing hulls have been per- formed. Brizzolara and Serra [3] compared CFD simulations of planing hulls with the Savitsky method. The flow was modelled with the RANS equations and the standard k-ε model was used to model the turbulence. The VOF method was used to resolve the free surface. It was found that the hull resistance from the CFD simulations differed in average 10% from the experimental results. This was lower than the error of the Savitsky method, and the study showed the potential of CFD simulations. In order to find the running attitude of a vessel, its motion must be included in the simulation. For this purpose, methods for coupling the fluid flow and body motions have been developed. Azcueta [4] implemented a method in 2001 where the inter- action between the hull and the fluid was modelled. The turbulence was modelled with the RANS equations and the standard k-ε model, and the free surface was cap- tured with a VOF method. The method was validated for a Series 60 Hull, which is a well-known displacement hull used for simulation benchmarking, and the results were compared with experimental data. The total resistance was underpredicted by 5.9% and the sinkage and trim were 8.2% and 6.0% lower than the experimental results respectively. Recent studies show that CFD simulations of displacement hulls have been made with a precision that start to approach that of towing tank tests. At the Workshop on Numerical Ship Hydrodynamics in Gothenburg 2010 [5, pp. 1–16], 33 partici- pating groups performed simulations of three large displacement ships. The results from all simulations show that the mean error of the resistance predictions, in com- parison to towing tank experiments, was only 0.1% with a standard deviation of 2.1%. The sinkage and trim predictions were less accurate, for the higher speeds the mean errors were around 4%. CFD simulation of planing hulls is more challenging in comparison to displacement hulls, and in general the accuracy of the predictions is lower. The viscous forces from the flow around a planing hull are strongly dependent on the wet surface of the hull, which in turn depends on the position in the water. It is therefore crucial to predict how a planing hull behaves in water before adequate resistance estimations can be made. [6] Also, in the higher speed range, nonlinear effects such as breaking waves and water sprays become more important. However, as the computer power steadily increases and better models are developed, the results from numerical simulations of planing hulls are continuously improved. [7] The results from recent studies [8, 9] have shown that CFD simulations of planing hulls can yield results that deviate from experimental data by less than 10%. 2 1. Introduction 1.2 Purpose and objectives The purpose of this study is to evaluate how accurate CFD is using commercially available software, and if these can fully or partly replace towing tank tests in planing hull design. This is done by calculating the resistance, sinkage and trim angle of a high speed planing hull, whereupon the results are compared with experimental data. In this study, the simulations have been performed using two of the most common general purpose CFD software, ANSYS FLUENT, developed by ANSYS, Inc., and STAR-CCM+, developed by CD-adapco. 1.3 Demarcations The focus of this study is on the stationary behaviour of planing hulls moving in calm water at a constant speed. The vessel is limited to heave and pitch motions. An economical evaluation of the method is not included. 1.4 Research questions The following questions are addressed in this report. • How accurately can the resistance and running attitude of a planing hull be predicted using CFD simulations? • Do the CFD simulations yield any additional information that is not obtained from model testing? • Are there any differences in the results obtained from FLUENT and STAR-CCM+? 1.5 Method A flow chart outlining the methodology of the CFD simulations is presented in Figure 1.1. The strategies differ in FLUENT and STAR-CCM+, and are based on guidelines for hull simulations provided by each software developer. The work starts with defining the computational domain required to perform the analysis, and then the computational meshes are generated. When appropriate models have been chosen in the model definitions step and the boundary conditions have been set, the solution can be initiated. The computational domain is the same for FLUENT and STAR-CCM+, while the mesh generation, model definitions, boundary conditions and solution steps are done separately. When the results are obtained, they are analysed in the post processing step. These steps are made in an iterative process where the results are evaluated and further simulations are performed. 3 1. Introduction Start Computational domain definition Mesh generation Model definitions Boundary conditions Solution Post processing End FLUENT STAR-CCM+ Figure 1.1: Work flow of the CFD simulations. 1.6 Report structure The report is organized in the following way. In Chapter 2, the theory behind planing hulls and the different types of resistance forces that are acting on the hull are explained. In the same chapter the mathematical models governing the fluids, the fluid interface and the vessel motion are presented. The numerical methods used for treating these equations are described in Chapter 3. Chapter 4 describes how the CFD simulations were performed and in Chapter 5, the results are presented and discussed. Finally, the conclusions are presented in Chapter 6 and recommendations for future work are given in Chapter 7. 4 2 Theory This chapter aims to provide the theoretical basis for understanding the remainder of the report. It starts with a section introducing the basics of naval architecture related to planing hulls, and continues with describing the mathematical models used in the simulation of a hull in motion. 2.1 Planing hull theory In Figure 2.1, some of the characteristic parts of a planing hull are highlighted. The bow and stern are the front and back of the vessel and the transom is the vertical surface located at the stern. The spray rails and lifting rails are used to guide the water flow along the hull in order to obtain the desired behaviour of the vessel when moving through the water. Lifting rail Spray rail Transom Bow Stern Figure 2.1: Characteristic parts of a hull. In naval architecture, the properties of a hull can be characterized by certain mea- surements shown in Figure 2.2. The measurements used are the overall length, LOA, length between perpendiculars, LPP , and the draughts at the forward and aft per- pendiculars, TFP and TAP . The aft perpendicular (AP) is a point on the vertical transom, and the forward perpendicular (FP) is the point where the waterline in- tersects the keel in the front. When the vessel is in motion, its position can be related to that of zero speed by measuring the change in vertical position of these perpendiculars. This change is known as the sinkage, s, which is defined as positive if the elevation of the perpen- dicular has increased. The trim angle, θ, is the angle of rotation around the y-axis 5 2. Theory in comparison to the position at zero speed, and a positive trim angle means that the bow has moved up in relation to the aft. Moreover, LW and AW denote the wetted length and the wetted area of the hull. The volume of the water displaced by the hull is denoted ∇. Waterline LP P LOA Base lineTF PTAP x z ⊗ y Figure 2.2: Characteristic hull measurements. 2.1.1 Vessel resistance The resistance of a vessel moving through calm, unrestricted water can be decom- posed into components. One common way is to divide the resistance into a pressure resistance and a friction resistance. Another frequently used decomposition is to use viscous resistance and wave resistance. [10, pp. 13–16] Viscous resistance arises due to shear forces when the vessel moves through the surrounding water and air. The contribution to the viscous resistance from the air is strongly dependent on the aerodynamic properties of the vessel, but it normally constitutes a minor part of the total resistance. Wave resistance is due to generation of water waves and water spray. A more intense wave generation and spray means that more energy is transferred from the hull to the water and thereby the resistance increases. The total resistance coefficient is a dimensionless quantity defined as CT = Rtot 1 2ρU 2 hullAW0 , (2.1) where Rtot is the total resistance, ρ is the density of the water, Uhull is the hull speed and AW0 is the wetted area at rest. This coefficient is used to characterize the total resistance and to compare the performance of different hulls. 2.1.2 Hull induced waves A small object moving through calm water will induce a wake pattern, known as the Kelvin wake pattern, which is illustrated in Figure 2.3. The intensity of the waves reflects the amount of energy that is continuously transferred from the object to the water. It can be shown that the angle between the object’s trajectory and 6 2. Theory the wave fronts, the Kelvin angle, is around 19.47◦ regardless of the speed of the object. The wake pattern behind a vessel can be approximated with the Kelvin wake pattern, even though a hull has many wave-generating features causing overlapping and interfering wave systems. [10, pp. 24–30] Moving object 19.47◦ Wave crest Wave trough Figure 2.3: Kelvin wake pattern behind a moving object. [10, pp. 24–30] 2.1.3 Planing hull characteristics At zero speed, the weight of a vessel is supported only by hydrostatic forces acting on the hull. If the speed increases, hydrodynamic forces will arise due to water that is accelerated away from the hull causing a reactionary force. When a hull is said to be planing, the hydrodynamic forces play a major role in supporting the weight of the vessel while the hydrostatic forces are of little importance. This requires that the vessel reaches a certain speed where the amount of water that is accelerated downwards by the hull is large enough. [11, pp. 342–389] A common measure used to characterize the resistance and wave pattern around a hull is the dimensionless Froude number, which relates inertial forces to external forces. It is defined as Fn = Uhull√ Lg , (2.2) where L is a characteristic hull length and g is the gravitational acceleration. In this report, all Froude numbers are based on the length between perpendiculars at zero speed. A common value of the Froude number which marks the transition region from non-planing to planing conditions is Fn = 1.0 [11, p. 342], but other values are used in the literature since the phenomenon has no strict definition. Planing hulls are designed to enhance the hydrodynamic forces acting in the ver- tical direction in order to reach planing. Theoretically, a flat underside provides a larger hydrodynamic lift compared to a so called v-shaped underside. However, planing hulls usually have v-shaped undersides with purpose to reduce the slamming in rough seas. To prevent the flow along the hull from bending outwards due to the 7 2. Theory v-shape, lifting rails are used to guide the flow backwards and thereby increase the hydrodynamic pressure under the hull. Spray rails can be used on the hull sides to increase the lift force by bending the water spray downwards. [10, pp. 175–181] Furthermore, it is important that the water separates from the hull at the transom in order to avoid sub-atmospheric pressures that could lead to instabilities. [12] This can be achieved by using a sharp edge at the transom where the water is supposed to detach from the hull. The same hull as in Figure 2.1 and 2.2 is shown in a front view in Figure 2.4. It can be seen that it has a v-shaped underside and is equipped with lifting rails and spray rails which are characteristic features of a planing hull. y z � x Figure 2.4: A planing hull in front view. 2.1.4 Towing tank tests When designing a vessel, the properties of the hull are traditionally tested at an early stage by constructing a smaller model which can be used in experiments. These experiments can be performed in open water, but to reduce the sources of errors, they are often carried out in large tanks known as towing tanks. In a towing tank, the hull is towed through the water while measurements are conducted. 2.1.4.1 Model test scaling When towing tank tests have been carried out on a model, the results must be scaled properly in order to apply to a hull in full scale. To achieve adequate scaling results, the model and the full scale hull must be geometrically similar and the relevant dimensionless numbers should be preserved. For a planing hull, the dimensionless numbers characterising the physics of the flow are the Froude number, the Reynolds number and the Weber number. The Froude number, as defined in equation (2.2), describes the effect of gravity on the water surface. The Reynolds number relates inertial forces to viscous forces, and is defined as Rn = UhullL ν , (2.3) 8 2. Theory where ν is the kinematic viscosity. The Weber number relates inertial forces to forces due to surface tension, and is defined as Wn = ρU2 hullL σ , (2.4) where σ is the surface tension. A dilemma related to scaling is the inability to satisfy the conservation of these three dimensionless numbers simultaneously. As can be seen from the definitions of the three dimensionless numbers, a model scale hull with shorter length must have a lower speed to preserve the Froude number but a higher speed to preserve the Reynolds and Weber numbers. The established solution to this problem is to use the same Froude number in the model tests and to account for the different Reynolds numbers in the scaling process. The difference in Weber numbers between model scale and full scale causes errors in the predictions of water spray, foaming and wave pattern, but it has been concluded that this has little influence on the resistance prediction of the full scale hull. [10, pp. 10–13] 2.2 Turbulence When a hull is moving through water at high speed, the flow around the hull is turbulent. In this section, the governing equations of turbulent flows are presented and turbulence modelling is explained. 2.2.1 Turbulent flow Turbulence has no physical definition, but it is characterized as a three-dimensional, irregular flow where turbulent kinetic energy is dissipated from the largest to the smallest turbulent scales. On the smallest turbulent scales, known as the Kol- mogorov scales, the energy is dissipated into heat due to viscous forces. Since turbulence is a dissipative phenomenon, energy must be continuously supplied in order to maintain a turbulent flow. The motion of a viscous fluid is governed by the Navier-Stokes equations, which are valid both for turbulent and laminar flow. For an incompressible, Newtonian fluid in three dimensions under the influence of an external gravitational field, the Navier-Stokes equations read ∂Ui ∂xi = 0, ∂Ui ∂t + Uj ∂Ui ∂xj = −1 ρ ∂P ∂xi + ν ∂2Ui ∂xj∂xj + gi. (2.5) Here, the equations are formulated using tensor notation. The indices i and j in the Navier-Stokes equations run over the spatial coordinates x, y and z. In these equations, Ui is the velocity in dimension i, xi is the spatial coordinate in dimension 9 2. Theory i, t is time, ρ is the density, P is the pressure, ν is the kinematic viscosity and gi is the gravitational acceleration. Analytical solutions to the Navier-Stokes equations only exist for a limited num- ber of simple cases such as laminar flow between flat plates. For turbulent flows in engineering applications, analytical solutions do not exist and the Navier-Stokes equations must be treated numerically. If they are solved using direct numerical simulation (DNS), the velocity field of the flow is obtained. However, since tur- bulence occurs on a wide range of time and length scales, DNS requires very high temporal and spatial resolutions to capture all the details of the flow. Thus, DNS is very computationally expensive and time consuming which limits the method to special applications such as academic research or simulation of simple flows. 2.2.2 Turbulence modelling The most common way of treating turbulence is to use turbulence models in which the turbulent features of the flow are not resolved in time. By performing Reynolds decomposition, the instantaneous velocity and pressure can be decomposed as Ui = Ui + ui, P = P + p, (2.6) where Ui and P denote the time averaged quantities while ui and p are the fluc- tuating components of the velocities and the pressure. By inserting the Reynolds decomposition into the Navier-Stokes equations given in equation 2.5, the Reynolds averaged Navier-Stokes (RANS) equations are obtained. These are written as ∂Ui ∂xi = 0, ∂Ui ∂t + Uj ∂Ui ∂xj = −1 ρ ∂P ∂xi + ν ∂2Ui ∂xj∂xj − ∂uiuj ∂xj + gi. (2.7) It can be noted that the RANS equations are very similar to the Navier-Stokes equations except for the additional term including uiuj, referred to as the Reynolds stress tensor. If the Reynolds stress term is modelled, the RANS equations describe the time-averaged flow quantities which requires substantially less computational resources in comparison to DNS. A common approach for modelling the Reynolds stress tensor of the RANS equa- tions is to use the Boussinesq approximation. In this assumption, the Reynolds stress tensor is modelled as a diffusion term by introducing a turbulent viscosity, νt, according to − uiuj = νt ( ∂Ui ∂xj + ∂Uj ∂xi ) − 2 3kδij. (2.8) In this equation, δij is the Kronecker delta which assumes a value of 1 if i = j and 0 otherwise, and k is the turbulent kinetic energy defined as k = 1 2uiui. (2.9) 10 2. Theory By using a model to describe how the turbulent viscosity depends on the flow, the RANS equations can be solved. The so called two-equation turbulence models, such as the k-ε model and the k-ω model, use two additional transport equations to describe the turbulent viscosity. They are referred to as complete models since they allow the turbulent velocity and length scales to be described independently. [13] 2.2.2.1 The standard k-ε model In the standard k-ε model, the transport equations for the turbulent kinetic en- ergy and its dissipation, ε, are used to obtain the turbulent viscosity. It has been described by Launder et al. [14], and the model equations are ∂k ∂t + ∂ ∂xi ( kUi ) = ∂ ∂xj [( ν + νt σk ) ∂k ∂xj ] + Pk − ε, ∂ε ∂t + ∂ ∂xi ( εUi ) = ∂ ∂xj [( ν + νt σε ) ∂ε ∂xj ] + ε k (Cε1Pk − Cε2ε) , νt = Cµ k2 ε , (2.10) where σk, σε, Cε1, Cε2 and Cµ are model constants and Pk is the production of turbulent kinetic energy. The latter is defined as Pk = −uiuj ∂Ui ∂xj (2.11) and is modelled using the Boussinesq approximation. The standard k-ε model is robust and gives good predictions for free flows with small pressure gradients. It is based on the assumption that the flow is fully turbulent which limits its applicability to high Reynolds number flows. [15] Over time, it has been observed that the standard k-ε model cannot be used to describe the wake behind a moving hull in a satisfactory manner. [10, p. 135] 2.2.2.2 The k-ω model In the k-ω model described by Wilcox [15], the transport equations for the turbulent kinetic energy and its specific dissipation, ω, are used in a similar way as for the standard k-ε model. The specific dissipation is related to the dissipation according to ω ∝ ε k . (2.12) The model equations for k and ω are ∂k ∂t + ∂ ∂xi ( kUi ) = ∂ ∂xj [( ν + νt σωk ) ∂k ∂xj ] + Pk − β∗kω, ∂ω ∂t + ∂ ∂xi ( ωUi ) = ∂ ∂xj [( ν + νt σω ) ∂ω ∂xj ] + ω k (Cω1Pk − Cω2kω) , νt = k ω , (2.13) 11 2. Theory where β∗, σωk , σω, Cω1 and Cω2 are model constants. The k-ω model has the advantage that it is also valid close to walls and in regions of low turbulence. [16, pp. 121–122] Thus, it is valid also in the low turbulent Reynolds number region close to walls, meaning that the transport equations can be used in the whole flow domain. A drawback with the k-ω model is that the results are sensitive to the choice of boundary conditions and initial conditions. 2.2.2.3 The SST k-ω model In order to make use of the collective advantages of the k-ε and k-ω models, Menter [17] developed the shear stress transport (SST) model by combining the two models into one using blending functions. In this hybrid model, the k-ω model is used in the boundary layer while the k-ε model, formulated on k-ω form, is used in the free flow. The SST k-ω model has shown good performance for many types of complex flows, such as in flows with adverse pressure gradients and separating flows, where the k-ε or the k-ω models have given results that differ significantly from experimental data. It has been recognized for its good overall performance [18] and it is the most commonly used turbulence model for simulations of ship hydrodynamics. [19, p. 7] 2.2.3 Boundary layers When a fluid flows along a surface, shear stresses give rise to a boundary layer in the vicinity of the surface. The structure of a boundary layer near the edge of a flat plate is illustrated in Figure 2.5, where the incident flow has a uniform velocity profile with velocity U0. When the flow reaches the plate, a laminar boundary layer starts to grow at the surface. After some distance, the boundary layer goes into a transition region after which a turbulent boundary layer is developed. The flow in the inner part of the turbulent boundary layer is laminar, and the turbulence increases further away from the wall. [20, pp. 464–475] U0 U(y) Laminar boundary layer Transition region Turbulent boundary layer Viscous sub-layer Buffer sub-layer Fully turbulent sub-layer y Figure 2.5: Schematic illustration of a boundary layer at a flat plate. [20, pp. 464–475] 12 2. Theory To characterize the flow near a wall, a dimensionless wall distance is often intro- duced. It is defined as y+ = u∗y ν , (2.14) where y is the distance to the wall and u∗ is a friction velocity. The friction velocity is defined as u∗ = √ τw ρ , (2.15) where τw is the wall shear stress, τw = ρν ∂U ∂y ∣∣∣∣∣ y=0 . (2.16) In the boundary layer, the gradients of the flow variables in the wall-normal direc- tion are generally very large in comparison to those of the free flow. This implies that a high spatial resolution is required by the solution method in order to capture the effects near the wall. A common alternative method used to circumvent the requirement of a high spatial resolution is to use wall functions, which are empirical models used to estimate the flow variables near walls. Wall functions can also be applied when the turbulence model used in a simulation is not valid close to the wall, which for example is the case for the standard k-ε model. [16, pp. 128–140] Although wall functions are undesired in computational ship hydrodynamics due to deviations for some types of flow, they are often used for numerical reasons. [10, pp. 152–154] Standard wall functions are based on the assumption that the boundary layer can be described as a flat plate boundary layer. This means that the time-averaged velocity can be expressed as a function of the dimensionless wall distance. In the viscous sub-layer, it can be shown that the velocity parallel to the wall is proportional to y+. In the fully turbulent sub-layer, the velocity follows the logarithmic law of the wall, meaning that the velocity is proportional to the natural logarithm of y+. Between these sub-layers, in the buffer sub-layer, there is a transition from linear to logarithmic y+ dependence. In order for a wall function to work properly, it should be used all the way to the fully turbulent layer which corresponds to a value of y+ above 30. The wall functions also estimate the turbulence quantities near the walls. [16, pp. 128–140] 2.3 Free water surface In order to simulate a hull moving in water, models are needed to resolve the interface between the water and air. There are different two-phase models available that either tracks the surface directly or tracks the different phases and then reconstructs the interface. One example is the level-set method, where all molecules of one phase are marked and then tracked in the fluid flow. The most frequently used method to capture the free surface in ship hydrodynamics is the volume of fluid (VOF) method [5, pp. 5–9]. In the VOF method, the different phases are tracked. [10, pp. 151–152] 13 2. Theory 2.3.1 The volume of fluid method In the VOF method, each phase is marked with a colour function, γ, which is the volume fraction of one of the phases. If only one phase is present, meaning that γ is either 0 or 1, the ordinary Navier-Stokes equations, see equation (2.5), are solved. If 0 < γ < 1, there is an interface present and the properties of the phases are averaged in order to get a single set of equations. The average density and viscosity are ρ = γρ1 + (1− γ)ρ2, ν = γν1 + (1− γ)ν2. (2.17) Then, a modified set of the Navier-Stokes equations can be used for the averaged fluid properties, ∂Ui ∂xi = 0, ∂Ui ∂t + Uj ∂Ui ∂xj = −1 ρ ∂P ∂xi + ν ∂2Ui ∂xj∂xj + gi + Si,s. (2.18) This formulation of the Navier-Stokes equations contains an additional source term, Si,s, accounting for the momentum exchange across the interface due to surface tension forces. This surface tension force has to be modelled correctly which can be a problem. The surface is captured by solving a transport equation for the colour function, ∂γ ∂t + Ui ∂γ ∂xi = 0. (2.19) As mentioned in Section 2.1.1, the free surface waves affect the forces on the hull. Therefore, it is important to get an accurate and stable solution of equation (2.19). The VOF method is conceptually simple and relatively accurate but the solution techniques may be diffusive. It is not as robust and accurate as the level-set method, but it is less computationally demanding. [21] 2.4 Fluid-structure interaction In order to simulate the dynamic behaviour of a hull before reaching the equilibrium position, the fluid-structure interaction (FSI), between the hull and the fluids has to be taken into account. This is done by solving the equations of motion and rotation of the vessel under the influence of the forces and moments from the surrounding fluids and gravity. The number of directions the body is allowed to translate and rotate in is called the number of degrees of freedom (DOF). 2.4.1 Rigid body motion A vessel can be approximated as a rigid body which can move in three dimensions and rotate around three axes, see Figure 2.6. The translations of a vessel along the x, y and z axes are often referred to as surge, sway and heave motions, respectively, 14 2. Theory while the rotations around the same axes are termed roll, pitch and yaw motions. Accordingly, the sinkage is affected by heave motion while the trim angle is related to pitch motion. x θx yθy z θz Figure 2.6: Coordinate system showing the 6 degrees of freedom of a rigid body. These consist of translation along and rotation around three axes in the x, y and z directions. For a rigid body, the translational motion of the centre of gravity is described by Newton’s second law, m d ~Ub dt = ~F , (2.20) where m is the mass, ~Ub is the velocity and ~F is the sum of forces acting on the body. The rotation of the body, expressed in body coordinates, is described by Euler’s equations, Md~Ω dt + ~Ω× (M · ~Ω) = ~τ , (2.21) where ~Ω is the angular velocity of the body and ~τ is the resultant torque acting on the body. Furthermore,M is a tensor of the moments of inertia and it is expanded into M = Mxx Mxy Mxz Myx Myy Myz Mzx Mzy Mzz  . (2.22) 2.4.2 Dynamic hull simulations Under most circumstances, a hull moving with constant speed will reach a steady position and orientation with respect to the free surface. [10, pp. 152–154] When such an equilibrium position is expected in a simulation, a 6 DOF solver can be implemented in the solution process as shown in Figure 2.7. When the motions and rotations have ended and the final position is reached, the net forces and moments acting on the hull are zero. 15 2. Theory Initial position Solve equa- tions of flow Solve equations of motion (6 DOF) Rigid body motion converged? Final position No Yes Figure 2.7: Iterative procedure of a hull simulation used to describe the fluid- structure interaction. 2.5 Non-dimensional coefficients For a more relevant comparison of the pressure and shear stress on a surface, these quantities can be scaled by the dynamic pressure of the free flow to obtain dimen- sionless numbers. The pressure coefficient is defined as CP = P − P∞ 1 2ρU 2 hull , (2.23) where P∞ is the undisturbed free stream pressure. Similarly, the skin friction coef- ficient is defined as Cf = τw 1 2ρU 2 hull , (2.24) where τw is the wall shear stress. 16 3 Numerical methods In this chapter, the numerical methods used for treating the mathematical models introduced in Chapter 2 are described. In some cases, FLUENT and STAR-CCM+ use different formulations and in these cases both methods are presented. 3.1 The finite volume method The finite volume method (FVM), is a numerical method of discretizing a continu- ous partial differential equation (PDE), into a set of algebraic equations. The first step of the discretization is to divide the computational domain into a finite number of volumes, forming what is called a mesh or a grid. Next, the PDE is integrated in each volume by using the divergence theorem, yielding an algebraic equation for each cell. In the centres of the cells, cell-averaged values of the flow variables are stored in so called nodes. This implies that the spatial resolution of the solution is limited by the cell size since the flow variables do not vary inside a cell. [22, pp. 115–118] The FVM is conservative, meaning that the flux leaving a cell through one of its boundaries is equal to the flux entering the adjacent cell through the same boundary. This property makes it advantageous for problems in fluid dynamics. [16, pp. 32–33] A stationary transport equation involving diffusion and convection of a general flow variable, φ, can be written as ρUi ∂φ ∂xi = ∂ ∂xi ( Γ ∂φ ∂xi ) + S(φ), (3.1) where Γ is the diffusivity and S is a source term which may depend on φ. It can be noted that the equations in Chapter 2 governing the transport of Ui, k, ε, ω and γ are all written on this form. By using the FVM, this equation can be written on discrete form as aPφP = ∑ nb anbφnb + SU , (3.2) where aP = ∑ nb anb − SP . (3.3) In these equations, where the summations run over all the nearest neighbours of each cell, φP is the value of the flow variable in the present cell and φnb are the values of the flow variable in the neighbouring cells. SU and SP are the constant and flow variable depending parts of the source term, respectively. Furthermore, aP 17 3. Numerical methods is the discretization coefficient associated to the present cell, anb are discretization coefficients describing the interaction with its neighbouring cells. The discretization coefficients depend on the discretization schemes used to approximate the values of the flow variables on the cell boundaries, also known as cell faces. By using appropriate discretization schemes to determine the coefficients of equation (3.2), a set of algebraic equations for the cell values is obtained. [22, pp. 115–118] 3.2 Spatial discretization schemes The convection and diffusion terms in equation (3.1) are discretized using different numerical schemes that estimate the face values of the flow variables. Most often, diffusion terms are discretized by using a central differencing scheme where the face values are calculated by interpolation between the closest cells. In order to dis- cretize the convection terms, the flow direction has to be taken into account. The simplest way is to let the face value between two cells be equal to the value of the first upstream cell which is done in the first order upwind scheme. In the second order upwind scheme, the face value is calculated from the two closest upwind cells. [22, pp. 134–178] It is usually recommended to start a numerical solution process with lower order schemes, such as the first order upwind scheme, since they are very stable. However, the low accuracy of these schemes can lead to a high degree of unphysical diffusion in the solution [23], known as numerical diffusion. When the flow field has started to settle, higher order schemes should therefore be used to obtain a more physically correct result. The second order upwind scheme is often considered as a suitable discretization scheme since it exhibits a good balance between numerical accuracy and stability. [24] 3.2.1 VOF discretization schemes The main problem related to the VOF method is to discretize the convection term in the transport equation for the colour function in order to get a sharp interface. The scheme has to be accurate and at the same time bounded because the colour function, γ, has to be between 0 and 1. Lower order numerical schemes are bounded but will smear out the interface due to numerical diffusion while higher order schemes are more accurate but less stable. A combination of higher and lower order schemes is often used like in HRIC and the Compressive schemes used in STAR-CCM+ and FLUENT respectively. 3.2.1.1 HRIC scheme The high resolution interface capturing scheme (HRIC), described by Muzaferija et al. [25], uses a combination of upwind and downwind interpolation. The blending of the schemes in each cell is a function of the volume fraction distribution over the neighbouring cells. The value of the flow variable is then corrected by the local value of the Courant number which is a measure of how much of one fluid that is 18 3. Numerical methods available in the donor cell. This is done in order to prevent that no more fluid flows out of a cell in one time step than what was available in the previous time step. In order to prevent that the interface becomes aligned with the numerical grid, another correction is introduced to account for the relative position of the free surface to the cell face. This is done by calculating the angle between the normal to the interface and the cell face normal. 3.2.1.2 Compressive scheme The Compressive scheme is a discretization scheme where the numerical order of accuracy can be varied by using a so called slope limiter in the range between 0 and 2. For low values of the slope limiter, first and second order schemes are used. For values above 1, higher order schemes are incorporated. [26, pp. 467–468] 3.3 Temporal discretization schemes For transient problems, the transport equation must also be discretized in time. This is done by integrating the PDE over a time step ∆t in addition to the spatial discretization. In order to solve this integrated equation, the cell values of the flow variables must be evaluated at a certain time. Implicit time integration means that the flow variables are evaluated at the future time, t + ∆t. Since these are not known in the current time step, implicit time integration requires iteration. In comparison to explicit time integration, where the flow variables are evaluated at the current time so that iteration is avoided, implicit time integration is more computationally expensive. On the other hand, implicit time integration is unconditionally stable, meaning that it is stable for all time step sizes. [22, pp. 243–248] 3.4 Pressure-velocity coupling The Navier-Stokes equations, as written in equation (2.5), contain one continuity equation and three momentum equations, if a three dimensional system is consid- ered. There are four unknown variables in these equations, the pressure and the three velocity components. The problem is that there is no equation for the pressure, so the continuity equation must be used as an indirect equation for the pressure. This is achieved by using a pressure-velocity coupling, which can be either segregated or coupled. The properties of these two groups of algorithms will be described briefly, a more thorough explanation has been given by Versteeg and Malalasekera [22, pp. 179–211]. The semi-implicit method for pressure linked equations (SIMPLE) is a segregated algorithm that solves each equation separately. First, a pressure is assumed and the velocities are calculated from the momentum equations. If the continuity equation is not satisfied by these velocities, the pressure is modified and the velocities are 19 3. Numerical methods calculated again. With a coupled algorithm, the Navier-Stokes equations are solved together. Coupled algorithms are suitable when the computational mesh has a poor quality, and they allow larger time steps than segregated algorithms. 3.5 Dynamic meshing To be able to handle motion, the mesh structure has to change dynamically with the moving object. There are different methods for the dynamic movement of the mesh and it is done differently in FLUENT and STAR-CCM+. The two that are most suitable for hull simulations are the diffusion-based smoothing method in FLUENT and the overset mesh with mesh rotation and translation in STAR-CCM+. In Figure 3.1, the diffusion-based smoothing method and the overset method are illustrated around a tilted square in a two-dimensional domain. The concepts of the two methods are described in the following two sections. (a) Diffusion-based smoothing. (b) Overset method. Figure 3.1: Schematic illustrations of dynamic meshes around a tilted square. 3.5.1 Diffusion-based smoothing In FLUENT, dynamic meshing can be incorporated using smoothing methdos where the cells are moved with a deforming boundary while the number of cells and their connectivity remain unchanged. Smoothing is suitable for relatively small boundary deformations, while larger deformations may require generation of new cells in order to maintain a high quality mesh. One smoothing method is the diffusion-based smoothing, where the motion of the cells is modelled as a diffusive process. The diffusion equation for the mesh motion 20 3. Numerical methods is ∂ ∂xj ( Γm ∂Ui ∂xj ) = 0, (3.4) where Γm is a mesh diffusivity and ~U is the velocity field of the mesh. The mesh diffusivity is calculated from Γm = 1 dα , (3.5) where d is a normalized boundary distance and α is a diffusion parameter ranging from 0 to 2. With a low diffusion parameter, the boundary motion diffuses uniformly throughout the mesh, while a high value means that most of the mesh motion takes place far away from the moving boundaries. [27, pp. 608–613] 3.5.2 Overset method The overset method in STAR-CCM+ uses two overlapping meshes, one for the mov- ing part and one for the stationary background. The moving part, referred to as the overset mesh, uses the mesh rotation and translation method where the fluid mesh is replaced with a rigid body mesh. All cells maintain their shape and the mesh motion is described by a displacement vector and rotation angles. In the case when having a solid that interacts with the fluid, the position of the mesh is determined by solving the equations of the motion and rotation of the body. The background mesh is stationary and exchanges information with the moving mesh in the following way. First, the cells around the interface of the overset mesh are identified and labelled as donor cells. Then the cells in the background closest to the donor cells are identified and set as acceptor cells. These cells have to form a continuous layer of cells around the overset mesh. The background cells that are completely covered by the overset region are inactivated. The donor and acceptor cells transfer information between the meshes. Each acceptor cell has one or more donor cells. Choosing the donor cells can be done differently, the method used in this study is linear interpolation. The advantage with the overset method is that only a certain part of the mesh is moving without requirement for altering the grid topology. A drawback is that the interpolation between the meshes can cause numerical errors. [28, pp. 2662–2704] 3.6 Convergence criteria To be able to decide if a solution has reached a desirable level of convergence it is useful to monitor the residuals of the flow variables in each iteration. A residual is a measure of the imbalance between the left and right hand sides of a discretized transport equation. An unscaled residual, Rφ, of a solution can thereby be obtained by calculating the sum of the residuals in all cells according to Rφ = N∑ i=1 ∣∣∣∣∣∑ nb anbφnb + SU − aPφP ∣∣∣∣∣ , (3.6) 21 3. Numerical methods where N is the total number of cells. The value of this residual can vary between different variables. In order to compare the residuals and judge the convergence, they must be related to something. This can be done by scaling them with an appropriate factor to make the residuals into dimensionless numbers. This scaling is done differently in FLUENT and STAR- CCM+. [16] In FLUENT [27, pp. 1540–1541], the residuals are scaled by a characteristic flow rate of the variable φ in the domain. This yields a globally scaled residual, defined as Rφ,s = Rφ N∑ i=1 |aPφP | . (3.7) In STAR-CCM+ [28, pp. 7178–7179], the root mean square residual is calculated, Rrms = √√√√ 1 N N∑ i=1 R2 φ, (3.8) where the square root of the squared residuals in all cells are summed up and divided by the number of cells. To be able to compare the residuals, Rrms is scaled with a normalization value to get the present value, Rpres = Rrms Rnorm . (3.9) The normalization value, Rnorm, is often taken as the maximum of Rrms of the first m iterations. By default, m is equal to five. As seen in the equations above the initial guess affects how much the residuals can decrease. Besides looking at the residuals, the mass, momentum and energy imbalance in each cell can be checked. It is also important to monitor the solution of important variables to see if they reach a steady value. In hull simulations, the resistance, sinkage and trim angle are monitored until they oscillate with low amplitude around a steady value. 3.7 Grid dependence When a continuous equation is discretized using a finite number of computational cells and the set of equations are solved iteratively, numerical errors may arise. In order to draw any reliable conclusions from the results of a CFD computation, this numerical error has to be estimated. The numerical uncertainty is defined as EN = √ E2 I + E2 D, (3.10) 22 3. Numerical methods where EI is the iterative uncertainty and ED is the discretization uncertainty. If the equations have converged below a properly chosen tolerance, it may be assumed that the iterative uncertainty is negligible. This means that the discretization uncertainty remains to be estimated, which can be done with a grid dependence study. 3.7.1 The grid convergence index procedure A common method used for grid dependence studies is described in detail by Eça and Hoekstra [29], who used the grid convergence index (GCI) method. In the GCI method, a set of systematically varied grids is used to estimate how the numerical error depends on the grid resolution. The first step is to use generalized Richardson extrapolation in order to estimate the discretization error according to δRE = Si − S0 = qhri + ... ≈ qhri , (3.11) where higher order terms are neglected and their effects are assumed to be incorpo- rated in the leading term. In equation (3.11), Si is the solution on grid i, S0 is the extrapolated solution for an infinite grid density, q is a constant, hi is a characteristic cell length of grid i and r is the observed order of accuracy of the numerical method used in the calculations. The values of S0, q and r can be determined using the least squares method for three or more solutions. The least squares method is used to minimise the square root of the squares of the residuals of equation (3.11), in other words the method is used to find values of S0, q and r which minimise the function f(S0, q, r) = √∑ i (Si − (S0 + qhri )) 2. (3.12) The next step in the GCI method is to describe the discretization uncertainty based on the results from the Richardson extrapolation. The discretization uncertainty is expressed as ED = FS |δRE| , (3.13) where FS is a factor of safety. Depending on the nature of the convergence of the solutions, different expressions for the discretization uncertainty and thereby the numerical uncertainty are used. The numerical uncertainty is then normalized with the asymptotic value, S0, which makes it easy to compare different properties. A minimum of three different meshes are required in order to estimate the parameters of equation (3.12). It is recommended [30, p. 4] that the step size hi of the grids is varied with a factor of √ 2 to get an appropriate variation. The refinement ratio in relation to the finest grid can be expressed as rr = hi h1 , (3.14) where h1 is the step size of the finest mesh. For an unstructured mesh [31], the refinement ratio can instead be obtained from the cube root of the fraction of the 23 3. Numerical methods cell counts, rr = 3 √ N1 Ni , (3.15) where N1 and Ni are the total number of cells in grid 1 and grid i, respectively. 24 4 CFD simulations In this chapter, the methodology used in the CFD analysis is described. The pro- cedure follows the steps in Figure 1.1, including computational domain definition, mesh generation, model definitions and properties, boundary conditions, solution and post-processing. At the end, a section presenting the simulated operating con- ditions is included. 4.1 Computational domain definition A large domain was created in order to avoid effects from the domain boundaries to affect the flow near the hull. The Kelvin wake pattern, presented in Section 2.1.2, was used to estimate the required width of the domain to avoid wave reflections from the boundaries. Since only the dynamic sinkage and trim angle were studied, a symmetry plane was defined along the longitudinal axis of the hull and only half the hull was included in the domain. In Figure 4.1, the computational domain is illustrated and its dimensions are expressed in terms of the overall hull length, LOA. These dimensions agree well with the minimum recommendations of ITTC [24]. 5 2LOA 1LOA 2LOA4LOA 2LOA Figure 4.1: Dimensions of the computational domain. The geometry of the hull was provided as a computer-aided design (CAD) model. After the hull geometry had been imported to the software, it was rotated and translated in order to obtain suitable initial values of the sinkage and trim angle. This reduced the simulation time substantially since an initial position far from the equilibrium position would require a long simulation time before the hull reaches equilibrium. By subtracting the hull geometry from the rest of the domain, the do- main of the fluid flow was established. The origin was positioned at the undisturbed free surface level and horizontally aligned with the centre of gravity of the hull 25 4. CFD simulations 4.2 Mesh generation To capture the important flow phenomena in the simulations, the mesh density was focused on certain regions of the domain. First, a fine mesh was used on the hull sur- face and prism layers were created along the hull to resolve the boundary layer and to obtain correct shear stresses on the hull. The prism layers were constructed with a total thickness corresponding to the estimated boundary layer thickness, and the height of the first cell layer was set to get a proper value of y+. Secondly, to resolve the free surface accurately, a high mesh density was used in the region around the free surface where the hull induced wakes were expected to be present. To prevent smearing of the free surface in front of the hull, a uniform cell height was used at the free surface. Outside these regions, the mesh was coarser. The mesh used in FLUENT was created using ANSYS Meshing. It was divided into two regions; one inner region of unstructured tetrahedral cells in a cuboid around the hull, and one outer region of structured hexahedral cells. This mesh structure utilizes the flexibility of unstructured cells near the hull, while the number of cells could be kept low with a structured mesh in the periphery of the domain. The interface between the structured and the unstructured regions was made non- conformal, meaning that the cell faces do not match at the interface, to keep the cell count low. In Figure 4.2, the structure of the FLUENT mesh is illustrated. (a) Overview. (b) Symmetry plane. (c) Prism layers. (d) Rail. Figure 4.2: Schematic illustrations of the FLUENT mesh structure. (a) shows an overview of the mesh and (b) shows the mesh on the symmetry plane with the inner region of tetrahedral elements and the outer region of hexahedral elements. (c) shows the prism layers at the hull on the symmetry plane and (d) shows the triangular surface mesh at a rail on the hull. 26 4. CFD simulations The mesh constructed in STAR-CCM+ was divided into one stationary background region and one moving overset region close to the hull. Both parts were meshed with trimmed, predominantly hexahedral cells with local refinements at the free surface and the wake. Since the whole overset region sinks and trims, thicker refinement zones around the free surface at the overlap were needed in order to maintain a uniform cell height in front of the hull, which is seen in Figure 4.3b. Care had to be taken so that the cells in the overlapping region were of the same size and that they formed a continuous layer around the overset region. In Figure 4.3, the structure of the STAR-CCM+ mesh is illustrated. (a) Overview. Overset regionOverlap (b) Symmetry plane. (c) Prism layers. (d) Rail. Figure 4.3: Schematic illustrations of the STAR-CCM+ mesh structure. (a) shows an overview of the mesh and (b) shows the mesh on the symmetry plane with the overset region and the background mesh. (c) shows the prism layers at the hull on the symmetry plane and (d) shows the surface mesh at a rail on the hull. 4.3 Model definitions and properties This section describes the mathematical models and numerical methods used in the CFD simulations. 4.3.1 Mathematical models Based on the theoretical background presented in Chapter 2, appropriate models were chosen for the simulations. These models were used for describing the two- phase flow and the interaction between the flow and the hull. 27 4. CFD simulations 4.3.1.1 Two-phase flow models The turbulence was modelled using the RANS equations with the SST k-ω tur- bulence model. Wall functions were used in order to avoid resolving the whole boundary layer. In FLUENT, blended wall functions were used, and the all y+ wall treatment was used in STAR-CCM+. Both these methods use a blend of low and high y+ wall functions and are suitable for a wide range of y+ values. The first layer thickness was set to obtain a y+ value around 100. The free surface was modelled and resolved with the VOF method. The fluid prop- erties were set to be the same as in the model experiments that were used to validate the results. Both the air and water were assumed to be incompressible. 4.3.1.2 Hull motion To enable the 6 DOF solver to solve the equations of motion and rotation, the mass and the moments of the inertia were specified. The moments of inertia were esti- mated using a flat plate with the same outer dimensions as the hull. Since only the heave and pitch motions of the hull were simulated, the 6 DOF solver of the software was limited to 2 DOF. This was done by only allowing translational motion along the z-axis and rotational motion around the y-axis. To account for the pulling mechanism used in a towing tank, an additional force had to be included in the 6 DOF solver as this towing force induces an additional torque around the centre of mass. It was assumed that the pulling mechanism exerted a horizontal force on the hull, meaning that the towing force had no components in the y- and z-directions. When the hull is towed at constant speed, the towing force magnitude is thus equal to the total resistance of the hull. Therefore, a force with the same magnitude as the resistance but with opposite direction was applied in the towing point of the hull in each time step. In FLUENT, the towing force was incorporated using a user defined function (UDF) which looped through all cell faces on the hull and summarized the friction and pressure forces acting in the horizontal direction. When the total force was known, the resulting torque could be calculated from the trim angle and the distance between the towing point and the centre of gravity. This torque was applied in the UDF of the 6 DOF solver. In STAR-CCM+, the horizontal component of the resistance was reported using a field function that calculates the friction and pressure forces on the hull surface. This force was then applied in the towing point as an external force. 4.3.2 Numerical methods After the mathematical models had been chosen, numerical methods were selected in the software. When choosing spatial discretization schemes, the second order upwind scheme was chosen for all convection terms except for the volume fraction equation, where the compressive scheme was used in FLUENT and the HRIC scheme was used in STAR-CCM+. The diffusion terms were discretized with the central differencing scheme. Since the steady behaviour of the hull was simulated and tem- 28 4. CFD simulations poral accuracy was not important, the first order implicit scheme was chosen in the temporal discretization. For the pressure-velocity coupling, the SIMPLE algorithm was used. The mesh motion in FLUENT was modelled using diffusion-based smoothing with a diffusion parameter of 1.5. To stabilize the fluid-structure interaction, the solution stabilization option was enabled. In STAR-CCM+, the overset mesh method was used. 4.4 Boundary conditions The boundaries of the computational domain are illustrated in Figure 4.4. The symmetry plane as well as the top, side and bottom of the domain were prescribed with symmetry boundary conditions, and the hull was set to a wall with no-slip. Air inlet Water inlet Outlet Top Bottom Symmetry plane Side Figure 4.4: Boundaries of the computational domain. At the inlet, located in front of the hull, the velocity of the incident air and water was set to the hull speed that was simulated. The turbulent flow variables were set by specifying values of the turbulent intensity and turbulent viscosity ratio. The outlet, located behind the hull, was set to a pressure outlet. The location of the free surface was defined using the open channel flow model in FLUENT, and in STAR-CCM+ the level of the free surface was set with a coordi- nate and a normal vector. The static pressure at the inlet and outlet boundaries could thereby be obtained from the free surface location by specifying an atmo- spheric reference pressure in a point located above the free surface. An undesired phenomenon that is often observed in hull simulations is a thin layer of air in the first cell layers at the hull surface. This air layer, which affects the shear stresses on the hull, emerges from the region where the free surface hits the hull. It is not observed in experiments, and any air layer forming in the first cell layers at the hull should therefore be removed. In FLUENT, a UDF was used to loop through the prism layer adjacent to the hull while the volume fraction of air was set to 0 if it was below 0.3. In STAR-CCM+ a sink term was added to the transport equation for the volume fraction. The sink term was applied in the cells closest to the hull where the volume fraction of air was below 0.3. 29 4. CFD simulations 4.5 Solution The inlet boundary conditions were used to initialize the flow field. Although a steady state solution was expected, transient simulations were run to enable mesh motion and increase the robustness of the solution. The iterative solution procedure was first run with a fixed sinkage and trim in order to prevent large motions of the hull when the flow field was developed from the initial condition. When the hull was released, the solution was run until the sinkage, trim angle and forces acting on the hull had stabilized around constant values. Since these measured values kept oscillating with relatively small amplitudes around the equilibrium values, the sim- ulation was continued for some time in order to sample enough data for accurate predictions. The sampled data was processed by calculating the average values and the corresponding standard errors. By doing this, confidence intervals for the esti- mations could be determined. The methods used for the data analysis are described in Appendix A. 4.6 Post-processing When results had been obtained from a simulation, they were analysed and com- pared and new simulations were set up. In this iterative procedure, appropriate settings for the simulations were found. Grid dependence studies were conducted in order to estimate the numerical uncertainty. 4.7 Simulated operating conditions The CFD simulations were performed on two different hulls – the hull R/V ATHENA, hereafter referred to as the Athena hull, and a hull designed by Swede Ship, referred to as the Swede Ship hull. 4.7.1 Athena hull The Athena hull was chosen because experimental data from model tests was avail- able. Therefore, the models and setup used in the CFD simulations could be vali- dated. The Athena hull is a large semi-planing hull with a full scale LPP of 46.9m. In the simulations, which were only performed using STAR-CCM+, a 1:8.556 scale model was used. Since the hull does not trim a lot the tow force was neglected. A grid dependence study was conducted at a Froude number of 0.545 which corre- sponds to a velocity of 7.78 knots and a length between perpendiculars of 5.48m. There were some uncertainties in the exact position of the centre of gravity of the hull, both for the Athena and the Swede Ship hull. Therefore two cases, CG1 and CG2, with different positions of the centre of gravity were simulated for Athena in order to get an idea of how sensitive the results are to the position of centre of grav- ity. Figure 4.5 shows the outline of the hull with the location of the two different centres of gravity used for the simulations. The dimensions of the Athena hull are 30 4. CFD simulations found in Table 4.1 and the fluid properties used in the simulations are presented in Table 4.2. • CG2 • CG1 Waterline x z ⊗ y Figure 4.5: Sketch of the Athena hull showing the two different centres of gravity, CG1 and CG2, used in the simulations. Table 4.1: Hull dimensions of R/V Athena. A model scale with scaling factor 1:8.556 was used in the simulations. Two different centres of gravity were simulated, CG1 and CG2. Description Symbol Value Length between perpendiculars LPP 46.9m Volume displacement ∇ 257.5m3 Longitudinal CG from AP, CG1 LCG 20.16m Longitudinal CG from AP, CG2 LCG 17.97m Vertical CG, CG1 V CG 2.396m Vertical CG, CG2 V CG 0.684m Table 4.2: Fluid properties used in the simulations of the Athena hull. Description Symbol Value Atmospheric pressure at water surface Patm 101325Pa Density of water ρw 998.83 kg/m3 Viscosity of water νw 8.8871×10−4 m2/s Density of air ρa 1.18415 kg/m3 Viscosity of air νa 1.85508×10−5 m2/s 4.7.2 Swede Ship hull The Swede Ship hull is a planing hull with an LPP of 13.5m. The simulations were conducted in full scale, and the hull dimensions and fluid properties used in the simulations are found in Table 4.3 and Table 4.4. A more detailed table with properties is presented in Appendix C. 31 4. CFD simulations Table 4.3: Properties and dimensions of the Swede Ship hull at a water density of 1025.87 kg/m3 and a shell plating thickness of 8.00mm. Description Symbol Value Length between perpendiculars LPP 13.500m Breadth waterline BWL 4.157m Volume displacement ∇ 22.8m3 Longitudinal CG from AP LCB 4.921m Vertical CG V CG 1.9m Vertical CB V CB 0.6785m Table 4.4: Fluid properties used in the simulations of the Swede Ship hull. Description Symbol Value Atmospheric pressure at water surface Patm 101325Pa Density of water ρw 1025.87 kg/m3 Viscosity of water νw 1.188×10−6 m2/s Density of air ρa 1.225 kg/m3 Viscosity of air νa 1.461×10−5 m2/s The simulations of the Swede Ship hull were performed both with fixed and free sinkage and trim. The purpose of the fixed sinkage and trim simulations was to ob- tain a detailed comparison of the results from FLUENT and STAR-CCM+. These simulations were conducted at a Froude number of 1.68 with a sinkage at LPP/2 of 0.685m and a trim angle of 3.19 ◦ that were determined from a free sinkage and trim simulation in STAR-CCM+. In the simulations with free sinkage and trim, the Swede Ship hull was simulated with a centre of gravity and a centre of buoyancy as shown in Figure 4.6. The centre of gravity was provided in the hull data, and the towing point was stated to be in the centre of buoyancy at zero speed. In order to find the centre of buoyancy, the hull data was used to position the hull in its zero speed position, and then the centre of mass of the displaced volume was calculated. Due to some uncertainties in the exact position of the centre of gravity in the vertical direction, one simulation was performed where it was shifted 0.1m up. The results for the two cases, CG1 and CG2, were compared in order to evaluate if, and how much, the uncertain position of the centre of gravity affects the results. •CG • CB Waterline x z ⊗ y Figure 4.6: Sketch of the Swede Ship hull showing the centre of gravity (CG) and the centre of buoyancy (CB) at zero speed. 32 4. CFD simulations The free sinkage and trim simulations were conducted at Froude numbers ranging from 0.447 to 1.79, corresponding to velocities of 10 knots to 40 knots. A summary of the simulated cases is shown in Table 4.5 and Table 4.6. After the final results were obtained, they were compared with experimental data. Table 4.5: Simulated cases for the Swede Ship hull with fixed sinkage and trim. Fn Velocity [knots] Mesh study, FLUENT 1.68 37.5 FUENT 1.68 37.5 STAR-CCM+ 1.68 37.5 Table 4.6: Simulated cases for the Swede Ship hull with free sinkage and trim. Fn Velocity [knots] Mesh Study, STAR-CCM+ 1.68 37.5 STAR-CCM+, CG1 1.68 37.5 STAR-CCM+, CG2 1.68 37.5 STAR-CCM+ 1.79 40.0 STAR-CCM+ 1.68 37.5 STAR-CCM+ 1.57 35.0 STAR-CCM+ 1.34 30.0 STAR-CCM+ 1.12 25.0 STAR-CCM+ 0.894 20.0 STAR-CCM+ 0.671 15.0 STAR-CCM+ 0.447 10.0 FLUENT 1.68 37.5 33 4. CFD simulations 34 5 Results and discussion This chapter presents the results from the CFD simulations performed as described in Chapter 4. The results are divided into three parts: free sinkage and trim sim- ulations of the Athena hull, fixed sinkage and trim simulations of the Swede Ship hull and free sinkage and trim simulations of the Swede Ship hull. 5.1 Athena hull In the following sections the results from the simulations in STAR-CCM+ of the Athena hull with free sinkage and trim are presented and compared with experimen- tal data. Two cases, CG1 and CG2, with different locations of the centre of gravity were simulated. The simulations were performed at a Froude number of 0.545. First the grid dependence study is presented and then the comparison between experi- mental data and the calculated properties is presented and discussed. 5.1.1 Grid dependence study The grid dependence study was carried out for CG1. The grid was systematically refined in steps resulting in six different meshes. The number of cells for the coarsest mesh was 0.287 million and the finest mesh had 11.9 million cells. The number of cells for all meshes is presented in Table D.1 in Appendix D. The convergence of the total resistance coefficient was studied and the results are presented in Figure 5.1. It shows that the second finest grid with 7.10 million cells is a good compromise between numerical accuracy and computational effort. This mesh has a numerical uncertainty of 1.3% and was used to simulate CG2. 35 5. Results and discussion 0 0.5 1 1.5 2 2.5 3 3.5 4 6.2 6.3 6.4 6.5 6.6 6.7 x 10 −3 3 √ N1/N i C T Computed Fitted curve Figure 5.1: Convergence of the total resistance coefficient with grid refinement for the free sinkage and trim simulations of the Athena hull performed in STAR-CCM+ at Fn = 0.545. The full circle represents the grid used for simulating CG2. 5.1.2 Calculated properties The computed properties for the two cases are presented and compared with exper- imental data in Table 5.1. It is seen that the results from the simulation of CG1 give a lower sinkage and trim angle and underpredict the resistance compared to the experimental results. The results from CG2 give a sinkage that is closer to the experimental value but the trim angle is larger and the resistance is overpredicted. The results indicate that the chosen models and setup are appropriate. Table 5.1: Experimental and calculated properties of the Athena hull with free sinkage and trim at Fn = 0.545. Experiment CG1 CG2 Resistance [N] 244.52 229.14 ± 0.01 255.39 ± 0.22 Error -6.31% 4.25% Sinkage [m] -0.0105 -0.00874 ± 0.00010 -0.00968 ± 0.00017 Error -16.56% -4.44% Trim angle [◦] 1.35 1.26 ± 0.00 1.50 ± 0.00 Error -6.76% 10.05% Athena is a semi-planing hull, and as seen in Table 5.1 the sinkage and trim angle are very small. A much higher mesh density is needed than what was used in these simulations to be able to obtain the exact position of the hull. Figure 5.2 shows the volume fraction of air around the hull at the symmetry plane. 36 5. Results and discussion Volume fraction of air, γ 0 0.2 0.4 0.6 0.8 1.0 x z ⊗ y Figure 5.2: Volume fraction of air at the symmetry plane of the Athena hull at Fn = 0.545. It is seen in the results that the position of the centre of gravity affects the calculated resistance, sinkage and trim. Another reason to why the calculated values differ from the experimental values is that the towing force was neglected. This implies that it is important to know exactly how the experiments were performed in order to obtain accurate and reliable results. 37 5. Results and discussion 5.2 Swede Ship hull, fixed sinkage and trim The following sections present the results from the simulations with fixed sinkage and trim for the Swede Ship hull. The hull speed was set to 37.5 knots, correspond- ing to a Froude number of 1.68, and the sinkage at LPP/2 and trim angle were 0.685m and 3.19◦, respectively. After a presentation of the grid dependence study, the calculated volume fraction distribution, pressure coefficient, skin friction coeffi- cient, dimensionless wall distance and free surface pattern are visualized. At the end, the calculated forces acting on the hull and the wetted area and length are presented. When cross sections are used to present the results, they have been constructed from planes as illustrated in Figure 5.3. According to this figure, two longitudinal cross sections run along the hull surface at y/LPP = 0.0037 and y/LPP = 0.037, two longitudinal cross sections are placed outside the hull at y/LPP = 0.30 and y/LPP = 0.37, and two transverse cross sections run across the hull at x/LPP = −0.22 and x/LPP = −0.15. Symmetry y/LP P = 0.0037 y/LP P = 0.037 y/LP P = 0.30 y/LP P = 0.37 x/LP P = −0.22 x/LP P = −0.15 • O x y ⊗ z Figure 5.3: Cross sections at different y-coordinates used in the illustrations of the results. The origin of the computational domain is marked by O, and its position is in level with the undisturbed free surface. 5.2.1 Grid dependence study Four meshes with a cell count ranging from 4.28 million to 24.2 million cells were used in the grid dependence study conducted in Fluent. In Table D.2 in Appendix D, the number of cells in each mesh are presented. The calculated total resistance coefficients for these meshes are presented in Figure 5.4. The second finest mesh with 16.9 million cells and a numerical uncertainty of 1.7% was chosen for further analysis. In STAR-CCM+, a mesh with 18.4 million cells was selected based on a grid dependence study for free sinkage and trim. Consequently, the meshes used in the following comparison of the results in FLUENT and STAR-CCM+ have similar cell counts. 38 5. Results and discussion 0 0.5 1 1.5 2 2.85 2.9 2.95 3 3.05 3.1 x 10 −3 3 √ N1/Ni C T Computed Fitted curve Figure 5.4: Convergence of the total resistance coefficient with grid refinement for the fixed sinkage and trim simulations of the Swede Ship hull performed in FLUENT at Fn = 1.68. The full circle represents the grid used in the further analysis. 5.2.2 Volume fraction of air at the hull In Figure 5.5, contour plots of the volume fraction of air on the hull are presented. The differences between the results obtained in FLUENT and STAR-CCM+ are found where the free surface hits the hull. In particular, the region around the symmetry plane shows large differences. The air at the symmetry in both cases and at the rails in STAR-CCM+ is due to that the air that hits the hull at the free surface is dragged down along the hull. This phenomenon is prevented by removing air from cells which has a volume fraction of water higher than 0.7. As seen in the contour plots, more air could be removed, but then some realistic air might be removed and the sides of the hull can get unrealistically wet. The areas with volume fraction of air around 0.7 in the left picture might be due to splashing water. The area that is completely covered by water has the same shape in the two cases. Volume fraction of air, γair 0 0.2 0.4 0.6 0.8 1.0 (a) FLUENT (b) STAR-CCM+ x y ⊗ z Figure 5.5: Contour plots of the volume fraction of air on the hull for the fixed sinkage and trim simulations at Fn = 1.68. 39 5. Results and discussion 5.2.3 Pressure coefficient Figure 5.6 shows contour plots of the pressure coefficient on the hull. The results from FLUENT and STAR-CCM+ agree very well, and no apparent differences can be identified. A spot of significantly higher pressure is observed where the free surface hits the second innermost rail on the hull. Pressure coefficient, CP 0 0.08 0.16 0.24 0.32 0.40 (a) FLUENT (b) STAR-CCM+ x y ⊗ z Figure 5.6: Contour plots of the pressure coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. In Figure 5.7, longitudinal cross sections of the pressure coefficient on the hull are shown. In these plots, the similarity between the results is very clear – the curves barely differ from each other in the two cross sections. −0.4 −0.2 0 0.2 0.4 0 0.02 0.04 0.06 0.08 0.1 0.12 C P x/LPP FLUENT STAR-CCM+ (a) y/LPP = 0.0037. −0.4 −0.2 0 0.2 0.4 0 0.02 0.04 0.06 0.08 0.1 0.12 C P x/LPP FLUENT STAR-CCM+ (b) y/LPP = 0.037. Figure 5.7: Longitudinal cross sections of the pressure coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. In Figure 5.8, transverse cross sections of the pressure coefficient on the hull are shown. It can be observed that the pressure predicted by FLUENT is higher in the outer region of the hull around the lifting rails. The origin of this difference is not known, but it indicates that the regions around the lifting rails require special attention. 40 5. Results and discussion 0 0.05 0.1 0.15 0.2 0 0.02 0.04 0.06 0.08 0.1 0.12 C P y/LPP FLUENT STAR-CCM+ (a) x/LPP = −0.22. 0 0.05 0.1 0.15 0.2 0 0.02 0.04 0.06 0.08 0.1 0.12 C P y/LPP FLUENT STAR-CCM+ (b) x/LPP = −0.15. Figure 5.8: Transverse cross sections of the pressure coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. 5.2.4 Skin friction coefficient In Figure 5.9, contour plots of the skin friction coefficient on the hull are shown. In general, the skin friction in the region where the free surface hits the hull is more diffuse in FLUENT and the results from STAR-CCM+ show a more sharp increase in the shear stress. When comparing these contour plots to those of the volume fraction in Figure 5.5, it can be noted that the shear stress is highly correlated to the volume fraction distribution on the hull. Skin friction coefficient, Cf 0 0.0008 0.0016 0.0024 0.0032 0.0040 (a) FLUENT (b) STAR-CCM+ x y ⊗ z Figure 5.9: Contour plots of the skin friction coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. Figure 5.10 shows the skin friction coefficient on two longitudinal cross sections of the hull. The cross section at y/LPP = 0.0037, just at the symmetry plane, shows the relatively large difference between FLUENT and STAR-CCM+ in the front of the hull. Further away from the symmetry plane, at y/LPP =0.037m, the results are similar. 41 5. Results and discussion −0.4 −0.2 0 0.2 0.4 0 1 2 3 4 x 10 −3 C f x/LPP FLUENT STAR-CCM+ (a) y/LPP = 0.0037. −0.4 −0.2 0 0.2 0.4 0 1 2 3 4 x 10 −3 C f x/LPP FLUENT STAR-CCM+ (b) y/LPP = 0.037. Figure 5.10: Longitudinal cross sections of the skin friction coefficient on the hull for the fixed sinkage and trim simulations at Fn = 1.68. 5.2.5 Dimensionless wall distance Figure 5.11 shows contour plots of the dimensionless wall distance at the hull. Since the regions of air were not taken into account when choosing the first layer thickness of the prism layers, the y+ values are not close to the target value of 100 where air flows at the hull. However, since the air resistance constitutes a very small part of the total friction resistance, this deviation is of little importance. Dimensionless wall distance, y+ 0 30 60 90 120 150 (a) FLUENT (b) STAR-CCM+ x y ⊗ z Figure 5.11: Contour plots of the dimensionless wall distance on the hull for the fixed sinkage and trim simulations at Fn = 1.68. 42 5. Results and discussion 5.2.6 Free surface wave pattern In Figure 5.12, the dimensionless height of the free surface is shown. The free surface is defined as the interface between water and air where the volume fractions are 0.5. The wake patterns predicted by FLUENT and STAR-CCM+ are similar, but the result obtained in STAR-CCM+ shows more details near the hull. This is probably because of a higher mesh density around the hull in STAR-CCM+, where more effort was put into the construction of an appropriate mesh in this region. Dimensionless free surface height, z/LP P -0.060 -0.036 -0.012 0.012 0.036 0.060 (a) FLUENT (b) STAR-CCM+ x y � z Figure 5.12: Contour plots of the dimensionless free surface height for the fixed sinkage and trim simulations at Fn = 1.68. z/LPP = 0 is the height of the undis- turbed free surface. Figure 5.13 shows longitudinal cross sections of the free surface height outside the hull. In these figures it is clearer that the mesh created in STAR-CCM+ captures more of the waves and the splashing close to the hull. Further away from the hull, the results are similar. 43 5. Results and discussion −4 −3 −2 −1 0 1 2 −0.05 0 0.05 z/ L P P x/LPP FLUENT STAR-CCM+ (a) y/LPP = 0.30. −4 −3 −2 −1 0 1 2 −0.05 0 0.05 z/ L P P x/LPP FLUENT STAR-CCM+ (b) y/LPP = 0.37. Figure 5.13: Cross sections of the dimensionless free surface height at Fn = 1